Design prototypes for stable, efficient, parametric IIR filters.

In the last article, I outlined the process of creating a parametric filter. The steps are

  1. Decide on the filter type.
  2. Design an analog prototype.
  3. Digitize the analog prototype using the bilinear transform.
  4. Implement the digital filter.

Here’s how the process looks:

Figure 1. Parametric filter design workflow.

In this article, we’ll discuss the second step of the process: designing the analog prototype.

Figure 2. In this article, we discuss analog prototype design.

Recap

As you remember from the previous article, parametric filters must have [VälimäkiReiss16]

  • interpretable, real-time-adjustable controls and
  • low processing delay.

This led us to choose infinite-impulse response (IIR) filters. To streamline the process of their design and to ensure that they remain stable, we said that the easiest way to come up with these filters is to design them in the analog domain and then digitize them.

How to design them in the analog domain, then?

Our Goal

Designing a filter in the analog domain is traditionally done by designing a low-pass filter with some of the desired characteristics and then transforming it to the desired filter type, for example, high-pass. This can be done with transformations like lowpass-to-bandpass transformation or lowpass-to-highpass transformation.

What is more, we can set the cutoff frequency of the low-pass filter to 1, because this frequency will be altered by the bilinear transform anyway.

So our first goal is to design a low-pass filter with the cutoff frequency equal to 1.

It’s all downhill from there. 😉

Analog Filters Design Methods

Filter design in the analog or digital domain is the process of approximating the desired frequency response with a certain set of constraints. [Smith07]

As such, it may be considered a form of constrained optimization.

There are many methods to achieve this, as there are many optimization methods. There are, however, 4 basic filter approximations considered as standard. Each of them is optimal in a different sense [OppenheimSchafer10].

In Short
Standard Analog Filter Design Methods
Method What is optimal?
Butterworth The amplitude response is maximally flat in the passband.
Chebyshev type I The amplitude response is equiripple (has ripples on the curve of a fixed width) in the passband and monotonic in the stopband.
Chebyshev type II The amplitude response is monotonic in the passband and equiripple in the stopband.
Elliptic functions The amplitude response has equiripple error in the passband and the stopband.

But in equalizer filters mostly Butterworth responses are used, because the amplitude response is monotonic (without any ripples) and the higher the frequency above the cutoff frequency, the bigger the filter’s attenuation [Zölzer08].

Additionally, it is easy to control the slope of the roll-off above the cutoff frequency with the filter order. If the filter order is NN, its attenuation in the stopband is N6N \cdot 6 dB per octave (each doubling of frequency) [Zölzer08].

Analog Prototype Butterworth Low-pass

We know that we want to design an analog prototype low-pass using the Butterworth approximation. What do we want to approximate exactly?

Approximation Goal

The goal of the approximation is the ideal low-pass filter.

Figure 3. Amplitude response of the ideal low-pass filter with cutoff frequency equal to 1.

Frequency ωa\omega_\text{a} is the analog cutoff frequency in radians per second. We assume that ωa=1\omega_\text{a} = 1, i.e., the filter’s transfer function is normalized.

Our only constraint is the filter order. According to [Zölzer08] the most commonly used orders are N=2N = 2 and N=4N = 4.

Butterworth Filter Derivation

WARNING: This part is math-heavy. It is intended for those who want to fully understand the derivation of analog prototypes. If you don’t want to get that deep, just use tabularized, ready-made formulas. You can skip to their examples here.

Note: This part is based on the great explanation from [ParksBurrus87].

The frequency response of an analog filter is found by evaluating its transfer function H(s)H(s) along the imaginary axis, i.e., for s=jωs=j\omega.

We will formulate the problem of approximating the ideal filter in terms of the squared magnitude response H(jω)2|H(j\omega)|^2. That is because H(jω)2|H(j\omega)|^2 is an analytic, real-valued function of a real variable, i.e., high-school mathematics apply. Another reason is that H(jω)2|H(j\omega)|^2 is proportional to the energy or power of the signal, which may be useful depending on the context.

