# IIR Filter Implementation in C results in unstable system response

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:

• 1e-15 precision on float will never work. Also, see this for workarounds when the poles are ridiculously low compared to the sampling frequency. Aug 19, 2022 at 20:34

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.

• You can ask MATLAB to compute the SOS representation of the IIR for you: MATLAB mathworks.com/help/dsp/ref/dsp.iirfilter.sos.html Aug 19, 2022 at 22:31
• @LucaCiti: the example in the documentation would not work here. Once you are in Numerator/Denominator form, it's too late and you can't go back. The trick is to design the filter as poles and zeros and/or second order sections. Don't use [b,a] = butter(N,Fc); use [z,p,k] = butter(N,Fc); instead. Aug 20, 2022 at 9:50
• @Hilmar I agree designing the filter as state space or zeros and poles (then use ss2sos or zp2sos) may be a better option. Aug 22, 2022 at 10:38