Finite Impulse Response (FIR) Filters

Time and frequency response of a lowpass FIR filter designed with Parks-McClellan algorithm for N=33

We learned in the concept of frequency that most signals of practical interest can be considered as a sum of complex sinusoids oscillating at different frequencies. The amplitudes and phases of these sinusoids shape the frequency contents of that signal and are drawn through magnitude response and phase response, respectively. In DSP, a regular goal is to modify these frequency contents of an input signal to obtain a desired representation at the output. This operation is called filtering and it is the most fundamental function in the whole field of DSP.

From this Eq, one can observe that it is possible to design a system, or filter, such that some desired frequency contents of the signal can be passed through by specifying H(k) \approx 1 for those values of k, while other frequency contents can be suppressed by specifying H(k) \approx 0. The terms lowpass, highpass and bandpass then make intuitive sense that low, high and band refer to desired spectral contents that need to be passed while blocking the rest.

Filter Frequency Response


Figure below shows the magnitude response |H(F)| (as a function of continuous frequency) of an ideal lowpass filter. A lowpass filter passes frequencies near 0 while blocks the remaining frequencies.

An ideal lowpass filter and aliases

As explained in the discussion about sampling, in a continuous frequency world, the middle filter is all that exists. However, in a sampled world, the frequency response of the filter — just like a sampled signal — repeats at intervals of F_S (sampling frequency) and there are infinite spectral replicas on both sides shown with dashed lines. The primary zone of the sampled frequency response is from -0.5F_S to +0.5F_S shown within dotted red lines with its spectral contents in solid line. From this Eq on sampling, we know that filter response near 0 is the same as that near F_S and the filter response around +0.5F_S is the same as around -0.5F_S.

On the other hand, a highpass filter passes both positive and negative frequencies near 0.5F_S while blocking everything else. Observe the symmetry around \pm 0.5F_S in Figure below. Similarly, bandpass filters only pass a specific portion of the spectrum.

An ideal highpass filter and aliases

Filter Impulse Response


Since a filter is a system that modifies the spectral contents of a signal, it naturally has an impulse response h[n] just like a system. When this impulse response is finite, the filter is called Finite Impulse Response (FIR). Impulse response h[n] of an FIR filter, i.e., the sequence of its coefficients, can be found by taking the iDFT of the frequency response H[k].

FIR Filter Design


Understanding of FIR filter design starts from the following argument for continuous frequency domain. Since ideal filters have unity passband and zero stopband and a direct transition between these two, this finite frequency domain support produces an impulse response h[n] in time domain which is infinite and hence non-realizable, as explained in the discussion about frequency. When this h[n] is truncated on the two sides, it is equivalent to multiplication of the ideal infinite h[n] by a rectangular window w_{\text{Rect}}[n] as

    \begin{equation*}         h_{\text{Truncated}}[n] = h[n] \cdot w_{\text{Rect}}[n]     \end{equation*}

As mentioned previously, windowing is a process of pruning an infinite sequence to a finite sequence. We know that multiplication in time domain is convolution in frequency domain, and that the spectrum of a rectangular window is a sinc signal. Thus, convolution between the spectrum of the ideal filter and this sinc signal generates the actual frequency response H(F) of the filter.

    \begin{equation*}         H(F) = \text{Ideal brickwall response} * sinc(F)     \end{equation*}

This sinc signal has decaying oscillations that extend from -\infty to \infty. Convolution with these decaying oscillations is the reason why a realizable FIR filter always exhibits ripple in the passband and in the stopband as well as a transition bandwidth between the two. Naturally, this transition bandwidth is equal to the width of the sinc mainlobe, which in turn is a function of the length of the rectangular window in time. Truncation with a rectangular window causes the FIR filters to have an impulse response h[n] of finite duration.

Sidelobes appearing in the stopband of the filter, transition bandwidth and ripple in the passband characterize the performance of a filter. These parameters are usually specified in terms of decibels (or dBs). The filter magnitude in dB is given by

    \begin{equation*}         |H(F)|_{dB} = 20 \log_{10} |H(F)|     \end{equation*}

It should be remembered that the multiplicative factor 20 above is replaced by 10 for power conversions (because power is a function of voltage or current squared). The advantage of converting magnitudes to dBs for plotting purpose is that very small and very large quantities can be displayed clearly on a single plot.

