Design Finite Impulse Response Digital Filters

Interpolation, as was noted last month, is an effective means of achieving a desired filter response (such as an equiripple response) with a reduced number of multipliers. The approach was applied to the implementation of interpolation finite-impulse-response (IFIR) filters. This month, in Part 3, the concept of interpolation will be further explored and applied to the design of various FIR filter types, including Nyquist filters and more complex filter banks.

In the context of multirate signal processing, interpolation usually refers to band-limited interpolation. Band-limited interpolation is based on the notion of an underlying band-limited continuous-time signal that is being sampled. Ideal band-limited interpolation will take a digital (sampled) signal and produced an interpolated signal that will be identical to the signal that would be obtained by sampling the underlying continuous-time signal at a higher rate. Ideal band-limited interpolation can be accomplished by means of upsampling and using an ideal lowpass filter. Especially interesting is a time-domain interpretation of the ideal interpolator, which leads naturally to polyphase implementations.

The key concept of band-limited-interpolation is that a signal to be interpolated is a sampled version of a band-limited continuous time signal. Denote the continuous-time signal by x_c(t) and suppose its spectrum is zero for all |f| > f_max. Figure 25 shows the frequency spectrum X(2πjf). If the signal is sampled at a sampling frequency f_s = 2f_max, the following signal results:

x_T = { x_c(nT) }, T = 1/f_s

Figure 26 shows its spectrum, xT (e^2πjf/fs). Now suppose the continuous-time signal was sampled at a rate of f′s = Lf_s = 2Lf_max. The sampled signal at the higher rate can be represented as:

x_T′ = {xc(mT′)}, T′ = 1/ f′s = T/L,

where:

m = Ln + k , k = 0 ... (L − 1)

will have a spectrum x_T′ (e^2πjf/fs) as shown in Fig. 27 for the case L = 2.

The job of the ideal interpolation filter should now be clear from the frequency-domain standpoint. Take the discrete-time signal with spectrum xT (e^2πjf/fs) and digitally produce the discrete-time signal x_T′ (e^2πjf/f′s) that would have been obtained from sampling the original continuous-time signal at a rate of f′_s = Lf′_s.

Figure 28 shows the response of an ideal interpolation filter. Clearly, it is a lowpass filter with periodicity of f′_s , i.e. it must be operating at the high sampling rate. For this reason, it is necessary to upsample the input signal (although this is not necessary in practice where efficient algorithms are used) by inserting an appropriate amount of zeros between samples in order to feed the interpolation filter a signal at the correct rate. More precisely, the response of the ideal filter H_D(e^2πjf/f′s) for the general case of interpolation by a factor of L is given by:

The impulse response of the ideal interpolation filter can be found from the inverse DTFT (ref. 1):

Using the fact that Lf′_s = 1/ T′ and T′ = 1/Lf_s, the following results:

As expected for an ideal lowpass filter, an infinite impulse response is required to realize it. Further insight for the ideal interpolation filter will be given in the next section as part of a time-domain analysis.

Once again, the key idea of ideal band-limited interpolation is to digitally produce a signal that would be exactly the same as a signal obtained by sampling a band-limited continuous-time signal at the higher sampling rate. Figure 29 depicts the situation in the time domain.

Assuming the Nyquist sampling criterion has been satisfied (i.e., the continuous-time signal is band-limited and has been sampled at a rate of f_s = 1/T ¾ 2f_max), no information has been lost from the continuous-time signal x_c(t). Therefore, it should be possible to somehow recreate any instantaneous value x_c(t₀) of the continuous-time signal from the sampled signal x_T. Looking at Fig. 29, it can be seen that the job of the fivefold interpolator is to take every input sample x_T and produce five output samples {x_T′}, m = 5n + k, k = 0...4 as follows (note that T = 0.5 and T′ = T/5 = 0.1):.

• x_T′ = x_T
• x_T′ = x_T
• x_T′ = x_T
• x_T′ = x_T
• x_T′ = x_T

