FIR Filtering Demonstration

Some initial setup

The code immediately below just sets up a few things that will be used for the demo, including generating the input signal that we will be working with. These lines do not directly relate to the filtering aspect of the demonstration.

In [1]:
from IPython.display import Audio
%run DemoSignals.ipynb
rcParams['figure.figsize'] = 10, 8

Foreword

This demonstration presents a very contrived scenario that is unlikely to crop up in the real world. I would have demonstrated a more relevant application of filtering (perhaps selecting a radio station, for example), but I instead chose to work with audio signals for a few reasons: - We can directly observe these signals using our sense of hearing - Audio deals with real-valued samples - Audio is already at baseband; no need to worry about mixing down higher frequency signals - Allows us to concentrate on the FIR filtering instead of getting hung up on these other aspects - (...and audio is of personal interest to me)

Despite the contrived nature of this demonstration, there are still analogies that can be drawn to more common filtering applications.

Problem Statement

Given an audio signal containing three melodies of non-overlapping frequencies, use FIR filtering to isolate each of the melodies.

Knowns

  • Sampling rate = 8kHz
  • Melody bands:
  • Channel 1: ~146Hz - ~785Hz
  • Channel 2: ~932Hz - ~1569Hz
  • Channel 3: ~1864Hz - ~3730Hz

First, we will look at using the Window Design Method to generate low pass, bandpass, and high pass filters to filter each of these signals. Then, we will take a quick look at using the Parks-McClellan method to filter the signals.

The Input Signal

In [2]:
fs = 8000
specgram(input_signal, Fs=fs)
Audio(input_signal, rate=fs)
Out[2]:

Window Design Method

Low Pass Filter

First we'll take a look at creating a low pass filter to pass only the "channel 1" audio.

From above the above list of knowns (and the spectrogram), we see that channel 1 contains frequencies up to around \(785\text{Hz}\). Channel 2 starts at around \(932\text{Hz}\). This leaves a very narrow transition band of \(932-785=147\text{Hz}\), or around \(0.018375f_s\). We will need a steep rolloff to accomplish this.

The first step to the window design method is to design an ideal "brick wall" frequency response in the frequency domain. What size FFT should we start with?

  • Consider a 64-point FFT
    • Each frequency bin would be \(8000/64=125\text{Hz}\)
    • Frequency bins around the transition region include \(\{\dots, 750, 875, 1000, \dots\}\)
    • If we made an ideal passband that includes the \(875\text{Hz}\) bin, we would ensure that the desired \(785\text{Hz}\) signal is in the passband, but since the stopband doesn't start until $1000, this means that channel 2's signal is not completely in the stopband
  • Consider a 128-point FFT
    • Each frequency bin would be \(8000/128=62.5\text{Hz}\)
    • Frequency bins around the transition region include \(\{\dots, 750, 812.5, 875, 937.5, 1000, \dots\}\)
    • Great! We'll make our ideal passband extend up through the \(812.5\text{Hz}\) bin and begin our ideal stopband at \(875\text{Hz}\). This unambiguously places all of the channel 1 data in the ideal passband and all of the channel 2 data in the ideal stopband
    • To create the ideal filter response, then, set \(812.5/62.5=13\) bins on either side of DC to 1, and the rest of the bins to 0

Let's give this a try...

In [3]:
N = 128
response = zeros(N)
response[range(-13, 14)] = 1

stem(response);

Now to get the filter coefficients, take the inverse FFT of the ideal filter response above

In [4]:
coefficients = real(fft.ifft(response)) # Numerical error results in a tiny imaginary part while taking IFFT; discard it
stem(coefficients);

We're not quite there. We're interested in the coefficients centered around time \(t=0\). Since the IFFT considers the time domain periodic, we simply wrap the signal (see Figure 5-18 in the text).

In [5]:
coefficients = roll(coefficients, N/2-1)
stem(coefficients);

We also want to truncate the coefficients so that they're symmetric. Lop off the unmatched 128th coefficient:

In [6]:
coefficients = coefficients[:-1]
stem(coefficients);

We've got our filter coefficients. Now we just need to convolve them with the input signal to get our filtered output. From Equation 5-6 in the text:

\[ y(n)=\sum_{k=0}^{M-1} h(k)x(n-k) \]

