11.2. Using IIR filters

In the previous section, we saw the general definition of IIR filters. In this section, we’ll see how to use them in practice, and begin to probe at understanding their behavior by testing their outputs in response to synthetic input signals.

As in the previous chapter, we’ll focus on low-pass filters, but the ideas will readily generalize to other settings.

11.2.1. IIR filters in Python

In the previous section, we defined a function filter that implements (11.2). Given an input signal \(x\) and filter coefficients b (feed-forward) and a (feed-back), we construct the output signal \(y\) one sample at a time, using nested loops in a manner not too different from (time-domain) convolution.

In practice, we would do better to use an existing implementation and not write our own from scratch. Luckily, the scipy.signal package provides exactly this function: lfilter (linear filter).

For example, the exponential moving average filter could be implemented as follows:

import scipy.signal

b = [1/2]  # Feed-forward
a = [1, -1/2]  # Feed-back

# Apply the filter to input signal x
y = scipy.signal.lfilter(b, a, x)

We can use this, for example, to compute the first few samples of the impulse response:

# Compute the first few samples of the impulse response
# of an exponential moving average

x = np.array([1, 0, 0, 0, 0, 0, 0])
y = scipy.signal.lfilter([1/2], [1, -1/2], x)
print(y)

[0.5       0.25      0.125     0.0625    0.03125   0.015625  0.0078125]

11.2.2. Example: Butterworth filters

The first type of filter that we’ll look at is called the Butterworth filter, after Stephen Butterworth [B+30]. We’ll not get into the details of how the filter coefficients are defined, but instead rely on the scipy.signal.butter function to construct them for us.

Using scipy.signal.butter is not too different from using the window-method function firwin to design an FIR low-pass filter. We’ll need to supply the order of the filter (the number of coefficients), as well as the cutoff frequency \(f_c\) and the sampling rate \(f_s\). The example code below constructs an order-10 filter with \(f_c = 500\) Hz and \(f_S = 44100\).

fs = 44100  # Sampling rate
fc = 500  # Cutoff frequency

# Build the low-pass filter
b, a = scipy.signal.butter(10, fc, fs=fs)

# Print the coefficients
print('Feed-forward: ', b)
print('Feed-back:    ', a)

Feed-forward:  [2.62867578e-15 2.62867578e-14 1.18290410e-13 3.15441093e-13
 5.52021913e-13 6.62426296e-13 5.52021913e-13 3.15441093e-13
 1.18290410e-13 2.62867578e-14 2.62867578e-15]
Feed-back:     [   1.           -9.54462136   41.00487909 -104.41737082  174.53480697
 -200.09486144  159.34094444  -87.02835847   31.2004603    -6.6300023
    0.6341236 ]

To demonstrate the filter’s behavior, let’s apply it to a 250-sample delayed impulse, padded with zeros out to 1000 samples. There’s nothing special about the particular delay that we’re using, but it will allow us to compare the magnitude and phase spectra before and after applying the filter. Remember that impulses (including delayed impulses) have energy at all frequencies, so this should give us a sense of how well the filter works at attenuating high frequencies. By default, the lfilter function produces an output \(y\) of the same length as the input, which might not capture all of the response behavior that we want to look at. Padding the input with trailing zeros gives us time to observe more of the response without truncating prematurely.

Warning

In this section, we’re investigating the behavior of an IIR filter by analyzing a finite bit of its output in response to a particular input signal. This should not be confused with analyzing the filter directly. We’ll see how to do that in the next chapter.

# Create the delayed impulse signal
x = np.zeros(1000)
x[249] = 1

# Apply the Butterworth filter from above
y = scipy.signal.lfilter(b, a, x)

We can now inspect the DFT of the input \(x\) and output \(y\) to see what the filter has done: the results are illustrated in Fig. 11.3.

Fig. 11.3 Top-left: a 250-sample delay signal \(\blue{x}\), padded with zeros to 1000 samples. Top-center: the magnitude spectrum \(|\darkblue{X}|\), measured in decibels, zoomed into the frequency range \([0, 1000]\). Top-right: the unwrapped phase spectrum \(\phi'\) of \(\blue{X}\). Bottom-left: the output \(\purple{y}\) after applying an order-10 low-pass Butterworth filter with cutoff \(f_c=500\). Bottom-center: the magnitude spectrum \(|\magenta{Y}|\), with stop-band marked by the shaded region. Bottom-right: the unwrapped phase spectrum \(\phi'\) of \(\magenta{Y}\).

Fig. 11.4 A more careful look at the group delay of \(\purple{y}\) within the passband shows that the delay is not constant over all frequencies, especially near the cutoff \(f_c=500\).

We can immediately observe a couple of interesting things in Fig. 11.3, especially when compared to the FIR filters demonstrated in Fig. 10.7 and Fig. 10.8. First, the time-domain output \(y\) is delayed relative to the input—it peaks much later—with asymmetric ripples that decay to 0.

Second, even though the filter has only order 10, its stop-band attenuation is comparable to the FIR filters, which have much higher order (in the hundreds).