Page Title

In general, the ideal interpolator consists of a bank of L filters which will fractionally advance the input signal by a factor of k/L, k = 0...(L − 1). The outputs of the filters are then interleaved (i.e., only one filter needs to operate per high rate output sample) to produce the interpolated signal. The L filters that comprise the filter bank are the fractional advance filters H_k(z):

H_k(z) = z ^k/L, k = 0...(L − 1)

Evaluating this on the unit circle results in

H_k(e^jω) = e^jωk/L, k = 0...(L − 1)

so that each filter H_k(e^jω) is allpass, i.e., |H_k(e^jω)| = 1 and has linear phase, arg{H_k(e^jω)} = ωk/L.Herein lies the impossibility of designing these filters. They cannot be designed as FIR filters because no FIR filter can be allpass (except for a pure delay). They cannot be designed as IIR filters, because no stable IIR filter can have linear phase. However, it is clear how to approximate the ideal interpolation filter bank. FIR approximations can produce the exact linear phase, while approximating an allpass response the best way possible. On the other hand, IIR approximations will be exactly allpass in response, while trying to produce the required phase. It is insightful to realize that the filters comprising the filter bank are the polyphase components of the ideal interpolation filter derived in Eq. 5. Thus, this view of the ideal interpolator has the efficient polyphase structure "built-in."

Indeed, the impulse response of each fractional advance filter in the filter bank is given by the inverse DTFT,

which corresponds to the L decimated sequences of the ideal impulse response by again writing uniquely m = L_n + k, k = 0 ... (L − 1) in Eq. 5.

While interpolation filters are simply lowpass filters that can be designed with the various techniques outlined previously, the polyphase filters that compose the ideal interpolation filter give some insight on things to be looking for when designing interpolation filters. Consider an interpolation by a factor of L. The ideal L polyphase filters will have a group-delay given by:

−k/L, k = 0 ... (L − 1)

For simplicity, consider an FIR approximation to the ideal interpolation filter where the order is of the form N = 2LM. Them each polyphase filter will have the order N/L = 2M. Note that the ideal interpolation filter is infinitely noncausal. After finite length truncation, it is possible to make the approximation causal by delaying by one-half the filter order, N/2. However, because the filter will be implemented in efficient polyphase form, it is possible to make each polyphase component causal by delaying it by M samples.

The delay will mean the introduction of a phase component in the response of each polyphase component. So that instead of approximating the ideal fractional advance e^jωk/L the polypase components will approximate e^{jω(k/L − M)}. Consequently, the group delay will have the form

A problem that arises is that even though the FIR approximation to the ideal interpolation filter is symmetric and thus has linear phase, the polyphase components are not necessarily symmetric and thus will not necessarily have exact linear phase. However, for each non symmetric polyphase filter, there is a mirror image polyphase filter which will have the exact same magnitude response with a mirror image group delay that will compensate any phase distortion.

When we analyzed the behavior of the ideal interpolation filter in the time-domain, we saw that for every input sample, L samples are produced including one that is exactly the same as the input sample. This exact copy is "produced" by the polyphase filter that has allpass magnitude and zero phase (i.e., the case where k = 0). In practice, this is the only polyphase filter that can be designed exactly, albeit with a group delay of M rather than 0. Roughly speaking, a Nyquist filter is one for which one of its polyphase components is a pure delay and thus leaves the input signal unchanged (except for a possible delay). When designing an interpolation filter, it is desirable for it to be a Nyquist filter since this will ensure that even a nonideal filter will allow the input samples to pass through unchanged. It can also be computationally advantageous since one of the polyphase subfilters will have no multipliers.

