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;Thanks Nigel.
Using the above calculation gives good results.
However, I noticed the Frequecy setting of these LS/HS filters is a bit different then accepted.
For example, MATLAB notes that for a LS/HS shelf, the Frequency will set the ‘half gain’ point, while by using this filter the frequency notes the positive peak frequency.
Please see the graph at bottom part fo this example:
https://www.mathworks.com/help/audio/ref/designshelvingeq.html
Any comments on this topic ?
Hi Malcolm,
Shelf definitions vary quite a bit. An analog low shelf EQ is often a lowpass filter summed with the input (that’s why you get that little unexpected dip where the shelf meats unity, that’s characteristic of vintage shelving EQs), so will behave like the definition I use. Checking quickly, I see the MOTU’s MW Equalizer works this way, as does the Waves REQ6—yet their EQ6 uses the center-point definition of frequency. Neutron 2 is more surprising: It uses the center-point definition for both frequency and gain when using the plain “Shelf” Band Shape—so a +5 dB setting gives a +10 dB shelf-top! But if you choose their Vintage setting, gain is the highest point of the shelf, and frequency is not the center point. But I suppose this is why most EQs display their resulting frequency response.
Nigel
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
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
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” 😉
Hy Nigel,
Could you help me to change the code in HP / LP / BP / etc. to adjust the gain?
Thank you in advance!
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.
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.
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.
Dear Nigel,
can you tell me, where have you found these formulas (web, book)?
Thanks
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.
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!
Hi Allex—sorry for the slow reply, I thought I had already answered this…
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.
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.
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.
What is the formula for doing a one-pole filter as shown in your biquad calculator?
A one-pole filter
Hi Nigel. Thanks for a valuable site here.
I need to ‘upsample’ a biquad from FS 22.05 to 44.1. Having the biquad-coefficients from the 22.05 version, do you know a recipe to get the same response from a FS=44.1 biquad?
I’ve had to this in the past, probably for the same reason you need to: you inherit a project where someone baked coefficient into ROM, or a table in code with no comments. If you know how it was calculated, you need to put the coefficients into the formulas and solve backwards. If you have simple filters like second order lowpass, it may be easier to look at the output with an impulse, white noise, or sweep, see where the -3 dB point is, and make a new filter. It’s a harder issue when you have something like a peaking filter that can be done many different ways, and when it’s not a fixed filter, so you need the same response from controls. In any case, the coefficients don’t change linearly with frequency, and don’t even move in the same direction or rate, so it’s not straight forward to upsample IIRs coefficients.
The motivation here is to modernise, and hopefully achieve a more forward-compatible version, a certain musical work, using live-electronics processing of voice. The very distinct response of the biquad in question makes up a special quality in the piece.
The biquad in question is (a’s in numerator, b’s in denominator):
a0: -0.5 a1: 0.5 a2: 0.0 b1: -0.67 b2: 0.74
Doing the maths backways to Z-plane, Freq seems to come out
right, but using the derived radius the resulting biquad for 44.1
misses important acoustical qualities.
Tweaking the Radius up a bit recreates more of the ‘resonant’
sound from the 22.05-version, but audibly (and visually) it’s
clearly not the same response.
Possibly it’s the ff-coeffs? I’ve left the ff-coeffs for now,
hoping these are only related to the gain. Would you know a
recipe for getting these through a similar rate-change?
Here’s the path i’m following in this attempt:
Radius:
b2 = 0.74 = radius^2, radius = sqrt(0.74) = .86
Freq (SR=22050) :
b1 = -0.67 = -2 * radius * cos(frequency * 2 * pi / 22050)
=> -0.67 / (2 * 0.86) = cos(freq * 2 * pi / 22050)
=> acos(-0.67 / (-2 * 0.86)) = freq * 2 * pi / 22050
=> freq = (22050/2pi) * acos(-0.67 / (-2 * 0.86)) = 4108
Back to time-domain coeffs, adjusted for SR=44100:
b2 = 0.74
b1 = (* -2 .74 (cos (/ (* 2 pi 4108) 44100))) = -1.2337
The ultimate goal here is to find a general solution for this problem. Any links to articles, or discussions about similar issue?
-anders:
The bilinear transform (https://en.wikipedia.org/wiki/Bilinear_transform) is a one-to-one mapping between continuous and discrete time. So you could map z back to s, assuming 22050 Hz sampling frequency, and then map s to z assuming 44100 Hz sampling frequency.
The whole thing can be done in one step, if z’ is for your new sampling rate (double) and z is for your old sampling rate, you get
(z-1)/(z+1) = 2*(z’-1)/(z’+1)
==>
z = (1-3*z’)/(z’-3)
which you should be able to substitute z in your biquad to express it in z’.
I haven’t studied your biquad though, if the interesting features are too close to the nyquist frequency, this might not be enough. Bilinear transform with prewarping might get you further in that case.