Third, the phase response looks approximately linear in the pass-band, but it turns out to not be exactly linear. As a result, different frequencies present in \(x\) have been delayed by slightly different amounts to create \(y\), which we perhaps could have anticipated from the fact that \(x\) is symmetric in time (after padding) but \(y\) is asymmetric.

These observations point to general features of IIR filters (not just Butterworth): they can be much more efficient than FIR filters, but they rarely have linear phase. Applying them naively can lead to phase distortion artifacts.

11.2.3. Selecting the order

In the example above, the choice of order 10 might seem arbitrary, but it was in fact chosen to satisfy certain criteria: no more than 3dB pass-band ripple, 60 dB of stop-band attenuation, and a transition band from 500 to 1000 Hz. The scipy helper function scipy.signal.buttord takes these parameters, and provides the minimal filter order which satisfies the constraints:

fc = 500  # Cutoff at 500 Hz
fstop = 1000  # End the transition band at 1000 Hz
passband_ripple = 3  # 3dB ripple
stopband_atten = 60  # 60 dB attenuation of the stop-band

# Get the smallest order for the filter
order = scipy.signal.buttord(fc, fstop, passband_ripple, stopband_atten, fs=fs)

# Now construct the filter
b, a = scipy.signal.butter(order, fc, fs=fs)

If you recall how The Parks-McClellan method optimizes the filter coefficients subject to constraints on ripple, attenuation, and transition bandwidth, the idea is similar here (although the underlying optimization is quite different). In general, when using IIR filters, it’s a good idea to look for helper functions (like order selection) which can help you determine the best settings to satisfy your constraints.

11.2.4. Compensating for phase distortion

The example above illustrates two things: 1) that IIR filters induce delay, just like FIR filters, and 2) IIR filters are often not linear-phase, resulting in phase distortion (see margin Fig. 11.4).

There is a nice trick for dealing with both of these issues: filter the signal twice, but in opposite directions. More specifically, we do the following:

1. Apply the filter once: y1 = filter(b, a, x)

2. Let y1_rev be y1 in reverse order: y1_rev = y1[::-1]

3. Apply the filter again: y2_rev = filter(b, a, y1_rev)

4. Reverse again to obtain the final output: y = y2_rev[::-1]

The key idea here is that if filtering the signal once (forward in time) induces some delay, then filtering it again in the opposite direction (backward in time) should undo that delay exactly. This even works if the delay is different for each frequency, so even if a filter has non-linear phase, we can undo any distortion artifacts.

There are two caveats to be aware of.

Tip

If you plan to use double-filtering, design your filters to have half the target gain. That is, if you want your total stop-band attenuation to be 80dB, design your filter to achieve 40dB.

1. Any gain that would be applied by the filter will be applied twice. If a frequency is attenuated by 40 dB in one pass, it will be attenuated by 80 dB in two passes. This also goes for gain and pass-band ripple, not just stop-band attenuation.

2. Double-filtering can only be applied if we can observe the entire input signal in advance: it is non-causal. This is fine for pre-recorded signals, but cannot work for real-time audio streams.

Caveats aside, this technique is so common that it is also provided by most signal processing frameworks. In scipy, it is called scipy.signal.filtfilt, and we can use it just like we would use lfilter. The example below illustrates how to apply this method to our example signal and low-pass filter, with results illustrated in Fig. 11.5.

fs = 44100  # Sampling rate
fc = 500  # Cutoff frequency

# Build the low-pass filter for double-filtering
# This means we'll cut the pass-band ripple and stop-band attenuation both in half
# We only need the first output from buttord
order, _ = scipy.signal.buttord(fc, 1000, 3/2, 60/2, fs=fs)

b2, a2 = scipy.signal.butter(order, fc, fs=fs)

# Apply the double filter
y2 = scipy.signal.filtfilt(b2, a2, x)

Fig. 11.5 Top: the input signal \(\blue{x}\) (an impulse delayed by 250 samples), along with the single-pass output \(\purple{y}\) (order=10) and the double-pass output \(\magenta{y_2}\) (order=6). Middle: the DFT magnitudes of \(\darkblue{X}\) (constant), \(\purple{Y}\), and \(\magenta{Y_2}\). Bottom: the unwrapped phase of \(\darkblue{X}\), \(\purple{Y}\), and \(\magenta{Y_2}\). Note that the phase of \(\magenta{Y_2}\) agrees with \(\darkblue{X}\) over the pass-band, indicating a total delay of 0.

Fig. 11.5 illustrates that the double-filtering approach does effectively align the final output \(y_2\) to the input \(x\): their time-domain representations peak in exactly the same position. Note, however, that reverse-filtering does introduce pre-echo: the oscillations in \(y_2\) preceding the peak at \(n=249\). In this respect, the result is similar to the FIR filter outputs in Fig. 10.8, and is to some extent unavoidable if we’re applying low-pass filtering to signals with strong transients. Remember: although the outputs are similar, we got there with much less work: an order-6 IIR filter applied twice, compares well to FIR filters with order in the hundreds.

Butterworth filters are just the beginning: in the next section, we’ll meet a few more types of IIR filters.