Allpole IIR filters
October 2021: I originally wrote this for a different website. I have adapted the article to this blog.
Part 1 of this article describes generalized formulas for any 2pole polynomial, nozero, lowpass or highpass, infinite impulse response (IIR) filter. We then extend the 2pole filter to a generalization for any evenorder allpole polynomial filter. For those who want answers without effort, Part 2 of this article shows recipes for constructing some traditional lowpass and highpass IIR filters, as well as tests of these filters.
I don't show the mathematics I went through in going from one step to the next. The intent here is to provide formulas for reference without having to know the details of how to derive them.
Part 1. Generalized calculations for any evenpole polynomial lowpass/highpass filter, having no zeros
The filters described in this article are allpole filters with no zeros. That is, an allpole filter has a frequency response function that goes infinite (poles) at specific frequencies, but there are no frequencies where the response function is zero. Basically the filter function (also called the transfer function) is a ratio with a constant in the numerator and a polynomial in the denominator. The poles (roots of the polynomial) are located in the negative frequency space, so the actual frequency response (the positive space) is well behaved.
Here we show how to calculate the IIR filter coefficients for the family of filters having two poles (i.e. a second order polynomial in the denominator of the transfer function). Next, we demonstrate how to break down a 4thorder polynomial filter into two different 2pole stages. This is important if, for example, you want a 4pole filter of a particular type; for example the properties of a 4pole Bessel cannot be realized simply by stacking two 2pole Bessels.
Generalized 2pole lowpass filter
First, define the frequency response of the filter as an inverted 2nddegree polynomial function \(H(s)\). This is called the Transfer Function in the \(s\) plane, where \(s\) is the imaginary axis. The amplitude response versus frequency is a function along this axis, analogous to the frequency ranging from −∞ to +∞, with DC at the origin. A generalized 2pole lowpass transfer function has the form:
$$H(s) = \frac{g}{s^2+ps+g}$$where \(g\) is a gain constant causing \(H(0)=1\), and \(p\) is a polynomial coefficient. If you have a coefficient on \(s^2\), you can divide the numerator and denominator by that coefficient and absorb it into \(p\) and \(g\).
For example, for the specific filters shown in Part 2 below, the transfer functions are:
 Butterworth: \(H(s)=\frac{1}{s^2+\sqrt{2}s+1}\)
 Criticallydamped: \(H(s)=\frac{1}{s^2+2s+1}\)
 Bessel: \(H(s)=\frac{3}{s^2+3s+3}\)