To be able to get from the squared magnitude response to the transfer function, we introduce an intermediate, complex-valued function of the complex variable ss

Ha(s)=Ha(s)Ha(s).(1)\mathcal{H}_\text{a} (s) = H_\text{a} (s) H_\text{a} (-s). \quad (1)

It can be easily shown that

Ha(s)s=jω=Ha(jω)2.(2)\mathcal{H}_\text{a}(s) \Bigr \rvert_{s=j\omega} = |H_\text{a}(j\omega)|^2. \quad (2)

The Butterworth squared magnitude response Ha(ω)=Ha(jω)2\mathcal{H}_\text{a}(\omega) = |H_\text{a}(j \omega)|^2 is a Taylor series approximation of the ideal squared magnitude response around ω=0\omega=0.

Taylor Series

Taylor series around ω=0\omega=0 is

Ha(jω)=K0+K1ω+K2ω2+=k=0Kkωk,(3)H_\text{a}(j\omega) = K_0 + K_1 \omega + K_2 \omega^2 + \dots = \sum_{k=0}^{\infty} K_k \omega^k, \quad (3)

where

Kk=1k!dkHa(ω)dωkω=0,(4)K_k = \frac{1}{k!} \frac{d^k H_\text{a}(\omega)}{d\omega^k} \Bigr\rvert_{\omega=0}, \quad (4)

with K0=Ha(0)K_0 = H_\text{a}(0).

General Squared Magnitude Response

The squared magnitude response Ha(jω)\mathcal{H}_\text{a}(j\omega) is an even function (symmetric with respect to the value axis) so it may be written in a general form as a function of ω2\omega^2, i.e., neglecting the odd powers of ω\omega because they are odd functions. The general form reads

Ha(jω)=d0+d2ω2+d4ω4++d2Mω2Mc0+c2ω2+c4ω4++c2Nω2N.(5)\mathcal{H}_\text{a}(j\omega) = \frac{d_0 + d_2 \omega^2 + d_4 \omega^4 + \dots + d_{2M} \omega^{2M}}{c_0 + c_2 \omega^2 + c_4 \omega^4 + \dots + c_{2N} \omega^{2N}}. \quad (5)

The following observations can already be made with regard to Equation 5.

  • We require that Ha(0)=1\mathcal{H}_\text{a}(0) = 1 so we can readily set c0=d0c_0 = d_0.
  • We require that Ha(j)=0\mathcal{H}_\text{a}(j\infty) = 0 which leads to the conclusion that the denominator must have a greater order than the numerator, i.e., N>MN > M and c2N0c_{2N} \neq 0.

Error Formulation

We can write Ha(jω)\mathcal{H}_\text{a}(j\omega) in terms of the sum of the desired value at 00 and approximation error E(ω)E(\omega)

Ha(jω)=1+E(ω).(6)\mathcal{H}_\text{a}(j\omega) = 1 + E(\omega). \quad (6)

We can insert Equation 6 into Equation 5 and obtain

d0+d2ω2++d2Mω2M=c0+c2ω2++c2Nω2N+E(ω)[c0+c2ω2++c2Nω2N].(7)d_0 + d_2 \omega^2 + \dots + d_{2M} \omega^{2M} = c_0 + c_2 \omega^2 + \dots + c_{2N} \omega^{2N} \newline + E(\omega) [c_0 + c_2 \omega^2 + \dots + c_{2N} \omega^{2N}]. \quad (7)

Error Minimization

To achieve the closest Taylor approximation possible (minimize error E(ω)E(\omega)), we want to eliminate as many low order terms as possible so we set

c0=d0(as before),(8)c_0 = d_0 \text{(as before)}, \quad (8)
c2=d2,(9)c_2 = d_2, \quad (9)
\vdots
c2M=d2M,(10)c_{2M} = d_{2M}, \quad (10)
c2M+2=0,(11)c_{2M+2} = 0, \quad (11)
\vdots
c2N2=0,(12)c_{2N - 2} = 0, \quad (12)
c2N0.(13)c_{2N} \neq 0. \quad (13)

