FIR Filters ----------- A **finite impulse response filter** is a filter whose impulse response function $h[n]$ is non zero on a bounded subset of $\setR$, i.e. $h[n]=0$ for $|n|>N$. The digital filter then is implemented as a convolution: .. math:: y[n] &= x[n] \ast h[n]\\ &= \sum_{k=-\infty}^{\infty} x[n-k] h[k] In the Z-domain a FIR filter is characterized as: .. math:: H(z) = \sum_{n=-\infty}^{\infty} h[n] z^{-n} Designing a FIR filter then boils down to the choice of the impulse response function $h[n]$. .. rubric:: Simple Filter Consider the following causal FIR filter: .. math:: y[n] = x[n-2] + 2 x[n-1] + x[n] The impulse response is: .. math:: h[n] = (\underline{1}\;2\;1) The DTFT is then given as: .. math:: H(\W) &= 1 + 2 e^{-j\W} + e^{-2j\W}\\ &= e^{-j\W}\left( e^{j\W} + 2 + e^{-j\W}\right)\\ &= e^{-j\W}\left(2+2\cos(\W)\right) We thus see that: .. math:: |H(\W)| &= 2+2\cos(\W)\\ \angle\W &= -\W Now consider a cosine inout signal: .. math:: x[n] = \cos(\W_0 n) The output then is .. math:: y[n] = |H(\W_0)| \cos\left(\W_0 n + \angle H(\W_0)\right) Below we check this in Python: .. exec_python:: FIRexample FIRfilters :linenumbers: :code: shutter :Code_label: Show code for figure :results: hide import numpy as np import matplotlib.pyplot as plt n = np.arange(32) t = np.linspace(0,32,1000) hn = np.array([0,0,1,2,1]) # to make the first 1 the origin when convolving W0 = np.pi/6 xn = np.cos(W0 * n) xt = np.cos(W0 * t) fig, axs = plt.subplots(2) axs[0].plot(t, xt, '--'); axs[0].stem(n, xn, use_line_collection=True); H = lambda W: np.exp(-1j*W) * (2 + 2*np.cos(W)) absW0 = np.abs(H(W0)) print("Gain: %5.2f" % absW0) angleW0 = np.angle(H(W0)) print("Angle: %5.2f" % angleW0) yt_theory = absW0 * np.cos(W0*t + angleW0) yn = np.convolve(xn, hn, mode='same') axs[1].plot(t, yt_theory, '--'); axs[1].stem(n, yn, use_line_collection=True); plt.savefig('source/figures/discretetimeconvolution121.png') .. figure:: /figures/discretetimeconvolution121.png :figwidth: 80% :align: center **FIR convolution filtering of cosine with kernel** $h=\{\underline{1}\;1\;1\}$ Convince yourself using the figure above that theory indeed predicts practice. Except for the first two samples of $y[n]$, why is that? .. rubric:: Running Average Filter Another simple FIR filter is to calculate the average of $K$ consecutive samples, for $K=5$ we get: .. math:: y[n] &= \frac{1}{5}\left( x[n-2] + x[n-1] + x[n] + x[n+1] + x[n+2]\right)\\ &= x[n] \ast h[n] where .. math:: h[n] = (1\;1\;\underline{1}\;1\;1) The discrete time Fourier transform is: .. math:: H(\W) &= \sum_{n=-\infty}^{\infty} h[n] e^{-j \W n}\\ &= \frac{1}{5} \left( e^{j 2 \W} + e^{j \W} + 1 + e^{-j\W} + e^{-j2\W} \right)\\ &= \frac{1}{5} \left( 1 + 2\cos(\W) + 2\cos(2\W) \right) On the interval $\W\in[-\pi,\pi]$ the frequency response looks like: .. exec_python:: fivepointaverage FIRfilters :linenumbers: :code: shutter :Code_label: Show code for figure :results: hide plt.close('all') W = np.linspace(-np.pi, np.pi, 1000) H = 1/5*(1 + 2*np.cos(W) + 2*np.cos(2*W)) plt.plot(W,H) plt.savefig('source/figures/freqresp5pointaverage.png') .. figure:: /figures/freqresp5pointaverage.png :figwidth: 80% :align: center **Frequency response of 5 point average filter** From the plot it is evident that there are 4 zeros in the frequency response. Of course we could set $H\W)=0$ and calculate those zeros. However we do so by analyzing this filter in the Z-domain. .. math:: H(z) &= \sum_{n=-\infty}^{\infty} x[n] z^{-n}\\ &= \frac{1}{5}\left( z^2 + z + 1 + z^{-1} + z^{-2}\right)\\ &= \frac{ \frac{1}{5} \left( 1 + z + z^2 + z^3 + z^4\right)}{z^2} Using the sum formula for a geometric sequence we get .. math:: H(z) &= \frac{\frac{1}{5} \frac{1-z^5}{1-z}}{z^2}\\ &= \frac{1}{5} \frac{1-z^5}{z^2(1-z)} Note that the zero in $z=1$ cancels the pole in $z=1$ so we are left with four zeros (obtained as solutions of $z^5=1$ excluding the solution $z=1$). The four zeros are at $z=\pm\pi/5$ and $z=\pm 2\pi/5$. .. figure:: /figures/poleszeros5taprunningaverage.png **A 5 tap running average filter.** Poles and Zeros plot and Frequency Resppnse plot. .. rubric:: Design in the Frequency Domain In this section we discuss the design based on frequency domain considerations. Suppose we want to approximate an ideal band-pass filter with unit gain for frequencies between $\Omega_\text{low}$ and $\Omega_\text{high}$. Note that due to the symmetry $H(-\Omega)=H(\Omega)$ the mathemetical formulation is: .. math:: H(\Omega) = u(\Omega+\Omega_\text{high}) - u(\Omega+\Omega_\text{low}) + u(\Omega-\Omega_\text{low}) - u(\Omega-\Omega_\text{high}) The inverse Fourier transform then gives: .. math:: h[n] &= \frac{1}{2\pi} \int_{\langle 2\pi \rangle} X(\Omega)e^{j\Omega n}d\Omega\\ &= \frac{1}{2\pi} \left( \int_{\Omega_\text{low}}^{\Omega_\text{high}} e^{j\Omega n} d\Omega + \int_{-\Omega_\text{high}}^{-\Omega_\text{low}} e^{j\Omega n} d\Omega \right)\\ &= \frac{1}{2\pi} \left( \frac{1}{jt} \left[e^{j\Omega n}\right]_{\Omega_\text{low}}^{\Omega_\text{high}} + \frac{1}{jn} \left[e^{j\Omega n}\right]_{-\Omega_\text{high}}^{-\Omega_\text{low}} \right)\\ &= \frac{1}{2\pi} \frac{1}{jn} \left( e^{j\Omega_\text{high}n} - e^{j\Omega_\text{low}n} + e^{-j\Omega_\text{low}n} - e^{-j\Omega_\text{high}n} \right)\\ &= \frac{1}{\pi n}\left( \frac{e^{j\Omega_\text{high}n} - e^{j\Omega_\text{high}n}}{2j} - \frac{e^{j\Omega_\text{low}n} - e^{-j\Omega_\text{low}n}}{2j}\right)\\ &= \frac{\sin(\Omega_\text{high}n)}{\pi n} - \frac{\sin(\Omega_\text{low}n)}{\pi n}\\ &= \frac{\Omega_\text{high}}{\pi}\sinc\left(\frac{\Omega_\text{high}}{\pi} n\right) - \frac{\Omega_\text{low}}{\pi}\sinc\left(\frac{\Omega_\text{low}}{\pi} n\right) Note that the sinc function is defined as $\sinc(t) = \sin(\pi t)/(\pi t)$, this is also the definition used in Numpy. In the code and plot below the formula above is used to generate $h[n]$ for a bandpass filter. The goal is to design a bandpass filter working on discrete signals for which the sampling frequency is 44100 Hz. The band start at low frequency 100 Hz and end at high frequency 4000 Hz. Be sure you understand how to get from frequency $f$ in Hz to radial frequency $\Omega$. Refer to the section on the DTFT. .. exec_python:: bandpassfreqdomain FIRfilters :linenumbers: :code: shutter :Code_label: Show code for figure :results: show plt.close('all') fs = 44100 flow = 100 fhigh = 4000 N = 32 n = np.arange(-N, N+1) Omega_low = 2*np.pi*flow/fs Omega_high = 2*np.pi*fhigh/fs h = Omega_high/np.pi * np.sinc(Omega_high/np.pi * n) - \ Omega_low/np.pi * np.sinc(Omega_low/np.pi * n) plt.stem(n, h, use_line_collection=True); # to check our code we compare it with firwin (with the right set of parameters) numtaps = 2*N+1 from scipy.signal import firwin hfirwin = firwin(numtaps, [flow, fhigh], pass_zero=False, fs=fs, window='boxcar', scale=False) print("Calculated h[n] equal to result from 'firwin' ? ", np.allclose(h, hfirwin)) plt.savefig('source/figures/discrete_bandpass_impulse_response.png') .. figure:: /figures/discrete_bandpass_impulse_response.png **Bandpass Filter Impulse Response.** In the code above there is a correctness check. The impulse response is compared with the impulse response calculated by the function scipy.signal.firwin. It should also be noted that the number of taps (65 in the example above) is too small: the sinc functions are still fluctuating with relative large amplitude for $|n|=32$. The frequency response therefore will not be that for an ideal bandpass filter. .. exec_python:: bandpassfreqdomain FIRfilters :linenumbers: :code: shutter :Code_label: Show code for figure :results: hide plt.close('all') from FIRfilters import test_bandpass test_bandpass() plt.savefig('source/figures/bandpass_sweep.png') .. figure:: /figures/bandpass_sweep.png **Bandpass Filter Response on a Chirp Signal.** In blue the chirp signal is shown. The time sampling is 44100 Hz and the total duration of the signal is 3s, therefore the entire plot is blue for the interval from -1 to 1. The instantaneous frequency in the signal ranges linearly from 1000 to 3000 Hz. Instead of the time the corresponding frequency is shown on the x-axis. In orange the response $x[n]\ast h[n]$ is shown. From the figure it is evident that the filter works as a bandpass filter. Also note that it is not an ideal bandpass filter as we had to truncate the impulse response. The code for this ezample is given below .. literalinclude:: /Python/FIRfilters.py