1
$\begingroup$

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}$$

Step Response

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: enter image description here Also, in both cases the filter response did not seem to settle toward a final value of 1 and did show some whiggling: enter image description here

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

1 Answer 1

7
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ You can ask MATLAB to compute the SOS representation of the IIR for you: MATLAB mathworks.com/help/dsp/ref/dsp.iirfilter.sos.html $\endgroup$
    – Luca Citi
    Aug 19 at 22:31
  • 1
    $\begingroup$ @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. $\endgroup$
    – Hilmar
    Aug 20 at 9:50
  • $\begingroup$ @Hilmar I agree designing the filter as state space or zeros and poles (then use ss2sos or zp2sos) may be a better option. $\endgroup$
    – Luca Citi
    Aug 22 at 10:38

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.