# Using IIR filters

## Contents

# 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.

```
# 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.

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:

Apply the filter once:

`y1 = filter(b, a, x)`

Let

`y1_rev`

be`y1`

in reverse order:`y1_rev = y1[::-1]`

Apply the filter again:

`y2_rev = filter(b, a, y1_rev)`

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.

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.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 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.