Use a new low-pass filter, based on the I3DL2 spec

Many thanks to Christopher Fitzgerald, for helping with it

2 files changed, 40 insertions, 479 deletions

diff --git a/Alc/ALu.c b/Alc/ALu.c
index 9baf6cee..3b8be4a2 100644
--- a/Alc/ALu.c
+++ b/Alc/ALu.c
@@ -157,31 +157,18 @@ __inline ALuint aluChannelsFromFormat(ALenum format)
 
 static __inline ALfloat lpFilter(FILTER *iir, ALfloat input)
 {
-    float *hist1_ptr,*hist2_ptr,*coef_ptr;
-    ALfloat 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++);
-
-    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++;
+    ALfloat *history = iir->history;
+    ALfloat a = iir->coeff;
+    ALfloat output = input;
+
+    output = output + (history[0]-output)*a;
+    history[0] = output;
+    output = output + (history[1]-output)*a;
+    history[1] = output;
+    output = output + (history[2]-output)*a;
+    history[2] = output;
+    output = output + (history[3]-output)*a;
+    history[3] = output;
 
     return output;
 }
@@ -234,10 +221,6 @@ static __inline ALvoid aluMatrixVector(ALfloat *vector,ALfloat matrix[3][3])
     memcpy(vector, result, sizeof(result));
 }
 
