1
$\begingroup$

I can't find anything on this topic, so either I'm in the wrong direction or else there is just nothing about it on the internet.

So let's say I have 3 $b_i$ coefficients ($b_0,b_1,$ and $b_2$) and 2 $a$ coefficients ($a_1$ and $a_2$) in a difference equation. How do I code the difference equation in C, so that I get the desired effect from the filter coefficients?

$\endgroup$
  • $\begingroup$ i presume your "b coefficients" are the feed-forward coefficients and your "a coefficients" are the feedback coefficients? dnopas, would you mind taking a quick second to learn a little about $\LaTeX$ so you can write nice looking equations to illustrate what you're asking? i'll post you some very clean and simple C code, if you can do that. $\endgroup$ – robert bristow-johnson Jan 6 '17 at 21:50
  • $\begingroup$ I can't believe you can't find anything on the internet. What have you searched for? Googling "IIR C implementation" yields copious amounts of examples. $\endgroup$ – Marcus Müller Jan 6 '17 at 21:50
  • $\begingroup$ @robertbristow-johnson imgur.com/a/FGLd7 This is basically the difference equation. My coefficients: a1=-1.42984, a2=0.4128056 b0=0.067455, b1=0.134911, b2=0.0674555 Holy shit sorry for the formatting $\endgroup$ – dnopas Jan 6 '17 at 21:57
  • $\begingroup$ okay, that is not factored into biquads, which i think is a good thing to do. $\endgroup$ – robert bristow-johnson Jan 6 '17 at 21:59
  • $\begingroup$ @robertbristow-johnson So basically the b constants are related to the input and previous inputs, and a for past outputs. So I should factor it into biquads, so into the Z domain? $\endgroup$ – dnopas Jan 6 '17 at 22:02
3
$\begingroup$
/**********************************************************************
*                                                                     *
*                                                                     *
*                                                                     *
*                        Copyright (c) 2005                           *
*                                                                     *
*                      Robert Bristow-Johnson                         *
*                                                                     *
*                     rbj@audioimagination.com                        *
*                                                                     *
*                                                                     *
*                                 EQ.c                                *
*                                                                     *
*                                                                     *
*   Utilities for computing EQ filter coefficients, plotting          *
*   frequency response, and EQ filtering sampled audio data.          *
*                                                                     *
*                                                                     *
**********************************************************************/






#include <stdio.h>
#include "__functions.h"

enum EQ_type
{
LPF,
HPF,
BPF,
BPF0,
notch,
APF,
peakingEQ,
lowShelf,
highShelf,
lowShelf1,
highShelf1
};


typedef struct
{
float Fs;
enum EQ_type type;
float f0;
float dBgain;
float Q;
float filterCoefficients[5];    /* b0, b1, b2, a1, a2   (a0 is normallized to 1) */
float graphCoefficients[6];     /* B2, B1, B0, A2, A1, A0 */
} sectionData;



int main(void);
int computeSectionCofficients(float Fs, enum EQ_type type, float f0, float dBgain, float Q, sectionData *data);
int clearArray(long N, float *array);
int plotSectionFrequencyResponse(int Npoints, float Fmin, float Fmax, float *frequencyValues, float *dBvalues, sectionData *data);
int filterBuffer(long Nsamples, int Nsections, float *input, float *output, float *filterStates, sectionData **sections);



#define BUFFER_SIZE 1024
#define GRAPH_SIZE 589
#define NUM_SECTIONS 4
#define FS 44100.0
#define FMIN 19.99955848414893
#define FMAX 19999.55848414893

#define F1 50.0
#define G1 0.0
#define Q1 1.0

#define F2 200.0
#define G2 0.0
#define Q2 2.0

#define F3 1000.0
#define G3 10.0
#define Q3 2.0

#define F4 5000.0
#define G4 0.0
#define Q4 1.0

sectionData section1, section2, section3, section4;         /* this has to be global */

void plotTheFrequencyResponse(void);
void runTheFilter(void);

int main(void)
{
computeSectionCofficients(FS, lowShelf1,  F1, G1, Q1, &section1);
computeSectionCofficients(FS, peakingEQ,  F2, G2, Q2, &section2);
computeSectionCofficients(FS, peakingEQ,  F3, G3, Q3, &section3);
computeSectionCofficients(FS, highShelf1, F4, G4, Q4, &section4);

plotTheFrequencyResponse();
runTheFilter();
}