Here's a straightforward implementation of the convolution operation (of course, we could just use numpy's convolve() function, or do this in some other more "Pythonic" manner, but let's be explicit for now):

In [7]:
def my_convolution(x, h):
    """ Naive convolution implementation
    
    Args:
        x: the input signal
        h: the filter coefficients        
    Returns:
        the convolution of x and h
        
    """
    M = len(h)
    y = []
    for n in arange(0, M + len(x) - 1):
        y.append(0)
        for k in arange(0, M):
            if (k < len(h) and (n -k) < len(x)):
                y[n] += h[k] * x[n - k]
    
    return y

Cross our fingers and try out our filter:

In [8]:
output_signal = my_convolution(input_signal, coefficients)

specgram(output_signal, Fs=fs)
Audio(output_signal, rate=fs)
Out[8]:

Results

That filter performed just "okay". We can see from the spectrogram (and we can hear if we listen carefully) that the other channels are still apparent in the signal.

To understand why, let's take a look at the frequency response of the filter that we created. Remember, we can obtain the frequency response of an FIR filter simply by taking the FFT of the filter coefficients.

In [9]:
pad_to = 2048
freq_axis = linspace(0, fs, pad_to, False)

plot(freq_axis, 20 * log10(abs(fft.fft(coefficients, pad_to))));

Windowing

The frequency response above shows us that although we achieved a quick rolloff, the side lobes in the stop band are still relatively high.

Let's see what happens if we apply a windowing function to our filter. We'll use the Blackman window because of it's good side lobe attenuation.

In [10]:
w = blackman(len(coefficients))
coefficients_windowed = w * coefficients

output_signal = my_convolution(input_signal, coefficients_windowed)
specgram(output_signal, Fs=fs)
Audio(output_signal, rate=fs)
Out[10]:

This looks better. We did a great job at rejecting most of the signal from the other channels, but the lowest tones from channel 2 come through pretty strong. What happened?

Let's compare the frequency response with (green) and without (blue) the windowing function.

In [11]:
# Compare frequency response of non-filtered and filtered window
subplot(211)
plot(freq_axis, 20 * log10(abs(fft.fft(coefficients, pad_to))))
plot(freq_axis, 20 * log10(abs(fft.fft(coefficients_windowed, pad_to))))

# Zoom in on region around end of channel 1 and start of channel 2
subplot(212)
plot(freq_axis, 20 * log10(abs(fft.fft(coefficients, pad_to))))
plot(freq_axis, 20 * log10(abs(fft.fft(coefficients_windowed, pad_to))))
xlim([700, 1000])
ylim([-100, 10]);

The Blackman window provides much better stopband rejection, but it has a much slower rolloff. The 932Hz tone in channel 2 is barely attenuated.

What can we do? We combatted the sideband attenuation by using a windowing function. We can get a steeper rolloff by increasing the number of filter coefficients.

For example, let's check out the effect of doubling our filter length.

In [12]:
# Repeat everything we did above, but using a 256-point FFT this time
# No need to plot every intermediate step along the way this time
N = 256
response = zeros(N)
response[range(-25, 26)] = 1

coefficients = real(fft.ifft(response)) # Numerical error results in a tiny imaginary part while taking IFFT; discard it
coefficients = roll(coefficients, N/2-1)
coefficients = coefficients[:-1]

w = blackman(len(coefficients))
coefficients_windowed = w * coefficients

output_signal = my_convolution(input_signal, coefficients_windowed)
specgram(output_signal, Fs=fs)
Audio(output_signal, rate=fs)
Out[12]:

This is what I call "good enough". Let's take a look at the frequency response of the new filter.

In [13]:
subplot(211)
plot(freq_axis, 20 * log10(abs(fft.fft(coefficients_windowed, pad_to))))

subplot(212)
plot(freq_axis, 20 * log10(abs(fft.fft(coefficients_windowed, pad_to))))
xlim([700, 1000])
ylim([-100, 10]);

Now we have good rolloff and good stopband attenuation, but it comes at the cost of a very long filter (computationally expensive) and a high latency/group delay. Of course, this was just one possible filter. We very well could have experimented with other windows and number of filter taps and come up more computationally efficient filter that still meets our needs.

