One of the most-used filter forms is the biquad. A biquad is a second order (two poles and two zeros) IIR filter. It is high enough order to be useful on its own, and—because of coefficient sensitivities in higher order filters—the biquad is often used as the basic building block for more complex filters. For instance, a biquad lowpass filter has a cutoff slope of 12 dB/octave, useful for tone controls; if you need a 24 dB/octave slope, you can cascade two biquads, and it will have less coefficient-sensitivity problems than a single fourth-order design.
Biquads come in several forms. The most obvious, a direct implementation of the second order difference equation (y[n] = a0*x[n] + a1*x[n-1] + a2*x[n-2] – b1*y[n-1] – b2*y[n-2]), called direct form I:
Direct form I
Direct form I is the best choice for implementation in a fixed point processor because it has a single summation point (fixed point DSPs usually have an extended accumulator that allows for intermediate overflows).
We can take direct form I and split it at the summation point like this:
We then take the two halves and swap them, so that the feedback half (the poles) comes first:
Now, notice that one pair of the z delays is redundant, storing the same information as the other. So, we can merge the two pairs, yielding the direct form II configuration:
Direct form II
In floating point, direct form II is better because it saves two memory locations, and floating point is not sensitive to overflow in the way fixed point math is. We can improve this a little by transposing the filter. To transpose a filter, reverse the signal flow direction—output becomes input, distribution nodes become summers, and summers become nodes. The characteristics of the filter are unchanged, but in this case it happens that the floating point characteristics are a little better. Floating point has better accuracy when intermediate sums are with closer values (adding small numbers to large number in floating point is less precise than with similar values). Here is the transposed direct form II:
Transposed direct form II
Notes and recommendations
Again, direct form I is usually the best choice for fixed point, and transposed direct form II for floating point.
At low frequency settings, biquads are more susceptible to quantization error, mainly from the feedback coefficients (b1 and b2) and the delay memory. Lack of resolution in the coefficients makes precise positioning of the poles difficult, which is particularly a problem when the poles are positioned near the unit circle. The second problem, delay memory, is because multiplication generates more bits, and the bits are truncated when stored to memory. This quantization error is fed back in the filter, causing instability. 32-bit floating point is usually good enough for audio filters, but you may need to use double precision, especially at very low frequencies (for control filtering) and at high sample rates.
For fixed point filters, 24-bit coefficients and memory work well for most filters, but start to become unstable below about 300 Hz at 48 kHz sample rate (or twice that at 96 kHz). Double precision is always costly on a fixed point processor, but fortunately there is a simple technique for improving stability. Looking at the direct form I drawing, the quantization occurs when the higher precision accumulator is stored in the lower precision delay memory on the right side. By taking the quantization error (the difference between the full accumulated value and its value after storing it to memory) and adding it back in for the next sample calculation, the filter performs nearly as well as using full double precision calculations, but at a much lower computational cost. This technique is called first order noise shaping. There are higher order noise shapers, but this one works well enough to handle almost all audio needs, even at high sample rates.
Direct form I with first-order noise shaping
In general, 16-bit fixed point processing is not suitable for audio without double precision coefficients and computation.
Finally, biquads are just one of a DSP programmers tools—they aren’t always the best filter form. There are other filters that don’t share the biquad’s low-frequency sensitivities (in general, biquad coefficient precision is very good at high frequencies, and poor at low ones; there are other filter forms that spread the precision out more evenly, or trade off reduced high frequency performance for better low frequency performance). However, biquads are well known and design tools are plentiful, so they are usually the first choice for an IIR filter unless you find a reason to use another.
There are too many filter forms to cover, but one other filter form that is popular for synthesizers is the state variable filter. It has very excellent low frequency performance, and limitations in the high frequencies that have to be worked around, but most importantly frequency and Q coefficients are separate and easy to change for dynamic filtering. It also make a great low frequency sine wave generator.
Thanks for the great article. This has helped me so much. I still don’t understand much of the process of generating biquad coefficients, but I’ve managed to build a cascading filter. Below is the code I’m using. Is this the correct way of cascading 3 biquads?
biquad = { a0, a1, a2, b1, b2 };
for (var i = 0; i < input.length; ++i) {
x = input[i];
y = x * biquad[0] + z1;
z1 = z2 – biquad[3] * y;
z2 = x * biquad[2] – biquad[4] * y;
x = y;
y = x * biquad[0] + z11;
z11 = z22 – biquad[3] * y;
z22 = x * biquad[2] – biquad[4] * y;
x = y;
y = x * biquad[0] + z111;
z111 = z222 – biquad[3] * y;
z222 = x * biquad[2] – biquad[4] * y;
output[i] = (float) y;
}
(Sorry for the slow approval—getting buried in spam posts lately…)
Yes, that’s pretty much it, except that you’re missing the a1 factor when calculating z1 (z11, z111):
z1 = x * biquad[1] + z2 – biquad[3] * y;
Thank you very much. for your reply. I forgot the reason for omitting that one, the coefficient in my case was probably zero. Cheers.
Well, now I see that it is missing for all three z1s. Thank you for pointing it out.
Hi, I just discovered this blog. Its great!
I am interested in implementing biquad filters in resource-limited fixed-point microcontrollers. I am having trouble finding a good resource that discusses the issues related to fixed-point coefficient calculation and scaling considerations that come with fixed-point numbers.
Right now, I’m relying on cascading single-pole LP IIR filters like the following:
y[n] = (1-a)y[n-1] + ax[n] = y[n-1] + a(x[n] - y[n-1])
where a = 1/2^k, k is some positive integer to be able to implement the multiplication as a right bitshift.
Are there any advantages to using a biquad instead of cascading two of the previously mentioned filters?
If you use two single pole filters, poles and zeros are limited to the x-axis. Try the pole-zero placement applet with the two “Pair” checkboxes unchecked, and you’ll see exactly what you get with two first-order poles and zeros.
Hi Nigel,
Excellent blog, very useful and good written !,
I’m trying to implement a fixed point biquad filter cascade for FPGA as
a part of an auditory filter channel simulating the cochlea for my Msc thesis.
My question is regarding the cascade of direct I forms, what considerations must be followed? Can i directly cascade all the biquad stages or there will be scaling issues ?
Thanks,
Hi Daniel,
There are no scaling issues, except that you may need to adjust the corners. That is, the cut-off is normally spec’d as the -3 dB point of the roll-off. When cascading two identical lowpass filters, the result is -6 dB at that point (the second filter cutting another 3 dB). So, depending on your needs, you might want to move the cutoff frequency up to compensate, or sharpen the corner by adjusting the Q of one of the filters, for instance. A lesser consideration is that if you’re cascading multiple filters of different types and different Q settings, optimal signal-to-noise depends on the order. Search the web for tips on cascading biquads for more info.
Nigel
Hi Nigel, great article and I *love* your site, it’s been very helpful to me. Question on your biquad diagrams, why do you call the feedback coefficients ‘b’ and the feedforward coefficients ‘a’? This seems to be backwards from what gnu octave uses (and wikipedia too, ref this article http://en.wikipedia.org/wiki/Digital_biquad_filter). The convention within octave seems to be that ‘b’ is the numerator coefficients, and ‘a’ are the denominator coefficients for a filter transfer function. The signal processing texts I’ve had also seem to follow this convention. Is it just a matter of personal preferance, or is there something else going on here that I’m missing?
A popular question. Unfortunately, there is no convention. Having ‘a’ on top of the transfer function makes sense to me, and when I went to put things in print, I checked my bookshelf and slightly more than half of my books agreed, so I went with it. Since these filters are almost always normalized to the output coefficient, it’s easy enough to see whether a0 or b0 is missing—if b0 is missing, then the b’s are on the bottom, as I have it.
Thanks, makes sense. One other question, what filter design algorithm are you using? As I understand it, a biquad strictly speaking is just a way to implement any 2nd order IIR filter so you could use it to implement a butterworth filter, chebyshev, or elliptical filter, etc. So which one do you use? Your frequency plots look like a butterworth (I can’t see any ripple in the passband), so I’m guessing that?
Yes, Butterworth (when Q is set to the reciprocal of the sure root of 2, ~0.7071).
This is a very useful set of articles and I’d like to thank Nigel for writing them. I’d also like to provide a little feedback – after spending a fair amount of time reading about BiQuads from various sources, it seems that the usual representation (“canonical”) of the transfer and difference equation uses the a and b coefficient labels in the swapped sense. For example, the transfer equation is almost always written as H(x) = [b0 + b1*z^-1 + b2*z^-2] / [ 1 + a1*z^-1 + a2*z^-2 ], and y[n] = b0*x[n] + b1*x[n-1] + b2*x[n-2] – a1*y[n-1] – a2*y[n-2] for the corresponding difference equation. Notice the labels a and b have been swapped.
I know this is a pretty trivial thing, but I must admit it confused me considerably early on. This is because there are some interesting modifications to the BiQuad structure that modify or add coefficients, and it was a little confusing to continually mentally transpose a and b in these articles. Until I’d read through a fair number of articles, I wasn’t sure where the discrepancy lay. Now I feel confident enough to say that this site uses the unusual labelling.
Other than that, great stuff! Thank you.
If you check the comments on the various biquad pages of this site, you’ll find that’s a popular comment, David. And I admit that it becoming more true as the web is perhaps encouraging more uniformity (wikipedia’s orientation is reversed from mine, for example). But historically it’s been a mix as to whether a is on top of the transfer function (as I show—the zeros) or bottom (poles). Note also that if you start with an FIR (no feedback), and call the coefficients a, then it’s logical to call feedback coefficients b. But as someone pointed out once, those who start with LPC (for speech), coefficients (feedback) are naturally a, and if they need to add input coefficients, they would be b. I also suspect that spelling gets a boost because the canonical, simples (non-transposed) form has a-b from left to right (even though the most useful version, direct form I and direct form II transposed swap sides coefficients).
As I mentioned to another commenter, more than a decade ago when I wrote this article, more than half of the texts on my bookshelf favored a on top. But note also that some include the minus sign as part of the coefficient (so that everything is straight addition in the summation), and other include it as a part of the operation, so you need to pay attention even when the coefficients are named as expected.
Anyway, while I might agree that the way I presented is less popular, it’s far from unusual. I could show you examples from Zölzer, Rabiner and Gold, Ifeachor and Jervis, Pohlmann, Datta, and Julius O. Smith (books I have on my shelf). I know that you can find it the other way from some of these same people too (Rabiner and Crochiere, for instance).
The popular, and freely readable book, The Scientist and Engineer’s Guide to Digital Signal Processing, has the same orientation that I use:
http://www.dspguide.com/ch19/1.htm
Hi Nigel,
I’ve a dubt: does the coefficients (a0, a1 …) useful for 1st or 2nd direct form, or we need specific coefficients for each type of direct form?
I think that the coefficients can be the same, regardless the direct form: don’t so?
Thank you in advance for your reply.
Regards
Hi Marcus. The direct forms are just different mathematical manipulations of the same equation. So, yes, the same coefficients work in all the forms.
Hi Nigel,
I have bit confusion in bandpass filter design through the biquad filter code that you have posted. How should i design a filter of desired bandwidth for bandpass filter. I get to know through internet search that bandwidth is related to Q (B(bandwidth) = W(cutoff freq)/Q). Could you please explain it ?
Many Thanks
Nitesh
Hi Nigel – your blog has helped me immensely – thanks.
I was wondering if there is there is a paper in which the
“Direct form I with first-order noise shaping” topology above is discussed in more detail. I am working on a dithering/noise-shaping algorithm and need some help.
Cheers,
Judd
I think I first saw it in Digital Audio Signal Processing (Zölzer, 1997), which has good detail on quantization effects, and presents first- and second-order noise shaping. I was working on TDM plugins and 56K-based embedded effects, and had run into the problem with biquads in the lower audio range (and exacerbated at higher sample rates). Double precision was awkward (meaning you had to code it) and a big hit for the 56K. I coded my at-risk biquads with first- and second-order noise shaping, and found that first order easily did the job, was a profound improvement over a straight biquad, and a very small hit in cycles.
Noise shaping for dither is a little different, since it’s not for stability within a filter but for plain truncation, and various shaping choices have different merits. Zölzer covers noise shaping associated with dither in the book as well, and you can probably find some of his papers on the web. (Quickly, I found a pdf called “Quantization”, only a slide presentation, but I recognize some of the figures from the book.)
Hi Nigel,
Like your blog very much.
A basic question about the biquad filter: Can we loop a biquad filter 2 times to get a 4th order IIR filter? Or we have to cascade 2 biquad filters to get the 4th order one. What’s the difference and benefit?
Thank you very much.
Finspoo
Good question, Finspoo. You won’t want to literally loop multiple biquads, if you mean run the output back the input and process it again. The reason is that you’ll effectively double the sample rate. You’ll want to run two biquads in series instead. That way they maintain independent internal states (delay elements). But then you need to adjust the Q settings separately to maintain the overall Q. Maybe you’ve run into my recent article on the subject, Cascading filters.
Hi Nigel!
Thank you for your article. I have a little question : you say at the end of your article that other filter forms exist such as SVF. But can’t a SVF be implemented by a biquad? Isn’t it what you do in your biquad design applet?
Thanks in advance for taking the time!
Florent
An SVF’s response can be duplicated by a direct form biquad (either can make a 12 dB/oct lowpass, for instance), yes. The advantages of the SVF are multiple simultaneous outputs and “parameterization”—change a frequency variable to change frequency without affecting Q, change a Q variable to alter it without affecting frequency. It’s all math implementing a transfer function, but there can be advantages to rearranging the math.
The two most popular forms for synths are the state variable and Moog ladder (which is four series one-pole filters plus a feedback path from output to input for resonance, since you can’t get Q from a one-pole. You can get the different filter forms—highpass, bandpass, etc., from taking outputs of the different one-pole, and combinations, with 12 or 24 dB/oct).
Hi Nigel. Thanks for maintaining this awesome resource for DSP. I had a question regarding the graphic here http://www.earlevel.com/DigitalAudio/images/BiquadDFI.gif
you have the two ‘z^-1’ blocks which refer to either input or output samples. When the equation calls for either x[n-2] or y[n-2], why doesn’t the block diagram say ‘z^-2’ ? is it convention or something within block diagrams that when you see multiple z^-n, each successive z^-n is a higher n in the actual equation? I don’t know much about the math behind all of this, so please excuse me if this is a newbish question.
Hi Charles,
These blocks are unit delays, so each is a -1 power. Remember, this is a power, as a superscript, not a locational designation, like a subscript or equivalent like x[n-2]. Make sense? Think of it as a mathematical operation. When you run through two, the product would be z raised to the -2 power. If you had two 2x gain stages, you wouldn’t write the second as 4x, for instance.
Thanks for this great web-site. Very useful indeed (thumbs up) !
I might have missed it, but I wonder about limiting the coefficients values.
For instance I know that some DSPs limit the coefficients to +/ 2.000, so a value of (say) -2.00132… is not legal.
There is some relation between the dynamic range of the signal itself, and the coefficients representation.
I wonder if limiting the coefficients is required for fixed point as well as for floating point?
Also, is this limitation needed more for “direct form 1” than for “direct form II” and “direct form 1 with feedback” ?
Yes, in fixed point DSPs like the 56k family, where coefficients range to 0.999…, it’s common to use half a coefficient value and MAC (multiply-accumulate) twice, for instance.
Great answer.
So, can one claim:
“With ’32 bit floating point’ implementation – no need for biquad limitation.
with ‘fixed point’ a limitation is needed and is based on word length (number of bits) for coefficients and that of the signal itself'”.
would you agree to the above statement ?
It depends on the filter type and settings. But I’d say floating point is easier than fixed point. 😉
Hi Nigel,
Given that the output of the filter is different from the input (else it wouldn’t be filtering anything!) and therefore that the two lots of z values will be different, (one is a history of what went in and the other a history of what went out) how is it possible to “split” the summer? There must be a piece of insight I’m missing.
Hi Chris. I think you’re basically asking how we can arrive at DF2 from DF1. I understand it can be confusing at first look, and I think that’s mostly from looking at the big picture, whereas it’s probably best to focus just on the parts that are changing until you can convince yourself that step is correct, before moving on to the next. The first step is splitting the summer into two summers. Let me describe the sum of products being fed into the summer, in the first diagram, sum = pa0 + pa1 + pa2 + -pb1 + -pb2. In the second diagram, sum1 = pa0 + pa1 + pa2. That sum flows into the second summer, where sum2 = sum1 + -pb1 + -pb2. You can see that sum2 could also be described as pa0 + pa1 + pa2 + -pb1 + -pb2, if you substitute for sum1. So, nothing has changed.
The second step is understanding how the two halfs can be reversed in order, in the third diagram. But the fact is, the first half is a linear filter on its own. And the second half is another linear filter. The order doesn’t matter, so we can reverse them.
Now, maybe the more confusing thing is how the two pairs of delays can become one. but if you can forget about the other stuff for a moment, you’ll see the second pair must contain the same values as the second pair, since they are both connected to the output of the first summer.
Ah, the penny just dropped!
When you show the two halves split, what you don’t mention is that the Z values on the right hand side will now change because there is no left hand side.. or at least there shouldn’t be one inferred – even though the left hand side is still hanging around.
So in that second image, where the summer is split, at that point the values in the right hand side z change.. after all, they won’t have any left hand side contributions – it’s been disconnected from the left hand side even though it is still shown.
Or, to put it another way, the z values stay until the order of the two halves swaps and then the Z values change.
Once again, thank you
OK, I understand the part you were getting hung up on. Yes, that’s right, the z-1 (and z-2) values don’t match, until the “equation” is rearranged. The z values, individually, aren’t important. What is important is the location of the resulting poles and zeros. And those don’t change by re-ordering the evaluation.