void plotTheFrequencyResponse(void)
{
float freq[GRAPH_SIZE];
float graph[GRAPH_SIZE];

clearArray(GRAPH_SIZE, graph);
plotSectionFrequencyResponse(GRAPH_SIZE, FMIN, FMAX, freq, graph, &section1);
plotSectionFrequencyResponse(GRAPH_SIZE, FMIN, FMAX, freq, graph, &section2);
plotSectionFrequencyResponse(GRAPH_SIZE, FMIN, FMAX, freq, graph, &section3);
plotSectionFrequencyResponse(GRAPH_SIZE, FMIN, FMAX, freq, graph, &section4);
}



void runTheFilter(void)
{
float aFilterBuffer[BUFFER_SIZE];
float aFilterStates[(2*NUM_SECTIONS+2)];
sectionData *sectionPtrArray[NUM_SECTIONS];

sectionPtrArray[0] = &section1;
sectionPtrArray[1] = &section2;
sectionPtrArray[2] = &section3;
sectionPtrArray[3] = &section4;

clearArray(2*NUM_SECTIONS+2, aFilterStates);
clearArray(BUFFER_SIZE, aFilterBuffer);
aFilterBuffer[0] = 1.0;                         /* input an impulse to get impulse response */
filterBuffer(BUFFER_SIZE, NUM_SECTIONS, aFilterBuffer, aFilterBuffer, aFilterStates, sectionPtrArray);
}





int computeSectionCofficients(float Fs, enum EQ_type type, float f0, float dBgain, float Q, sectionData *data)
{
register float temp, *coef_ptr;
float A = __exp2(0.08304820237218*dBgain);
float A2 = A*A;
float Ap1, Am1;
float cos_w0 = __cos_pi(2.0*f0/Fs);
float sin_w0 = __sin_pi(2.0*f0/Fs);
float alpha = 0.5*sin_w0/Q;
float a0, a1, a2, b0, b1, b2;

switch (type)
    {
    case LPF:
        b0 =  0.5*(1.0 - cos_w0)*A2;
        b1 =  2.0*b0;
        b2 =      b0;
        a0 =   1.0 + alpha;
        a1 =  -2.0*cos_w0;
        a2 =   1.0 - alpha;
        break;

    case HPF:
        b0 =  0.5*(1.0 + cos_w0)*A2;
        b1 = -2.0*b0;
        b2 =      b0;
        a0 =   1.0 + alpha;
        a1 =  -2.0*cos_w0;
        a2 =   1.0 - alpha;
        break;

    case BPF:
        b0 =   0.5*sin_w0*A2;
        b1 =   0.0;
        b2 =  -b0;
        a0 =   1.0 + alpha;
        a1 =  -2.0*cos_w0;
        a2 =   1.0 - alpha;
        break;

    case BPF0:
        b0 =   alpha*A2;
        b1 =   0.0;
        b2 =  -b0;
        a0 =   1.0 + alpha;
        a1 =  -2.0*cos_w0;
        a2 =   1.0 - alpha;
        break;

    case notch:
        b0 =   A2;
        b1 =  -2.0*cos_w0*A2;
        b2 =   A2;
        a0 =   1.0 + alpha;
        a1 =  -2.0*cos_w0;
        a2 =   1.0 - alpha;
        break;

    case APF:
        b0 =  (1.0 - alpha)*A2;
        b1 =  -2.0*cos_w0*A2;
        b2 =  (1.0 + alpha)*A2;
        a0 =   1.0 + alpha;
        a1 =  -2.0*cos_w0;
        a2 =   1.0 - alpha;
        break;

    case peakingEQ:
        if (A > 1.0)
            {
            alpha *= A;
            }
         else
            {
            alpha /= A;
            }
        b0 =   1.0 + alpha*A;
        b1 =  -2.0*cos_w0;
        b2 =   1.0 - alpha*A;
        a0 =   1.0 + alpha/A;
        a1 =  -2.0*cos_w0;
        a2 =   1.0 - alpha/A;
        break;

    case lowShelf:
        temp = 2.0*__sqrt(A)*alpha;
        Ap1  = A + 1.0;
        Am1  = A - 1.0;
        b0 =  A*(Ap1 - Am1*cos_w0 + temp);
        b1 =  2.0*A*(Am1 - Ap1*cos_w0);
        b2 =  A*(Ap1 - Am1*cos_w0 - temp);
        a0 =  Ap1 + Am1*cos_w0 + temp;
        a1 =  -2.0*( Am1 + Ap1*cos_w0);
        a2 =  Ap1 + Am1*cos_w0 - temp;
        break;

    case highShelf:
        temp = 2.0*__sqrt(A)*alpha;
        Ap1  = A + 1.0;
        Am1  = A - 1.0;
        b0 =  A*(Ap1 + Am1*cos_w0 + temp);
        b1 = -2.0*A*(Am1 + Ap1*cos_w0);
        b2 =  A*(Ap1 + Am1*cos_w0 - temp);
        a0 =  Ap1 - Am1*cos_w0 + temp;
        a1 =  2.0*( Am1 - Ap1*cos_w0);
        a2 =  Ap1 - Am1*cos_w0 - temp;
        break;

    case lowShelf1:
        temp = cos_w0 + 1.0;
        b0 =  A2*sin_w0 +    temp;
        b1 =  A2*sin_w0 -    temp;
        b2 =  0.0;
        a0 =     sin_w0 +    temp;
        a1 =     sin_w0 -    temp;
        a2 =  0.0;
        break;

    case highShelf1:
        temp = cos_w0 + 1.0;
        b0 =     sin_w0 + A2*temp;
        b1 =     sin_w0 - A2*temp;
        b2 =  0.0;
        a0 =     sin_w0 +    temp;
        a1 =     sin_w0 -    temp;
        a2 =  0.0;
        break;
    }

coef_ptr = &(data->filterCoefficients[0]);
temp = 1/a0;
*coef_ptr++ = temp*b0;                      /* b0/a0 */
*coef_ptr++ = temp*b1;                      /* b1/a0 */
*coef_ptr++ = temp*b2;                      /* b2/a0 */
*coef_ptr++ = temp*a1;                      /* a1/a0 */
*coef_ptr   = temp*a2;                      /* a2/a0 */

coef_ptr = &(data->graphCoefficients[0]);
temp = b0*b2;
*coef_ptr++ = temp;                         /* B2 */
*coef_ptr++ = -0.25*b1*(b0 + b2) - temp;    /* B1 */
temp = b0 + b1 + b2;
temp = 0.0625*temp*temp;
*coef_ptr++ = temp;                         /* B0 */
temp = a0*a2;
*coef_ptr++ = temp;                         /* A2 */
*coef_ptr++ = -0.25*a1*(a0 + a2) - temp;    /* A1 */
temp = a0 + a1 + a2;
temp = 0.0625*temp*temp;
*coef_ptr   = temp;                         /* A0 */

data->Fs = Fs;
data->type = type;
data->f0 = f0;
data->dBgain = dBgain;
data->Q = Q;

return 0;
}


