Lab 6: Filter Design

Lab 6: Filter Design

We will turn our attention to filter design in this lab and consider some practical design trade-offs in creating effective filters. In Exercises 1 and 2, we will try out two popular methods of FIR filter design and apply one of them in Exercise 3 to clean up a noisy audio signal. Exercise 4 will explore the benefit of linear phase and Exercise 5 will give us a look at more sophisticated, modern filtering techniques.

FIR and IIR Filters
LSI digital filters can be separated into two types: Finite Impulse Response (FIR) and Infinite Impulse Response (IIR) filters. FIR filters have finitely many terms, only use input signal values, and are always BIBO stable since they have no poles other than at $z=0$. Conversely, IIR filters have infinitely many terms in their impulse response, use both input signal and previous output signal values, and can be BIBO unstable if their ROC does not contain the unit circle in the complex z-domain. In the z-domain, we see that FIR filters have a constant denominator of one in their transfer function, while IIR filters will have some polynomial in their transfer function's denominator according to the system feedback in the LCCDE. We will primarily focus on FIR filter design in this lab.

Consider the elements of a filter's transfer function in the below figure:



There are three key aspects to consider:

Passband ripple - The passband ripple shows how much the frequency resposne varies in the frequencies we would like to keep. Ideally, the passband should be flat so each frequency in the passband is emphasized equally.
Stopband Attenuation - The stopband attenuation tells us how much rejected frequencies will be suppressed. Clearly, we would like the stopband attenuation to be as strong as possible to removed unwanted frequency content.
Transition bandwidth - Transition bandwidth is given by the bandwidth between the passband edge and stopband edge. We would like a narrow transition band to make sure we clearly accept and reject the appropriate frequencies.
Like any other design problem, there is a tradeoff between these parameters with how expensive (how large) the filter is or the method by which we create the filter. Exercise 1 will explore some of these tradeoffs. Now let's consider how we may go about designing an FIR filter.

Linear Phase FIR Filter Design
Consider the problem of designing an ideal low-pass filter. The frequency response of this ideal filter would be

$$ H_d(\omega) = \begin{cases} 1,\quad |\omega| < \omega_c\\ 0,\quad \omega_c \leq |\omega| \leq \pi \end{cases}, $$
where $\omega_c \in [0,\pi]$ is the cut-off frequency of our filter. Note that this ideal frequency response has a constant phase response ($H_d(\omega)$ is real-valued in this case)! Following the DTFT pair of interest, we would obtain the following impulse response for our ideal low-pass filter:

$$ h_d[n] = \begin{cases} \frac{\sin \omega_c n}{\pi n},\quad n\neq 0\\ \frac{\omega_c}{\pi},\quad n = 0\\ \end{cases}. $$
Notice that this impulse response is even symmetric. This even symmetry guarantees that our frequency response will have constant phase. Now, there are two problems with this representation of our filter:

It is of infinite length.
It is non-causal since we have non-zero impulse response values for $n < 0$.
What is one way we can resolve these two issues? With the windowing method!

Windowing Method
We may limit the length of our filter by choosing to keep only the middle $N$ samples and guarantee causality by shifting our impulse response by $M = \frac{N-1}{2}$. This way, our impulse response will be for $0\leq n \leq N-1$. Recall from our DTFT properties that a shift in the time domain will give linear phase in the frequency domain:

$$ x[n-M] \leftrightarrow e^{-jM\omega}X_d(\omega). $$
This "causal delay" $M$ is what gives us our linear phase in FIR filter design since the ideal response had constant phase. If we put all these steps together, we arrive at the following design process for an FIR filter using the windowing method:

Let $D(\omega)$ be your ideal frequency response like $H_d(\omega)$ above.

Apply linear phase to compensate for the necessary causal shift in the time domain to obtain $G(\omega)$, where $M = \frac{N-1}{2}$ and $N$ is the filter length.