A filter is designed in accordance with the parameters in Table below and shown graphically in the next Figure.

Lowpass filter specifications table

A lowpass filter marked with specs frequencies and their amplitudes

FIR filter design basically requires finding the values of filter taps (or coefficients) that translate into a desired frequency response. Many software routines are available to accomplish this task. A standard method for FIR filter design is the Parks-McClellan algorithm. The Parks-McClellan algorithm is an iterative algorithm for finding the optimal FIR filter such that the maximum error between the desired frequency response and the actual frequency response is minimized. Filters designed this way exhibit an equiripple behavior in their frequency responses and are sometimes called equiripple filters (equiripple implies equal ripple within the passband and within the stopband, but not necessarily equal in both bands).

Example


We design a lowpass FIR filter through Parks-McClellan algorithm in Matlab that meets the following specifications.

    \begin{align*}             &\text{Sample rate}, \quad F_S = 12\: \text{kHz} \\             &\text{Passband edge} \quad F_{\text{pass}} = 2\: \text{kHz} \\             &\text{Stopband edge} \quad F_{\text{stop}} = 3\: \text{kHz} \\             &\text{Passband ripple} \quad \delta_{\text{pass,dB}} = 0.1\: \text{dB} \\             &\text{Stopband ripple} \quad \delta_{\text{stop,dB}} = -60\: \text{dB}           \end{align*}

First, noting from above Figure that the maximum amplitude in the passband of frequency domain is 1+\delta_{\text{pass}}, we can convert dB units into linear terms as

    \begin{align*}       20 \log_{10} \big(1+\delta_{\text{pass}}\big) &= \delta_{\text{pass,dB}} \\       1+\delta_{\text{pass}} &= 10^{\delta_{\text{pass,dB}}/20} \\         \delta_{\text{pass}} &= 0.012     \end{align*}

For the stopband,

    \begin{align*}       20 \log_{10} \delta_{\text{stop}} &= \delta_{\text{stop,dB}} \\       \delta_{\text{stop}} &= 10^{-60/20} = 0.001     \end{align*}

A Matlab function firpm() returns the coefficients of a length N+1 linear phase (which implies real and symmetric coefficients, see the concept of phase) FIR filter. For finding an approximate N, the Matlab command firpmord() — which stands for FIR Parks-McClellan order — can be used.

    \begin{equation*}         N = \text{firpmord} ([F_{\text{pass}}\quad F_{\text{stop}}],[1 \quad 0],[\delta_{\text{pass}} \quad \delta_{\text{stop}}],F_S)     \end{equation*}

where the vector [1\quad 0] specifies the desired amplitude pointing towards passband and stopband, respectively. For our example,

    \begin{equation*}         N = \text{firpmord} ([2000\quad 3000],[1 \quad 0],[0.012 \quad 0.001],12000) =29     \end{equation*}

A word of caution: it can raise confusion in the sense that firpmord() uses the band widths for input amplitudes while firpm() uses the band edges. For example, firpmord() uses the desired amplitude vector [1 \quad 0] to denote passband and stopband. On the other hand, firpm() uses the amplitude vector [1\quad  1\quad 0 \quad 0] to denote band edges.

For a trial value of N=29,

    \begin{align*}       h &= \text{firpm}(N-1,[0 \quad F_{\text{pass}} \quad F_{\text{stop}} \quad 0.5F_S]/(0.5F_S),[1\quad 1\quad 0\quad 0]) \\         &= \text{firpm}(28,[0 \quad 2000 \quad 3000 \quad 6000]/6000,[1\quad 1\quad 0\quad 0])     \end{align*}

where the vector [1\quad 1 \quad 0 \quad 0] specifies the desired amplitude at band edges of passband and stopband, respectively. The normalization by 0.5 F_S is of no significance as it is just a requirement of this function. Furthermore, the filter is designed with N+1 coefficients by Matlab and that is why we used N-1 as its argument.

If you try this little code, you will see that the result thus obtained had stopband levels of only around -45 dB and a passband ripple of 0.05 dB. The sidelobes do not meet the 60 dB attenuation while the passband ripple exceeds the requirement of 0.1 dB. Therefore, N can be increased until 60 dB attenuation is achieved for N = 41. To avoid excessive filter length increase which increases the processor workload, different weighting factors can be employed for passband and stopband ripples as well. For example, a weighting factor of 1:10 for passband to stopband ([1 \quad 10] in vector form) yields,

    \begin{equation*}       h = \text{firpm}(28,[0 \quad 2000 \quad 3000 \quad 6000]/6000,[1\quad 1\quad 0\quad 0],[1 \quad 10])     \end{equation*}

