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