$$ G(\omega) = D(\omega)e^{-jM\omega} $$
Take the inverse DTFT of $G(\omega)$ to obtain your shifted, infinite-length filter
$$ g[n] = \textrm{DTFT}^{-1}\{G(\omega)\} = \frac{1}{2\pi}\int_{-\pi}^{\pi}G(\omega)e^{j\omega n}d\omega $$
Apply a window function centered at $n=M$ and of length $N$ to guarantee the causality and finiteness of your filter.
$$ h[n] = g[n]w[n], $$
where $w[n]$ is the window function centered at $n=M$. As we discussed in Lab 5, these windowing functions may take on a number of shapes like the rectangular window, Hamming window, and so on. The choice of window will affect important elements of our frequency response like the passband ripple, stopband attenuation, and transition bandwidth.

Linear and Generalized Linear Phase
Before we get started with the lab activities, we should briefly discuss what identifies linear and generalized linear phase responses. Suppose we have the following frequency response:

$$ H_d(\omega) = R(\omega)e^{j(\alpha - \omega M)}, $$
where $R(\omega)$ is real-valued, $\alpha$ is a real-valued constant, and $M=\frac{N-1}{2}$ is our causal delay. Clearly, the phase response of this DTFT has some kind of linear form since we see a linear function in the complex exponential's exponent. So, how do we decide if this phase is linear or generalized linear? A phase response is purely linear if the constant $\alpha$ is zero. Thus, we have generalized linear phase when $\alpha \neq 0$. Still, in practice, generalized linear phase is typically acceptable since it will offer the important benefits of linear phase. We will revisit the importance of linear phase in Exercise 3.

Now, let's get started designing some filters!

#Import libraries
import numpy as np
import matplotlib.pyplot as plt

from scipy import signal
from skimage.io import imread
from numpy.random import randn
from IPython.display import Audio
from scipy.io import wavfile

#Utility function for dB scaling of magnitude spectra
def sig2db(mag_spec):
    return 20*np.log10(mag_spec)

%matplotlib inline
Exercise 1: Window Method
Suppose we have a generic frequency response for a low-pass filter like described earlier:

$$ H_d(\omega) = \begin{cases} 1,\quad |\omega| < \omega_c\\ 0,\quad \omega_c \leq |\omega| \leq \pi \end{cases}. $$
a. Fill in the function $\textrm{windowed_lpf()}$ that computes the impulse response of a lowpass filter of length $N$ with cutoff frequency $\omega_c$ using the windowing method. (Hint: there is a closed-form expression for this impulse response, try finding it by hand!)

$\textbf{Programming Note}$: The third argument of this function is an optional window function that can be used to window the impulse response. You may specify a different window function by passing a different Python function. For example, we could use a Hamming window by calling the function as $\textrm{windowed_lpf(omega_c, N, window=np.hamming)}$. This way, typing $\textrm{window(N)}$ will return a length $N$ Hamming window. Take note that the default argument would be for a rectangular window.

b. Use your function to generate low-pass filters with cutoff frequency $\omega_c=\frac{\pi}{3}$ and length $N=11,25,$ and $101$ using a rectangular window. Plot the magnitude response of each of the three filters on a dB scale using $\textrm{signal.freqz()}$ and the provided $\textrm{sig2db()}$ function. Comment on the differences between the three filters. Consider elements of the frequency response like passband ripple, transition bandwidth, etc. You may use the default number of points (512) when using $\textrm{signal.freqz()}$.

c. Now use your function to generate low-pass filters with cutoff frequency $\omega_c=\frac{\pi}{3}$ and length $N=25$ for Hamming, Hanning, and Bartlett windows (Hint: use corresponing Numpy functions). Plot the magnitude response of each of the three filters on a dB scale as you did in the previous part. Comment on the differences between the filters including your length 25 filter for a rectangular window in part 1.(b). Please refer to the elements of the frequency response like in part 1.(b).

d. Create a half-band high-pass filter with the following ideal frequency response