Bandpass Filter

Now that we've demonstrated our ability to isolate the "channel 1" signal via a low-pass filter, let's design a band-pass filter to isolate channel 2. Recall, channel 2 occupies around \(932\text{Hz}\) to \(1569\text{Hz}\).

Bandpass filter design begins with low-pass filter design. Thus, we first need to develop a low pass filter with a passband of \(1569-932=637\text{Hz}\) centered around DC. The closest neighbor to channel 2 is channel 1, which, again, is only \(932-785=147\text{Hz}\) away. Thus, for the same reasons as the first filter, we'll need a high-order filter to get the steep cutoff required to isolate these two channels. We know from above that a 256-point FFT yielded good results for this narrow transition band, so we'll use that again here.

Our passband this time must span \(637\text{Hz}\), which correponds to 11 bins on either side of DC.

In [14]:
N = 256
response = zeros(N)
response[range(-11, 12)] = 1

stem(response);
In [15]:
coefficients = real(fft.ifft(response)) # Numerical error results in a tiny imaginary part while taking IFFT; discard it
coefficients = roll(coefficients, N/2-1)
coefficients = coefficients[:-1]

w = blackman(len(coefficients))
coefficients_windowed = w * coefficients

Our windowed filter coefficients for a low-pass filter with a \(637\text{Hz}\) passband look like this:

In [16]:
stem(coefficients_windowed);

Which gives us a low pass filter with the following frequency response:

In [17]:
# Frequency response over entire spectrum
subplot(211)
plot(freq_axis, 20 * log10(abs(fft.fft(coefficients_windowed, pad_to))))

# Zoomed in around transition band
subplot(212)
plot(freq_axis, 20 * log10(abs(fft.fft(coefficients_windowed, pad_to))))
xlim([0, 500])
ylim([-100, 10]);

Now, here's the part where we transform this low-pass filter into a bandpass filter.

From Equation 5-20 in the book:

\[ h_{bp} \left(k\right) = h_{lp}\left(k\right) * s_{shift}\left(k\right) \]

where, \(h_{bp}\) are our desired bandpass filter coefficients, \(h_{lp}\) are our low-pass filter coefficients, and \(s_{shift}\) is a sinusoidal sequence having our desired bandcenter frequency.

We want our bandpass filter centered at the middle of channel 2, or \(\frac{932 + 1569}{2}=1250.5\text{Hz}\). Our \(s_{shift}\) sequence is therefore:

\[ s_{shift}\left(k\right) = \cos{\left(2 \pi k \frac{1250.5}{8000} \right)} \]

In code:

In [18]:
s_shift = cos(2 * pi * arange(len(coefficients_windowed)) * 1250.5 / fs)
h_bp = coefficients_windowed * s_shift

# Plot the filter coefficients
subplot(211)
stem(h_bp)

# Plot the frequency response
subplot(212)
plot(freq_axis, 20 * log10(abs(fft.fft(h_bp, pad_to))));

Result

Those coefficients look a little bit strange, but let's see how they fare.

In [19]:
output_signal = my_convolution(input_signal, h_bp)

specgram(output_signal, Fs=fs)
Audio(output_signal, rate=fs)
Out[19]:

Looks good! Now on to the high pass filter.

High Pass Filter

The last thing that we want to try with the Window design method approach is to design a high pass filter to isolate channel 3. The process is very similar to that of the bandpass filter, except that we will shift the band center all the way up to \(f_s/2\).

Because there is a wider gap (about \(1864-1569=295\text{Hz}\)) between channel 3 and channel 2 compared to the other channels, we do not need as steep a rolloff for this filter. Let's try going back to a 128-point FFT for this one.

A high pass filter for Channel 3 would extend from \(f_s/2=4000\text{Hz}\) down to around \(1864\text{Hz}\), so we need a filter that passes \(4000-1864=2136\text{Hz}\) to the left of \(f_s/2\). We therefore start by designing a low-pass filter that has a passband extending \(2136\text{Hz}\) from DC. This corresponds to 35 bins on either side of DC for a 128-point FFT.

Let's go through the motions for the low-pass filter one last time.

In [20]:
N = 128
response = zeros(N)
response[range(-35, 36)] = 1