Nyquist filters are also called Lth-band filters because the passband of their magnitude response occupies roughly 1/L of the Nyquist interval. In the special case of an interpolation by a factor of 2, the filters are known as half-band filters. Halfband filters are commonly used when interpolating (or decimating) by a factor of 2. The cutoff frequency for a halfband filter is always 0.5π. Moreover, the passband and stopband ripples are identical, limiting the degrees of freedom in the design. The function "firhalfband" designs FIR halfband filters. The specifications set still follows the triangle metaphor shown in Fig. 2, taking into account the limitations just described. The following three function calls design three equiripple linear-phase halfband filters using a different pair of specifications in each case from the three available-order (N), transition-width (TW), and peak passband/stopband ripple (R):

b1 = firhalfband(102, .47); %N and TW
b2 = firhalfband(102, .01, 'dev'); %N and TW
b3 = firhalfband('minorder', .47, .01); %TW and R

Page Title

It is possible to analyze how the design compares to the ideal interpolation filter by creating an FIR interpolator object and looking at its polyphase subfilters, for example, using the third filter, b3,

h = mfilt.firinterp(2,2*b3); polyphase(h)

Figures 30 and 31 show the magnitude and group-delay responses for the polyphase components of this filter. Note that M = N/2L is 16.5 in this case, so that the group delays are exactly M − k/L, k = 0, 1. The only deviation from an ideal filter (ignoring an overall delay of M samples) comes from the fact that one of the polyphase subfilters is not a perfectly allpass response.

Nyquist filters are characterized in the time-domain by their impulse response being exactly equal to zero every L samples (except the exact middle sample of the impulse response). This is precisely why we get a polyphase sub filter that is a perfect allpass delay and allows the samples to be interpolated to pass through the filter unchanged. Designing a filter that is both a lowpass and simultaneously satisfies the just mentioned time-domain characteristic is not a trivial task except for the case of window-based designs.^{13, 14}

Nevertheless, the advantage of conventional optimal equiripple designs over a Nyquist window-based design is not as clear in this case as it is with any conventional lowpass filter. As an example, consider a Kaiser window Nyquist filter design with a stopband attenuation of 40 dB. Nyquist filters are often designed in terms of their rolloff factor, ρ, due to their applications in communications. The rolloff factor is related to transition-width (TW) supply by TW = ρπ/L. In this example, ρ = 0.1 and L = 4; thus, the transition width is 0.0025π:

b1 = firnyquist('minorder', 4, .1, .01); L = 4

The resulting filter is of 90th order. If an equiripple filter of the same order and same attenuation is designed, it will result in a filter with a smaller transition width, but on that does not satisfy the time-domain requirement:

b2 = fircequip(90, .25, );

Figure 32 shows the magnitude responses of the polyphase subfilters for the Nyquist window-based design. For comparison, Fig. 33 shows the magnitude responses for the optimal equiripple design. Note the better approximation to allpass filters in the Nyquist design compared to the equiripple design (albeit for a slightly smaller interval, which is the tradeoff). Similarly, if we compare the group-delay response of the polyphase subfilters, the Nyquist design once again better approximates the ideal constant group-delay as compared to the equiripple design. Figure 34 shows the group-delay responses for the polyphase subfilters of the Nyquist design. Figure 35 shows the group-delay responses for the polyphase subfilters of the equiripple design.

Figure 36 shows a two-channel filter bank. Filters H₀(z) and H₁(z) are called the analysis filters while G₀(z) and G₁(z) are the synthesis filters. The filter bank is called perfect reconstruction if the end-to-end system acts as a delay, i.e., if the output signal is simply a delayed (and possibly scaled) version of the input.

It is well known^{10, 15, 16} that perfect reconstruction can be achieved if

0.5G₀(z)H₀(-z) + 0.5G₁(z)H₁(-z) = 0

and

0.5G₀(z)H₀(z) + 0.5G₁(z)H₁(-z) = z^-k

Starting with a prototype lowpass filter H(z) of odd order N, the following selection for the filters results in perfect reconstruction using solely FIR filters¹⁵:

H₀(z) = H(z)
H₁(z) = z^-NH₀(-z^-1)
G₀(z) = 2z^-NH₀(z^-1)
G₁(z) = 2z^-NH₁(z^-1)