In this case, the sidelobe level reached -55 dB, a 10 dB improvement over the same N=29 when there were no band weights. Next, the filter order is gradually increased until it meets the desired specifications at N = 33, a difference of 8 taps from N=41 and no weighting factors. Such trial and error approach is often necessary to design a filter with given specifications.

The impulse response h[n] and frequency response 20 \log_{10} H(F) of the final filter are drawn in Figure below. Note the stopband attenuation of 60 dB and passband ripple within 0.1 dB.

Time and frequency response of a lowpass FIR filter designed with Parks-McClellan algorithm for N=33

Filtering as Convolution


We also learned above that the higher the number of taps, the more closely the filter follows an expected response (e.g., lowpass, highpass, bandpass, etc.) at a cost of increased computational complexity. Therefore, the defining feature of an FIR filter is the number of coefficients that determines the response length, the number of multipliers and the delay in samples required to compute each output. An FIR filter is built of multipliers and adders. A delay element, which is just a clocked register, is used between coefficients. A single delay element is represented with a D_1 symbol. Figure below illustrates a 5-tap FIR filter.

FIR filter implementation

Taking into account the structure in Figure above, this 5-tap filter output can be mathematically written as

    \begin{equation*}         r[n] = h_0  s[n] + h_1  s[n-1] + h_2  s[n-2] + h[3]  s[n-3] + h[4]  s[n-4]     \end{equation*}

For a general FIR filter of length N,

(1)   \begin{equation*}         r[n] = \sum \limits _{m=0} ^{N-1} h[m] s[n-m]     \end{equation*}

If this equation sounds familiar to you, you guessed it right. Compare it with convolution Eq, and this is nothing but a convolution between the input signal and the filter coefficients. Just like convolution, at each clock cycle, samples from input data and filter taps are cross-multiplied and the outputs of all multipliers are added to form a single output for that clock cycle. On the next clock cycle, the input data samples are shifted right by 1 relative to filter taps, and the process continues.

Filter Delay


Due to convolution, the length of the output signal r[n] is given as in the article on convolution, which is repeated below.

    \begin{equation*}         \textmd{Length}\{r[n]\} = \textmd{Length}\{s[n]\} + \textmd{Length}\{h[n]\} - 1     \end{equation*}

which is longer than the input by \textmd{length}\{h[n]\} - 1 samples. We can say that every FIR filter causes the output signal to be delayed by a number of samples given by

(2)   \begin{equation*}         D = \frac{\textmd{Length}\{h[n]\} - 1}{2}     \end{equation*}

As a result, D samples of the output at the start — when the filter response is moving into the input sequence — can be discarded. Similarly, the last D samples — when the filter is moving out of the input sequence — can be discarded as well.

Filter Phase Response


As discussed in the concept of phase, the phase response of filters is also important because the shape of the signal in time domain is determined by phase alignment of its constituent complex sinusoids. A filter with unity magnitude across the spectrum does not affect the amplitude of those sinusoids but if its phase response is not linear, it can completely change the shape of the filtered signal due to different sinusoids adding up with different phases at the output. A linear phase on the other hand just delays the output without distorting it. In the post on convolution, we discussed that the output response of an LTI system to a complex sinusoid at its input is nothing but the same complex sinusoid at the same frequency but with different amplitude and phase. If this phase of the filter is a linear function of frequency k, then rearranging Eq on time and phase shifts,

    \begin{flalign*}         \begin{aligned}               \textmd{Phase shift} \quad - \textmd{constant} \cdot k \cdot \textmd{delay} ~& \xrightarrow{\hspace*{1cm}}~ s[n-\textmd{delay}] \quad \textmd{Time shift}         \end{aligned}     \end{flalign*}

for all complex sinusoids in its spectrum. Fortunately, maintaining a linear phase response in FIR filters is straightforward by making the coefficients symmetric, see this Figure as an example. This is one of the main reasons for widespread use of FIR filters in DSP applications.

Leave a Comment

Your email address will not be published. Required fields are marked *