Part 1 of this article first 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.
Good news 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.
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.
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 the s axis, where s is analogous to the frequency, ranging from −∞ to +∞, with DC at the origin. A generalized 2pole lowpass transfer function has the form:
H(s) = 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) = 1 ⁄ (s^{2} + √2 s + 1) 
Criticallydamped:  H(s) = 1 ⁄ (s^{2} + 2s + 1) 
Bessel:  H(s) = 3 ⁄ (s^{2} + 3s + 3) 
If you only know the complex roots (α ± jβ) of your 2nddegree polynomial in the denominator (which usually means the constants g and p aren't nice clean integers), then use
p = −2α
g = α^{2} + β^{2}
Next, we calculate the cutoff frequency correction factor c. To do this, we substitute jω for s in H(s), where j=√−1. Then, we multiply this by its complex conjugate (i.e. substituting −j for j) to get the squared amplitude response:
H(ω)^{2} = H(jω) H(−jω) = g^{2} ⁄ (ω^{4} − (2g−p^{2})ω^{2} + g^{2})
Cascading the filters n times is equivalent to raising H(ω)^{2} to the nth 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(ω)^{2n} = ½ and solve for ω. Because H(ω)^{2} is a quadratic function of ω^{2}, we use the quadratic formula to solve for ω^{2}. The correction factor c is simply the inverse of the square root of the result. In other words:
c^{2} = 1 ⁄ ω^{2} = 2 ⁄ { 2g − p^{2} + [ (2g − p^{2})^{2} − 4g^{2}(1 − 2^{1⁄n}) ]^{½} }
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 will use the adjusted digital cutoff frequency ω_{0}:
ω_{0} = tan(π c f_{0} ⁄ f_{s})
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 = (z−1) ⁄ [ω_{0}(z+1)] in H(s), collecting the powers of 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) will be in the numerator and the b_{i} coefficients (which are multiplied by the most recent y_{i} outputs) will be their negative in the denominator. After simplifying, the coefficients are:
a_{0} = g ω_{0}^{2} ⁄ (1 + p ω_{0} + g ω_{0}^{2})
a_{1} = 2a_{0}
a_{2} = a_{0}
b_{1} = 2a_{0} (1 ⁄ (gω_{0}^{2}) − 1)
b_{2} = 1 − (a_{0} + a_{1} + a_{2} + b_{1})
Finally, calculate the filter:
y_{k} = a_{0} x_{k} + a_{1} x_{k}_{−1} + a_{2} x_{k}_{−2} + b_{1} y_{k}_{−1} + b_{2} y_{k}_{−2}
To get a highpass filter, use ω_{0} = 1 ⁄ tan[πf_{0} ⁄ (cf_{s})] as the adjusted digital cutoff frequency (that is, invert c before applying to f_{0} and invert the tangent to get ω_{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, will result in a highpass filter with a 3 dB cutoff at f_{0}.
A 4thorder polynomial filter has a transfer function of the form:
H(s) = 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 1pass 4th order 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), which will result 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) = g ⁄ [ (s − R_{1}) (s − R_{2}) (s − R_{3}) (s − R_{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, SymbMath by Dr. Weiguang Huang is a free webbased math solver, implemented as a Java applet, that will solve for the roots of a 4thorder polynomial (in the Run menu, select "Solve Polynomial").
Let's use the 4pole Bessel transfer function as an example:
n pole Bessel: H(s) = g ⁄ (p_{0} + p_{1}s + p_{2}s^{2} + ... + p_{n}s^{n}), p_{i} = (2n−i)! ⁄ [ 2^{n−i} i! (n−i)! ], g = p_{0}
4 pole Bessel: H(s) = 105 ⁄ (s^{4} + 10s^{3} + 45s^{2} + 105s + 105)
SymbMath shows that the four roots of the denominator are:
R_{1}, R_{2} = α_{1} ± jβ_{1} = −2.1037893971796278 ± 2.657418041856753j
R_{3}, R_{4} = α_{2} ± jβ_{2} = −2.8962106028203722 ± 0.8672341289345028j
Setting g_{1} = R_{1}R_{2} = α_{1}^{2} + β_{1}^{2} and g_{2} = R_{3}R_{4} = α_{2}^{2} + β_{2}^{2} and combining the corresponding root pairs in the denominator of H(s), we get a product of two 2ndorder transfer functions:
H(s) = g_{1} ⁄ (s^{2} − 2α_{1}s + g_{1}) × g_{2} ⁄ (s^{2} − 2α_{2}s + g_{2})
We set p = −2α 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 ω in H(jω)H(−jω) = ½. This is a quartic function of ω^{2}, and unlike quadratics, quartics are nontrivial to solve. SymbMath can solve them. For a 4pole Bessel lowpass filter, the correction factor c = 1 ⁄ ω = 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 2ndorder filter stages making the whole 4thorder filter are different. That's why cascading two identical 2pole filters won't result in a 4pole filter of the same type. That is, if you want a 4thorder Bessel, you can't start out by cascading two 2ndorder Bessels; you must instead decompose the 4thorder Bessel into two different 2ndorder 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 will likely have to resort to numerical approximation techniques to solve for the polynomial coefficients and for the frequency correction factor.
If you print this, you may want to do so in landscape mode so that the width of the table fits on the page.
You can also download an Excel spreadsheet (74 kilobytes) 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.
Definitions 
n = number of filter passes 

Filter type 
Low Pass 
High Pass 

Butterworth 
Critically Damped 
Bessel 
Butterworth 
Critically Damped 
Bessel 

3 dB cutoff 
c = (2^{1⁄n} − 1)^{−¼} 
c = (2^{1/(2n)} − 1)^{−½} 
c = 
c = (2^{1⁄n} − 1)^{¼} 
c = (2^{1/(2n)} − 1)^{½} 
c = 
polynomial 
g = 1 
g = 1 
g = 3 
g = 1 
g = 1 
g = 3 
corrected 
f* = c f_{0} ⁄ f_{s} 
f* = ½ − c f_{0} ⁄ f_{s} 

ensure f_{s} is OK for 
0 < f* < ^{1}⁄_{8} 
^{3}⁄_{8} < f* < ^{1}⁄_{2} 

warp cutoff 
ω_{0} = tan(π f*) 

filter 
K_{1} = p ω_{0} A_{0} = K_{2} ⁄ (1 + K_{1} + K_{2}) 

filter 
a_{0} = A_{0} 
a_{0} = A_{0} 

filter (1 pass) 
y_{k} = a_{0} x_{k} + a_{1} x_{k}_{−1} + a_{2} x_{k}_{−2} + b_{1} y_{k}_{−1} + b_{2} y_{k}_{−2} 
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.
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:
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."