Designing IIR Digital Filters with the Bilinear Transform

As we've seen the implementation of digital filters into code is not too hard. The complicated bit is calculating the coefficients and having suitable equations for setting those coefficients so that the filter implements the amplitude and phase characteristics as the frequency changes. Getting this right needs some understanding of the underlying theory and builds a platform from which a number of other filters can be created.

I'm going to walk through the design for a general 2nd order low pass filter and use this as the basis for understanding the different concepts and poking around a bit with the theory. This is not intended in any way to be comprehensive, there are many great posts and books out there with the full detail. However, I wanted to join up the code and theory a bit and have some of fun with the equations.

There are a number of different ways of designing digital filters. The approach I'm going to use here is looking at IIR filters using the bilinear transform method. There are other transforms for designing IIR filters and Fourier methods for FIR design. I'm going to use this approach as the maths isn't too bad and it's enough to get the kind of audio filters that are interesting to play about a bit with.

The bilinear approach starts off looking at the a 'classic' analogue design as these are well established and described, then maps this from the s-domain into the z-domain to create a digital implementation. The idea is to imitate the frequency response in some way and the approach is to find a way to mathematically imitate this in the discrete domain. (I'm very thankful for that course in discrete maths many moons ago).

Starting off then, the generalised transfer function in the s-plane describing the frequency response for a second-order low-pass filter is:

$$ H(s) = \frac {g} {s^2 + ps + g} $$

where s corresponds to the complex frequency and g is the gain to make H(0) = 1.

This is the equation for a classic damped oscillating system.

What we have here is a quadratic polynomial in the denominator (the bit of the fraction below the line). The trick is to calculate the values of g and p that make this function do what we want. First of all, let's see how this function represents functions we want to use. The frequency response is described by setting /( s = j \omega /) and then to see the amplitude response we calculate the square-root of the complex-conjugates of the transfer function. What this means is putting \( s = j \omega \) and \( s = -j \omega \) into H(s) and doing the simplifying maths.

$$ |H(\omega)| = \sqrt { H(j\omega) H(-j\omega) } $$

Which for our transfer function comes out as follows:

$$ |H(\omega)|^2 = \frac {g^2}{(\omega^4+(p^2-2g)\omega^2+g^2)} $$
 
\(\omega\) can now be plugged into the equation to look at the behaviour of the function over different values. In this case \(\omega\) is in the value of radians, or \(0-2\pi\). We can also tweak the parameter p to get some different responses. Which looks like this:



Before moving on we're going to simplify this and normalise the gain g=1:

$$ H(s) = \frac {1} {s^2 + ps + 1} $$

$$ |H(\omega)|^2 = \frac {1}{(\omega^4+(p^2-2)\omega^2+1)} $$

Usually p is set to \(p=\frac {1}{Q}\) giving:

$$ |H(\omega)|^2 = \frac {1}{(\omega^4+(\frac {1}{Q^2}-2)\omega^2+1)} $$

All good. Now the bilinear transform method of getting this function into the digital domain is to substitute \( s \leftarrow \frac {1}{K} \frac {(z-1)}{(z+1)} = \frac {1}{K} \frac {(1-z^{-1})}{(1+z^{-1})}\) using whichever from suites for the maths manipulation where \( K = tan \frac {\omega}{2} \). I've covered the maths for this juggling in a different post, so we can skip it here. What it leads to for our 2nd order filter is a bi-quadratic (e.g. quadratic function in the numerator and denominator), which represents our old friend the BiQuad digital filter:

$$ H(z) = \frac {b_2 z^{-2} + b_1 z^{-1} + b_0}{a_2 z^{-2} + a_1 z^{-1} + a_0} $$

Fantastic, just what we wanted. We now have some numbers for the coefficients for our BiQuad filter implementation, remembering that \( a_0\) is usually normalised to 1.

Now we've got the basic process there's a couple of questions to consider, how do we tune this and tweak the frequency?

In the analog case this is a whole different business, however, it's simplified a little bit in the digital filter design world as the process we've used here is to design a 'prototype' filter response normalised in gain and frequency and we're then mapping this into the digital domain using the bilinear transform using frequency warping. To map our prototype onto different frequency cut-offs we simply adjust the frequency warping term as follows:

$$ K = tan \frac {\omega}{2} \leftarrow tan \left( \pi \frac {f_c}{f_s} \right)$$

Where \(f_c\) is the corner frequency and \(f_s\) is the sampling frequency, so  {f_c}{f_s} can look something like \( \frac {10,000}{44,100} \) for a 10kHz corner frequency and a 44.1kHz sample rate.

Note that the bilinear transform mapping maps the frequency in the s-domain at \( \omega=1 \leftarrow f_c \) in the z-domain for our normalised prototype filter.