If you know only the complex roots \(\alpha \pm j\beta\) of your 2nddegree polynomial in the denominator (which usually means the constants \(g\) and \(p\) aren't nice clean integers), then use
$$p= 2\alpha$$ $$g = \alpha^2+\beta^2$$Next, we calculate the cutoff frequency correction factor \(c\). To do this, we substitute \(j\omega\) for \(s\) in \(H(s)\), where \(j=\sqrt{1}\). Then, multiply this by its complex conjugate (i.e. substituting \(j\) for \(j\) to get the squared amplitude response:
$$H(\omega)^2 = H(j\omega)H(j\omega) = \frac{g^2}{\omega^4  (2g  p^2)\omega^2 + g^2}$$Cascading the filters \(n\) times is equivalent to raising \(H(\omega)^2\) to the \(n\)th power. Because we want the 3 dB cutoff frequency (where the squared amplitude is 0.5) to be the same for any number of cascades, simply set \(H(\omega)^{2n}=\frac{1}{2}\) and solve for \(\omega\), using the quadratic formula to solve for \(\omega^2\). The correction factor \(c\) is simply the inverse of the square root of the result. In other words:
$$c^2 = \frac{1}{\omega^2} = \frac{2}{2g  p^2 +\sqrt{ (2gp^2)^2  4g^2 (1\sqrt[n]{2}) }}$$Take the square root of the above to get the correction factor \(c\), which we apply to the cutoff frequency \(f_0\) in the next step, to preserve the 3 dB cutoff position for any number of \(n\) filter passes.
Because this is an analog filter, we need to make another adjustment to the cutoff frequency, to "warp" it from the analog world to the digital world. For the rest of the calculations, we use the adjusted digital cutoff frequency \(\omega_0\):
$$\omega_0 = \tan\left(\pi c \frac{f_0}{f_s}\right)$$where \(f_0\) is the analog cutoff frequency and \(f_s\) is the sampling frequency.
Now we need the coefficients of the filter. To do this, we perform a bilinear transformation from the \(s\) domain to the \(z\) domain. This isn't as hard as it sounds; it's just algebraic manipulation. Basically it involves substituting \(s = \frac{z1}{\omega_0(z+1)}\) in \(H(s)\), collecting the powers of \(\frac{1}{z}\) together in the numerator and denominator, normalizing things so that the coefficient for \(z^0\) in the denominator is 1, and reading off the coefficients. The \(a_i\) coefficients (which are multiplied by the most recent \(x_i\) inputs) are in the numerator and the \(b_i\) coefficients (which are multiplied by the most recent \(y_i\) outputs) appear as their negative in the denominator. After simplifying, the coefficients are:
$$a_0 = \frac{g\omega_0^2}{1+p\omega_0+g\omega_0^2}$$ $$a_1 = 2a_0$$ $$a_2 = a_0$$ $$b_1 = 2a_0\left(\frac{1}{g\omega_0^2}1\right)$$ $$b_2 = 1(a_0+a_1+a_2+b_1)$$Finally, calculate the filter:
$$y_k = a_0 x_k + a_1 x_{k1} + a_2 x_{k2} + b_1 y_{k1} + b_2 y_{k2}$$Generalized 2pole highpass filter
To get a highpass filter, use \(\omega_0 = 1/\tan\left(\frac{\pi f_0}{c f_s}\right)\) as the adjusted digital cutoff frequency (that is, invert \(c\) before applying to \(f_0\) and invert the tangent to get \(\omega_0\)), and calculate the lowpass coefficients. Then, negate the coefficients \(a_1\) and \(b_1\) — but be sure you calculate \(b_2\) before negating those coefficients! Applying these coefficients in the final filter formula above results in a highpass filter with a 3 dB cutoff at \(f_0\).
Generalized 4pole filter
A 4thorder polynomial filter has a transfer function of the form:
$$H(s) = \frac{g}{s^4+ps^3+qs^2+rs+g}$$Although we can perform the same bilinear transformation here as with the 2ndorder filter, and obtain a onepass fourthorder filter, we'll run into trouble if we do this. We'd find that the filter coefficients may have magnitudes different by several orders (say, one coefficient might be 0.013 and another might be 128.2), resulting in instability due to the finite precision of the computer. Using double precision helps, but the real solution is to decompose the 4thorder filter into a cascade of two 2ndorder filters.
Usually, a 4thorder polynomial cannot be factored readily into two clean 2ndorder polynomials. However, we know that the polynomal can be rewritten in terms of its roots; i.e. the four values \(s= \{R_1, R_2, R_3, R_4\}\) that cause the polynomial to equal zero. Therefore the transfer function becomes:
$$H(s) = \frac{g}{(sR_1) (sR_2) (sR_3) (sR_4)}$$...and we know that \(g\) is the product of all the roots; i.e. \(g = R_1 R_2 R_3 R_4\).
We can calculate the complex roots of a 4thorder polynomial using a number of online tools. For example, Wolfram Alpha is a free webbased math solver that can solve for the roots of a 4thorder polynomial.
Let's use the 4pole Bessel transfer function as an example. The general form of an npole Bessel transfer function is:
$$H(s) = \frac{ \frac{(2n)!}{2^{n} n!}} {\sum_{i=0}^n \frac{(2ni)!}{2^{ni} i! (ni)! } s^n}$$In our case, \(g\) is the numerator. The 4pole Bessel transfer function is then:
$$H(s) = \frac{105}{s^4+10s^3+45s^2+105s+105}$$Wolfram Alpha shows that the four roots of the denominator are:
$$R_1, R_2 = \alpha_1 \pm j\beta_1 = 2.10378939717963 \pm 2.65741804185675 j$$ $$R_3, R_4 = \alpha_2 \pm j\beta_2 = 2.89621060282037 \pm 0.86723412893450 j$$Setting \(g_1 = R_1 R_2 = \alpha_1^2 + \beta_1^2\) and \(g_2 = R_3 R_4 = \alpha_2^2 + \beta_2^2\) and combining the corresponding root pairs in the denominator of \(H(s)\), we get a product of two secondorder transfer functions:
$$H(s) = \frac{g_1}{s^2  2\alpha_1 s + g_1} \cdot \frac{g_2}{s^2  2\alpha_2 s + g_2}$$We set \(p = 2\alpha\) and solve for the coefficients of each 2pole filter stage separately, as described in the previous section.
Our remaining problem involves the frequency correction factor \(c\), which must be the same for each stage. This must be calculated, as with the 2pole filter, by solving for \(\omega\) in \(H{j\omega) H(j\omega) = \frac{1}{2}. This is a quartic function of \(\omega^2\), and unlike quadratics, quartics are nontrivial to solve. Wolfram Alpha can solve them. For a 4pole Bessel lowpass filter, the correction factor \(c = \frac{1}{\omega} = \frac{1}{2.1139176749042}\).
After all this, the 4pole filter is then realized by feeding the input data into the first stage, and feeding that output into the second stage.
Note that the two secondorder filter stages making the whole 4thorder filter are different. That's why cascading two identical 2pole filters doesn't result in a 4pole filter of the same type. That is, if you want a fourthorder Bessel, you can't start out by cascading two secondorder Bessels; you must instead decompose the fourthorder Bessel into two different secondorder filters (which aren't themselves Bessels).
Any evenorder filter can be decomposed into a cascade of 2pole filters in the same manner described above, although for higherorder filters you likely must resort to numerical approximation techniques to solve for the polynomial coefficients and for the frequency correction factor.
Part 2. Recipes for 2Pole IIR filters
You can download an Excel spreadsheet (43 kilobytes, rightclick and "save as") that demonstrates the generalized 2pole polynomial filter. You can play with the g and p polynomial coefficients to see how the filter reacts to an input signal consisting of a square wave impulse followed by Gaussian noise.
Calculation steps
Definitions: \(n\) = number of filter passes (cascaded filters)
 \(\lambda_0\) = 3 dB cutoff wavelength
 \(f_0\ = 1 / \lambda_0\) = 3 dB cutoff frequency
 \(f_s\) = sample frequency
 \(f_c\) = corrected cutoff frequency
 \(\omega_0\) = digital domain cutoff frequency
 \(x_k\) = input value at index \(k\)
 \(y_k\) = filter output value at index \(k\)
Filter type 
Low Pass 
High Pass 


Butterworth 
Critically 
Bessel 
Butterworth 
Critically 
Bessel 

3 dB cutoff correction \(c\) 
$$\frac{1}{ \sqrt[4]{\sqrt[n]{2}1} }$$  $$\frac{1}{\sqrt{\sqrt[2n]{2}1}} $$  $$\frac{1}{\sqrt{3}\sqrt{ \sqrt{\sqrt[n]{2}\frac{3}{4}}\frac{1}{2}} } $$  $$\sqrt[4]{\sqrt[n]{2}1}$$  $$\sqrt{\sqrt[2n]{2}1}$$  $$\sqrt{3}\sqrt{ \sqrt{\sqrt[n]{2}\frac{3}{4}}\frac{1}{2}}$$ 
Polynomial coefficients  $$g=1$$ $$p=\sqrt{2}$$  $$g=1$$ $$p=2$$  $$g=3$$ $$p=3$$  $$g=1$$ $$p=\sqrt{2}$$  $$g=1$$ $$p=2$$  $$g=3$$ $$p=3$$ 
Corrected cutoff frequency \(f_c\)  $$f_c = c \frac{f_0}{f_s}$$  $$f_c = \frac{1}{2}  c \frac{f_0}{f_s}$$  
Ensure \(f_c\) OK for stability constraint  $$0 \lt f_c \lt \frac{1}{8}$$ (unstable when \(f_c \ge \frac{1}{4}\))  $$\frac{3}{8} \lt f_c \lt \frac{1}{2}$$ (unstable when \(f_c \le \frac{1}{4}\))  
Warp cutoff freq from analog to digital domain 
$$\omega_0 = \tan \pi f_c $$  
Filter coefficient calculations  $$K_1 = p \omega_0$$ $$K_2 = g \omega_0^2$$ $$A_0 = \frac{K_2}{1+K_1+K_2}$$ $$A_1 = 2 A_0$$ $$A_2 = A_0$$ $$B_1 = 2 A_0 \left(\frac{1}{K_2}1\right)$$ $$B_2 = 1  (A_0+A_1+A_2+B_1)$$  
Filter coefficients  $$a_0 = A_0$$ $$a_1 = A_1$$ $$a_2 = A_2$$ $$b_1 = B_1$$ $$b_2 = B_2$$  $$a_0 = A_0$$ $$a_1 = A_1$$ $$a_2 = A_2$$ $$b_1 = B_1$$ $$b_2 = B_2$$  
Filter (1 pass)  $$y_k = a_0 x_k + a_1 x_{k1} + a_2 x_{k2} + b_1 y_{k1} + b_2 y_{k2}$$ 
Notes on the three filters
 As shown in the table above, a highpass filter is constructed from a lowpass filter by inverting \(c\) and \(\omega_0\), and negating the lowpass coefficients \(a_1\) and \(b_1\). Note that the identity \(\frac{1}{\tan \pi x} = \tan[\pi(\frac{1}{2}x)]\) is used above to accomplish the inversion of \(\omega_0\) for the highpass filter.
 The Bessel has the most linear phase response, with a stopband frequency response midway between the Butterworth and criticallydamped filters. See the figure below.
 Cascading n 2pole Butterworth filters does not yield a 2npole Butterworth filter. The same is true for Bessel. Cascading two Butterworths, for example, results in what is known as a "LinkwitzRiley" filter.
 On the other hand, cascading criticallydamped filters results in another criticallydamped filter.
 You can cancel the overshoot of a Butterworth with the lag of the criticallydamped filter to obtain another criticallydamped filter that has the positive features of both. For example, a Butterworth of cutoff frequency \(2f_0\) cascaded with a two criticallydamped filters of cutoff \(f_0\) results in a criticallydamped filter that has less lag than the original criticallydamped filter cascaded three times, yet greater stopband attenuation (steeper falloff). Similarly, an arithmetic average of a 1pass \(2f_0\) Butterworth and a 2pass criticallydamped filter of cutoff \(f_0\) results in a criticallydamped filter having a frequency response similar to the 1pass Butterworth (i.e. better than a 1pass critical), and with a lag in between the Butterworth and the 2pass criticallydamped.
Tests
Frequency Response
The following frequency response curves show the performance of each filter above, for different combinations of cutoff frequency and number of filters cascaded together (passes). These plots also demonstrate that the cutoff frequency correction \(c\) above does indeed cause the 3 dB cutoff frequency to occur in the same place for each filter, regardless of the number of passes.
The frequency responses shown resulted from using a chirp as the input signal (i.e., a sinewave that increases in frequency over time). The chirp contained 32,000 data points and covered over 5 octaves, with the wavelength \(1/f\) decreasing linearly to allow each frequency to have nearly the same number of cycles approximating that frequency. The plotted curves are maximum squaredamplitude output levels, updated each cycle. The curves look noisy at the short wavelengths (high frequencies) due to the chirp wavelength nearing the Nyquist limit of 2, resulting in inconsistent wave amplitudes in each cycle. The sampling frequency \(f_s = 1\) in all cases.
Temporal response to a step function
The following timedomain responses to a step function show the performance of each filter, for the same cutoff frequency (1/80 in this case), plotted on the same time scale, for different number of passes. We make the following observations:
 All highpass filters oscillate, even when critically damped. For a step response one can generally expect a criticallydamped filter to exhibit one less zerocrossing than the number of poles in the filter (and the more one cascades filters, the more poles one has). Underdamped filters like the Butterworth oscillate more.
 The criticallydamped filter has the fastest rise time, yet the Bessel appears to converge faster.
 The Bessel actually has a microscopic amount of overshoot, so that it appears to settle faster than the criticallydamped filter; however, the amplitude of the Bessel's oscillations in the highpass case are higher than the criticallydamped filter.
 Although the Butterworth filter has the flattest passband response, sharpest rolloff, and steepest falloff of all the filters shown, in the time domain it suffers from overshoot. The lowpass Butterworth also exhibits the most lag.
References
 YoungHoo Kwon, Butterworth Digital Filters
 D. Gordon E. Robertson and James J. Dowling, "Design and responses of Butterworth and critically damped digital filters," Journal of Electromyography and Kinesiology 13 (2003) pp 569573.
 S.D. Murphy and D.G.E. Robertson, "Construction of a Highpass Digital Filter," presentation to the North American Conference on Biomechanics II (CSB VII), Chicago, USA, 1992.
 Peter Nachtwey, Bessel High Pass
 R.W. Erickson, Filter Circuits, 10 March 1999.
 C.R. Bond, Filter Algorithms, Parameters and Software for Electronic Circuits
 Usenet thread in comp.dsp, "Wanted: criticallydamped highpass IIR filter."
Comments
Post a Comment