IIR Filter Implementation in C results in unstable system response

Question

I need to implement a basic 6th order Butterworth lowpass filter on a DSP that sampled data at 50kHz. I calculate the coefficients using matlab for a filter with the following step response: $$\begin{align} y[n] {} = & \frac{1}{a_0}(b_0 x[n] + b_1 x[n-1] + \cdots + b_P x[n-P] \\ & {} - a_1 y[n-1] - a_2 y[n-2] - \cdots - a_Q y[n-Q]) \end{align}$$

The coefficients of this filter are given as:

a =

    1.0000   -5.9757   14.8789  -19.7584   14.7590   -5.8798    0.9760

>> b

b =

   1.0e-13 *

    0.0089    0.0533    0.1332    0.1776    0.1332    0.0533    0.0089

I tried implementing & simulating this filter in C but I keep getting an unstable response, I am not sure where the error is in my code:

#include <stdio.h>

#define N_LOWPASS 7

struct lowpass{
    float bs[N_LOWPASS];
    float as[N_LOWPASS];
    float xs[N_LOWPASS];
    float ys[N_LOWPASS];
    int xpointer;
    int ypointer;
};

struct lowpass lp={ {8.881784197001252e-16,5.329070518200751e-15,1.332267629550188e-14,1.776356839400251e-14,1.332267629550188e-14,5.329070518200751e-15,8.881784197001252e-16},
                    {1,-5.975723643994615,14.878912715306338,-19.758412157954773,14.758996633559045,-5.879789434468497,0.976015887552560},
                    {0},
                    {0},
                    0,
                    0
};
                    
float update_lowpass(struct lowpass* lp,float xn){
    //note : at the beginning of the function call, the pointers into the circular buffers point to the values that will be replaced
    //update x[n] position in circular buffer
    lp->xs[lp->xpointer]=xn;
    //compute y[n]
    unsigned int i;
    float sum=lp->bs[0]*xn;
    for(i=1; i<N_LOWPASS; i++){
        sum+=lp->bs[i]*lp->xs[(lp->xpointer-i+N_LOWPASS)%N_LOWPASS]-lp->as[i]*lp->ys[(lp->ypointer-i+N_LOWPASS)%N_LOWPASS];
    }
    //update y[n] in circular buffer
    float yn=sum/lp->as[0];
    lp->ys[lp->ypointer]=yn;

    lp->xpointer=(lp->xpointer+1)%N_LOWPASS;
    lp->ypointer=(lp->ypointer+1)%N_LOWPASS;
    return yn;
}

int main()
{
    //simulate the step response of the filter
    unsigned int N=50000;
    float yout[50000];
    unsigned int i=0;
    for(i=0; i<N; i++){
        float yn=update_lowpass(&lp,1.0);
        yout[i]=yn;
        printf("%f\n",yn);
    }
    return 0;
}

Update After changing the data type of the floats to double the step response did show stable behaviour but the two solutions from Matlab and C did differ somewhat: Also, in both cases the filter response did not seem to settle toward a final value of 1 and did show some whiggling:

Answer 1

Using 64-bit floats just papers over an underlying problem. As a polynomial increases in order, the positions of the zeros gets ever more sensitive to the values of the lowest-order terms. Moreover, in filtering, the sensitivity of the filter to quantization noise goes up similarly.

Moreover, the narrower the filter bandwidth or the sharper the filter's skirts, the more precision you need in both coefficients and data paths.

So while this measure worked for your particular 6th-order filter, at some point as you increase filter order or decrease bandwidth, it will stop doing the job for you.

It is standard practice -- and for good reasons -- to split IIR filters into 2nd-order filter sections (and possibly one 1st-order section) and cascade them. The biggest two advantages are that it minimizes the consequences of coefficient rounding and it minimizes the consequences of truncation in filtering.

There's even a standard form for these filters: it is the biquad filter. Search for that term and you'll get all sorts of useful information on how to choose which form to use and how to code them up.