$$ H_d(\omega) = \begin{cases} 1,\quad \frac{\pi}{2} < |\omega| \leq \pi\\ 0,\quad |\omega| \leq \frac{\pi}{2} \end{cases} $$
for a rectangular window and $N=25$. Plot the magnitude response of your high-pass filter as you did in parts 1.(b) and 1.(c). Hint: Use your low-pass filter function from part 1.(a) and use the modulation property to your advantage!

#Fill in this function for part 1.a:
def windowed_lpf(omega_c,N,window=np.ones):
    #create ideal impulse response in time domain before windowing
    M = (N-1)/2
    g_n = np.array([1./(np.pi*(n-M))*np.sin(omega_c*(n-M)) if n != M else omega_c/np.pi for n in range(N)])
    #apply window function 
    lpf = g_n * window(N)
    return lpf
 
#Code for part 1.b:

lpf_rect11 = windowed_lpf(np.pi/3, 11, window = np.ones) #initiate all LPFs
lpf_rect25 = windowed_lpf(np.pi/3, 25, window = np.ones)
lpf_rect101 = windowed_lpf(np.pi/3, 101, window = np.ones)

plt.figure(figsize=(30,30))

b = lpf_rect11
a = [1,0]
w,h = signal.freqz(b,a,512) #w = omega/digital frequencies, h = frequency response
plt.subplot(431)
plt.title('Magnitude Response of Low-Pass Filter, Rectangular Window (N=11) ')
plt.plot(w,sig2db(h)) #plot magnitude of frequency response with db-scaling on y-axis
plt.xlabel('$\omega$')
plt.ylabel('Magnitude Response (dB)')

b = lpf_rect25
a = [1,0]
w,h = signal.freqz(b,a,512) #w = omega/digital frequencies, h = frequency response
plt.subplot(432)
plt.title('Magnitude Response of Low-Pass Filter, Rectangular Window (N=25) ')
plt.plot(w,sig2db(h)) #plot magnitude of frequency response with db-scaling on y-axis
plt.xlabel('$\omega$')
plt.ylabel('Magnitude Response (dB)')

b = lpf_rect101
a = [1,0]
w,h = signal.freqz(b,a,512) #w = omega/digital frequencies, h = frequency response
plt.subplot(433)
plt.title('Magnitude Response of Low-Pass Filter, Rectangular Window (N=101) ')
plt.plot(w,sig2db(h)) #plot magnitude of frequency response with db-scaling on y-axis
plt.xlabel('$\omega$')
plt.ylabel('Magnitude Response (dB)')

#Code for part 1.c:
lpf_ham = windowed_lpf(np.pi/3, 25, np.hamming)
lpf_han = windowed_lpf(np.pi/3, 25, np.hanning)
lpf_bart = windowed_lpf(np.pi/3, 25, np.bartlett)

b = lpf_ham
a = [1,0]
w,h = signal.freqz(b,a,512) #w = omega/digital frequencies, h = frequency response
plt.subplot(434)
plt.title('Magnitude Response of Low-Pass Filter, Hamming Window (N=25) ')
plt.plot(w,sig2db(h)) #plot magnitude of frequency response with db-scaling on y-axis
plt.xlabel('$\omega$')
plt.ylabel('Magnitude Response (dB)')

b = lpf_han
a = [1,0]
w,h = signal.freqz(b,a,512) #w = omega/digital frequencies, h = frequency response
plt.subplot(435)
plt.title('Magnitude Response of Low-Pass Filter, Hanning Window (N=25) ')
plt.plot(w,sig2db(h)) #plot magnitude of frequency response with db-scaling on y-axis
plt.xlabel('$\omega$')
plt.ylabel('Magnitude Response (dB)')

b = lpf_bart
a = [1,0]
w,h = signal.freqz(b,a,512) #w = omega/digital frequencies, h = frequency response
plt.subplot(436)
plt.title('Magnitude Response of Low-Pass Filter, Bartlett Window (N=25) ')
plt.plot(w,sig2db(h)) #plot magnitude of frequency response with db-scaling on y-axis
plt.xlabel('$\omega$')
plt.ylabel('Magnitude Response (dB)')

