Finite Impulse Response (FIR) Filters

Digital filters that have an impulse response which reaches zero in a finite number of steps are (appropriately enough) called Finite Impulse Response (FIR) filters. An FIR filter can be implemented non-recursively by convolving its impulse response (which is often used to define an FIR filter) with the time data sequence it is filtering. FIR filters are somewhat simpler than Infinite Impulse Response (IIR) filters, which contain one or more feedback terms and must be implemented with difference equations or some other recursive technique.

IDL's DIGITAL_FILTER function computes the impulse response of an FIR filter based on Kaiser's window, which in turn is based on the modified Bessel function. The Kaiser filter is "nearly optimum in the sense of having the largest energy in the mainlobe for a given peak sidelobe level" [Jackson, Leland B., Digital Filters and Signal Processing]. The DIGITAL_FILTER function constructs lowpass, highpass, bandpass, or bandstop filters.

The figure below plots a bandstop filter which suppresses frequencies between 7 cycles per second and 15 cycles per second for data sampled every 0.02 seconds.

; Sampling period in seconds: 
delt = 0.02 
; Frequencies above f_low will be passed: 
f_low = 15. 
; Frequencies below f_high will be passed: 
f_high = 7. 
; Ripple amplitude will be less than -50 dB: 
a_ripple = 50. 
; The order of the filter: 
nterms = 40 
; Compute the impulse response = the filter coefficients: 
bs_ir_k = DIGITAL_FILTER(f_low*2*delt, f_high*2*delt, $ 
   a_ripple, nterms) 
; The frequency response of the filter is the FFT of its 
; impulse response: 
nfilt = N_ELEMENTS(bs_ir_k) 
; where nfilt = number of points in impulse response. 
; Scale frequency response by number of points: 
bs_fr_k = FFT(bs_ir_k) * nfilt 
; Create a log plot of magnitude in decibels: 
f_filt = FINDGEN(nfilt/2+1) / (nfilt*delt) 
; Magnitude of bandstop filter transfer function: 
mag = ABS(bs_fr_k(0:nfilt/2)) 
PLOT, f_filt, 20*ALOG10(mag), YTITLE='Magnitude in dB', $ 
   XTITLE='Frequency in cycles / second', /XLOG, $ 
   XRANGE=[1.0,1.0/(2.0*delt)], XSTYLE=1, $ 
   TITLE='Frequency Response for Bandstop!CFIR Filter (Kaiser)'

Other FIR filters can be designed based on the Hanning and Hamming windows (see Using Windows), or any other user-defined window function. The design procedure is simple:

The figure below shows the plot using the same sampling period and frequency cutoffs as above, and the corresponding ideal filter is constructed in the frequency domain using the Hanning window.

; Sampling period in seconds: 
delt = 0.02 
; Frequencies above f_low will be passed: 
f_low = 15. 
; Frequencies below f_high will be passed: 
f_high = 7. 
; The length of the filter: 
nfilt = 81 
f_filt = FINDGEN(nfilt/2+1) / (nfilt*delt) 
; Pass frequencies greater than f_low and less than f_high: 
ideal_fr = (f_filt GT f_low) OR (f_filt LT F_high) 
; Convert from byte to floating point: 
ideal_fr = FLOAT(ideal_fr) 
; Replicate to obtain values for negative frequencies: 
ideal_fr = [ideal_fr, REVERSE(ideal_fr[1:*])] 
; Now use an inverse FFT to get the impulse response 
; of the ideal filter: 
ideal_ir = FLOAT(FFT(ideal_fr, /INVERSE)) 
; Ideal_fr is an even function, so the result is real. 
; Scale by the # of points: 
ideal_ir = ideal_ir / nfilt 
; Shift it before applying the window: 
ideal_ir = SHIFT(ideal_ir, nfilt/2) 
; Apply a Hanning window to the shifted ideal impulse response. 
; These are the coefficients of the filter: 
bs_ir_n = ideal_ir*HANNING(nfilt) 
; The frequency response of the filter is the FFT 
; of its impulse response. Scale by the number of points: 
bs_fr_n = FFT(bs_ir_n) * nfilt 
; Create a log plot of magnitude in decibels 
; Magnitude of Hanning bandstop filter transfer function: 
mag = ABS(bs_fr_n(0:nfilt/2)) 
PLOT, f_filt, 20*ALOG10(mag), YTITLE='Magnitude in dB', $ 
   XTITLE='Frequency in cycles / second', /XLOG, $ 
   XRANGE=[1.0,1.0/(2.0*delt)], XSTYLE=1, $ 
   TITLE='Frequency Response for Bandstop!CFIR Filter (Hanning)'