Alc/lpfilter.c

   1 /* ----------------- file filterIIR00.c begin ----------------- */
   2 /*
   3 Resonant low pass filter source code.
   4 By baltrax@hotmail.com (Zxform)
   5 */
   6
   7 #include <stdlib.h>
   8 #include <stdio.h>
   9 #include <math.h>
  10
  11 #include "alMain.h"
  12 #include "alFilter.h"
  13
  14
  15 static void szxform(
  16     double *a0, double *a1, double *a2,   /* numerator coefficients */
  17     double *b0, double *b1, double *b2,   /* denominator coefficients */
  18     double fc,           /* Filter cutoff frequency */
  19     double fs,           /* sampling rate */
  20     double *k,           /* overall gain factor */
  21     float *coef);        /* pointer to 4 iir coefficients */
  22
  23 /*
  24  * --------------------------------------------------------------------
  25  *
  26  * lpFilter - Perform IIR filtering sample by sample on floats
  27  *
  28  * Implements cascaded direct form II second order sections.
  29  * Requires FILTER structure for history and coefficients.
  30  * The size of the history array is 2*FILTER_SECTIONS.
  31  * The size of the coefficient array is 4*FILTER_SECTIONS + 1 because
  32  * the first coefficient is the overall scale factor for the filter.
  33  * Returns one output sample for each input sample.
  34  *
  35  * float lpFilter(FILTER *iir,float input)
  36  *
  37  *     FILTER *iir        pointer to FILTER structure
  38  *     float input        new float input sample
  39  *
  40  * Returns float value giving the current output.
  41  * --------------------------------------------------------------------
  42  */
  43
  44 /*** moved to ALu.c ***/
  45
  46 /*
  47  * --------------------------------------------------------------------
  48  *
  49  * InitLowPassFilter()
  50  *
  51  * Initialize filter coefficients.
  52  * We create a 4th order filter (24 db/oct rolloff), consisting
  53  * of two second order sections.
  54  * --------------------------------------------------------------------
  55  */
  56 int InitLowPassFilter(ALCcontext *Context, FILTER *iir)
  57 {
  58     float    *coef;
  59     double   fs, fc; /* Sampling frequency, cutoff frequency */
  60     double   Q;      /* Resonance > 1.0 < 1000 */
  61     unsigned nInd;
  62     double   a0, a1, a2, b0, b1, b2;
  63     double   k;      /* overall gain factor */
  64     struct {
  65         double a0, a1, a2;       /* numerator coefficients */
  66         double b0, b1, b2;       /* denominator coefficients */
  67     } ProtoCoef[FILTER_SECTIONS];      /* Filter prototype coefficients,
  68                                           1 for each filter section */
  69
  70
  71     /*
  72      * Setup filter s-domain coefficients
  73      */
  74     /* Section 1 */
  75     ProtoCoef[0].a0 = 1.0;
  76     ProtoCoef[0].a1 = 0;
  77     ProtoCoef[0].a2 = 0;
  78     ProtoCoef[0].b0 = 1.0;
  79     ProtoCoef[0].b1 = 0.765367;
  80     ProtoCoef[0].b2 = 1.0;
  81
  82     /* Section 2 */
  83     ProtoCoef[1].a0 = 1.0;
  84     ProtoCoef[1].a1 = 0;
  85     ProtoCoef[1].a2 = 0;
  86     ProtoCoef[1].b0 = 1.0;
  87     ProtoCoef[1].b1 = 1.847759;
  88     ProtoCoef[1].b2 = 1.0;
  89
  90     /* Clear the coefficient and history arrays */
  91     memset(iir->coef, 0, sizeof(iir->coef));
  92     memset(iir->history, 0, sizeof(iir->history));
  93
  94     k = 1.0;              /* Set overall filter gain */
  95     coef = iir->coef + 1; /* Skip k, or gain */
  96
  97     Q = 1;                   /* Resonance */
  98     fc = LOWPASSFREQCUTOFF;  /* Filter cutoff (Hz) */
  99     fs = Context->Frequency; /* Sampling frequency (Hz) */
 100
 101     /*
 102      * Compute z-domain coefficients for each biquad section
 103      * for new Cutoff Frequency and Resonance
 104      */
 105     for (nInd = 0; nInd < FILTER_SECTIONS; nInd++)
 106     {
 107         a0 = ProtoCoef[nInd].a0;
 108         a1 = ProtoCoef[nInd].a1;
 109         a2 = ProtoCoef[nInd].a2;
 110
 111         b0 = ProtoCoef[nInd].b0;
 112         b1 = ProtoCoef[nInd].b1 / Q;      /* Divide by resonance or Q */
 113         b2 = ProtoCoef[nInd].b2;
 114         szxform(&a0, &a1, &a2, &b0, &b1, &b2, fc, fs, &k, coef);
 115         coef += 4;                       /* Point to next filter section */
 116     }
 117
 118     /* Update overall filter gain in coef array */
 119     iir->coef[0] = k;
 120
 121     return 0;
 122 }
 123 /* ----------------- file filterIIR00.c end ----------------- */
 124
 125
 126 /* ----------------- file bilinear.c begin ----------------- */
 127 /*
 128  * ----------------------------------------------------------
 129  *      bilinear.c
 130  *
 131  *      Perform bilinear transformation on s-domain coefficients
 132  *      of 2nd order biquad section.
 133  *      First design an analog filter and use s-domain coefficients
 134  *      as input to szxform() to convert them to z-domain.
 135  *
 136  * Here's the butterworth polinomials for 2nd, 4th and 6th order sections.
 137  *      When we construct a 24 db/oct filter, we take to 2nd order
 138  *      sections and compute the coefficients separately for each section.
 139  *
 140  *      n       Polinomials
 141  * --------------------------------------------------------------------
 142  *      2       s^2 + 1.4142s +1
 143  *      4       (s^2 + 0.765367s + 1) (s^2 + 1.847759s + 1)
 144  *      6       (s^2 + 0.5176387s + 1) (s^2 + 1.414214 + 1) (s^2 + 1.931852s + 1)
 145  *
 146  *      Where n is a filter order.
 147  *      For n=4, or two second order sections, we have following equasions for each
 148  *      2nd order stage:
 149  *
 150  *      (1 / (s^2 + (1/Q) * 0.765367s + 1)) * (1 / (s^2 + (1/Q) * 1.847759s + 1))
 151  *
 152  *      Where Q is filter quality factor in the range of
 153  *      1 to 1000. The overall filter Q is a product of all
 154  *      2nd order stages. For example, the 6th order filter
 155  *      (3 stages, or biquads) with individual Q of 2 will
 156  *      have filter Q = 2 * 2 * 2 = 8.
 157  *
 158  *      The nominator part is just 1.
 159  *      The denominator coefficients for stage 1 of filter are:
 160  *      b2 = 1; b1 = 0.765367; b0 = 1;
 161  *      numerator is
 162  *      a2 = 0; a1 = 0; a0 = 1;
 163  *
 164  *      The denominator coefficients for stage 1 of filter are:
 165  *      b2 = 1; b1 = 1.847759; b0 = 1;
 166  *      numerator is
 167  *      a2 = 0; a1 = 0; a0 = 1;
 168  *
 169  *      These coefficients are used directly by the szxform()
 170  *      and bilinear() functions. For all stages the numerator
 171  *      is the same and the only thing that is different between
 172  *      different stages is 1st order coefficient. The rest of
 173  *      coefficients are the same for any stage and equal to 1.
 174  *
 175  *      Any filter could be constructed using this approach.
 176  *
 177  *      References:
 178  *             Van Valkenburg, "Analog Filter Design"
 179  *             Oxford University Press 1982
 180  *             ISBN 0-19-510734-9
 181  *
 182  *             C Language Algorithms for Digital Signal Processing
 183  *             Paul Embree, Bruce Kimble
 184  *             Prentice Hall, 1991
 185  *             ISBN 0-13-133406-9
 186  *
 187  *             Digital Filter Designer's Handbook
 188  *             With C++ Algorithms
 189  *             Britton Rorabaugh
 190  *             McGraw Hill, 1997
 191  *             ISBN 0-07-053806-9
 192  * ----------------------------------------------------------
 193  */
 194
 195 static void prewarp(double *a0, double *a1, double *a2, double fc, double fs);
 196 static void bilinear(
 197     double a0, double a1, double a2,    /* numerator coefficients */
 198     double b0, double b1, double b2,    /* denominator coefficients */
 199     double *k,                          /* overall gain factor */
 200     double fs,                          /* sampling rate */
 201     float *coef);                       /* pointer to 4 iir coefficients */
 202
 203
 204 /*
 205  * ----------------------------------------------------------
 206  *      Pre-warp the coefficients of a numerator or denominator.
 207  *      Note that a0 is assumed to be 1, so there is no wrapping
 208  *      of it.
 209  * ----------------------------------------------------------
 210  */
 211 static void prewarp(
 212     double *a0, double *a1, double *a2,
 213     double fc, double fs)
 214 {
 215     double wp, pi;
 216
 217     pi = 4.0 * atan(1.0);
 218     wp = 2.0 * fs * tan(pi * fc / fs);
 219
 220     *a2 = (*a2) / (wp * wp);
 221     *a1 = (*a1) / wp;
 222     (void)a0;
 223 }
 224
 225
 226 /*
 227  * ----------------------------------------------------------
 228  * bilinear()
 229  *
 230  * Transform the numerator and denominator coefficients
 231  * of s-domain biquad section into corresponding
 232  * z-domain coefficients.
 233  *
 234  *      Store the 4 IIR coefficients in array pointed by coef
 235  *      in following order:
 236  *             beta1, beta2    (denominator)
 237  *             alpha1, alpha2  (numerator)
 238  *
 239  * Arguments:
 240  *             a0-a2   - s-domain numerator coefficients
 241  *             b0-b2   - s-domain denominator coefficients
 242  *             k       - filter gain factor. initially set to 1
 243  *                       and modified by each biquad section in such
 244  *                       a way, as to make it the coefficient by
 245  *                       which to multiply the overall filter gain
 246  *                       in order to achieve a desired overall filter gain,
 247  *                       specified in initial value of k.
 248  *             fs      - sampling rate (Hz)
 249  *             coef    - array of z-domain coefficients to be filled in.
 250  *
 251  * Return:
 252  *             On return, set coef z-domain coefficients
 253  * ----------------------------------------------------------
 254  */
 255 static void bilinear(
 256     double a0, double a1, double a2,    /* numerator coefficients */
 257     double b0, double b1, double b2,    /* denominator coefficients */
 258     double *k,                          /* overall gain factor */
 259     double fs,                          /* sampling rate */
 260     float *coef                         /* pointer to 4 iir coefficients */
 261 )
 262 {
 263     double ad, bd;
 264
 265     /* alpha (Numerator in s-domain) */
 266     ad = 4. * a2 * fs * fs + 2. * a1 * fs + a0;
 267     /* beta (Denominator in s-domain) */
 268     bd = 4. * b2 * fs * fs + 2. * b1* fs + b0;
 269
 270     /* update gain constant for this section */
 271     *k *= ad/bd;
 272
 273     /* Denominator */
 274     *coef++ = (2.*b0 - 8.*b2*fs*fs) / bd;         /* beta1 */
 275     *coef++ = (4.*b2*fs*fs - 2.*b1*fs + b0) / bd; /* beta2 */
 276
 277     /* Nominator */
 278     *coef++ = (2.*a0 - 8.*a2*fs*fs) / ad;         /* alpha1 */
 279     *coef = (4.*a2*fs*fs - 2.*a1*fs + a0) / ad;   /* alpha2 */
 280 }
 281
 282
 283 /*
 284  * ----------------------------------------------------------
 285  * Transform from s to z domain using bilinear transform
 286  * with prewarp.
 287  *
 288  * Arguments:
 289  *      For argument description look at bilinear()
 290  *
 291  *      coef - pointer to array of floating point coefficients,
 292  *                     corresponding to output of bilinear transofrm
 293  *                     (z domain).
 294  *
 295  * Note: frequencies are in Hz.
 296  * ----------------------------------------------------------
 297  */
 298 static void szxform(
 299     double *a0, double *a1, double *a2, /* numerator coefficients */
 300     double *b0, double *b1, double *b2, /* denominator coefficients */
 301     double fc,                          /* Filter cutoff frequency */
 302     double fs,                          /* sampling rate */
 303     double *k,                          /* overall gain factor */
 304     float *coef)                        /* pointer to 4 iir coefficients */
 305 {
 306     /* Calculate a1 and a2 and overwrite the original values */
 307     prewarp(a0, a1, a2, fc, fs);
 308     prewarp(b0, b1, b2, fc, fs);
 309     bilinear(*a0, *a1, *a2, *b0, *b1, *b2, k, fs, coef);
 310 }
 311
 312
 313 /* ----------------- file bilinear.c end ----------------- */
 314
 315 /* ----------------- file filter.txt begin -----------------
 316 How to construct a kewl low pass resonant filter?
 317
 318 Lets assume we want to create a filter for analog synth.
 319 The filter rolloff is 24 db/oct, which corresponds to 4th
 320 order filter. Filter of first order is equivalent to RC circuit
 321 and has max rolloff of 6 db/oct.
 322
 323 We will use classical Butterworth IIR filter design, as it
 324 exactly corresponds to our requirements.
 325
 326 A common practice is to chain several 2nd order sections,
 327 or biquads, as they commonly called, in order to achive a higher
 328 order filter. Each 2nd order section is a 2nd order filter, which
 329 has 12 db/oct roloff. So, we need 2 of those sections in series.
 330
 331 To compute those sections, we use standard Butterworth polinomials,
 332 or so called s-domain representation and convert it into z-domain,
 333 or digital domain. The reason we need to do this is because
 334 the filter theory exists for analog filters for a long time
 335 and there exist no theory of working in digital domain directly.
 336 So the common practice is to take standard analog filter design
 337 and use so called bilinear transform to convert the butterworth
 338 equasion coefficients into z-domain.
 339
 340 Once we compute the z-domain coefficients, we can use them in
 341 a very simple transfer function, such as iir_filter() in our
 342 C source code, in order to perform the filtering function.
 343 The filter itself is the simpliest thing in the world.
 344 The most complicated thing is computing the coefficients
 345 for z-domain.
 346
 347 Ok, lets look at butterworth polynomials, arranged as a series
 348 of 2nd order sections:
 349
 350  * Note: n is filter order.
 351  *
 352  *      n       Polynomials
 353  *      --------------------------------------------------------------------
 354  *      2       s^2 + 1.4142s +1
 355  *      4       (s^2 + 0.765367s + 1) * (s^2 + 1.847759s + 1)
 356  *      6       (s^2 + 0.5176387s + 1) * (s^2 + 1.414214 + 1) * (s^2 + 1.931852s + 1)
 357  *
 358  * For n=4 we have following equasion for the filter transfer function:
 359  *
 360  *                     1                              1
 361  * T(s) = --------------------------- * ----------------------------
 362  *        s^2 + (1/Q) * 0.765367s + 1   s^2 + (1/Q) * 1.847759s + 1
 363  *
 364
 365 The filter consists of two 2nd order secions since highest s power is 2.
 366 Now we can take the coefficients, or the numbers by which s is multiplied
 367 and plug them into a standard formula to be used by bilinear transform.
 368
 369 Our standard form for each 2nd order secion is:
 370
 371        a2 * s^2 + a1 * s + a0
 372 H(s) = ----------------------
 373        b2 * s^2 + b1 * s + b0
 374
 375 Note that butterworth nominator is 1 for all filter sections,
 376 which means s^2 = 0 and s^1 = 0
 377
 378 Lets convert standard butterworth polinomials into this form:
 379
 380        0 + 0 + 1                  0 + 0 + 1
 381 -------------------------- * --------------------------
 382 1 + ((1/Q) * 0.765367) + 1   1 + ((1/Q) * 1.847759) + 1
 383
 384 Section 1:
 385 a2 = 0; a1 = 0; a0 = 1;
 386 b2 = 1; b1 = 0.5176387; b0 = 1;
 387
 388 Section 2:
 389 a2 = 0; a1 = 0; a0 = 1;
 390 b2 = 1; b1 = 1.847759; b0 = 1;
 391
 392 That Q is filter quality factor or resonance, in the range of
 393 1 to 1000. The overall filter Q is a product of all 2nd order stages.
 394 For example, the 6th order filter (3 stages, or biquads)
 395 with individual Q of 2 will have filter Q = 2 * 2 * 2 = 8.
 396
 397 These a and b coefficients are used directly by the szxform()
 398 and bilinear() functions.
 399
 400 The transfer function for z-domain is:
 401
 402        1 + alpha1 * z^(-1) + alpha2 * z^(-2)
 403 H(z) = -------------------------------------
 404        1 + beta1 * z^(-1) + beta2 * z^(-2)
 405
 406 When you need to change the filter frequency cutoff or resonance,
 407 or Q, you call the szxform() function with proper a and b
 408 coefficients and the new filter cutoff frequency or resonance.
 409 You also need to supply the sampling rate and filter gain you want
 410 to achive. For our purposes the gain = 1.
 411
 412 We call szxform() function 2 times becase we have 2 filter sections.
 413 Each call provides different coefficients.
 414
 415 The gain argument to szxform() is a pointer to desired filter
 416 gain variable.
 417
 418 double k = 1.0;      // overall gain factor
 419
 420 Upon return from each call, the k argument will be set to a value,
 421 by which to multiply our actual signal in order for the gain
 422 to be one. On second call to szxform() we provide k that was
 423 changed by the previous section. During actual audio filtering
 424 function iir_filter() will use this k
 425
 426 Summary:
 427
 428 Our filter is pretty close to ideal in terms of all relevant
 429 parameters and filter stability even with extremely large values
 430 of resonance. This filter design has been verified under all
 431 variations of parameters and it all appears to work as advertized.
 432
 433 Good luck with it.
 434 If you ever make a directX wrapper for it, post it to comp.dsp.
 435
 436
 437  *
 438  * ----------------------------------------------------------
 439  *References:
 440  *Van Valkenburg, "Analog Filter Design"
 441  *Oxford University Press 1982
 442  *ISBN 0-19-510734-9
 443  *
 444  *C Language Algorithms for Digital Signal Processing
 445  *Paul Embree, Bruce Kimble
 446  *Prentice Hall, 1991
 447  *ISBN 0-13-133406-9
 448  *
 449  *Digital Filter Designer's Handbook
 450  *With C++ Algorithms
 451  *Britton Rorabaugh
 452  *McGraw Hill, 1997
 453  *ISBN 0-07-053806-9
 454  * ----------------------------------------------------------
 455
 456
 457
 458 // ----------------- file filter.txt end ----------------- */