How to implement IIR filter in C?

Question

I am trying to implement an IIR filter.

I tried to implement the following filter but

FIR : y(n) = b0(x[n]) + ... +bM-1(x[n-M+1])

IIR : y(n) = {b0(x[n]) + ... +bM-1(x[n-M+1])} - {a1(y[n-1]) + ... +aN(y[n-N}

I am confused on how to implement y[n-1]......

Here is my code.

void IIRFloat(double *coeffs_B, double *coeffs_A, double *input, double *output, int length, int filterLength)
{
    double bcc, acc;
    double *coeffa, *coeffb;
    double *inputp;
    double *outputp;
    int n,k;

    //filter length =7
    for (n = 0; n < length; n++) {
        coeffa = coeffs_A;
        coeffb = coeffs_B;
        inputp = &insamp[filterLength - 1 + n]; //insamp[6]~insamp[85]

        acc = 0;
        bcc = 0;

        for (k = 0; k < filterLength; k++) {
            bcc += (*coeffb++) * (*inputp--); //b[0] * x[6] + b[1] * x[5]...+b[6] * x[0]
        }

        for (k = 1; k < filterLength; k++) {
            acc += (*coeffa++) * (*output--); //a[1] * y[5] + a[2] * y[4]...+a[6] * y[0]
        }
        output[n] = bcc-acc;

    }
}

I will not copy here the code for memove function for seek of brevity.

Answer 1

An IIR filter is recursive which means it takes the past values of the output (equivalently it could be expressed as infinite sequence of the input).

Suppose you have the filter

y[n] = b*x[n]-a*y[n-1]

Generally the first output is initialized to a given value, for example 0.

Suppose you want to filter the signal x of length N then you should do something like:

double out[N]; //output vector
double a = 0.3, b=0.5; //assign the coefficients a value
out[0] = 0; //initialize the first element

for (int i=1; i<N; i++)
{
    out[i] = b*x[i] -a*[i-1];  
}

In the case of your code I cannot know what you are doing in the line inputp = &insamp[filterLength - 1 + n]; and that might be a problem.

I will assume inputp is the signal you want to filter.

Another issue: you used filter filterLength to indicate the lengths of both the input elements and the output elements: that is not generally in a IIR filter.

The output elements from 0 to filterlength should be initialized somehow, suppose to 0. Then in the second loop you made the loop index starting from 1 but the coefficient array should start from 0.

(DOUBT: If the oldest element is y[5] how can the length be 7?)

Using index instead of dereferencing the array your code should be something like this:

void IIRFloat(double *coeffs_B, double *coeffs_A, double *input, double *output, int length, int filterLength)
{
    double bcc, acc;
    double *inputp;
    int n,k;


    for (int ii=0; ii<filterLength; ii++)
    {
        output[ii] = 0;
    }

    //filter length =7
    for (n = 0; n < length; n++) {
        inputp = &insamp[filterLength - 1 + n]; //insamp[6]~insamp[85]

        acc = 0;
        bcc = 0;

        for (k = 0; k < filterLength; k++) 
        {
            bcc += coeffb[k] * inputp[filterLength-k-1]; //b[0] * x[6] + b[1] * x[5]...+b[6] * x[0]
        }

        for (k = 0; k < filterLength; k++) 
        {
            acc += coeffa[k] * output[filterLength-k-1]; //a[1] * y[5] + a[2] * y[4]...+a[6] * y[0]
        }
        output[n] = bcc-acc;

    }
}

EDIT I think the two for being the same length can be merged together:

for (k = 0; k < filterLength; k++) 
{
    output[n] += (coeffb[k] * inputp[filterLength-k-1] - coeffa[k] * output[filterLength-k-1]); 
}