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

This entry was posted in Biquads, Digital Audio, Filters, IIR Filters. Bookmark the permalink.

1. Miles Egan says:

Your entire series of posts on filters and resampling has been very helpful. Thanks!

2. weili says:

useful

3. weili says:

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?

• Nigel Redmon says:

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.

4. Marika Picardal says:

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

5. Malcolm Rest says:

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?

• Nigel Redmon says:

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.

• Nigel Redmon says:

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…

• Malcolm Rest says:

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

• Nigel Redmon says:

Thanks Malcolm!

• Nigel Redmon says:

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;

6. Malcolm rest says:

Nigel hi,

I believe there is a typo at “High Shelf/cut” calculation.
“Norm” seems to appear twice.

Thanks for the post !

Malcolm

• Nigel Redmon says:

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

7. Malcolm rest says:

My mistake – just double checked and there is a need for two different “norm” params in the code.

8. Malcolm rest says:

Dear Nigel,

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 !

• Nigel Redmon says:

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…

9. Malcolm rest says:

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

• Nigel Redmon says:

I’ll keep that in mind, Malcolm. I have a few other things in the works first, though…

10. Kermes says:

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?

• Nigel Redmon says:

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

• Kermes says:

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!

• Nigel Redmon says:

Great! Yes, I meant to mention the b0 normalization too…

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

• Nigel Redmon says:

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

12. Vermeille says:

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.

• Nigel Redmon says:

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

• Matthias says:

Hello,
I am basically trying to get the same thing done. Could you please explain a bit more detailed how the coefficients can be scaled? As I understand it, a0, a1, a2 can be scaled freely, as this will result just in a different volume. However, with b0,b1 this doesn’t seem possible so easy.

Basically, I am trying to do:

[code omitted]

Maybe a more complicated approach is required?
Thank you

• Nigel Redmon says:

I encourage you to post these kind of questions, which are not of general need or interest, in a more interactive DSP forum such as kvraudio. I assume you’re doing this because of limitations in whatever environment you’re implementing them in. It’s just basic math, and the exactly answer depends on you needs, so an interactive discussion is what you want if you need help.

BTW, the lack of new articles here is not from lack to stuff to put up or motivation, just my development obligations over the past months keeping me from finishing and posting. But check back because I’ll have some new items posted “soon” 😉

13. Gábor Kakuk says:

Hy Nigel,

Could you help me to change the code in HP / LP / BP / etc. to adjust the gain?

• Nigel Redmon says:

That part’s easy—just multiply the result by your gain factor. Or, if you don’t want to add another operation, multiply the coefficients.

• Gábor Kakuk says:

Wow!
Thank youso much!!!

Another question…or it is an interesting idea I think.

Is it possible to change the code of the shelf filter (Low & High), to get “asymmetric” slopes?
Ok, so. The Q could be positive but also negative.
If it is negative, the only the negative slope will appears on the filters curve, if it is positive, only the positive slope will appears. if it is 0, then no slope.

• Nigel Redmon says:

I’m not sure if I follow exactly what you are looking for, but as far as biquads go, there is a limit on possible shapes with only a pair of complex conjugate poles, and zeros. Point me to a response plot if you can.

14. Niro says:

Dear Nigel,
can you tell me, where have you found these formulas (web, book)?
Thanks

• Nigel Redmon says:

All of the formulas were derived from analog (s-domain) prototypes, which describe the frequency responses of the filters, via the bilinear transform, as shown here. Additionally, I got help with the derivation of the shelves from Digital Audio Signal Processing, by Udo Zölzer. Robert Bristow-Johnson’s cookbook (search: rbj cookbook) has some good information on the analog prototypes as well. He uses a equivalent sin/cos substitution instead of the one involving tan that I use, but although the formulas may look different, for the most part they give identical results. However, note that his shelves aren’t symmetrical, so you’ll find the ones here more practical.

15. Allex says:

Hello Nigel,

Why are the ‘else’ statements for the peak and shelf in your code?
I build the filters in Matlab and the ‘else’ statements result in coefficients that are very similar for negative gains and positive gains and result in the same frequency response.

If i leave the ‘else’ out and use the positive part of your code for all gain values the results are good. (I’m using the ‘dfilt.df1()’ structure and the ‘fvtool()’ to verify the filters)

I’m also trying to determine the coefficient formulas by hand and was wondering if you used the prototypes used in Robert Bristow-Johnson’s cookbook?

Your posts helped me out a lot, thanks!

• Nigel Redmon says:

The else statements are to reverse the roles of the poles and zeros. If you take only the positive half, it will look correct when you give it negative gain, but if you start playing with it at different gain values you’ll see that the negative gains aren’t symmetrical with the corresponding positive ones. There are other ways to make the poles and zeros move the way we want. In this case I took my lead from Zölzer, but rbj’s are equally valid.

16. Rob says:

Hi Nigel,

Let me first say thanks for this great site.

I’ve used your implementation of the Biquad in a filter module for Synthedit, it works very well. I’ve also added the option to cascade up to 4 times.
I want to implement gain compensation, because when someone cascades the Q peak goes really high. How can I do that? Should I subtract a value from the Q at each cascading stage? Or is there a nice formulea I could use to determine what the max Q should be?

Thanks for any help.

• Nigel Redmon says:

Yes, cascading filters gives you issues to consider, not just for the amount of resonance, but even for Butterworth filters the -3dB point changes for the composite filter. (BTW, the wikipedia page for “Butterworth filter” gives a table of how to combine filters yield the proper -3dB cutoff, if that’s the goal.) Some synthesizer filters compensate the resonance gain to reduce the overall gain, lowering the passband gain. (Note this is often in a Moog-style feedback loop, where the compensation factor is obvious.) Completely reducing the overall gain by the amount of peak gain due to Q is not necessarily the best thing to do, especially in more general filters. Just because you have +15dB gain due to Q doesn’t mean you’ll get that must gain on the signal—it depends on the signal content at that frequency. Boosting 12kHz on voice, for instance. So, for those kinds of filters, it’s usually best to have a clipping indicator and let the user adjust accordingly.