It is often the case in signal processing algorithm design that there are linear phase filtering applications wherein we’re interested in designing very high performance (very narrow passband widths, and/or very high attenuation ) nonrecursive FIR filters.
But consider the possibility that you’ve used one of your favorite algoritms to determine that you need, say, a 2000-tap linear phase FIR filter. But when you try to design such a filter using your trusty Remez-Exchange-based (Parks-McClellan) filter design software, you obtain unusable design results.
This occurs because it happens that some software incarnations of the Remez-Exchange algorithm have convergence problems (inaccurate results) when the number of filter taps, or filter order, exceeds four to five hundred.
There’s a slick way around this high-order FIR filter design problem using a frequency-domain zero stuffing technique.
If our FIR filter design software cannot generate FIR coefficient sets whose lengths are in the thousands, then we can design a shorter length set of coefficients and interpolate those coefficients (time-domain impulse response ) to whatever length we desire.
Rather than use time-domain interpolation schemes and account for their inaccuracies, we can simplify the process by performing time-domain interpolation by means of frequency-domain zero stuffing.
An example of the process is as follows: assume that we have a signal sampled at a rate of fs = 1000 Hz. We want a lowpass filter whose cutoff frequency is 20 Hz with 60 dB of stopband attenuation. Compounding the problem are the requirements for linear phase and removal of any DC (zero Hz) component from the signal. (Those last two requirements preclude using the DC-removal schemes discussed elsewhere. )
First, we design a prototype nonrecursive FIR filter having, say, N = 168 coefficients whose desired frequency response magnitude is shown in Figure 13–72(a) below , and its h p (k ) coefficients depicted in Figure 13–72(b). Next, we compute a 168-point DFT of the coefficients to obtain the frequency-domain samples H p (m ) whose magnitudes are shown in Figure 13–72(c).
Figure 13-72. Prototype FIR filter: (a) magnitude response; (b) h p( k ) coefficients; (c) | H p( m )| magnitudes of the 168-point DFT of h p( k ).
Under the assumption that our final desired filter required roughly 1600 taps, we’ll interpolate the h p (k ) prototype impulse response by a factor of M = 10. We perform the interpolation by inserting (M –1)N zeros in the center of the H p (m ) frequency-domain samples, yielding a 1680-point H (m ) frequency-domain sequence whose magnitudes are shown in Figure 13–73(a) below.
Figure 13-73. Desired FIR filter: (a) magnitude of zero-stuffed H p( m ); (b) interpolated h ( k ) coefficients; (c) magnitude of desired frequency response.
Finally, we perform a 1680-point inverse DFT on H (m ) to obtain the interpolated h (k ) impulse response (coefficients), shown in Figure 13–73(b), for our desired filter. (The ten-fold compression of the H p (m ) passband samples results in a ten-fold expansion of the h p (k ) impulse response samples.)
The frequency magnitude response of our final very high-order FIR filter, over the frequency range of –30 -to- 30 Hz, is provided in Figure 13–73(c).
With this process, the prototype filter’s h p (k ) coefficients are preserved within the interpolated filter’s coefficients if the H p (N /2) sample (fs /2) is zero. That condition ensures that H (m ) exhibits conjugate symmetry and forces the h (k ) coefficients to be real-only.
The design steps for this high-order filter design scheme are:
1) With the desired filter requiring MN taps, set the number of prototype filter coefficients, N , to an integer value small enough so your FIR filter design software provides useable results. The integer interpolation factor M equals the number of desired taps divided by N .
2) Design the N -tap prototype FIR filter accounting for the M -fold frequency compression in the final filter. (That is, cutoff frequencies for the prototype filter are M times the desired final cutoff frequencies.)
3) Perform an N -point DFT on the prototype filter’s h p (k ) coefficients to obtain H p (m ).
4) Insert M –1 zero-valued samples just before the H p (N /2) sample of H p (m ) to obtain the new MN -point H (m ) frequency response.
5) Compute the MN -point inverse DFT of H ( m ) yielding an MN -length interpolated h ( k ) coefficient set. (Due to computational errors, discard theimaginary part of h ( k ), making it real-only. )
6) Multiply h ( k ) by M to compensate for the 1/ M amplitude loss induced by interpolation.
7) Test the h ( k ) coefficient set to determine its actual frequency response using standard filter analysis methods. (One method: append thousands of zeros to h ( k ) and perform a very large FFT on the expanded sequence. )
An example application of this filter design is when you’re building a high-performance lowpass polyphase filter where the structures of the high-performance interpolated FIR and frequency sampling lowpass filters don’t permit their decomposition into polyphase sub-filters for such an application.
Used with the permission of the publisher, Prentice Hall, this on-going series of articles on Embedded.com is based on copyrighted material from “Understanding Digital Signal Processing, Second Edition” by Richard G. Lyons. The book can be purchased on line.
Richard Lyons is a consulting systems engineer and lecturer with Besser Associates. As a lecturer with Besser and an instructor for the University of California Santa Cruz Extension, Lyons has delivered digitasl signal processing seminars and training course at technical conferences as well at companies such as Motorola, Freescale, Lockheed Martin, Texas Instruments, Conexant, Northrop Grumman, Lucent, Nokia, Qualcomm, Honeywell, National Semiconductor, General Dynamics and Infinion.