This type of perfect reconstruction filter bank is called an orthogonal filter bank or a power-symmetric filter bank.¹⁵

The function "firpr2chfb" designs equiripple FIR filters H0(z), H1(z), G0(z), G1(z) such that the filter bank is an orthogonal perfect reconstruction filter bank. The parameters to specify are simply the filter order N and the passband-edge frequency wp. A prototype lowpass filter is designed from which the four required filters are obtained.

For example,

= firpr2chfb(99,.45);

Page Title

Alternatively, the peak stopband ripple can be specified. As usual, we can obtain minimum-order designs by specifying both the passband-edge frequency and the peak stopband ripple. In all cases, the design specifications apply to the prototype filter H(z),

= firpr2chfb(99, le-3, 'dev');

= firpr2chfb('minorder', .45, le-3);

The power-symmetric term is used because for these designs, the following holds:

It is possible to look at the magnitude-squared responses of H₀(z), and H₁(z) using the "fvtool" function. Figure 37 shows the magnitude-squared responses for N = 19 and ω_p = 0.45π. Notice how where one filter's ripple rises the other filter's ripple declines to add up to one.

Increasing the filter order (and possibly the passband-edge frequency) improves the lowpass/highpass separation provided by the analysis filters but does not have an effect on the perfect reconstruction characteristic of the overall system.

Several factors have to be taken into account when implementing an FIR filter using fixed-point arithmetic. For one thing, the coefficients have to be quantized from double-precision floating point in which they are designed into fixed-point representation with usually a smaller number of bits. It is necessary to make sure to make the most of the limited number of bits available. Furthermore, performing the arithmetic in fixed-point will introduce further quantization errors when actually filtering with the quantized coefficients. These quantization errors should be minimized as much as the hardware at hand allows.

The Filter Design Toolbox of MATLAB employs certain notation to represent fixed-point numbers. Consider a register used to store a fixed-point number,

The register has B bits (it has a wordlength of B), the value of the kth bit is given by b_k which can obviously be only 0 or 1. A two's complement fixed-point number stored in such a register has a value given by

where L is a positive or negative integer to be described in the following section.

From Eq. 6, it can be seen that the value of a fixed-point number is determined by assigning weights of the form 2 − m to each bit. The leftmost bit, b₀ , has the largest weight, 2^{B - L - 1} , and is called the most-significant bit (MSB). The rightmost bit, b_{B - 1}, has the smallest weight, 2^-L , and is called the least-significant bit (LSB).

Given the bit values, b_k, the pair {B, L} completely characterizes a fixed-point number, i.e., suffices in determining the value that the bits represent. Such a pair is called the format of a given quantity, and is stored in a two-element vector, .

Consider the following filter

b = gremez('minorder',,...

,);

The filter has an attenuation of 80 dB and a its largest coefficient is 0.1206.

The first thing to do is check if there are enough bits available to represent the coefficients and provide the required dynamic range. A good rule of thumb¹⁷ is to assume 5 dB/b for the dynamic range (note that the usual 6 dB/b rule does not apply because the quantization error for the filter coefficients tends to be correlated, especially at the extremes of the impulse response). In this example, at least 16 b is needed to provide the required dynamic range (80 dB). But it is not sufficient to simply apply 16 b to the design. For example, the following code creates a fixed-point FIR filter using 16 b to represent the coefficients in fractional format:

Hq = qfilt('fir', {b},...'CoefficientFormat', );

Figure 38 shows the magnitude response of the quantized filter, with the nonquantized magnitude response also shown for comparison. Note that the stopband attenuation for the quantized response is significantly less than 80 dB at various frequency bands. The problem is the poor utilization of the available range for the format (Fig. 39).

