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

Aldin

This question appears to be off-topic because it is about signal processing theory and not programming – Silvasilvain

So, just to ask the obvious, have you checked the related wikipedia entry? – Barfield

The formula, as you've provided it, looks like a generic second-order differential equation and you haven't provided the boundary conditions (or equivalent), so currently the question appears unanswerable. Can you provide more context? – Lydia

@Silvasilvain I have used stack overflow before and received good responses, I was unsure of where else I would be able to ask – Aldin

@KarolyHorvath Yes, I have checked Wikipedia, first place I went – Aldin

@DanNissenbaum This is literally all the information I have, maths isn't my strong suite especially when I have no idea how to get numbers from that information. Once I know how to obtain the coefficients I can implement it in code and then work out Pole Zero positions and do it with varying sampling and cut-off frequencies – Aldin

This question appears to be off-topic because it is about signal processing theory and not programming. – Fayina

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

Mistymisunderstand

This helps, is there a way I can modify this to be able to calculate different cutoff and sampling freqs to get the coefficients for any entered? – Aldin

Just use the provided formula, ff is the ratio of cutoff and sampling freq: ff=f_cutoff / f_sampling – Mistymisunderstand

Here's a bit more about computing butterworth filter gains and coding butterworth filters. – Adlee

I'm not from a signals/EE background. Can someone please explain why the coefficients don't have to add to 1? It looks like the equation posted by OP would take in x[n], x[n-1], etc and output y[n] centered around a completely different "mean" signal path – San

The numbers do add up to 1, b0+b1+b2+a1+a2=1, so I'm not sure what you mean. – Mistymisunderstand

Can I upvote this again? Seems like every 15 months I need to type in a Butterworth filter into Excel and all I remember is that it is on StackOverflow somewhere (and not on Wikipedia). Having the exact answers right there in front of me makes it completely trivial that I haven't made a typo. Anyway, see you in 15 months. – Trinidad

I've tried your algorithm, but I don't get a 0 dB gain in the pass band. I asked another question here – Sick

You appear to be off by a factor of "2" in the definition of FF. I replied on your other post. – Trinidad

Does anyone have a reference / explainer about these formulas? – Retaliation

See post by djvg below and his link to robots.ox.ac.uk/~sjrob/Teaching/SP/l6.pdf – Trinidad

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

Apocynthion

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

Bobbobb

Note that Fnorm as used by Matlab and Wn (as used by SCIPY) is half the value of the FF used in Pentadecagon's answer. I just wasted half a day not noticing this.... – Trinidad

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

Simulate

Clarification: a1 and a2 remain the same. only b values change – Tafia

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

Acrosstheboard

Hi, your a's and b's are switched but also you have a negative sign on your b's compared to what, say, pentadecagon posted below. Actually, there seems to be a subtle difference that I cant work out. You have 2*lambda whereas in the math below sqrt(2)*lambda seems to be correct. I typed in your equation and it definitely smooths, but for the same value of lambda it is a bit more aggressive. – Trinidad

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.

Olericulture

Hey Bharat, thanks for your response. Are you able to edit the post and add some examples? – Endicott