Calculate Coefficients of 2nd Order Butterworth Low Pass Filter

Question

With the,

Sampling Freq: 10kHz
Cut-off Freq: 1kHz

How do I actually calculate the coefficients for the difference equation below?

I know the difference equation will be in this form, but do not know how to actually work out and come up with the numbers for the coefficients b0, b1, b2, a1, a2

y(n)  =  b0.x(n) + b1.x(n-1) + b2.x(n-2) + a1.y(n-1) + a2.y(n-2)

I will eventually be implementing this LPF in C++ but I need to know how to actually calculate the coefficients first before I can get anywhere with it

Answer 1

Here you go. ff is the frequency ratio, 0.1 in your case:

    const double ita =1.0/ tan(M_PI*ff);
    const double q=sqrt(2.0);
    b0 = 1.0 / (1.0 + q*ita + ita*ita);
    b1= 2*b0;
    b2= b0;
    a1 = 2.0 * (ita*ita - 1.0) * b0;
    a2 = -(1.0 - q*ita + ita*ita) * b0;

and the result is:

b0=0.0674553
b1=0.134911
b2=0.0674553
a1=1.14298
a2=-0.412802

Answer 2

For those wondering where those magical formulas from the other answers come from, here's a derivation following this example.

Starting with the transfer function for the Butterworth filter

G(s) = wc^2 / (s^2 + s*sqrt(2)*wc + wc^2)

where wc is the cutoff frequency, apply the bilinear z-transform, i.e. substitute s = 2/T*(1-z^-1)/(1+z^-1):

G(z) = wc^2 / ((2/T*(1-z^-1)/(1+z^-1))^2 + (2/T*(1-z^-1)/(1+z^-1))*sqrt(2)*wc + wc^2)

T is the sampling period [s].

The cutoff frequency needs to be pre-warped to compensate for the nonlinear relation between analog and digital frequency introduced by the z-transform:

wc = 2/T * tan(wd*T/2)

where wd is the desired cutoff frequency [rad/s].

Let C = tan(wd*T/2), for convenience, so that wc = 2/T*C.

Substituting this into the equation, the 2/T factors drop out:

G(z) = C^2 / ((1-z^-1)/(1+z^-1))^2 + (1-z^-1)/(1+z^-1)*sqrt(2)*C + C^2)

Multiply the numerator and denominator by (1+z^-1)^2 and expand, which yields:

G(z) = C^2*(1 + 2*z^-1 + z^-2) / (1 + sqrt(2)*C + C^2 + 2*(C^2-1)*z^-1 + (1-sqrt(2)*C+C^2)*z^-2')

Now, divide both numerator and denominator by the constant term from the denominator. For convenience, let D = 1 + sqrt(2)*C + C^2:

G(z) = C^2/D*(1 + 2*z^-1 + z^-2) / (1 + 2*(C^2-1)/D*z^-1 + (1-sqrt(2)*C+C^2)/D*z^-2')

This form is equivalent to the one we are looking for:

G(z) = (b0 + b1*z^-1 + b2*z^-1) / (1 + a1*z^-1 +a2*z^-2)

So we get the coefficients by equating them:

a0 = 1

a1 = 2*(C^2-1)/D

a2 = (1-sqrt(2)*C+C^2)/D

b0 = C^2/D

b1 = 2*b0

b2 = b0

where, again, D = 1 + sqrt(2)*C + C^2, C = tan(wd*T/2), wd is the desired cutoff frequency [rad/s], T is the sampling period [s].

Answer 3

You can use this link to get the coefficients of n-order Butterworth filter with specific sample rate and cut of frequency. In order to test the result. You can use MATLAB to obtain the coefficients and compare with program's output

http://www.exstrom.com/journal/sigproc

fnorm = f_cutoff/(f_sample_rate/2); % normalized cut off freq, http://www.exstrom.com/journal/sigproc
% Low pass Butterworth filter of order N
[b1, a1] = butter(nth_order, fnorm,'low');

Answer 4

FYI If you need a high pass filter coefs, all you have to do is use the same code:

const double ita =1.0/ tan(M_PI*ff);
const double q=sqrt(2.0);
b0 = 1.0 / (1.0 + q*ita + ita*ita);
b1= 2*b0;
b2= b0;
a1 = 2.0 * (ita*ita - 1.0) * b0;
a2 = -(1.0 - q*ita + ita*ita) * b0;

but then after multiply all your b terms by ita^2 and negate b1

b0 = b0*ita*ita;
b1 = -b1*ita*ita;
b2 = b2*ita*ita;

now you have a 2nd order high pass filter

Answer 5

my a's and b's are switched here, but my code looks like this:

//Boulanger and Lazzarini, The Audio Programming Book, pg 484
void calculateDifferenceEquation() {
    float lambda = 1.0 / (tan(M_PI * mFrequency / mSampleRate));
    a0 = 1.0 / (1.0 + (2.0 * lambda) + pow(lambda, 2.0));
    a1 = 2.0 * a0;
    a2 = a0;
    b1 = 2.0 * a0 * (1.0 - pow(lambda, 2.0));
    b2 = a0 * (1.0 - (2.0 * lambda) + pow(lambda, 2.0));
}

Answer 6

The best way would be to use something like lab view to simulate your filter and get the coefficients as per your fc and fs. And then use them in c. And finally burn in the code to your microcon. And compare the response with the ones in lab view or simulink.