coefficients = real(fft.ifft(response)) # Numerical error results in a tiny imaginary part while taking IFFT; discard it
coefficients = roll(coefficients, N/2-1)
coefficients = coefficients[:-1]

w = blackman(len(coefficients))
coefficients_windowed = w * coefficients

Now we shift the filter up to \(f_s/2\) to turn it into a high pass filter by multiplying the filter coefficients by a sinusoid, \(s_{shift}\) of frequency \(f_s/2\).

Notice that this can be simplified slightly:

\[ h_{hp}\left(k\right) = h_{lp}\left(k\right) * s_{shift}\left(k\right) \] where \[ \begin{align} s_{shift}\left(k\right) &= \cos\left(2 \pi k \frac{f_s/2}{f_s}\right) \\ &= \cos\left(2 \pi k \frac{1}{2}\right) \\ &= \cos\left(\pi k\right) \\ &= -1^k, k \in \mathbb{Z} \end{align} \] hence \[ \begin{align} h_{hp}\left(k\right) &= h_{lp}\left(k\right) * -1^k \\ &= h_{lp} * (1, -1, 1, -1,\dots) \end{align} \]

In [21]:
s_shift = (-1) ** arange(len(coefficients_windowed))
h_hp = coefficients_windowed * s_shift

# Plot the filter coefficients
subplot(211)
stem(h_hp)

# Plot the frequency response
subplot(212)
plot(freq_axis, 20 * log10(abs(fft.fft(h_hp, pad_to))));

Finally, let's see how it performs

In [22]:
output_signal = my_convolution(input_signal, h_hp)

specgram(output_signal, Fs=fs)
Audio(output_signal, rate=fs)
Out[22]:

That's the window design method. Now, let's see what we can do with the...

Parks-McClellan Method

Take the easy way out.

Feel free to see how changing the number of filter taps affects the filter's response.

Channel 1

In [ ]:
import scipy.signal

# The desired number of filter taps
n_taps = 128

# Specify the frequencies for pass, transition, and stop bands
bands = [0, 0, 146, 785, 932, 4000]

# Desired gains for the frequency bands between frequency pairs above
desired = [0, 1, 0]

# Leverage Parks-McClellan
h = scipy.signal.remez(n_taps, bands, desired, Hz=fs)

# Plot the filter coefficients
subplot(211)
stem(h_hp)

# Plot the frequency response
subplot(212)
plot(freq_axis, 20 * log10(abs(fft.fft(h, pad_to))));

-c:21: RuntimeWarning: divide by zero encountered in log10
In [ ]:
output_signal = my_convolution(input_signal, h)

specgram(output_signal, Fs=fs)
Audio(output_signal, rate=fs)

Channel 2

In [ ]:
# The desired number of filter taps
n_taps = 128

# Specify the frequencies for pass, transition, and stop bands
bands = [0, 785, 932, 1569, 1864, 4000]

# Desired gains for the frequency bands between frequency pairs above
desired = [0, 1, 0]

# Leverage Parks-McClellan
h = scipy.signal.remez(n_taps, bands, desired, Hz=fs)

# Plot the filter coefficients
subplot(211)
stem(h_hp)

# Plot the frequency response
subplot(212)
plot(freq_axis, 20 * log10(abs(fft.fft(h, pad_to))));

Notice that the algorithm is allowed to do whatever it wants in the transition bands

In [ ]:
output_signal = my_convolution(input_signal, h)

specgram(output_signal, Fs=fs)
Audio(output_signal, rate=fs)

Channel 3

In [ ]:
# The desired number of filter taps
n_taps = 128

# Specify the frequencies for pass, transition, and stop bands
bands = [0, 1569, 1864, 3730, 4000, 4000]

# Desired gains for the frequency bands between frequency pairs above
desired = [0, 1, 0]

# Leverage Parks-McClellan
h = scipy.signal.remez(n_taps, bands, desired, Hz=fs)

# Plot the filter coefficients
subplot(211)
stem(h_hp)

# Plot the frequency response
subplot(212)
plot(freq_axis, 20 * log10(abs(fft.fft(h, pad_to))));
In [ ]:
output_signal = my_convolution(input_signal, h)

specgram(output_signal, Fs=fs)
Audio(output_signal, rate=fs)