• 11/18/99 05:32 PM
    Sign in to follow this  
    Followers 0

    Resonant Lowpass Filter Soruce

    General and Gameplay Programming

    Myopic Rhino
    // ----------------- file filterIIR00.c begin -----------------
    /*
    Resonant low pass filter source code.
    By baltrax@hotmail.com (Zxform)
    */

    #include
    #include
    #include

    /**************************************************************************

    FILTER.C - Source code for filter functions

    iir_filter IIR filter floats sample by sample (real time)

    *************************************************************************/

    /* FILTER INFORMATION STRUCTURE FOR FILTER ROUTINES */

    typedef struct {
    unsigned int length; /* size of filter */
    float *history; /* pointer to history in filter */
    float *coef; /* pointer to coefficients of filter */
    } FILTER;

    #define FILTER_SECTIONS 2 /* 2 filter sections for 24 db/oct filter */

    typedef struct {
    double a0, a1, a2; /* numerator coefficients */
    double b0, b1, b2; /* denominator coefficients */
    } BIQUAD;

    BIQUAD ProtoCoef[FILTER_SECTIONS]; /* Filter prototype coefficients,
    1 for each filter section */

    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 */

    /*
    * --------------------------------------------------------------------
    *
    * iir_filter - Perform IIR filtering sample by sample on floats
    *
    * Implements cascaded direct form II second order sections.
    * Requires FILTER structure for history and coefficients.
    * The length in the filter structure specifies the number of sections.
    * The size of the history array is 2*iir->length.
    * The size of the coefficient array is 4*iir->length + 1 because
    * the first coefficient is the overall scale factor for the filter.
    * Returns one output sample for each input sample. Allocates history
    * array if not previously allocated.
    *
    * float iir_filter(float input,FILTER *iir)
    *
    * float input new float input sample
    * FILTER *iir pointer to FILTER structure
    *
    * Returns float value giving the current output.
    *
    * Allocation errors cause an error message and a call to exit.
    * --------------------------------------------------------------------
    */
    float iir_filter(input,iir)
    float input; /* new input sample */
    FILTER *iir; /* pointer to FILTER structure */
    {
    unsigned int i;
    float *hist1_ptr,*hist2_ptr,*coef_ptr;
    float output,new_hist,history1,history2;

    /* allocate history array if different size than last call */

    if(!iir->history) {
    iir->history = (float *) calloc(2*iir->length,sizeof(float));
    if(!iir->history) {
    printf("\nUnable to allocate history array in iir_filter\n");
    exit(1);
    }
    }

    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 < iir->length; 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);
    }


    /*
    * --------------------------------------------------------------------
    *
    * main()
    *
    * Example main function to show how to update filter coefficients.
    * We create a 4th order filter (24 db/oct roloff), consisting
    * of two second order sections.
    * --------------------------------------------------------------------
    */
    int main()
    {
    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 */

    /*
    * 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;

    iir.length = FILTER_SECTIONS; /* Number of filter sections */

    /*
    * Allocate array of z-domain coefficients for each filter section
    * plus filter gain variable
    */
    iir.coef = (float *) calloc(4 * iir.length + 1, sizeof(float));
    if (!iir.coef)
    {
    printf("Unable to allocate coef array, exiting\n");
    exit(1);
    }

    k = 1.0; /* Set overall filter gain */
    coef = iir.coef + 1; /* Skip k, or gain */

    Q = 1; /* Resonance */
    fc = 5000; /* Filter cutoff (Hz) */
    fs = 44100; /* Sampling frequency (Hz) */

    /*
    * Compute z-domain coefficients for each biquad section
    * for new Cutoff Frequency and Resonance
    */
    for (nInd = 0; nInd < iir.length; 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;

    /* Display filter coefficients */
    for (nInd = 0; nInd < (iir.length * 4 + 1); nInd++)
    printf("C[%d] = %15.10f\n", nInd, iir.coef[nInd]);
    /*
    * To process audio samples, call function iir_filter()
    * for each audio sample
    */
    return (0);
    }




    // ----------------- file filterIIR00.c end -----------------


    Reposting bilinear.c just in case the other one was not the latest version.

    // ----------------- 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
    * ----------------------------------------------------------
    */

    #include

    void prewarp(double *a0, double *a1, double *a2, double fc, double fs);
    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.
    * ----------------------------------------------------------
    */
    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;
    }


    /*
    * ----------------------------------------------------------
    * 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
    * ----------------------------------------------------------
    */
    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.
    * ----------------------------------------------------------
    */
    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 -----------------
    0


    Sign in to follow this  
    Followers 0


    User Feedback

    Create an account or sign in to leave a review

    You need to be a member in order to leave a review

    Create an account

    Sign up for a new account in our community. It's easy!


    Register a new account

    Sign in

    Already have an account? Sign in here.


    Sign In Now

    There are no reviews to display.