By Don Orofino, MathWorks and Paul Pacheco, MathWorks
Engineers commonly face the challenge of designing filters with narrow passbands and transition bands. As we try to pack more information into a finite allocation of the frequency spectrum, very tight constraints begin to percolate out of the product requirements phase and into the underlying filter specifications. At a minimum, this requires filters with high orders, i.e., longer lengths, which implies more storage for the coefficients and higher execution rates for the processor. Moreover, very high order filters can be difficult (if not impossible) to design.
This article introduces the concept of multirate filters as a practical approach to the design and implementation of finite impulse response (FIR) filters with narrow spectral constraints. Multirate filters change the input data rate at one or more intermediate points within the filter itself while maintaining an output rate that is identical to the input rate. These filters can achieve both greatly reduced filter lengths and computational rates as compared to standard single-rate filter designs, and thereby provide a practical solution to an otherwise difficult problem. Note that reference textbooks are available on the topic of multirate signal processing [1,2], and an excellent tutorial article was recently published [3].
Multirate FIR filters can leverage many standard FIR filter design methods; in this article, we will use an equiripple design technique (Remez) throughout. To illustrate the design process, we will give examples of multirate filters for lowpass, highpass, and bandpass filters with very narrow passbands using MATLAB for filter design and Simulink for simulation results. The concepts learned for very narrow filters will be extended to the design of filters with very wide passbands. But since multirate filters are most readily understood by first considering a lowpass filter with a very narrow passband, we will focus on this design first.
Consider the design of a lowpass filter that preserves all frequencies below 70 Hz, removes all frequencies above 80 Hz, and runs at an input sample rate of 8 kHz. The filter cutoff frequency is therefore 80 Hz, or in terms of normalized frequency, 80/4000 or 0.02, where the normalization is done with respect to half the sample rate. In order to completely specify an FIR filter, we must also specify the amount of desired maximum passband attenuation (for this example, 0.15 dB or 0.01 linear) and minimum stopband attenuation (80 dB or 10^{-4} linear) that we wish to achieve.
For a single-rate equiripple FIR filter, we can estimate the order of the required filter using MATLAB and the Signal Processing Toolbox function remezord
:
remezord([70 80], [1 0], [0.01 1E-4], 8000)
This design requires an order 2511 (i.e., 2512-coefficient) FIR filter, a very long filter indeed! The actual design of the filter proceeds as follows:
b = remez(2511, [0 70 80 4000]/4000,... [1 1 0 0], [1 50]);
Assuming that one multiply and accumulate (MAC) operation is required per pair of (symmetric) filter coefficients, the required computational rate is 2512/2 MACs/sample * 8000 samples/sec or just over 10 million MACs/sec (multiply-accumulate operations per second). This will serve as our reference as we explore the benefits of multirate filtering.
Note that the bandwidth of the signal output from the single-rate filter designed above (see Fig. 2) is quite narrow with respect to the sample rate. The minimum sample rate necessary to reconstruct a signal with 80 Hz bandwidth is 160 Hz, and this requires only 160/8000 or 1/50 the rate at which the output is presently sampled. This implies that we could decimate the output by 50 (i.e., remove 49 out of every 50 samples) without aliasing the output signal.
Pursuing this further, the cascade of an FIR filter with a decimator is most efficiently performed with a polyphase FIR decimation filter [1] which reduces the computational rate by the decimation factor, M. This achievement is most easily understood by considering that only one out of every M input samples is output by the decimator; all others are "thrown away." Thus, one could consider computing one out of every M output samples. The polyphase FIR decimation filter is structured in such a way as to balance the computation of this one output sample over M consecutive inputs, thereby obtaining the noted reduction in computational rate.
However, decimation is only half the story. For general multirate filtering, we assume that the output rate is equal to the input rate. Thus, we need to interpolate the output of the polyphase FIR decimation filter back up to the original input rate. This may be performed by the cascade of an upsample operation (in which L-1 zeros are inserted after every sample) followed by an FIR filter which is most efficiently performed with a polyphase FIR interpolation filter. Similar to the decimator, this implementation achieves the interpolation while reducing the computational rate by the factor L. This outcome is most easily understood by considering that only one out of every L interpolated samples is non-zero prior to entering the interpolation filter. This means we may effectively skip the filter computation for the zero-valued samples, achieving a reduction in the rate of computation.
Based on this development, we can see that the cascade of polyphase FIR decimation and interpolation filters will form an efficient, multirate filter. A single-stage multirate filter will do nothing to ease the filter length requirement, but it can go a long way toward reducing the computational load. Returning to the lowpass filter example, we can choose to decimate the filter output by a factor of 48 (slightly smaller than the maximum decimation factor of 8000/160, or 50 as a margin of safety).
In the Simulink model shown in Fig. 2, the multirate filter is constructed using the FIR Decimator and FIR Interpolator blocks from the DSP Blockset. The total number of filter coefficients is 2512 (i.e., the same number as for the single-rate FIR filter design). The computational rate is
(2512/2 MACs/sample) * (8000 samples/sec) * (1/48 decimation factor) * 2 (decimator plus interpolator stage),
or approximately 419 kMACs/sec (thousands of MACs/sec).
Additional design efficiency can be achieved by using several cascaded multirate stages, which reduces the total number of filter coefficients. For example, instead of implementing a design based on a single decimation factor of 48, we could consider two stages of decimation using 12 and 4, or 24 and 2, or any other combination of factors that multiply to 48. Choosing the first pair of factors, the corresponding implementation is shown in the Simulink model of Fig. 3.
The first stage decimates by a factor of 12; thus, the corresponding first-stage filter must cut off before 1/12 of the normalized sampling frequency (approximately 0.083) to prevent aliasing. Note that this requirement is less constraining than our overall specification, which requires an 80 Hz cutoff. (Recall that an 80 Hz cutoff translates to 80/4000 or 0.02 in normalized frequency, and is smaller and thus more constraining than 0.083.) This can lead to a shorter (first stage) filter design. The second stage filter is then responsible for attaining our final filter requirements (70 Hz passband, 80 Hz stopband), but now operates at 1/12 the original sample rate. Thus, the design constraints placed on the second stage are eased, resulting in a shorter second-stage filter. However, since we are now cascading two filters, the passband ripple requirement must be halved for each filter stage, and reducing the ripple specification tends to increase overall filter length. Choice of decimation factor for each stage allows a degree of optimization across these design factors, and can lead to a design that is more efficient than single-stage designs.
Using the same equiripple filter design technique as the previous examples, an estimate of the first-stage filter order is 225, and is computed using:
remezord([1/12-.03 1/12], [1 0],... [0.5e-2 1e-4], 2)
Note that the passband ripple specification has been reduced from 0.01 to 0.005, and the width of the transition band (.03) was chosen so as not to encroach on the desired final passband (80 Hz bandwidth). The second-stage filter requires (coincidentally) another order 225 filter:
remezord([70 80], [1 0],... [0.5e-2 1e-4], 8000/12)
This design therefore requires a total of 452 coefficients, which is significantly fewer than the single-stage design (2512 coefficients). The rate of computation for the decimator stage of the two-stage design is
(226/2 MACs/sample) * (8000 samples/sec)
* (1/12 decimation factor)
+ (226/2 MACs/sample) * (8000/12 samples/sec)
* (1/4 decimation factor).
The interpolator is again identical, leading to a total of 188 kMACs/sec.
A comparative summary of the single-rate, single-stage multirate, and two-stage, multirate design examples are given in Table 1. Note that the two-stage, multirate filter design can be further improved (e.g., shorter filters and reduced computational rate) by optimizing the choice of decimation factors.
Finally, frame-based implementation methods are always available to further increase the throughput of a multirate model. The source block (Sine Wave) in Fig. 4 produces 48-element frames of data; changes in sample rates are revealed by corresponding changes in frame size. Thus, decimation by 12 in the first filter stage yields an output frame size of 48/12, or four samples. Note that the frame size output by the source must be a multiple of the aggregate decimation factor (here, 48). All of the multirate techniques explored in this article can be implemented using sample-based or frame-based methods.
The preceding lowpass filter with narrow passband can be readily extended to bandpass and highpass filter designs. Multirate techniques can be leveraged to achieve a bandpass filter with a narrow passband by employing quadrature modulation. Basically, the center frequency of the passband of interest is modulated to baseband (i.e., to zero frequency), filtered using a narrow lowpass filter, and then modulated back to the original center frequency. Quadrature modulation requires a complex signal to enter the lowpass filter; this increases the amount of computation because both the real and imaginary parts of the signal must be processed. As shown in Fig. 5, complex data is supported by Simulink and is denoted by "(c)" on the signal lines in the model.
Note that the modulation and demodulation operations are realized by multiplying with a discrete complex exponential modulation function, m(n), sampled at Fs samples/sec and generated at the center frequency f of interest:
\[m(n) = \cos\left(\frac{2\pi f_n}{F_s}\right) + j \sin\left(\frac{2\pi f_n}{F_s}\right)\]
A highpass filter with narrow passband is an extension of the bandpass design above, with the center frequency set to half the sample rate (\(F_s/2\)). In effect, the narrow lowpass filter is modulated to half the sample rate and becomes a narrow highpass filter. In this case, the modulation function degenerates to
\[m(n) = \cos(\pi n) + j \sin(\pi n) = \cos(\pi n)\]
The modulation function is no longer complex, and the values of the cosine term are simply \((-1)^n\). Moreover, instead of multiplying the input signal by \(m(n)\), we can choose to multiply the corresponding FIR filter coefficients by the same alternating sequence of +1 and 0, effectively modulating the filter itself to the frequency band of interest. Thus, the multirate implementation of a narrow highpass filter is identical to that of a narrow lowpass filter (as shown in Figs. 3 and 4) with the sign of every other coefficient negated.
The previous multirate narrow passband filter designs can be used to develop highpass and lowpass filters with wide passbands. Assuming the transfer function of a lowpass filter is \(L(z)\), the transfer function of a corresponding highpass filter with a wide passband is \(H(z) = 1 - L(z)\). Essentially, we are subtracting the lowpass filter response from an allpass filter response, leaving ourselves with a wide highpass filter.
The wide highpass filter can be realized in the time domain by subtracting the multirate lowpass filter output from a delayed replica of the original input signal, as shown in Fig. 6. Care must be taken to select the delay that exactly equals the group delay of the multirate filter. For a single-rate FIR filter, the group delay is half the filter length. For a two-stage multirate lowpass design, the group delay of each filter must be calculated relative to the original input rate, i.e., multiplying the single-rate group delay values by the preceding decimation factors. Keeping in mind that there are both decimation and interpolation filter stages, the final group delay N of the earlier two-stage lowpass filter design is
N = 225 + 225 * 12 = 2925 samples
In much the same manner, a wideband lowpass filter may be designed from a narrow highpass filter.
Efficient filters with very narrow or very wide passbands can be designed and implemented using multistage multirate DSP filtering techniques. MATLAB and Simulink provide an excellent environment for performing the design, analysis, and simulation of these complex systems.
Published 2000
P. P. Vaidyanathan, Multirate Systems and Filter Banks, Prentice Hall, Englewood Cliffs, New Jersey, 1993.