Design of FIR Filters Using the Frequency Sampling Method
2022-09-12 07:18:38

This article will review the basics of the frequency sampling method for designing Finite Impulse Response (FIR) filters.

FIR filters have two main advantages over Infinite Impulse Response (IIR) filters: 

  • An FIR filter has no feedback loop and is inherently stable. (We always need to check the stability of an IIR filter.)
  • An FIR filter can easily provide a linear-phase response, which is crucial in phase-sensitive applications such as data communications, seismology, etc.

This article will teach you how to design an FIR filter using the frequency sampling method.

FIR Filter Design by Windowing

Let's assume that we want to design a Finite Impulse Response (FIR) filter with the desired frequency response Hd(ω). The window method applies the inverse discrete-time Fourier transform to find the corresponding time-domain representation, hd(n), as given by

 

hd(n)=12πππHd(ω)ejnωdω

Equation 1

 

Generally, hd(n) is not of finite length and since we were looking for a finite impulse response, we can truncate hd(n) and achieve an approximation of the desired response. This is equivalent to multiplying hd(n) by a rectangular window.

As discussed in a previous article, From Filter Specs to Window Parameters in FIR Filter Design, we can use other window functions to achieve a trade-off between ripples in the passband and the sharpness of the transition band.

For many classic filters, we can easily calculate the integral of Equation 1. For example, assume that we are designing an ideal lowpass filter with cutoff frequency of ωc. Then, we have

 

