The bilinear transform is the most popular method of converting analog filter prototypes in the s domain to the z domain so we can implement them as digital filters. The reason we are interested in these s domain filters is that analog filter theory has been around longer than digital filter theory, so we have available to us a number of popular and useful filters in the s domain. The standard biquads used heavily in audio processing are direct implementations from the s domain using the bilinear transform.
The theory of the bilinear transform is well documented in DSP texts and on the web, so rather than spend time going into the theory, we’ll cut to the chase and show how to do the transform, using a second order lowpass filter as an example.
The s domain transfer function of a second order lowpass filter is
To move to the z domain, we need to substitute for s in terms of z. At the same time, we “pre-warp” the response using the tan function in order to wrap the infinite s plane onto the circular z plane. The combined substitution is
Here’s our z domain transfer function after substituting for s:
In order to solve for our filter coefficients, we need to get the equation into the form of a second order polynomial in the numerator (the zeros) and another second order polynomial in the denominator (the poles). That is, we need it to look something like this:
It helps to use a math package to cut down on some of the tedious work. If we take the right side of the equation and ask MathCad to “simplify”, it gives
The first thing we see is that there are too many Q’s, so divide the top and bottom by Q, resulting in
Now combine terms in the numerator and rearrange in descending powers of z; do the same for the denominator:
The z operator means delay, and its power tell us how many samples. Positive powers of z are in the future, and negative powers are in the past. Our filter can’t see into the future (we call such filters “non-causal”), so we need to reduce the powers of z to zero and below (to make it causal). To do this we multiply the numerator and denominator by 
, which subtracts two from each power of z, and gives us our final form:
Our coefficients are now evident:
Because the b0 coefficient is associated with y[n], the filter output, we divide everything by b0 to normalize the filter. Here we also note some convenient relationships between the “a” coefficients:
As a further efficiency, you can compute the common numerator term once and use it where needed. When calculating the coefficients, Q is the filter Q value (for a Butterworth lowpass, where the passband is flat with no corner peaking, use 0.7071, which is the reciprocal of the square root of two). Recall that we defined K as
where Fc is the filter corner frequency you want, and Fs is the sample rate.
Flawless explanation Nigel, excellent tutorials!
The clearest explanations of the process I have ever read! I now get that it’s kinda mechanical in Mathematica or Maple, yet requires human assistance to get right.
Any advice on approaching polynomials where y[0] isn’t present? Do you timeshift the whole expression one sample and go from there or is there a different “correct thing” to do?
Thanks in advance.
Well, y[0] is the output, but yes, in general each additional increment in the index indicates another unit delay—so that if you skip from y[1] to y[3], there would be two unit delays.
Thanks for these very good explanations !
But I’m wondering about gain ? where is it in these equations ? implicit ?
I can’t figure out where it is…!
Hope it’s not to late ! thank you
There is gain, but we normalize to unity gain by dividing all the factors by the output coefficient. So, the lowpass filter has unity gain in the passband, and only corner frequency and Q adjustment. Something like a peaking and shelving filters are a little different, with positive and negative gain for the peaking and shelf portions, and unity gain elsewhere.
Okay, thanks. I understand for unity gain now. But, for instance, in the case of a peaking filter, how do you adjust the gain of the peak ?
Trying to understand that, I played with your calculator and noticed that only a0 and a2 change (I also noticed that a0 and a2 aren’t the same in the case of a peaking filter and that’s where I understood that your equation above is only for lowpass filter, isn’t it ?)
How can I reproduce your reasoning for the peaking filter ?
To clarify, my goal is to convert the coefficient used in the direct form 1 equation to those of the direct form 2 equation (as the DF1 is used in Max/MSP and the DF2 in Pure Data, and I need to “import” curve from Max/MSP to Pure data).
Tank you again !
OK, first, you don’t need to convert coefficients for the different biquad direct forms. The different forms are simply transpositions. My C++ source uses transposed direct form II, which works well for floating point, but I’ve used the same in direct form I, which suits the 56k DSP family. So, if you want filters that behave the same as in my Biquad calculator v2, then you can use the equations given in the source code.
Also, you noted that only a0 and a2 change for the peaking filter gain, but that’s only true for boost—cut requires a different arrangement to keep the response symmetrical (the poles are fixed and the zeros as moved for boost, zeros are fixed and poles are moves for cut).
Basically, a peaking filter is a bandpass response added to the input. But if you simply mix a bandpass with the direct signal, you’ll have phase issues where they come together (because the IIR filter is not phase linear). So, you actually have to work it into the equation then do the bilinear transform on the whole thing.
For instance, since we have the lowpass equation above, 1 / (s^2 + s/Q +1), I’ll use the example of a low shelf filter. It would be the lowpass, times a gain value (H, which is zero-based), plus the direct signal (1); Q is a complicated issue with shelving filters, so I’ll revert to Butterworth response, where 1/Q is sqrt(2): H / (s^2 + s*sqrt(2) +1) + 1. We put that into a common denominator: (H + (s^2 + s*sqrt(2) +1)) / (s^2 + s*sqrt(2) +1). We want unity gain to be a gain of 1, so we substitute V for H + 1, and end up with: (V + s^2 + s*sqrt(2)) / (s^2 + s*sqrt(2) +1). Now, you run that through the bilinear z transform, and that’s where you get the coefficient formulas. For the boost case anyway…
Wow, okay thanks a lot. I’m a bit lost in math (and english also!).
So, you say that I don’t need to convert coefficients between DFI and DFII but, when I try to apply those of DFI (from max msp’s biquad~) to a biquad~ using DFII (the one of puredata), I can’t ear a sound. It should mean that the two biquads don’t use the coefficients the same way (and It’s normal since they are not based on the same equation). But there should be a way to convert (or transpose) the coefficients between one and the other…
Your coefficients don’t work with puredata so It’s a pitty for me, I can’t use your sources…
I’ll have to be creative to do what I want to do ! 🙂
Thank you for your help anyway !
One key gotcha is that while “a” and “b” coefficient naming is common, there is no convention on which one is the feedback (pole coefficients) and which are the feed-forward (zero coefficients). Try using the a2 coefficient for b2, etc. I use the convention that “a” is on top of the transfer function, the left side of a direct form I diagram, as do many other classic text books.
Yes, I tried changing a for b, etc, but it still doesn’t work.
I think I’m going to calculate the coefficients, starting with frequency, gain and q. Now I’ve to find how…!
I tried your calculator again and some type of filter seem to work (I can ear sound through a biquad objet) but not as expected (for instance a “selective” bandpass filter – with Q of 10 – doesn’t filter the audio a lot as it should). Very strange behavior…
I looked up the Max biquad~ object, if that’s what you are using:
y[n] = a0 * x[n] + a1 * x[n-1] + a2 * x[n-2] – b1 * y[n-1] – b2 * y[n-2]
This is precisely the convention that I use, so there should be no change needed.
For a filter that’s not working right for you, give me one exact set of parameters for the calculator (v2), and the list of coefficients it produces. Describe the incorrect filter response that you are seeing as clearly as you can.
Yes that’s what I thought. But my goal is to create a filter in max msp and to import it into pure data. My issue is in pure data because it uses these equations :
y[n] = b0*w[n] + b1*w[n-1] + b2*w[n-2]
with :
w[n] = x[n] – a1*w[n-1] – a2*w[n-2]
And I’m pretty sure that it’s not the same set of coefficients between those two (I tried to put them into pure data, to change a and b etc…)
With those equations, a and b coefficients need to be swapped. However, if that still doesn’t work, it may be that the equations are misstated—another very common difference in implementations is whether the minus signs are rolled into the coefficients or not (some people like to use all additions). So after you swap a and b, and confirm that it still doesn’t work, flip the signs of the feedback coefficients. (If that doesn’t work, you can try other combinations, but the feedback coefficients are the usual issue.)
This is it !
So first, I had to invert a and b, then, to flip the signs of b1 and b2. In fact, Puredata equation is :
y[n] = b0*w[n] + b1*w[n-1] + b2*w[n-2]
with :
w[n] = x[n] + a1*w[n-1] + a2*w[n-2] (all additions, as you said. I made a mistake in my previous reply)
Thank you very much !
PS: May I post a link to this website to inform the Max/MSP community ? A post talks about this with no solutions, It could help.
Yes, of course you can post a link! Glad you got the problem solved!
Very good explanation! can it be used for a notch filter as well?
Yes—I derived all the filters, including notch, in the biquad calculator and biquad source code, exactly this way.
Hi Nigel,
I have a question: When doing the pre-warping and substitute s to Z at the same time , the K=tan(wT/2) , how does this tan(wT/2) came from ?
(My understanding is to pre-warp by Wc=2/T*tan(wc/2) first than substitute s to Z by s=2/T*(z-1)/(z+1) )
thanks
This shows the derivation of the tan warping: https://en.wikipedia.org/wiki/Bilinear_transform#Frequency_warping
I have to ask; how do you arrive at s = 1/tan(wT/2) ((z-1) / (z+1)). There are multiple sources that uses that particular formula for s, but neglects to explain how the omega in the numerator got lost? The wikipedia page you’re linking to states the formula for s with an omega in the numerator.
So, why s = 1/tan(wT/2) ((z-1) / (z+1)) instead of s = w/tan(wT/2) ((z-1) / (z+1)) ?
Ah, too broad of a question for a quick comment. There are so many places to find the derivation of prewarping—I’d rather not get into looking at why a particular one was stated a certain way.
Darn — most articles will simply state it, and neglect the derivation, leaving the reader wondering how the …. it got to be 1 instead of w, or vice versa.
Where else? I’m not sure what they are even getting at on wikipedia formula, as s would be in radians per second in that case, so it seems like a mistake. Anyway, there are different conventions—some would say that you scale s by dividing by 2/T * tan(ωT/2), then substitute 2/T * (z-1)/(z+1) for s. In that case, 2/T cancels, and it reduces to the same substitution I gave. So, either take my word for it, or someone else’s, or understand why you’re substituting (z-1)/(z+1) as a mapping approximation and scaling with tangent, and follow through with the logic. It you trust that derivation, it’s probably a mistake to spend much time trying to figure out why an isolated formula doesn’t fit to the scheme you’re trying to use.
hi. i know its been a long time since this topic started but i hoped you can help me with something.
I’m trying to find a connection between filter coeficients and filter design im getting from here: http://www.schwietering.com/jayduino/filtuino/index.php?characteristic=bu&passmode=lp&order=1&usesr=usesr&sr=100&frequencyLow=10¬eLow=¬eHigh=&pw=pw&calctype=float&run=Send
i wanted to use those formulas for digital filters but i cant find how are they calculated. its either ready application that gives me program like v[1]=(2.452372752527856026e-1 * x)+(0.50952544949442879485 * v[0]);
return(v[0] + v[1]);
or explanation thats hard to use in program (im not that advanced).
i want to find a way to calculate coeficcients and maybe formulas i can use in microcontroler for n-order filters
Can you help me?
Explore other articles here, I think you’ll find what you’re looking for. You gave an example of first order filter, you can find C++ code here: A one-pole filter. You can find biquad code here: Biquad C++ source code.
i tried to calculate those coefficients and it looks like a site I mentioned is omitting coefficients? is there some siplifcation at work or i dont see something? 1st order filter was just example. You can generate 2nd and 3rd order there but i cant get same coefficients or equation like them…
That calculator gives the same results as mine, basically. Go to Biquad calculator v3 the default is 100 Hz lowpass at 44.1k. Here is that same spec with the other calculator: schwietering 100 Hz. Compare my coefficient output with the generated code. For instance, notice that for a lowpass, a0 and a2 are the same, so the number appears just once. And a1 is twice that, and you have a “2 *” in the code. It works out the same.
but we are missing one coefficient (0.00010048283637747806). in your calculator we multiply those with with past inputs and outputs. in there we save not output but some inermediary value. thats why im confused. its close but different
That’s the one I was referring to about the “2 *” in the code (2 * 0.00005024141818873903 = 0.00010048283637747806).