int clearArray(long N, float *array)
{
register long n;

for (n=0; n<N; n++)
    {
    array[n] = 0.0;
    }

return 0;
}


int plotSectionFrequencyResponse(int Npoints, float Fmin, float Fmax, float *frequencyValues, float *dBvalues, sectionData *data)
{
register int n;
float logFreqStep = __log2(Fmax/Fmin)/(float)(Npoints-1);
register float phi, tau = 1.0/(data->Fs);
float *coef_ptr = &(data->graphCoefficients[0]);
register float B2, B1, B0, A2, A1, A0;

B2 = *coef_ptr++;
B1 = *coef_ptr++;
B0 = *coef_ptr++;
A2 = *coef_ptr++;
A1 = *coef_ptr++;
A0 = *coef_ptr;

for (n=0; n<Npoints; n++)
    {
    phi = Fmin*__exp2((float)n*logFreqStep);
    frequencyValues[n] = phi;
    phi = __sin_pi(tau*phi);
    phi = phi*phi;                                                                              /* sin(w/2)^2 */
    dBvalues[n] += DB_LOG2_ENERGY*(__log2((B2*phi + B1)*phi + B0) - __log2((A2*phi + A1)*phi + A0));
    }

return 0;
}


int filterBuffer(long Nsamples, int Nsections, float *input, float *output, float *filterStates, sectionData **sections)
{
register long n;
register int i;
register sectionData **section_ptr;
register float sampleValue, state1, state2, *state_ptr, *coef_ptr;

for (n=0; n<Nsamples; n++)
    {
    section_ptr = sections;                         /* reset section pointer */
    state_ptr = filterStates;                       /* reset state pointer */

    sampleValue = input[n];                         /* get input sample */

        state2 = *state_ptr++;                      /*   x[n-2]            */
        state1 = *state_ptr--;                      /*   x[n-1]            */
        *state_ptr++ = state1;                      /*   x[n-1] -> x[n-2]  */
        *state_ptr++ = sampleValue;                 /*   x[n]   -> x[n-1]  */

    i = Nsections;
    while (--i >= 0)
        {
        coef_ptr = &((*section_ptr++)->filterCoefficients[0]);  /* point to section filter coefficients */

        sampleValue =  *coef_ptr++ * sampleValue;   /*   b0*x[n]           */
        sampleValue += *coef_ptr++ * state1;        /*   b1*x[n-1]         */
        sampleValue += *coef_ptr++ * state2;        /*   b2*x[n-2]         */

        state2 = *state_ptr++;                      /*   y[n-2]            */
        state1 = *state_ptr--;                      /*   y[n-1]            */
        *state_ptr++ = state1;                      /*   y[n-1] -> y[n-2]  */

        sampleValue -= *coef_ptr++ * state1;        /*  -a1*y[n-1]         */
        sampleValue -= *coef_ptr++ * state2;        /*  -a2*y[n-2]         */

        *state_ptr++ = sampleValue;                 /*   y[n]   -> y[n-1]  */
/*
*      at this point sampleValue is the output of the section just completed 
*     and will be the input to the section about to be computed unless this 
*     is the last section, in which case sampleValue is the input.
*     state1 and state2 are the output states of the section just completed
*     and will be the input states of of the section about to be computed.
*/
        }

    output[n] = sampleValue;                        /* store output sample  */
    }

return 0;
}
$\endgroup$
  • 1
    $\begingroup$ Robert, just a suggestion: for such a non-trivial piece of code, you may want to consider attaching some kind of license to it, to give users some clarity about what they can do with it. $\endgroup$ – MBaz Jan 6 '17 at 22:17
  • $\begingroup$ Thanks for helping. I was just looking for some simple code to implement this low order iir filter with difference equations. I'll take a look at this to see if I can find the part I need to understand! $\endgroup$ – dnopas Jan 6 '17 at 22:18
  • $\begingroup$ Btw isn't it redundant to divide with a0 since its 1 anyway? $\endgroup$ – dnopas Jan 6 '17 at 22:24
  • $\begingroup$ @MBaz, i am not expecting to sell this code to anyone in the future. it's 12 years old (many parts are much older). i am not really concerned about what users do with it. $\endgroup$ – robert bristow-johnson Jan 7 '17 at 1:01