Equations 8-13 tell us that the numerator of Ha(jω)\mathcal{H}_\text{a}(j\omega) from Equation 5 can be arbitrary, because any setting of parameters d0,,d2Md_0, \dots, d_{2M} and subsequent setting of parameters c0,,c2Mc_0, \dots, c_{2M} will yield equally good Taylor approximation.

That allows us to pick the numerator as we wish. In order to have Ha(j)=0\mathcal{H}_\text{a}(j \infty) = 0, we set c0=d0=1c_0 = d_0 = 1 and d2=d4==d2M=0d_2 = d_4 = \dots = d_{2M} = 0.

We, thus, obtain

Ha(jω)=11+c2Nω2N.(14)\mathcal{H}_\text{a}(j \omega) = \frac{1}{1 + c_{2N} \omega^{2N}}. \quad (14)

Determining the Cutoff Frequency

Parameter c2Nc_{2N} determines the analog cutoff frequency ωa\omega_\text{a} so that Ha(ωa)=12\mathcal{H}_\text{a}(\omega_\text{a}) = \frac{1}{2} (-3 dB point).

We have already decided that for us ωa=1\omega_\text{a} = 1 so we need to set

c2N=1.(15)c_{2N} = 1. \quad (15)

Final Squared Magnitude Response

We obtained the final formula for the analog Butterworth low-pass filter of the NN-th order

Ha(jω)=11+ω2N.(16)\mathcal{H}_\text{a}(j\omega) = \frac{1}{1 + \omega^{2N}}. \quad (16)

Equation 16 is the Taylor approximation of the ideal low-pass filter’s squared magnitude response at ω=0\omega = 0. This means that Ha(jω)\mathcal{H}_\text{a}(j\omega) is maximally flat at ω=0\omega = 0.

It turns out that Equation 16 is at the same time the Taylor approximation at ω=\omega = \infty. So Ha(jω)\mathcal{H}_\text{a}(j\omega) is maximally flat at both ends: ω=0\omega=0 and ω=\omega = \infty. That is why, Butterworth filter is said to have maximally flat amplitude response at the endpoints [Smith07].

Transfer Function Derivation

You may wonder:

Since Equation 16 is the squared magnitude response, how do we obtain the transfer function over the ss-domain?

We can use the definition of Ha(s)\mathcal{H}_\text{a}(s) from Equation 1.

Ha(s)=Ha(s)Ha(s)=11+(s2)N,(17)\mathcal{H}_\text{a}(s) = H_\text{a}(s) H_\text{a}(-s) = \frac{1}{1 + (-s^2)^N}, \quad (17)

because if we substitute s=jωs = j\omega, we arrive back at Equation 16.

The fractional on the right side of Equation 17 has exactly 2N2N poles. What are they?