#Code for 1.d: 

def windowed_hpf(omega_c,N,window = np.ones):
    M = (N-1)/2
    g_n = np.array([1./(np.pi*(n-M))*np.sin(omega_c*(n-M))*np.exp(-1j*np.pi*n) if n != M else omega_c/np.pi for n in range(N)])
    hpf = g_n * window(N)
    return hpf

hpf_rect25 = windowed_hpf(np.pi/2, 25, window = np.ones)
b = hpf_rect25
a = [1,0]
w,h = signal.freqz(b,a,512) #w = omega/digital frequencies, h = frequency response
plt.subplot(437)
plt.title('Magnitude Response of High-Pass Filter, Rectangular Window (N=25) ')
plt.plot(w,sig2db(h)) #plot magnitude of frequency response with db-scaling on y-axis
plt.xlabel('$\omega$')
plt.ylabel('Magnitude Response (dB)')

hpf_rect25 = windowed_hpf(np.pi/2, 25, np.hamming)
b = hpf_rect25
a = [1,0]
w,h = signal.freqz(b,a,512) #w = omega/digital frequencies, h = frequency response
plt.subplot(438)
plt.title('Magnitude Response of High-Pass Filter, Hamming Window (N=25) ')
plt.plot(w,sig2db(h)) #plot magnitude of frequency response with db-scaling on y-axis
plt.xlabel('$\omega$')
plt.ylabel('Magnitude Response (dB)')

hpf_rect25 = windowed_hpf(np.pi/2, 25, np.hanning)
b = hpf_rect25
a = [1,0]
w,h = signal.freqz(b,a,512) #w = omega/digital frequencies, h = frequency response
plt.subplot(439)
plt.title('Magnitude Response of High-Pass Filter, Hanning Window (N=25) ')
plt.plot(w,sig2db(h)) #plot magnitude of frequency response with db-scaling on y-axis
plt.xlabel('$\omega$')
plt.ylabel('Magnitude Response (dB)')

hpf_rect25 = windowed_hpf(np.pi/2, 25, np.bartlett)
b = hpf_rect25
a = [1,0]
w,h = signal.freqz(b,a,512) #w = omega/digital frequencies, h = frequency response
plt.subplot(4,3,11)
plt.title('Magnitude Response of High-Pass Filter, Bartlett Window (N=25) ')
plt.plot(w,sig2db(h)) #plot magnitude of frequency response with db-scaling on y-axis
plt.xlabel('$\omega$')
plt.ylabel('Magnitude Response (dB)')
/Users/ryanyoseph/anaconda3/lib/python3.7/site-packages/numpy/core/_asarray.py:85: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)
/Users/ryanyoseph/anaconda3/lib/python3.7/site-packages/numpy/core/_asarray.py:85: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)
/Users/ryanyoseph/anaconda3/lib/python3.7/site-packages/numpy/core/_asarray.py:85: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)
/Users/ryanyoseph/anaconda3/lib/python3.7/site-packages/numpy/core/_asarray.py:85: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)
/Users/ryanyoseph/anaconda3/lib/python3.7/site-packages/numpy/core/_asarray.py:85: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)
/Users/ryanyoseph/anaconda3/lib/python3.7/site-packages/numpy/core/_asarray.py:85: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)
/Users/ryanyoseph/anaconda3/lib/python3.7/site-packages/numpy/core/_asarray.py:85: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)
/Users/ryanyoseph/anaconda3/lib/python3.7/site-packages/numpy/core/_asarray.py:85: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)
/Users/ryanyoseph/anaconda3/lib/python3.7/site-packages/numpy/core/_asarray.py:85: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)
/Users/ryanyoseph/anaconda3/lib/python3.7/site-packages/numpy/core/_asarray.py:85: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)
Text(0, 0.5, 'Magnitude Response (dB)')