CODEWITHSUNDEEP
c++signal-processing
Daniel
Daniel

Reputation: 379

Calculate Coefficients of 2nd Order Butterworth Low Pass Filter

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

Upvotes: 12

Views: 32318

Answers (6)

Clay Bailey
Clay Bailey

Reputation: 11

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));
}

Upvotes: 1

djvg
djvg

Reputation: 14255

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].

Upvotes: 12

wmohr
wmohr

Reputation: 41

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

Upvotes: 4

Bharat
Bharat

Reputation: 1

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.

Upvotes: -1

Richard Le
Richard Le

Reputation: 669

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');

Upvotes: 4

pentadecagon
pentadecagon

Reputation: 4847

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

Upvotes: 19

Related Questions