-static __inline ALfloat aluComputeSample(ALfloat GainHF, ALfloat sample, ALfloat LowSample)
-{
-    return LowSample + ((sample - LowSample) * GainHF);
-}
 
 static ALvoid CalcSourceParams(ALCcontext *ALContext, ALsource *ALSource,
                                ALenum isMono, ALenum OutputFormat,
@@ -258,6 +241,7 @@ static ALvoid CalcSourceParams(ALCcontext *ALContext, ALsource *ALSource,
     ALfloat RoomRolloff;
     ALfloat DryGainHF = 1.0f;
     ALfloat WetGainHF = 1.0f;
+    ALfloat cw, a, g;
 
     //Get context properties
     DopplerFactor   = ALContext->DopplerFactor * ALSource->DopplerFactor;
@@ -604,6 +588,25 @@ static ALvoid CalcSourceParams(ALCcontext *ALContext, ALsource *ALSource,
                 break;
         }
 
+        // Update filter coefficients. Calculations based on the I3DL2 spec.
+        cw = cos(2.0f*3.141592654f * LOWPASSFREQCUTOFF / ALContext->Frequency);
+        // We use four chained one-pole filters, so we need to take the fourth
+        // root of the squared gain, which is the same as the square root of
+        // the base gain.
+        // Be careful with gains < 0.0001, as that causes the coefficient to
+        // head towards 1, which will flatten the signal
+        g = aluSqrt(__max(DryGainHF, 0.0001f));
+        a = 0.0f;
+        if(g < 0.9999f) // 1-epsilon
+            a = (1 - g*cw - aluSqrt(2*g*(1-cw) - g*g*(1 - cw*cw))) / (1 - g);
+        ALSource->iirFilter.coeff = a;
+
+        g = aluSqrt(__max(WetGainHF, 0.0001f));
+        a = 0.0f;
+        if(g < 0.9999f) // 1-epsilon
+            a = (1 - g*cw - aluSqrt(2*g*(1-cw) - g*g*(1 - cw*cw))) / (1 - g);
+        ALSource->Send[0].iirFilter.coeff = a;
+
         *drygainhf = DryGainHF;
         *wetgainhf = WetGainHF;
     }
@@ -648,8 +651,6 @@ ALvoid aluMixData(ALCcontext *ALContext,ALvoid *buffer,ALsizei size,ALenum forma
     ALuint rampLength;
     ALfloat dryGainStep[OUTPUTCHANNELS];
     ALfloat wetGainStep[OUTPUTCHANNELS];
-    ALfloat dryGainHFStep;
-    ALfloat wetGainHFStep;
     ALuint BlockAlign,BufferSize;
     ALuint DataSize=0,DataPosInt=0,DataPosFrac=0;
     ALuint Channels,Frequency,ulExtraSamples;
@@ -667,7 +668,7 @@ ALvoid aluMixData(ALCcontext *ALContext,ALvoid *buffer,ALsizei size,ALenum forma
     ALbufferlistitem *BufferListItem;
     ALuint loop;
     ALint64 DataSize64,DataPos64;
-    FILTER *Filter;
+    FILTER *DryFilter, *WetFilter;
     int fpuState;
 
     SuspendContext(ALContext);
@@ -742,7 +743,8 @@ ALvoid aluMixData(ALCcontext *ALContext,ALvoid *buffer,ALsizei size,ALenum forma
                     //Get source info
                     DataPosInt = ALSource->position;
                     DataPosFrac = ALSource->position_fraction;
-                    Filter = &ALSource->iirFilter;
+                    DryFilter = &ALSource->iirFilter;
+                    WetFilter = &ALSource->Send[0].iirFilter;
                     DrySend = ALSource->DryGains;
                     WetSend = ALSource->WetGains;
 
@@ -752,10 +754,6 @@ ALvoid aluMixData(ALCcontext *ALContext,ALvoid *buffer,ALsizei size,ALenum forma
                         dryGainStep[i] = (newDrySend[i]-DrySend[i]) / rampLength;
                         wetGainStep[i] = (newWetSend[i]-WetSend[i]) / rampLength;
                     }
-                    dryGainHFStep = (DryGainHF-ALSource->DryGainHF) / rampLength;
-                    wetGainHFStep = (WetGainHF-ALSource->WetGainHF) / rampLength;
-                    DryGainHF = ALSource->DryGainHF;
-                    WetGainHF = ALSource->WetGainHF;
 
                     //Compute 18.14 fixed point step
                     increment = (ALint)(Pitch*(ALfloat)(1L<<FRACTIONBITS));
@@ -815,19 +813,16 @@ ALvoid aluMixData(ALCcontext *ALContext,ALvoid *buffer,ALsizei size,ALenum forma
                             DrySend[i] += dryGainStep[i];
                             WetSend[i] += wetGainStep[i];
                         }
-                        DryGainHF += dryGainHFStep;
-                        WetGainHF += wetGainHFStep;
 
                         if(Channels==1)
                         {
-                            ALfloat sample, lowsamp, outsamp;
+                            ALfloat sample, outsamp;
                             //First order interpolator
                             sample = (Data[k]*((1<<FRACTIONBITS)-fraction) +
                                       Data[k+1]*fraction) >> FRACTIONBITS;
-                            lowsamp = lpFilter(Filter, sample);
 
                             //Direct path final mix buffer and panning
-                            outsamp = aluComputeSample(DryGainHF, sample, lowsamp);
+                            outsamp = lpFilter(DryFilter, sample);
                             DryBuffer[j][FRONT_LEFT]  += outsamp*DrySend[FRONT_LEFT];
                             DryBuffer[j][FRONT_RIGHT] += outsamp*DrySend[FRONT_RIGHT];
                             DryBuffer[j][SIDE_LEFT]   += outsamp*DrySend[SIDE_LEFT];
@@ -835,7 +830,7 @@ ALvoid aluMixData(ALCcontext *ALContext,ALvoid *buffer,ALsizei size,ALenum forma
                             DryBuffer[j][BACK_LEFT]   += outsamp*DrySend[BACK_LEFT];
                             DryBuffer[j][BACK_RIGHT]  += outsamp*DrySend[BACK_RIGHT];
                             //Room path final mix buffer and panning
-                            outsamp = aluComputeSample(WetGainHF, sample, lowsamp);
+                            outsamp = lpFilter(WetFilter, sample);
                             WetBuffer[j][FRONT_LEFT]  += outsamp*WetSend[FRONT_LEFT];
                             WetBuffer[j][FRONT_RIGHT] += outsamp*WetSend[FRONT_RIGHT];
                             WetBuffer[j][SIDE_LEFT]   += outsamp*WetSend[SIDE_LEFT];
@@ -923,8 +918,6 @@ ALvoid aluMixData(ALCcontext *ALContext,ALvoid *buffer,ALsizei size,ALenum forma
                     //Update source info
                     ALSource->position = DataPosInt;
                     ALSource->position_fraction = DataPosFrac;
-                    ALSource->DryGainHF = DryGainHF;
-                    ALSource->WetGainHF = WetGainHF;
                 }
 
                 //Handle looping sources
@@ -1015,7 +1008,7 @@ ALvoid aluMixData(ALCcontext *ALContext,ALvoid *buffer,ALsizei size,ALenum forma
                 ALfloat LateReverbGain = ALEffectSlot->effect.Reverb.LateReverbGain;
                 ALfloat sample, lowsample;
 
-                Filter = &ALEffectSlot->iirFilter;
+                WetFilter = &ALEffectSlot->iirFilter;
                 for(i = 0;i < SamplesToDo;i++)
                 {
                     sample  = WetBuffer[i][FRONT_LEFT] +WetBuffer[i][SIDE_LEFT] +WetBuffer[i][BACK_LEFT];
@@ -1027,7 +1020,7 @@ ALvoid aluMixData(ALCcontext *ALContext,ALvoid *buffer,ALsizei size,ALenum forma
                     DelayBuffer[LatePos] *= LateReverbGain;
 
                     Pos = (Pos+1) % Length;
-                    lowsample = lpFilter(Filter, DelayBuffer[Pos]);
+                    lowsample = lpFilter(WetFilter, DelayBuffer[Pos]);
                     lowsample += (DelayBuffer[Pos]-lowsample) * DecayHFRatio;
 
                     DelayBuffer[LatePos] += lowsample * DecayGain;
diff --git a/Alc/lpfilter.c b/Alc/lpfilter.c
deleted file mode 100644
index 98019e2e..00000000
--- a/Alc/lpfilter.c
+++ /dev/null
@@ -1,432 +0,0 @@
-/* ----------------- 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"
-
-
-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.
- * --------------------------------------------------------------------
- */
-
-/*** moved to ALu.c ***/
-
-/*
- * --------------------------------------------------------------------
- *
- * InitLowPassFilter()
- *
- * Initialize filter coefficients.
- * We create a 2nd order filter (12 db/oct rolloff), consisting
- * of one second order section.
- * --------------------------------------------------------------------
- */
-int InitLowPassFilter(ALCcontext *Context, FILTER *iir)
-{
-    float    *coef;
-    double   fs, fc; /* Sampling frequency, cutoff frequency */
-    double   Q;      /* Resonance > 1.0 < 1000 */
-    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 prototype coefficients,
-                                          1 for each filter section */
-
-
-    /*
-     * Setup filter s-domain coefficients
-     */
-    /* Section 1 */
-    ProtoCoef.a0 = 1.0;
-    ProtoCoef.a1 = 0;
-    ProtoCoef.a2 = 0;
-    ProtoCoef.b0 = 1.0;
-    ProtoCoef.b1 = 1.4142;
-    ProtoCoef.b2 = 1.0;
-
-    /* Clear the coefficient and history arrays */
-    memset(iir->coef, 0, sizeof(iir->coef));
-    memset(iir->history, 0, sizeof(iir->history));
-
-    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
-     */
-    a0 = ProtoCoef.a0;
-    a1 = ProtoCoef.a1;
-    a2 = ProtoCoef.a2;
-
-    b0 = ProtoCoef.b0;
-    b1 = ProtoCoef.b1 / Q;      /* Divide by resonance or Q */
-    b2 = ProtoCoef.b2;
-    szxform(&a0, &a1, &a2, &b0, &b1, &b2, fc, fs, &k, coef);
-
-    /* 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 12 db/oct filter, we take a 2nd order
- *      section and compute the coefficients.
- *
- *      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=2, or one second order section, we have following equasion:
- *
- *      (1 / (s^2 + (1/Q) * 1.4142s + 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 = 1.4142; 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)
-{
-    static const double pi = 3.141592653589793238;
-
-    double 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 12 db/oct, which corresponds to 2nd
-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.
-
-Each 2nd order section is a 2nd order filter, which has 12 db/oct
-rolloff. So, we only need one section
-
-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=2 we have following equasion for the filter transfer function:
- *
- *                     1
- * T(s) = -------------------------
- *        s^2 + (1/Q) * 1.4142s + 1
- *
-
-The filter consists of one 2nd order section 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
-------------------------
-1 + ((1/Q) * 1.4142) + 1
-
-Section 1:
-a2 = 0; a1 = 0; a0 = 1;
-b2 = 1; b1 = 1.4142; 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 1 time becase we have 1 filter section.
-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 following calls 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 ----------------- */