(s2)N+1=0,(18)(-s^2)^N + 1 = 0, \quad (18)
(s2)N=1,(19)(-s^2)^N = -1, \quad (19)
(1)Ns2N=1,(20)(-1)^N s^{2N} = -1, \quad (20)
s2N=(1)N+1,(21)s^{2N} = (-1)^{N+1}, \quad (21)
s2N={1,if N is even,1,if N is odd.(22)s^{2N} = \begin{cases} -1, \quad \text{if } N \text{ is even},\newline 1, \quad \text{if } N \text{ is odd}.\end{cases} \quad (22)
sk={ei(π+2kπ)/2N,if N is even,ei2kπ/2N,if N is odd.,k=(N1),(N2),,0,1,2,,N1,N.(23)s_k = \begin{cases} e^{i(\pi + 2k\pi)/2N}, \quad \text{if } N \text{ is even},\newline e^{i2k\pi/2N}, \quad \text{if } N \text{ is odd.}\end{cases}, \newline \quad k=-(N-1), -(N-2), \dots, 0, 1, 2, \dots , N-1, N. \quad (23)

Note: This indexing of kk is chosen so as to simplify further derivations. Another but equivalent indexing is k=0,1,2,,2N1k=0,1,2,\dots,2N-1.

So Ha(s)\mathcal{H}_\text{a}(s) has 2N2N poles evenly spaced around the unit circle. As their number is even, they are placed symmetrically: NN on the left half-plane, NN on the right half-plane.

Since we want just Ha(s)H_\text{a}(s), we need to factorize Ha(s)\mathcal{H}_\text{a}(s) into Ha(s)H_\text{a}(s) and Ha(s)H_\text{a}(-s).

For Ha(s)H_\text{a}(s) to be stable, we need all of its poles to lie on the left half-plane of the ss-plane, i.e., have negative real parts.

We can obtain it by finding the poles of Ha(s)H_\text{a}(-s) (which lie on the right half-plane) and negating their real parts (because the poles are symmetrical with respect to the imaginary axis).

Poles of Ha(s)H_\text{a}(-s) are sks_k from Equation 23 with arguments in the (π/2,π/2)(-\pi/2,\pi/2) range, i.e., for k=0,±1,±2,,±(N1)/2k=0,\pm 1,\pm 2,\dots,\pm (N-1)/2 if NN is odd, and for k=0,±1,±2,,±(N/21),N/2k=0, \pm 1, \pm 2, \dots, \pm (N/2 -1),-N/2 if NN is even. After negating the real part of these sks_k, we obtain the pole locations of Ha(s)H_\text{a}(s)

skHa(s)={ei(π+2kπ)/2N,k=0,±1,,±(N/21),N/2if N is even,ei2kπ/2N,k=0,±1,±2,,±(N1)/2if N is odd.(24)s_k^{H_\text{a}(s)} = \begin{cases} -e^{-i(\pi + 2k\pi)/2N},k=0, \pm 1, \dots, \pm (N/2 -1),-N/2 \quad \text{if } N \text{ is even},\newline -e^{-i2k\pi/2N}, k=0,\pm 1,\pm 2,\dots,\pm (N-1)/2 \quad \text{if } N \text{ is odd.}\end{cases} \quad (24)

Note: Negating the real part of a complex number is equivalent to negating its complex conjugate.

Knowing all the poles, we can write out Ha(s)H_\text{a}(s) in the factorized version. Let’s consider NN odd first (so k=0,±1,±2,,±(N1)/2k=0,\pm 1,\pm 2,\dots,\pm (N-1)/2).

Ha(s)=k1ssk=k1s+ei2kπ/2N.(25)H_\text{a}(s) = \prod \limits_k \frac{1}{s - s_k} = \prod \limits_k \frac{1}{s + e^{-i 2 k \pi / 2N}}. \quad (25)

We did it! Now let’s just polish this formula.

Tidying Up the Product

The polynomial in the denominator of Equation 17 has real coefficients. Therefore, all roots occur in complex conjugate pairs apart from s0s_0, which is a real number (remember, we consider NN-odd case now).

Since the complex conjugate lies on the same half-plane, we can combine the roots with their conjugates to create a neat-looking real polynomial in the denominator. The single real pole (1-1) must be factored out of the product, because it doesn’t have a conjugate pair.

Ha(s)=1s+1k1(s+ei2kπ/2N)(s+ei2kπ/2N)=1s+1k1s2+2cos(kπ/N)s+1,(26)H_\text{a}(s) = \frac{1}{s+1} \prod \limits_k \frac{1}{(s + e^{-i 2 k \pi / 2N})(s + e^{i 2 k \pi / 2N})} \newline = \frac{1}{s+1} \prod \limits_k \frac{1}{s^2 + 2 \cos (k \pi / N) s + 1}, \quad (26)

with k=1,2,,(N1)/2k=1, 2,\dots, (N-1)/2 (note the lack of nonpositive integers).

For NN even, we obtain analogously (without any real roots)

Ha(s)=m1s2+2cos(mπ/2N)s+1,(27)H_\text{a}(s) = \prod \limits_m \frac{1}{s^2 + 2 \cos (m \pi / 2N) s + 1}, \quad (27)

where m=1,3,,N1m = 1, 3, \dots, N-1. This comes from having m=2k+1,k=0,1,,(N/21)m=2k +1, k=0, 1, \dots, (N/2 -1) to facilitate derivations.

According to [ParksBurrus87], Equations 26 and 27 are very convenient forms for implementation.

We did it! We obtained our analog prototype!

Now, let’s analyze it a little bit. 😉

Butterworth Low-pass Transfer Function

As an example, the low-pass transfer function of the second-order Butterworth low-pass is [Zölzer08]

H2(s)=1s2+2s+1.(28)H_2(s) = \frac{1}{s^2 + \sqrt{2} s + 1}. \quad (28)

The fourth-order Butterworth low-pass has the following transfer function

H4(s)=1(s2+1.848s+1)(s2+0.765s+1).(29)H_4(s) = \frac{1}{(s^2 + 1.848 s + 1)(s^2 + 0.765 s + 1)}. \quad (29)

Hint: To obtain an arbitrary analog cutoff frequency ωa\omega_\text{a}, simply replace ss with s/ωas/\omega_\text{a} in the above transfer functions.

Visualization

To see, how much the Butterworth low-pass filter deviates from the ideal response from Figure 3, let’s plot the amplitude responses of both filters against the ideal response.

Figure 4. Butterworth low-pass amplitude response of 2nd, 4th, and 11th order plotted against the ideal response.

The 11th order is shown for additional comparison.

We can observe that Butterworth filters cross the cutoff frequency with exactly the same gain, which is 1/21/\sqrt{2}. That means that our derivations are correct.

Additionally, we can observe that the higher the filter order, the more steep the slope of the transition band (between the passband and the stopband).

Figure 4 also shows that the Butterworth approximation is indeed maximally flat at frequencies ω=0\omega=0 and ω=\omega=\infty.

Figure 5 shows the same amplitude responses but this time the magnitude is expressed in decibels (20log10()20 \log_{10}(\cdot)). The -3 dB at cutoff frequency is clearly visible.

Figure 5. Butterworth low-pass amplitude response in decibels for 2nd, 4th, and 11th orders plotted against the ideal response.

To obtain the transfer function, I used the Polynomial class from the numpy.polynomial.polynomial module of the NumPy library.

SciPy has a ready-made function to obtain the transfer function of an arbitrary analog filter with one of 5 different design methods. Alternatively, you may check out my implementation for the Butterworth case.

To obtain plots in Figures 4 and 5, I used the scipy.signal.freqs function from the Python SciPy library.

All these SciPy functions have their equivalents in Matlab.

Summary

We did it! We obtained the transfer function of an analog low-pass filter which we can now digitize with the bilinear transform and then transform to the desired form (high-pass, band-pass, etc.).

I put a lot of effort into this article: if you find it useful, please, let me know in the comments!

If you have any questions, I would be happy to answer them in the comments as well.

Thank you for reading! 🙂

Bibliography

[OppenheimSchafer10] Alan V Oppenheim, Ronald W. Schafer, Discrete-Time Signal Processing, 3rd Edition, Pearson 2010.

[ParksBurrus87] T. W. Parks, C. S. Burrus, Digital Filter Design, John Wiley & Sons, Inc., 1987.

[Smith07] Julius O. Smith, Introduction to Digital Filters with Audio Applications, http://ccrma.stanford.edu/~jos/filters/, online book, 2007 edition, accessed November 26, 2021.

[VälimäkiReiss16] Vesa Välimäki, Joshua D. Reiss, All About Audio Equalization: Solutions and Frontiers [PDF], Applied Sciences, Vol. 6, Issue 5, May 6, 2016.

[Zölzer08] Zölzer Udo, Digital Audio Signal Processing, 2nd ed., Helmut Schmidt University, Hamburg, Germany, John Wiley & Sons Ltd, 2008.

Links above may be affiliate links. That means that I may earn a commission if you decide to make a purchase. This does not incur any cost for you. Thank you.