Two equivalent approaches can be used to make the most of the 16 b. To use the format, the coefficients can be scaled by multiplying them by a factor of 8 to make the largest coefficients as close to 1 as possible without overflowing. Alternatively, the format can be used so that the quantization range becomes Figure 40 shows the magnitude response using this format, with an obvious improvement over the first case.

Note that whether the coefficients are scaled and the format is used, or the format is used without scaling, the actual stored value (the binary bits) of each coefficient is the same. However, in the former case, the filter now has a gain of 18 dB due to the multiplication by eight. But this can be compensated at the end, by moving the binary point 3 b to the left, without changing the bits.

Page Title

To emphasize the point regarding the need to use both the right number of bits and use them wisely, Fig. 41 shows the magnitude responses of four different quantizations of the same filter. In all cases, the format has been selected to cover the range . If fewer than 16 b are available, it is probably wise to redesign the filter, since many more multipliers will be required than be applied with Eq. 6. (If the design specification is changed from 80 dB to 60 dB, 178 multipliers are required as opposed to 220. If it is reduced to 40 dB, 134 multipliers are required. Of course, it is not a given that the application can allow this change in specifications. The point is having less than 16 b makes it unfeasible to attain 80 dB.) On the other hand, increasing the precision to 24 b provides only modest improvements in this case.

In the April 2004 issue, the final installment of this four-part article series on FIR filter design with MATLAB will explore the use of fixed-point arithmetic and how to use a digital-signal-processor (DSP) accumulator to maximum precision in the implementation of digital filters. This final part will also provide a design example based on creating two FIR filters for a GSM digital downconverter.

REFERENCES

21. P.P. Vaidyanathan and T.Q. Nguyen, "Eigenfilters: a new approach to least squares FIR filter design and applications including Nyquist filters," IEEE Transactions on Circuits and Systems, Vol. CAS-34, pp. 11-23, 1987.
22. S.K. Mitra, Digital Signal Processing. A Computer-Based Approach, 2nd ed., McGraw-Hill, New York, 2001.
23. N.J. Fliege, Multirate Digital Signal Processing, Wiley, Chichester, 1994.
24. F. Harris, "Multirate Signal Processing in Transmitter & Receiver Designs," UCLA Extension course, November, 2000.
25. This spectrum shows a band-limited continuous-time signal.
26. This spectrum shows a sampled signal with sampling rate, f_s = 2f_max.
27. This spectrum shows a sampled signal with sampling rate of f
9_s = 4f_max.
28. The response of an ideal interpolation filter is overlaid here with the spectrum of a signal sampled at a rate of f_s = 2f_max.
29. This plot of an ideal band-limited interpolation is depicted in the time domain.
30. This plot shows the magnitude response for the polyphase subfilters of a halfband FIR filter. Ideally, both subfilters would be perfectly allpass.
31. This plot shows the group-delay response for the polyphase subfilters of a halfband FIR filter.
32. This plot shows the magnitude response for the polyphase subfilters of a Nyquist FIR filter designed with the window method. The polyphase subfilters better approximate allpass filters than a comparable equiripple design for the majority of the frequency band.
33. This plot shows the magnitude response for the polyphase subfilters of a optimal equiripple lowpass FIR filter. None of the subfilters behaves as a perfect allpass, an indication that this is not a Nyquist filter.
34. This group-delay response was generated for the polyphase subfilters of a Nyquist FIR filter of order 90 and L = 4.
35. This group-delay response was generated for the polyphase subfilters of a conventional equiripple lowpass design that could be used for interpolation with L = 4.
36. This block diagram shows the two-channel FIR subband filter bank.
37. This plot shows the magnitude-squared responses of the analysis filters in an FIR perfect reconstruction filter bank. The two filters are power-complementary.
38. This plot shows the magnitude response of the filter quantized with format.
39. This plot shows the impulse response of filter to be quantized shown relative to the available range for the coefficient format selected.
40. This plot shows the magnitude response of the filter quantized with format.
41. This plot shows the magnitude responses for various quantizations of a filter with 80-dB stopband attenuation.