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.