For fixed filters, we can plug biquad coefficients into our programs. But often, we need to calculate them on the fly, to user settings or changes in sample rate. As a companion to the biquad calculator, here are the formulas used, in JavaScript; refer to our article on biquads for implementation diagrams:
function calcBiquad(type, Fc, Fs, Q, peakGain) {
var a0,a1,a2,b1,b2,norm;
var V = Math.pow(10, Math.abs(peakGain) / 20);
var K = Math.tan(Math.PI * Fc / Fs);
switch (type) {
case "lowpass":
norm = 1 / (1 + K / Q + K * K);
a0 = K * K * norm;
a1 = 2 * a0;
a2 = a0;
b1 = 2 * (K * K - 1) * norm;
b2 = (1 - K / Q + K * K) * norm;
break;
case "highpass":
norm = 1 / (1 + K / Q + K * K);
a0 = 1 * norm;
a1 = -2 * a0;
a2 = a0;
b1 = 2 * (K * K - 1) * norm;
b2 = (1 - K / Q + K * K) * norm;
break;
case "bandpass":
norm = 1 / (1 + K / Q + K * K);
a0 = K / Q * norm;
a1 = 0;
a2 = -a0;
b1 = 2 * (K * K - 1) * norm;
b2 = (1 - K / Q + K * K) * norm;
break;
case "notch":
norm = 1 / (1 + K / Q + K * K);
a0 = (1 + K * K) * norm;
a1 = 2 * (K * K - 1) * norm;
a2 = a0;
b1 = a1;
b2 = (1 - K / Q + K * K) * norm;
break;
case "peak":
if (peakGain >= 0) { // boost
norm = 1 / (1 + 1/Q * K + K * K);
a0 = (1 + V/Q * K + K * K) * norm;
a1 = 2 * (K * K - 1) * norm;
a2 = (1 - V/Q * K + K * K) * norm;
b1 = a1;
b2 = (1 - 1/Q * K + K * K) * norm;
}
else { // cut
norm = 1 / (1 + V/Q * K + K * K);
a0 = (1 + 1/Q * K + K * K) * norm;
a1 = 2 * (K * K - 1) * norm;
a2 = (1 - 1/Q * K + K * K) * norm;
b1 = a1;
b2 = (1 - V/Q * K + K * K) * norm;
}
break;
case "lowShelf":
if (peakGain >= 0) { // boost
norm = 1 / (1 + Math.SQRT2 * K + K * K);
a0 = (1 + Math.sqrt(2*V) * K + V * K * K) * norm;
a1 = 2 * (V * K * K - 1) * norm;
a2 = (1 - Math.sqrt(2*V) * K + V * K * K) * norm;
b1 = 2 * (K * K - 1) * norm;
b2 = (1 - Math.SQRT2 * K + K * K) * norm;
}
else { // cut
norm = 1 / (1 + Math.sqrt(2*V) * K + V * K * K);
a0 = (1 + Math.SQRT2 * K + K * K) * norm;
a1 = 2 * (K * K - 1) * norm;
a2 = (1 - Math.SQRT2 * K + K * K) * norm;
b1 = 2 * (V * K * K - 1) * norm;
b2 = (1 - Math.sqrt(2*V) * K + V * K * K) * norm;
}
break;
case "highShelf":
if (peakGain >= 0) { // boost
norm = 1 / (1 + Math.SQRT2 * K + K * K);
a0 = (V + Math.sqrt(2*V) * K + K * K) * norm;
a1 = 2 * (K * K - V) * norm;
a2 = (V - Math.sqrt(2*V) * K + K * K) * norm;
b1 = 2 * (K * K - 1) * norm;
b2 = (1 - Math.SQRT2 * K + K * K) * norm;
}
else { // cut
norm = 1 / (V + Math.sqrt(2*V) * K + K * K);
a0 = (1 + Math.SQRT2 * K + K * K) * norm;
a1 = 2 * (K * K - 1) * norm;
a2 = (1 - Math.SQRT2 * K + K * K) * norm;
b1 = 2 * (K * K - V) * norm;
b2 = (V - Math.sqrt(2*V) * K + K * K) * norm;
}
break;
}
return [ a0, a1, a2, b1, b2 ];
}
Your entire series of posts on filters and resampling has been very helpful. Thanks!
useful
Hi Nigel,
Do you have plan to write a tutorial how to calculate the FIR coefficient based on cut-off frequency for lowpass or bandpass filters?
I’ll keep that in mind. I’ve been busy—I have a couple of partially written articles on other subjects that I need to get back to and finish up, but I’ll think about an article on FIRs. Thanks for the feedback.
Hello Nigel,
I am planning to implement a Biquad IIR filter on an FPGA platform. Right now, I’m searching for a simple yet useful Audio Application for the biquad filter. I would like to ask for your suggestions.
Thanks,
Marika
Nigel hi,
Thanks for the great site.
What is the applicational difference between a notch and negative peak with hi-Q?
Could you share/upload a calculation for Low/High Shelf EQ with “Q factor” for their sharpness?
Hi Malcolm,
The application difference between a notch and a negative peak is that the notch, which comes to a very sharp point, is usually used to target a specific frequency; most often it’s used to remove hum from the AC power line. Power line frequency (60 Hz in the US, 50 Hz in the UK) is regulated precisely (if not, power grids would have trouble combining power from multiple generators), so surgical precision can be used to remove it and its harmonics. Moving notches are also used as an effect (as in a phase-shifter).
Negative peak filters are more often used to de-emphasize a frequency area—perhaps a recording that was made in a room with a bad resonance.
Oh, and the Q for shelves…that’s a good idea, I’ll keep that in mind as I think about the next filter article…
Dear Nigel hi,
Let me contribute to this great web site, not just “nag” you all the time
I have calculated and tested the needed code change for “Q_Factor” for Low_Shelf & High_Shelf.
It is done by simply dividing the second mathematical expression of each coefficient, besides a1 & b1, by Q.
Here is an example based on your “Low Shelf”:
case “lowShelf”:
if (peakGain >= 0) { // boost
norm = 1 / (1 + Math.SQRT2 * K /Q+ K * K);
a0 = (1 + Math.sqrt(2*V) * K /Q+ V * K * K) * norm;
a1 = 2 * (V * K * K – 1) * norm;
a2 = (1 – Math.sqrt(2*V) * K /Q+ V * K * K) * norm;
b1 = 2 * (K * K – 1) * norm;
b2 = (1 – Math.SQRT2 * K /Q+ K * K) * norm;
}
else { // cut
norm = 1 / (1 + Math.sqrt(2*V) * K/Q + V * K * K);
a0 = (1 + Math.SQRT2 * K /Q+ K * K) * norm;
a1 = 2 * (K * K – 1) * norm;
a2 = (1 – Math.SQRT2 * K /Q+ K * K) * norm;
b1 = 2 * (V * K * K – 1) * norm;
b2 = (1 – Math.sqrt(2*V) * K /Q+ V * K * K) * norm;
}
Hope this is useful to some…
Malcolm
Thanks Malcolm!
Just getting back around to this…After my recent C++ article, it just occurred to me that I didn’t work in Q into the shelves. I did, today, then remembered your comment…yours was different than mine, so I tried it out, and see a problem. Your version won’t work right for “flat” (Q = 0.7071) response.
I did a quickie implementation, replacing all sqrt(2) with 1/Q. That should be correct, but I didn’t go through the BLT for a couple of reasons: first, lack of time, second, doing it this way (tweaking the pole/zero relationships on both sides of the shelf frequency) gives an overshoot on both sides. From what I’ve seen, the “good sounding” way to do it is to have the overshoot on the unity side of the shelf. So I need to look closer at the best “musical” implementation, but for now here’s the version with overshoot on both sides for Q > 0.7071 (in C++):
case bq_type_lowshelf: // with Q if (peakGain >= 0) { // boost norm = 1 / (1 + 1/Q * K + K * K); a0 = (1 + sqrt(V)/Q * K + V * K * K) * norm; a1 = 2 * (V * K * K - 1) * norm; a2 = (1 - sqrt(V)/Q * K + V * K * K) * norm; b1 = 2 * (K * K - 1) * norm; b2 = (1 - 1/Q * K + K * K) * norm; } else { // cut norm = 1 / (1 + sqrt(V)/Q * K + V * K * K); a0 = (1 + 1/Q * K + K * K) * norm; a1 = 2 * (K * K - 1) * norm; a2 = (1 - 1/Q * K + K * K) * norm; b1 = 2 * (V * K * K - 1) * norm; b2 = (1 - sqrt(V)/Q * K + V * K * K) * norm; } break;Nigel hi,
I believe there is a typo at “High Shelf/cut” calculation.
“Norm” seems to appear twice.
Thanks for the post !
Malcolm
Hi Malcolm. As you noted later, yes, “norm” was recalculated. In fact, I grabbed the calculation from Udo Zölzer’s Digital Audio Signal Processing, essentially, for the biquad calculator code, and reproduced that code here. But you made me think…why should it be? The “norm” variable is short for “normalize”—it’s the b0 term by which all the other components are normalized. I decided to rearrange the b1 and b2 calculations to make them more uniform with the rest of the filters—thanks for mentioning it. I suppose that the lack of uniformity in Zölzer’s book was just a minor oversight—either form produces the same results. I’ll list the original interpretation here, and make the change (which is to multiply numerators and denominators for b1 and b2 by V) in the article:
case "highShelf": ... else { // cut norm = 1 / (V + Math.sqrt(2*V) * K + K * K); a0 = (1 + Math.SQRT2 * K + K * K) * norm; a1 = 2 * (K * K - 1) * norm; a2 = (1 - Math.SQRT2 * K + K * K) * norm; norm = 1 / (1 + Math.sqrt(2/V) * K + K * K / V); b1 = 2 * (K * K / V - 1) * norm; b2 = (1 - Math.sqrt(2/V) * K + K * K / V) * norm; }My mistake – just double checked and there is a need for two different “norm” params in the code.
Dear Nigel,
Thanks for your prompt answers, must appreciated.
I was looking up some other sources for IIR EQ filtering and noticed some of them use trigonometric functions to shift from exp.(-jw) field for the calculation part itself.
It looks like trigonometric calculations make the code much shorter. Is there a side effect to using trigonometric calculations?
Are they less stable/precise?
Thanks again !
Hi Malcolm,
Well, the “j” in exp(-jw) is the square root of minus one, so you have to use trig identities to expand it out then…it’s no longer shorter
But I plan to get to some more filters, and look at something more synth-oriented, some time after I get through the last couple of installments on the oscillator…in progress…
Dear Nigel,
I wonder if you can find the time to upload an IIR based formula for “all pass” filter for some “phase correction games”.
Thanks !
Malcolm
I’ll keep that in mind, Malcolm. I have a few other things in the works first, though…
I implemented the code for calculating low pass filter coefficients and got the same coefficients in C code and in GNU Octave (MATLAB-alike software) as in your BiQuad calculator. For retro computing purposes of emulating an ancient sound card, I wish to emulate a fourth order Sallen-Key topology which should implement a Butterworth response to my understanding. Fs=49716, Fc=15392, Q=1.25 for the first BiQuad section, and Fs=49716, Fc=15392, Q=0.5405 for the second.
But when I tried the coefficients in GNU Octave and tried to plot frequency response of one BiQuad with freqz, it is all wrong. MATLAB (and Octave) seems to want the coefficients in Transposed Direct Form II, while there is no mention what form your calculator gives the coefficients, but I have to assume they are Direct Form I coefficients. In fact based on your Biquad tutorial, the coefficients should be exactly the same in all forms, the arrangement of delays just differ?
Can you give any pointers how to use those coefficients in MATLAB/Octave to plot a frequency response with freqz command? Do you also know how to plot the frequency response of both BiQuads combined?
Yes, the coefficients are the same in all of the direct forms—they are just computational rearrangements.
The problem you have is that for freqz, the “numerator and denominator coefficients are b and a, respectively” (from the freqz definition). If you look at my derivation here, you’ll see that I have a as the numerator, and b as the denominator (see the fourth equation down). Unfortunately, this is a typical problem—there is no standard “spelling”. I once surveyed all of my DSP books, and the way I spell the coefficients won out—and happened to be the one I prefer—but it just depends on whether you start with the transfer function, where it makes sense to have a on top, or the direct form I implementation, where it makes sense to have a on the left (which is equivalent to a on the bottom of the transfer function).
Sorry, the problem was that one needs to have also three b coefficients, b1 b2 from your calculator, and b0 is 1. After building a vector b with [1 b1 b2] it works. Of course b and a correctly placed as numerator and denominator. Thank you!
Great! Yes, I meant to mention the b0 normalization too…
Thanks for this article and the biquad calculator!
I just discovered this today. I’m trying to add a biquad filter to my own little dsp project. I started out from examples like http://peabody.sapp.org/class/350.838/lab/biquad/ and the pure data source code. I got a basic object ready, based on this algorithm:
for (Int i = 0; i < length; i++) {
output = *inPtr++ + fb1 * last + fb2 * previous;
*outPtr++ = ff1 * output + ff2 * last + ff3 * previous;
previous = last;
last = output;
}
which works with the example parameters provided. But it would be nice if I could use the object by providing frequency, gain and Q instead of coefs. After looking at this code, that seemed an easy task. Just map a0-2 to ff1-3 and b1-2 to fb1-2, right?
That didn't work. I haven't written an object yet to plot out what the results are, but at any rate there's no sound and there should be. I've looked at all code I could get my hands on, but I can't figure out what to do different. Do I need to map this differently, or is it not just about mapping?
I'm sorry if I'm asking a very dumb question here. I'm more an artist and math is not my strongest point.
Regards,
yvan
Hi Yvan,
You are correct about the mapping. The code you gave implements the direct form II diagram shown here… [And after further discussion privately, we determined that Yvan needed to negate the feedback coefficients (b1 and b2) as well. This is a common problem—some people use a for the feedforward terms and b for feedback, and others do the opposite; and others subtract the feedback terms and some add negated feedback terms. DSP texts vary on this choice, so you need to pay attention that the coefficient calculations agreed with the filter realization when implementing your biquads.]
Regards,
Nigel
Thank you for your amazing blog. It helped me a lot.
Following your formulas, I don’t see a scale factor, and for Q > 1 values, it may cause overflows. How to calculate a scale factor ?
Thank you again.
Vermeille.
Hi Vermeille—The calculations are floating point. If you want to implement a filter in fixed-point, you may have to scale and adjust when evaluating the filter. (For instance, coefficients can be in the range of ±2; you’d halve the relevant coefficients, then double products in the filter evaluation.)
Nigel