Hd(ω)={1|ω|ωc0else

Equation 2

 

In this particular example, some mathematical manipulations can easily lead us to the corresponding hd(n) (see Equation (8) in this article) because Hd(ω) is given by a simple equation. Then, we can apply a suitable window function to arrive at a finite-length response which is an approximation of the desired filter. For more details of this discussion, please see this article on FIR filter design by windowing.

However, if Hd(ω) has a complicated equation, we wouldn’t be able to easily calculate Equation 1. In these cases, we can use the frequency sampling method to design FIR filters with an arbitrary Hd(ω).

In this article, we will first review a practical example where the required FIR filter has a complicated frequency response and the frequency sampling method is quite helpful. Then, we will review the basics of FIR filter design using this method.

An FIR Filter to Compensate for the sin(x)x Distortion

We can think of a practical DAC as an ideal DAC, which produces an impulse train, followed by a zero-order hold block.

 

Figure 1. Practical DACs generally apply a zero-order hold to the output values. Image courtesy of Digital Signal Processing.

 

The hold circuit is required because, in practice, we cannot produce the narrow pulses of an ideal DAC. However, placing a zero-order hold at the output of an ideal DAC (as shown in Figure 1) means that the ideal output, ys(t), will be modified by the transfer function of the hold circuit. It can be shown that the transfer function of the zero-order hold is given by

 

H(f)=thsin(πfth)πfth

 

where th is the hold time which is generally equal to the sampling period. The following figure shows the normalized frequency response of the hold circuit, i.e., HDAC=H(f)th.

 

Figure 2. The normalized frequency response of the hold circuit versus ffs.

As the above figure shows, the frequency components near fs2 experience an attenuation of about 4 dB compared to the low-frequency components. We discussed in a previous article on multirate DSP in D/A conversion that it is possible to employ a multirate system to alleviate the problem. However, for relatively small values of oversampling ratio (e.g., L8), we may still need to compensate for the hold circuit distortion. To this end, we can apply a filter with the frequency response of 1HDAC to y(n) in Figure 1.

Assuming that we are using a multirate system (see Figure 3), we can modify the interpolation filter itself to avoid the use of an extra filter. In this case, the magnitude of the interpolation filter response will be proportional to 1HDAC in the filter passband (0ωπL for an L-fold interpolation) and zero otherwise.

 

Figure 3. H(z) can be designed to both eliminate the frequency components above πL and compensate for the hold circuit distortion. Image courtesy of Digital Signal Processing.

 

To summarize, we can use Equation 1 to design ideal filters such as the one given by Equation 2. However, this method is not easy when designing a filter with an arbitrary frequency response, e.g., when designing the magnitude of H(z) to be equal to 1HDAC in a particular band. In these cases, we can use the frequency sampling method.

The Frequency Sampling Method

The frequency sampling method is based on replacing the integral of Equation 1 by an approximate sum. Let's assume that we need to design a filter with the frequency response shown in Figure 4.

 

Figure 4. The frequency sampling method can be used to design a filter with an arbitrary frequency response such as H(ω).

 

To have simpler mathematical expressions, we can approximate H(ω) using the red curve in the following figure.

 

Figure 5. The red curve is an approximation of H(ω).

In this particular example, the desired frequency response is sampled at regular frequency intervals 2π13k for k=0,1,2,,12 (the samples are shown by dots in the figure). If we were to calculate the integral of H(ω) itself, we could simply approximate the integral calculating the area under the red curve which is sum of 13 rectangles of width 2π13 and length determined by the samples of H(ω). To perform the integral of Equation 1, we also need to take the term ejnω into account. Hence, for N=2M+1 samples of H(ω), we can approximate Equation 1 using

 

hd(n)=12π(k=MMH(ωk)ejnωk2π2M+1)

Equation 3

 

where ωk=2π2M+1k.

In the above equation, H(ωk)ejnωk and 2π2M+1 respectively represent the length and width of the rectangles used to approximate H(ω)ejnω. Interestingly, comparing Equation 3 with Equation 2 of this article on DFT analysis in DSP, we observe that the frequency sampling method is actually calculating the inverse discrete Fourier transform of samples of H(ω).

We can further simplify Equation 3 by taking into account that the frequency response is symmetrical around the origin, i.e., H(ωk)=H(ωk). Hence, we obtain

 

hd(n)=12M+1(H(ω0)+k=1MH(ωk)(ejnωk+ejnωk))

 

which gives

 

hd(n)=12M+1(H(ω0)+2k=1MH(ωk)cos(nωk))

Equation 4

 

Note that, after calculating hd(n), we can apply a window function in a way similar to the windowing method discussed above.

Frequency Sampling Method Example 1

Use the frequency sampling method to design a 9-tap lowpass FIR filter with a cutoff frequency of 0.25π radians/sample.

First, we need to find the value of the frequency response samples. Assuming an ideal response, the samples below 0.25π are equal to 1 and the other samples are zero. Then, we obtain

 

 

Hence, Equation 4 gives

hd(n)=19(1+2cos(2π9n))

 

Based on this equation, we obtain the corresponding time-domain representation as given by the following table:

 

 

The frequency response of the designed filter is shown in Figure 6.

 

Figure 6. The frequency response of the filter designed in Example 1.

 

This figure shows that there are some ripples in the passband and stopband of the designed filter. These ripples make the obtained response deviate from that of the target ideal filter. However, note that, at 2π9k, k=0,1,,8, the filter exhibits exactly the magnitude response that we were originally looking for (see Table 1). This can be explained by noting that the frequency sampling is based on the inverse DFT analysis. We applied the inverse DFT to the frequency domain samples of Table 1 and obtained the corresponding time-domain sequence, hd(n). Figure 6 is the frequency content of hd(n) and its equally-spaced samples at 2π9k, k=0,1,,8 give the nine-point DFT of hd(n). Clearly, these samples should match those given by Table 1.

We can increase the number of the frequency domain samples to get a better approximation of the target filter and achieve a sharper transition. Obviously, this will require a filter of a higher order. The next example explains the design of a 25-tap filter for the lowpass filter of example 1.

Frequency Sampling Method Example 2

Use the frequency sampling method to design a 25-tap lowpass FIR filter with a cutoff frequency of 0.25π radians/sample.

First, we find the value of the frequency response samples. Assuming an ideal response, the samples below 0.25π are equal to 1 and the other samples are zero. Then, we obtain

 

 

and

 

 

Hence, Equation 4 gives

hd(n)=125(1+2(cos(2π25n)+cos(4π25n)+cos(6π25n)))

 

Based on this equation, we obtain the corresponding time-domain representation as given by the following table:

 

 

and

 

 

The frequency response of the achieved 25-tap filter is shown in Figure 7. 

 

Figure 7. The frequency response of the 25-tap filter with cutoff frequency of 0.25 radians/sample.

 

To reduce the passband ripples, we can employ a window function; however, this will increase the transition band and we may need to further increase the filter length to achieve a sharper response.

Summary

  • When Hd(ω) is given by a simple equation, we can easily calculate the integral of Equation 1. However, for an arbitrary Hd(ω), the frequency sampling method will be easier to apply.

  • For relatively small values of oversampling ratio (e.g., L8), we may need to compensate for the hold circuit distortion. To this end, we can design the interpolation filter to have a frequency response of 1HDAC.

  • For N=2M+1 samples of H(ω), we can determine the filter coefficients of frequency sampling method using Equation 4.