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