Implement yet another low-pass filter

This one using the Butterworth IIR filter design

2 files changed, 518 insertions, 16 deletions

diff --git a/Alc/ALu.c b/Alc/ALu.c
index 43837123..ea0f0dcb 100644
--- a/Alc/ALu.c
+++ b/Alc/ALu.c
@@ -596,8 +596,6 @@ ALvoid aluMixData(ALCcontext *ALContext,ALvoid *buffer,ALsizei size,ALenum forma
     ALfloat Pitch;
     ALint Looping,State;
     ALint fraction,increment;
-    ALint LowFrac;
-    ALuint LowStep;
     ALuint Buffer;
     ALuint SamplesToDo;
     ALsource *ALSource;
@@ -609,6 +607,7 @@ ALvoid aluMixData(ALCcontext *ALContext,ALvoid *buffer,ALsizei size,ALenum forma
     ALbufferlistitem *BufferListItem;
     ALuint loop;
     ALint64 DataSize64,DataPos64;
+    FILTER *Filter;
 
     SuspendContext(ALContext);
 
@@ -667,6 +666,7 @@ ALvoid aluMixData(ALCcontext *ALContext,ALvoid *buffer,ALsizei size,ALenum forma
                     //Get source info
                     DataPosInt = ALSource->position;
                     DataPosFrac = ALSource->position_fraction;
+                    Filter = &ALSource->iirFilter;
 
                     //Compute 18.14 fixed point step
                     increment = (ALint)(Pitch*(ALfloat)(1L<<FRACTIONBITS));
@@ -680,8 +680,6 @@ ALvoid aluMixData(ALCcontext *ALContext,ALvoid *buffer,ALsizei size,ALenum forma
                     DataPos64 <<= FRACTIONBITS;
                     DataPos64 += DataPosFrac;
                     BufferSize = (ALuint)((DataSize64-DataPos64+(increment-1)) / increment);
-                    LowStep = Frequency/LOWPASSFREQCUTOFF;
-                    if(LowStep < 1) LowStep = 1;
 
                     BufferListItem = ALSource->queue;
                     for(loop = 0; loop < ALSource->BuffersPlayed; loop++)
@@ -718,16 +716,14 @@ ALvoid aluMixData(ALCcontext *ALContext,ALvoid *buffer,ALsizei size,ALenum forma
                     {
                         k = DataPosFrac>>FRACTIONBITS;
                         fraction = DataPosFrac&FRACTIONMASK;
-                        LowFrac = ((DataPosFrac+(DataPosInt<<FRACTIONBITS))/LowStep)&FRACTIONMASK;
+
                         if(Channels==1)
                         {
                             ALfloat sample, lowsamp, outsamp;
                             //First order interpolator
                             sample = (Data[k]*((1<<FRACTIONBITS)-fraction) +
                                       Data[k+1]*fraction) >> FRACTIONBITS;
-                            lowsamp = (Data[((k+DataPosInt)/LowStep    )*LowStep - DataPosInt]*((1<<FRACTIONBITS)-LowFrac) +
-                                       Data[((k+DataPosInt)/LowStep + 1)*LowStep - DataPosInt]*LowFrac) >>
-                                      FRACTIONBITS;
+                            lowsamp = lpFilter(Filter, sample);
 
                             //Direct path final mix buffer and panning
                             outsamp = aluComputeSample(DryGainHF, sample, lowsamp);
@@ -915,9 +911,9 @@ ALvoid aluMixData(ALCcontext *ALContext,ALvoid *buffer,ALsizei size,ALenum forma
                 ALfloat Gain = ALEffectSlot->effect.Reverb.Gain;
                 ALfloat ReflectGain = ALEffectSlot->effect.Reverb.ReflectionsGain;
                 ALfloat LateReverbGain = ALEffectSlot->effect.Reverb.LateReverbGain;
-                ALfloat LastDecaySample = ALEffectSlot->LastDecaySample;
-                ALfloat sample;
+                ALfloat sample, lowsample;
 
+                Filter = &ALEffectSlot->iirFilter;
                 for(i = 0;i < SamplesToDo;i++)
                 {
                     DelayBuffer[Pos] = ReverbBuffer[i] * Gain;
@@ -927,12 +923,10 @@ ALvoid aluMixData(ALCcontext *ALContext,ALvoid *buffer,ALsizei size,ALenum forma
                     DelayBuffer[LatePos] *= LateReverbGain;
 
                     Pos = (Pos+1) % Length;
-                    DelayBuffer[Pos] *= DecayHFRatio;
-                    DelayBuffer[Pos] += LastDecaySample * (1.0f-DecayHFRatio);
-                    LastDecaySample = DelayBuffer[Pos];
-                    DelayBuffer[Pos] *= DecayGain;
+                    lowsample = lpFilter(Filter, DelayBuffer[Pos]);
+                    lowsample += (DelayBuffer[Pos]-lowsample) * DecayHFRatio;
 
-                    DelayBuffer[LatePos] += DelayBuffer[Pos];
+                    DelayBuffer[LatePos] += lowsample * DecayGain;
 
                     sample += DelayBuffer[LatePos];
 
@@ -950,7 +944,6 @@ ALvoid aluMixData(ALCcontext *ALContext,ALvoid *buffer,ALsizei size,ALenum forma
                 ALEffectSlot->ReverbPos = Pos;
                 ALEffectSlot->ReverbLatePos = LatePos;
                 ALEffectSlot->ReverbReflectPos = ReflectPos;
-                ALEffectSlot->LastDecaySample = LastDecaySample;
             }
 
             ALEffectSlot = ALEffectSlot->next;
diff --git a/Alc/lpfilter.c b/Alc/lpfilter.c
new file mode 100644
index 00000000..cf94c9f6
--- /dev/null
+++ b/Alc/lpfilter.c
@@ -0,0 +1,509 @@
+/* ----------------- file filterIIR00.c begin ----------------- */
+/* 
+Resonant low pass filter source code.
+By [email protected] (Zxform)
+*/
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+
+#include "alMain.h"
+#include "alFilter.h"
+
+
+#define FILTER_SECTIONS   2   /* 2 filter sections for 24 db/oct filter */
+
+
+static void szxform(
+    double *a0, double *a1, double *a2,   /* numerator coefficients */
+    double *b0, double *b1, double *b2,   /* denominator coefficients */
+    double fc,           /* Filter cutoff frequency */
+    double fs,           /* sampling rate */
+    double *k,           /* overall gain factor */
+    float *coef);        /* pointer to 4 iir coefficients */
+
+/*
+ * --------------------------------------------------------------------
+ *
+ * lpFilter - Perform IIR filtering sample by sample on floats
+ *
+ * Implements cascaded direct form II second order sections.
+ * Requires FILTER structure for history and coefficients.
+ * The size of the history array is 2*FILTER_SECTIONS.
+ * The size of the coefficient array is 4*FILTER_SECTIONS + 1 because
+ * the first coefficient is the overall scale factor for the filter.
+ * Returns one output sample for each input sample.
+ *
+ * float lpFilter(FILTER *iir,float input)
+ *
+ *     FILTER *iir        pointer to FILTER structure
+ *     float input        new float input sample
+ *
+ * Returns float value giving the current output.
+ * --------------------------------------------------------------------
+ */
+float lpFilter(FILTER *iir, float input)
+{
+    unsigned int i;
+    float *hist1_ptr,*hist2_ptr,*coef_ptr;
+    float output,new_hist,history1,history2;
+
+    coef_ptr = iir->coef;                /* coefficient pointer */
+
+    hist1_ptr = iir->history;            /* first history */
+    hist2_ptr = hist1_ptr + 1;           /* next history */
+
+    /* 1st number of coefficients array is overall input scale factor,
+     * or filter gain */
+    output = input * (*coef_ptr++);
+
+    for(i = 0;i < FILTER_SECTIONS;i++)
+    {
+        history1 = *hist1_ptr;           /* history values */
+        history2 = *hist2_ptr;
+
+        output = output - history1 * (*coef_ptr++);
+        new_hist = output - history2 * (*coef_ptr++);    /* poles */
+
+        output = new_hist + history1 * (*coef_ptr++);
+        output = output + history2 * (*coef_ptr++);      /* zeros */
+
+        *hist2_ptr++ = *hist1_ptr;
+        *hist1_ptr++ = new_hist;
+        hist1_ptr++;
+        hist2_ptr++;
+    }
+
+    return output;
+}
+
+
+/*
+ * --------------------------------------------------------------------
+ *
+ * InitLowPassFilter()
+ *
+ * Initialize filter coefficients.
+ * We create a 4th order filter (24 db/oct rolloff), consisting
+ * of two second order sections.
+ * --------------------------------------------------------------------
+ */
+int InitLowPassFilter(ALCcontext *Context, FILTER *iir)
+{
+    float    *coef;
+    double   fs, fc; /* Sampling frequency, cutoff frequency */
+    double   Q;      /* Resonance > 1.0 < 1000 */
+    unsigned nInd;
+    double   a0, a1, a2, b0, b1, b2;
+    double   k;      /* overall gain factor */
+    struct {
+        double a0, a1, a2;       /* numerator coefficients */
+        double b0, b1, b2;       /* denominator coefficients */
+    } ProtoCoef[FILTER_SECTIONS];      /* Filter prototype coefficients,
+                                          1 for each filter section */
+
+
+    /*
+     * Setup filter s-domain coefficients
+     */
+    /* Section 1 */
+    ProtoCoef[0].a0 = 1.0;
+    ProtoCoef[0].a1 = 0;
+    ProtoCoef[0].a2 = 0;
+    ProtoCoef[0].b0 = 1.0;
+    ProtoCoef[0].b1 = 0.765367;
+    ProtoCoef[0].b2 = 1.0;
+
+    /* Section 2 */
+    ProtoCoef[1].a0 = 1.0;
+    ProtoCoef[1].a1 = 0;
+    ProtoCoef[1].a2 = 0;
+    ProtoCoef[1].b0 = 1.0;
+    ProtoCoef[1].b1 = 1.847759;
+    ProtoCoef[1].b2 = 1.0;
+
+    /*
+     * Allocate array of z-domain coefficients for each filter section
+     * plus filter gain variable
+     */
+    iir->coef = (float*)calloc(4*FILTER_SECTIONS + 1, sizeof(float));
+    if(!iir->coef)
+    {
+        AL_PRINT("Unable to allocate coef array\n");
+        return 1;
+    }
+
+    /* allocate history array */
+    iir->history = (float*)calloc(2*FILTER_SECTIONS, sizeof(float));
+    if(!iir->history) {
+        AL_PRINT("\nUnable to allocate history array\n");
+        return 1;
+    }
+
+
+    k = 1.0;              /* Set overall filter gain */
+    coef = iir->coef + 1; /* Skip k, or gain */
+
+    Q = 1;                   /* Resonance */
+    fc = LOWPASSFREQCUTOFF;  /* Filter cutoff (Hz) */
+    fs = Context->Frequency; /* Sampling frequency (Hz) */
+
+    /*
+     * Compute z-domain coefficients for each biquad section
+     * for new Cutoff Frequency and Resonance
+     */
+    for (nInd = 0; nInd < FILTER_SECTIONS; nInd++)
+    {
+        a0 = ProtoCoef[nInd].a0;
+        a1 = ProtoCoef[nInd].a1;
+        a2 = ProtoCoef[nInd].a2;
+
+        b0 = ProtoCoef[nInd].b0;
+        b1 = ProtoCoef[nInd].b1 / Q;      /* Divide by resonance or Q */
+        b2 = ProtoCoef[nInd].b2;
+        szxform(&a0, &a1, &a2, &b0, &b1, &b2, fc, fs, &k, coef);
+        coef += 4;                       /* Point to next filter section */
+    }
+
+    /* Update overall filter gain in coef array */
+    iir->coef[0] = k;
+
+    return 0;
+}
+/* ----------------- file filterIIR00.c end ----------------- */
+
+
+/* ----------------- file bilinear.c begin ----------------- */
+/*
+ * ----------------------------------------------------------
+ *      bilinear.c
+ *
+ *      Perform bilinear transformation on s-domain coefficients
+ *      of 2nd order biquad section.
+ *      First design an analog filter and use s-domain coefficients
+ *      as input to szxform() to convert them to z-domain.
+ *
+ * Here's the butterworth polinomials for 2nd, 4th and 6th order sections.
+ *      When we construct a 24 db/oct filter, we take to 2nd order
+ *      sections and compute the coefficients separately for each section.
+ *
+ *      n       Polinomials
+ * --------------------------------------------------------------------
+ *      2       s^2 + 1.4142s +1
+ *      4       (s^2 + 0.765367s + 1) (s^2 + 1.847759s + 1)
+ *      6       (s^2 + 0.5176387s + 1) (s^2 + 1.414214 + 1) (s^2 + 1.931852s + 1)
+ *
+ *      Where n is a filter order.
+ *      For n=4, or two second order sections, we have following equasions for each
+ *      2nd order stage:
+ *
+ *      (1 / (s^2 + (1/Q) * 0.765367s + 1)) * (1 / (s^2 + (1/Q) * 1.847759s + 1))
+ *
+ *      Where Q is filter quality factor in the range of
+ *      1 to 1000. The overall filter Q is a product of all
+ *      2nd order stages. For example, the 6th order filter
+ *      (3 stages, or biquads) with individual Q of 2 will
+ *      have filter Q = 2 * 2 * 2 = 8.
+ *
+ *      The nominator part is just 1.
+ *      The denominator coefficients for stage 1 of filter are:
+ *      b2 = 1; b1 = 0.765367; b0 = 1;
+ *      numerator is
+ *      a2 = 0; a1 = 0; a0 = 1;
+ *
+ *      The denominator coefficients for stage 1 of filter are:
+ *      b2 = 1; b1 = 1.847759; b0 = 1;
+ *      numerator is
+ *      a2 = 0; a1 = 0; a0 = 1;
+ *
+ *      These coefficients are used directly by the szxform()
+ *      and bilinear() functions. For all stages the numerator
+ *      is the same and the only thing that is different between
+ *      different stages is 1st order coefficient. The rest of
+ *      coefficients are the same for any stage and equal to 1.
+ *
+ *      Any filter could be constructed using this approach.
+ *
+ *      References:
+ *             Van Valkenburg, "Analog Filter Design"
+ *             Oxford University Press 1982
+ *             ISBN 0-19-510734-9
+ *
+ *             C Language Algorithms for Digital Signal Processing
+ *             Paul Embree, Bruce Kimble
+ *             Prentice Hall, 1991
+ *             ISBN 0-13-133406-9
+ *
+ *             Digital Filter Designer's Handbook
+ *             With C++ Algorithms
+ *             Britton Rorabaugh
+ *             McGraw Hill, 1997
+ *             ISBN 0-07-053806-9
+ * ----------------------------------------------------------
+ */
+
+static void prewarp(double *a0, double *a1, double *a2, double fc, double fs);
+static void bilinear(
+    double a0, double a1, double a2,    /* numerator coefficients */
+    double b0, double b1, double b2,    /* denominator coefficients */
+    double *k,                          /* overall gain factor */
+    double fs,                          /* sampling rate */
+    float *coef);                       /* pointer to 4 iir coefficients */
+
+
+/*
+ * ----------------------------------------------------------
+ *      Pre-warp the coefficients of a numerator or denominator.
+ *      Note that a0 is assumed to be 1, so there is no wrapping
+ *      of it.
+ * ----------------------------------------------------------
+ */
+static void prewarp(
+    double *a0, double *a1, double *a2,
+    double fc, double fs)
+{
+    double wp, pi;
+
+    pi = 4.0 * atan(1.0);
+    wp = 2.0 * fs * tan(pi * fc / fs);
+
+    *a2 = (*a2) / (wp * wp);
+    *a1 = (*a1) / wp;
+    (void)a0;
+}
+
+
+/*
+ * ----------------------------------------------------------
+ * bilinear()
+ *
+ * Transform the numerator and denominator coefficients
+ * of s-domain biquad section into corresponding
+ * z-domain coefficients.
+ *
+ *      Store the 4 IIR coefficients in array pointed by coef
+ *      in following order:
+ *             beta1, beta2    (denominator)
+ *             alpha1, alpha2  (numerator)
+ *
+ * Arguments:
+ *             a0-a2   - s-domain numerator coefficients
+ *             b0-b2   - s-domain denominator coefficients
+ *             k       - filter gain factor. initially set to 1
+ *                       and modified by each biquad section in such
+ *                       a way, as to make it the coefficient by
+ *                       which to multiply the overall filter gain
+ *                       in order to achieve a desired overall filter gain,
+ *                       specified in initial value of k.
+ *             fs      - sampling rate (Hz)
+ *             coef    - array of z-domain coefficients to be filled in.
+ *
+ * Return:
+ *             On return, set coef z-domain coefficients
+ * ----------------------------------------------------------
+ */
+static void bilinear(
+    double a0, double a1, double a2,    /* numerator coefficients */
+    double b0, double b1, double b2,    /* denominator coefficients */
+    double *k,                          /* overall gain factor */
+    double fs,                          /* sampling rate */
+    float *coef                         /* pointer to 4 iir coefficients */
+)
+{
+    double ad, bd;
+
+    /* alpha (Numerator in s-domain) */
+    ad = 4. * a2 * fs * fs + 2. * a1 * fs + a0;
+    /* beta (Denominator in s-domain) */
+    bd = 4. * b2 * fs * fs + 2. * b1* fs + b0;
+
+    /* update gain constant for this section */
+    *k *= ad/bd;
+
+    /* Denominator */
+    *coef++ = (2.*b0 - 8.*b2*fs*fs) / bd;         /* beta1 */
+    *coef++ = (4.*b2*fs*fs - 2.*b1*fs + b0) / bd; /* beta2 */
+
+    /* Nominator */
+    *coef++ = (2.*a0 - 8.*a2*fs*fs) / ad;         /* alpha1 */
+    *coef = (4.*a2*fs*fs - 2.*a1*fs + a0) / ad;   /* alpha2 */
+}
+
+
+/*
+ * ----------------------------------------------------------
+ * Transform from s to z domain using bilinear transform
+ * with prewarp.
+ *
+ * Arguments:
+ *      For argument description look at bilinear()
+ *
+ *      coef - pointer to array of floating point coefficients,
+ *                     corresponding to output of bilinear transofrm
+ *                     (z domain).
+ *
+ * Note: frequencies are in Hz.
+ * ----------------------------------------------------------
+ */
+static void szxform(
+    double *a0, double *a1, double *a2, /* numerator coefficients */
+    double *b0, double *b1, double *b2, /* denominator coefficients */
+    double fc,                          /* Filter cutoff frequency */
+    double fs,                          /* sampling rate */
+    double *k,                          /* overall gain factor */
+    float *coef)                        /* pointer to 4 iir coefficients */
+{
+    /* Calculate a1 and a2 and overwrite the original values */
+    prewarp(a0, a1, a2, fc, fs);
+    prewarp(b0, b1, b2, fc, fs);
+    bilinear(*a0, *a1, *a2, *b0, *b1, *b2, k, fs, coef);
+}
+
+
+/* ----------------- file bilinear.c end ----------------- */
+
+/* ----------------- file filter.txt begin -----------------
+How to construct a kewl low pass resonant filter?
+
+Lets assume we want to create a filter for analog synth.
+The filter rolloff is 24 db/oct, which corresponds to 4th
+order filter. Filter of first order is equivalent to RC circuit
+and has max rolloff of 6 db/oct.
+
+We will use classical Butterworth IIR filter design, as it
+exactly corresponds to our requirements.
+
+A common practice is to chain several 2nd order sections,
+or biquads, as they commonly called, in order to achive a higher
+order filter. Each 2nd order section is a 2nd order filter, which
+has 12 db/oct roloff. So, we need 2 of those sections in series.
+
+To compute those sections, we use standard Butterworth polinomials,
+or so called s-domain representation and convert it into z-domain,
+or digital domain. The reason we need to do this is because
+the filter theory exists for analog filters for a long time
+and there exist no theory of working in digital domain directly.
+So the common practice is to take standard analog filter design
+and use so called bilinear transform to convert the butterworth
+equasion coefficients into z-domain.
+
+Once we compute the z-domain coefficients, we can use them in
+a very simple transfer function, such as iir_filter() in our
+C source code, in order to perform the filtering function.
+The filter itself is the simpliest thing in the world.
+The most complicated thing is computing the coefficients
+for z-domain.
+
+Ok, lets look at butterworth polynomials, arranged as a series
+of 2nd order sections:
+
+ * Note: n is filter order.
+ *
+ *      n       Polynomials
+ *      --------------------------------------------------------------------
+ *      2       s^2 + 1.4142s +1
+ *      4       (s^2 + 0.765367s + 1) * (s^2 + 1.847759s + 1)
+ *      6       (s^2 + 0.5176387s + 1) * (s^2 + 1.414214 + 1) * (s^2 + 1.931852s + 1)
+ *
+ * For n=4 we have following equasion for the filter transfer function:
+ *
+ *                     1                              1
+ * T(s) = --------------------------- * ----------------------------
+ *        s^2 + (1/Q) * 0.765367s + 1   s^2 + (1/Q) * 1.847759s + 1
+ *
+
+The filter consists of two 2nd order secions since highest s power is 2.
+Now we can take the coefficients, or the numbers by which s is multiplied
+and plug them into a standard formula to be used by bilinear transform.
+
+Our standard form for each 2nd order secion is:
+
+       a2 * s^2 + a1 * s + a0
+H(s) = ----------------------
+       b2 * s^2 + b1 * s + b0
+
+Note that butterworth nominator is 1 for all filter sections,
+which means s^2 = 0 and s^1 = 0
+
+Lets convert standard butterworth polinomials into this form:
+
+       0 + 0 + 1                  0 + 0 + 1
+-------------------------- * --------------------------
+1 + ((1/Q) * 0.765367) + 1   1 + ((1/Q) * 1.847759) + 1
+
+Section 1:
+a2 = 0; a1 = 0; a0 = 1;
+b2 = 1; b1 = 0.5176387; b0 = 1;
+
+Section 2:
+a2 = 0; a1 = 0; a0 = 1;
+b2 = 1; b1 = 1.847759; b0 = 1;
+
+That Q is filter quality factor or resonance, in the range of
+1 to 1000. The overall filter Q is a product of all 2nd order stages.
+For example, the 6th order filter (3 stages, or biquads)
+with individual Q of 2 will have filter Q = 2 * 2 * 2 = 8.
+
+These a and b coefficients are used directly by the szxform()
+and bilinear() functions.
+
+The transfer function for z-domain is:
+
+       1 + alpha1 * z^(-1) + alpha2 * z^(-2)
+H(z) = -------------------------------------
+       1 + beta1 * z^(-1) + beta2 * z^(-2)
+
+When you need to change the filter frequency cutoff or resonance,
+or Q, you call the szxform() function with proper a and b
+coefficients and the new filter cutoff frequency or resonance.
+You also need to supply the sampling rate and filter gain you want
+to achive. For our purposes the gain = 1.
+
+We call szxform() function 2 times becase we have 2 filter sections.
+Each call provides different coefficients.
+
+The gain argument to szxform() is a pointer to desired filter
+gain variable.
+
+double k = 1.0;      // overall gain factor
+
+Upon return from each call, the k argument will be set to a value,
+by which to multiply our actual signal in order for the gain
+to be one. On second call to szxform() we provide k that was
+changed by the previous section. During actual audio filtering
+function iir_filter() will use this k
+
+Summary:
+
+Our filter is pretty close to ideal in terms of all relevant
+parameters and filter stability even with extremely large values
+of resonance. This filter design has been verified under all
+variations of parameters and it all appears to work as advertized.
+
+Good luck with it.
+If you ever make a directX wrapper for it, post it to comp.dsp.
+
+
+ *
+ * ----------------------------------------------------------
+ *References:
+ *Van Valkenburg, "Analog Filter Design"
+ *Oxford University Press 1982
+ *ISBN 0-19-510734-9
+ *
+ *C Language Algorithms for Digital Signal Processing
+ *Paul Embree, Bruce Kimble
+ *Prentice Hall, 1991
+ *ISBN 0-13-133406-9
+ *
+ *Digital Filter Designer's Handbook
+ *With C++ Algorithms
+ *Britton Rorabaugh
+ *McGraw Hill, 1997
+ *ISBN 0-07-053806-9
+ * ----------------------------------------------------------
+
+
+
+// ----------------- file filter.txt end ----------------- */