1
$\begingroup$

From you question, I'm assuming you have a transfer function of the form $H(z)$:

$$ H(z) = \frac{Y(z)}{X(z)} = \frac {b_0 + b_1 z^{-1} + b_2 z^{-2}} {1 + a_1 z^{-1} + a_2 z^{-2}} \tag{1} $$

Which has finite difference equation:

$$ y[n] = b_0 x[n] + b_1 x[n-1] + b_2 x[n-2] - a_1 y[n-1] - a_2 y[n-2] \tag{2} $$

You can boil this down to very simple C code by processing one sample at a time.

For convention let's assume input sample variables are x0 ($x[n]$), x1 ($x[n - 1]$), x2 ($x[n - 2]$).

Also lets assume the same for output sample variables y0 ($y[n]$), y1 ($y[n - 1]$), y2 ($y[n - 2]$).

Then you can write a simple structure like this:

struct difference_equation
{
   // coefficients
   double b0;
   double b1;
   double b2;
   double a1;
   double a2;
   // input history
   double x1;
   double x2;
   // output history
   double y1;
   double y2;

   double process(double x0)
   {
      double y0 = b0 * x0 + b1 * x1 + b2 * x2
                - a1 * y1 - a2 * y2;

      x2 = x1;
      x1 = x0;

      y2 = y1;
      y1 = y0;

      return y0;
   }
}

And you can use it like so:

difference_equation de;

//assign coefficients (I'm assuming you know how to calculate your own values):
de.b0 = 1;
de.b1 = 0;
de.b2 = 0;
de.a1 = 0;
de.a2 = 0;

// reset input/output history
de.x1 = 0;
de.x2 = 0;
de.y1 = 0;
de.y2 = 0;

double input_samples[3];
double output_samples[3];

for (int i = 0; i < 3; i++)
{
   output_samples[i] = de.process(input_samples[i]);
}

Now, this is direct form I. There are slightly more efficient ways to express this code using direct form II but the other forms aren't so obvious to express in code for someone learning.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.