Finding the Frequency Response of an FIR Filter, at Frequencies Other than the Original Sample Points

We would like to know the frequency response of an FIR filter we construct having $N$ filter coefficients by using the IDFT of the desired frequency response. We know it is correct at the sample points, but what about other frequencies. Maybe it is bouncing all over the place there. This is a way to find this information using the computer. There is an analytical way of finding it at every frequency too, but that requires a little more theoretical work to explain.

$H(n) = H(n+N)$ is a sampled version of the desired frequency response. It is sampled every $1/{NT}$ Hz for frequencies from $-1/{2T} < f < 1/{2T}$ to obtain $N$ points. Then the negative frequency samples are moved ahead by one period, or $N$ points, so that, for example, $H(-1)$ is moved to $H(N-1)$, and $H(-2)$ is moved to $H(N-2)$. This leaves us with $\{ H(0), H(1), ..., H(N-1)\}$. We then use the IDFT to obtain the FIR filter coefficients. $$h(k) = \sum_{n=0}^{N-1} H(n)e^{j2\pi nk/N}$$. A python example is shown below.

In [57]:
# Low Pass FIR Filter Design
import numpy as np
import matplotlib.pyplot as plt

T = 0.1
N = 32
f_c = 1.5  # Cutoff of low pass filter

f = np.arange(-1/2/T, 1/2/T,1/N/T)
H = np.piecewise(f, [(-f_c < f)&(f < f_c), (f > f_c)|(f < -f_c)],[1, 0])
plt.plot(f,H,'o')
plt.title('Desired Frequency Response Samples')
plt.xlabel('f (Hz)')
plt.show()

Notice that the desired frequency respones is only specified at the $N$ points as described above.

In [58]:
H = np.fft.fftshift(H)  # Now shift the frequency components for doing the IFFT.
plt.plot(H,'o')
plt.show()
h = np.fft.ifft(H)
t = np.arange(0, N*T, T)
plt.plot(t, h, 'o')
plt.xlabel('time, t (s)')
plt.title('$h(k)$ from the IFFT')
plt.show()

Notice that $h(k)$ as given above is large at first and at the end. Remember that $h(k) = h(k+N)$. If you really want to know what the impulse response looks like you need to shift the last half of the data back $N$ points. This makes it evident that this filter is noncausal, because it's impulse response starts before the impulse occurs. It is also important to make this shift if the data to be filtered is not periodic with period, $N$.

In [59]:
h = np.fft.fftshift(h)  # Comment this line to see what happens
# when you don't shift, but use this to filter a signal x(t) with
# longer duration than NT.

# h = h*np.blackman(N)  # Uncomment this line if you want to see how 
# a better response can be had by using the Blackman window function.
# After uncommenting it you have to re-execute all the python blocks.

ts = np.arange(-N*T/2, N*T/2, T)
plt.plot(ts, h, 'o')
plt.xlabel('time, t (s)')
plt.title('$h(kT)$ Shifted')
plt.show()

It really isn't possible to have a noncausal filter function such as the one you see above. The solution is to delay it by $NT/2$. This delays the output by the same amount, but at least you can do do the filtering. With that delay, the index $k$, runs from zero to $N-1$ so that $0 \le k <N$.

Now we would like to view the frequency spectra at points in between the sampled points to see how it really looks. To do this, recall that the spacing between FFT points is $1/{NT}$. If we added a bunch of zeros before and after the impulse response above, and took the FFT, we would have more resolution. The zeros basically are assuming that the impulse response before and after what you see above is zero, and that is actually exactly right, the way we are filtering the incoming A/D samples, $x(kt)$ where $k$ keeps advancing as long as the song lasts.

In [61]:
M = 100  # Half the number of zeros to add.
h_ext = np.concatenate((np.zeros(M), h), 0)
h_ext = np.concatenate((h_ext, np.zeros(M)), 0)
plt.plot(h_ext,'o', h_ext)
plt.title('Impulse Response with Attached Zeros')
plt.xlabel('Index, $k$')
plt.ylabel('FIR Filter Coefficients')
plt.show()
H_interp = np.fft.fft(h_ext)
f_interp = np.arange(0, 1/T-1/(2*M+N)/T, 1/(2*M+N)/T)
# print(np.shape(f_interp))  # Debugging
H_interp_mag = np.abs(H_interp)
H_interp_angle = np.angle(H_interp)
plt.plot(f_interp, H_interp_mag,'o',f_interp, H_interp_mag)
plt.title('Interpolated Frequency Response')
plt.xlabel('Frequency, $f$ (Hz)')
plt.ylabel('Magnitude of FIR Filter Response')
plt.show()

Note that the frequencies $ 5 < f < 10 $ are really the negative frequency components, $-5 < f < 0$. The periodic nature of the frequency response and the Gibbs phenomena (side lobes of the filter) are very evident. If you plot the response in decibels, you get this.

In [62]:
H_m_db = 20*np.log10(H_interp_mag)
plt.plot(f_interp, H_m_db,'o',f_interp, H_m_db)
plt.title('Interpolated Frequency Response')
plt.xlabel('Frequency, $f$ (Hz)')
plt.ylabel('Magnitude of FIR Filter Response (dB)')
Out[62]:
Text(0, 0.5, 'Magnitude of FIR Filter Response (dB)')

Note that the sidelobes of the filter are really pretty big. To take care of that you need to smooth out the transition at the end of the original filter coefficients. This is done with a windowing function. I leave you readers to try some windowing functions out for yourselves. Another thing that would be instructive is to construct what happens if you didn't delay the impulse response by $NT/2$. You can comment out a line of code above that does the fftshift and run the cells further down to see how terrible the result is.

In [ ]: