CS3291 DSP 4.15 06/03/2006 BMGC
CS3291: Digital Signal Processing '05-'06
Section 4: A design technique for FIR digital filters
4.1.Introduction:
An FIR digital filter of order M may be implemented by programming the signal flow graph shown below. Its difference equation is:
y[n] = a0x[n] + a1x[n-1] + a2x[n-2] + ... + aMx[n-M]
Fig. 4.1
Its impulse-response is {..., 0, ..., a0, a1, a2,..., aM, 0, ...} and its frequency-response is the DTFT of the impulse response, i.e.
M
H(e j W ) = å a n e - j W n
n=0
Now consider the problem of choosing the multiplier coefficients. a0, a1,..., aM such that H( ejW ) is close to some desired or target frequency response. We shall use the inverse DTFT:
To take an example, assume we require a low-pass filter whose gain-response approximates the ideal 'brick-wall' gain-response in Figure 4.2.
If we take the phase-response f(W) to be zero for all W, the required frequency-response is:-
and by the inverse DTFT,
= (1/3)sinc(n/3) for all n.
where
A digital filter whose impulse-response is {h[n]} with each sample h[n] thus defined would have exactly the required low-pass frequency-response. However {h[n]} would have non-zero samples extending from n=-¥ to n=¥ as illustrated in fig 4.3, and would therefore not be a finite impulse response. It would also not be causal.
To produce a realisable impulse-response:
(1) Truncate {h[n]} to a FIR by setting h[n] to zero for all values of n outside the range
-M/2 £ n £ M/2 ( assume M is the order of the required FIR filter and is even ).
(2) Delay resulting sequence by M/2 samples to ensure that the first non-zero sample occurs at n = 0.
The resulting causal impulse response may be realised by setting an = h[n] for n=0,1,2,...,M. Taking M=4, for example, the finite impulse response obtained for the p/3 cut-off low-pass specification is :
{..,0,..,0, 0.14, 0.28, 0.33 , 0.28 , 0.14 , 0 ,..,0,..}
The resulting FIR filter is as shown in Figure 4.1 with a0=0.14, a1=0.28, a2=0.33, a3=0.28, a4=0.14. ( Note: a 4th order FIR filter has 4 delays & 5 multiplier coefficients ).
The gain & phase responses of this FIR filter, calculated by computer (see later), are given approximately in Figure 4.4.
Fig. 4.4
Clearly, the effect of the truncation of {h[n]} to ±M/2 and the M/2 samples delay is to produce gain and phase responses which are different from those originally specified.
Considering the gain-response first, the cut-off rate is by no means sharp, and two "ripples" appear in the stop-band, the peak of the first one being at about -21dB.
The phase-response is not zero for all values of W as was originally specified, but is linear phase ( i.e. a straight line graph through the origin ) in the pass-band of the low-pass filter ( -p/3 to p/3 ) with slope arctan( M/2 ) with M = 4 in this case. This means that f( W ) = - ( M/2 )W for | W | £ p/3; i.e. we get a linear phase-response ( for | W | £ p/3 ) with a phase-delay of M/2 samples.
Question.Why does delaying the impulse-response by M/2 produce this effect on the phase-response, remembering that we originally asked for f( W ) to be zero for all W?
Answer: if the DTFT of {h[n]} is H(ejW) = G(W)ejf(W), the DTFT of the delayed impulse-response:
Therefore, delaying {h[n]} by M/2 samples multiplies the frequency-response corresponding to {h[n]} by e-MWj/2 which increases f(W) by -MW / 2. This increases the phase-delay -f(W) / W by M/2 sampling intervals, as expected.
Question. Can we improve the low-pass filter by increasing the order to say ten?
Answer: Taking 11 terms of { (1 / np) sin (np / 3) } we get, after delaying by 5 samples:
{...0,-0.055,-.069, 0,.138,.276,.333,.276,.138,0,-.069,-.055,0,...}.
If we analyse, by computer, a tenth order FIR filter with this impulse-response, we obtain the gain response shown in fig 4.5. The graph was produced by the following MATLAB statements:
for n=1:11; a(n)= 0.3333*sinc((n-6)/3); end; % In MATLAB, sinc(x)=sin(px) / (px)
H=freqz(a,1,1000); % 1000 points of A( ejW ) into complx array H
whitebg;
plot( [0:999]/1000,20*log10(abs(H) );
axis( 0,1,-60,10] );
grid on;
xlabel( ‘rel_freq / pi );
ylabel( ‘Gain (dB)’ )
An similar graph may be produced by the single statement:
freqz( [-0.055, -0.069, 0, 0.138, 0.276, 0.333, 0.276, 0.138, 0, -0.069, -0.055] );
In may be seen in fig 4.5 that the cut-off rate for the 10th order FIR filter is sharper than for the 4th order case, there are more stop-band ripples and, rather disappointingly, the gain at the peak of the first ripple after the cut-off remains at about -21 dB. This effect is due to a well known property of Fourier series approximations, known as Gibb's phenomenon. The phase-response is linear phase in the passband ( -p/3 to p/3 ) with a phase delay of 5 samples. As seen in fig 4.6, going to 20th order produces even faster cut-off rates and more stop-band ripples, but the main stop-band ripple remains at about -20dB. This trend continues with 40th and higher orders as may be easily verified by modifying the MATLAB program above and running it again. To improve matters we need to discuss "windowing ".
Fig 4.5: Gain response of tenth order lowpass FIR filter with WC = p/3
Fig 4.6: Gain response of 20th order lowpass FIR filter with WC = p/3
4.2. Windowing:
To design the FIR filters considered up to now we effectively multiplied {h[n]} as calculated by the inverse DTFT by a rectangular window sequence {r[n]} where
This causes a sudden transition to zero at the window edges and it is these transitions that produce the stop-band ripples. The levels of these ripples may be reduced if {r[n]} is replaced by a non-rectangular window sequence, denoted { w[n] } say, which produces a more gradual transition at the window edges. The simplest non-rectangular window sequence in common use is the Hann or von-Hann window which is essentially a " raised cosine ":-
Many other types of window sequence exist ( e.g. Hamming, Kaiser ). Perhaps the most well known window is the Hamming window whose formula is very similar to that of the Hann window. It is also a 'raised cosine' window and has a similar effect except that it is generally considered to be slightly better. In these notes, we will stay with the Hann window for simplicity. Multiplying {h[n]} by {w[n]} instead of {r[n]} gradually tapers the impulse-response towards zero at the window edges. Consider again the low-pass filter example with cut-off frequency p/3 radians per sampling interval. The ideal impulse-response was found to be:
{h[n]} = {... , 0.14, 0.28, 0.33, .28, .14, ...}
When M = 4, the Hann window {w[n]} = {..,0,..,0, 0.25, 0.75, 1, .75, .25, 0,..,0,..}
Multiplying term by term and delaying by M/2 = 2 samples we obtain the finite impulse-response:
{..,0,..,0, .04, .21,.33,.21,.04, 0,..,0,..}
The resulting "Hann-windowed" FIR filter of order 4 is as shown in Figure 4.1 with a0=0.04, a1=0.21, etc. Its gain-response is as shown in Figure 4.9.
Figure 4.9: 4th order FIR filter (Hann)
4.3 Effect of windowing on frequency response of FIR digital filter
The effect of the Hann window is to gradually reduce the amplitude of the ideal impulse-response towards zero at the edges of the window rather than to abruptly truncate the amplitude as does the rectangular window. The effect on the gain-response of the FIR filter obtained is:
i) to greatly reduce stop-band ripples ( good ).
ii) to reduce the cut-off rate ( bad ).
The phase-response is not affected in the passband. We can improve the cut-off rate by going to higher orders. The graphs below are for 10th and 20th order ( Hann windowed ):
Fig 4.10: Tenth order FIR filter with WC = p/3 ( Hann window )
%MATLAB program to design & graph 10th order FIR lowpass filter with Hann window
w=hanning(11);
for n=1:11; a(n) = 0.3333*sinc((n-6)/3)*w(n);
end;
h = freqz(a,1,1000);
whitebg;
plot([0:999]/1000,20*log10(abs(h)),'k');
axis([0,1,-50,0]);
grid on;
xlabel('Rel_freq / pi');
ylabel('Gain(dB)');
Fig 4.11: 20th order FIR filter with WC = p/3 (Hann window)
A more convenient way of designing such FIR filters is to use signal processing toolbox commands such as FIR1 (see later)
4.4. High-pass, band-pass & band-stop linear phase FIR filters
These can be designed almost as easily as low-pass as long as you remember to define the required gain-response GI(W) correctly from -p to +p; i.e. make GI(-W) = GI(W). For example, a band-pass filter with pass-band from fS/8 to fS/4 would have the following gain response ideally:-
Fig 4.12: Gain response of ideal band-pass filter
Taking f(W) = 0 for all W initially as before, we obtain:
and applying the inverse DTFT gives the following formula (not forgetting negative values of W:
Evaluating the integrals in this expression will produce a nice formula for h[n]. An alternative way of obtaining the same formula is to express
H(ejW) = H1(ejW) - H2(ejW) where:
Now, applying the inverse DTFT, we obtain:
and it then becomes obvious (referring back to the low-pass result earlier) that:
h[n] = (1/2)sinc(n/2) - (1/4)sinc(n/4) for -¥ < n < ¥
Truncating, windowing and delaying this sequence as for the low-pass case earlier produces the causal impulse-response of an FIR filter of the required order.
This will be linear phase for the same reasons as in the low-pass case. Its impulse-response will be symmetric about n=M/2 when M is the order.
Exercise: By a similar method, or otherwise, show that the impulse-response for an ideal zero phase 'brick-wall' high-pass filter with cut-off frequency p/6 radians per sample (i.e. one twelfth of the sampling frequency) is:
h[n] = sinc(n) - (1/6)sinc(n/6) for -¥ < n < ¥
Exercise: Derive the impulse-response for an ideal zero phase 'brick-wall' band-stop digital filter with cut-off frequencies WL and WU radians per sample.
The technique is not restricted to these “conventional” types of brick-wall gain-responses. For example, it is not difficult to design a linear phase filter whose gain-response approximates that below:
Fig 4.13 : Further example of an ideal gain-response.
4.5. Summary of design technique
To design an FIR digital filter of even order M, with gain response G I (W) and linear phase by the windowing method,
1) Set H(ejW) = G I (W) the required gain-response. This assumes f(W) = 0.
2) IDTFT to produce the ideal impulse-response {h[n]}.
3) Window to ±M/2 using chosen window.
4) Delay windowed impulse-response by M/2 samples.
5) Realise by setting multipliers of FIR filter.
Instead of obtaining H( ejW ) = G I ( W ), we get e-jWM/2G(W) with G(W) a distorted version of GI(W) the distortion being due to windowing.
The phase-response is therefore f(W) = -WM/2 which is a linear phase-response with phase-delay M/2 samples at all frequencies W in the range 0 to p. This is because -f(W) / W = M/2 for all W.
Notice that the filter coefficients, and hence the impulse-response of each of the digital filters we have designed so far are symmetric in that h[n] = h[M-n] for all n in the range 0 to M where M is the order. If M is even, this means that h[M/2 - n] = h[M/2 + n] for all n in the range 0 to M/2. The impulse response is then said to be 'symmetric' about sample M/2. The following example illustrates this for an example where M=6 and there are seven non-zero samples within {h[n]}:
{… 0, …, 0, 2, -3, 5, 7, 5, -3, 2, 0, …,0, … }
The most usual case is where M is even, but, for completeness, we should briefly consider the case where M is odd. In this case, we can still say that {h[n]} is 'symmetric about M/2' even though sample M/2 does not exist. The following example illustrates the point for an example where M=5 and {h[n]}therefore has six non-zero sample:
(…, 0,…, 0, 1, 3, 5, 5, 3, 1, 0, …, 0, …}
When M is odd, h[(M-1)/2 - n] = h[(M+1)/2 + n] for n = 0, 1, …, (M-1)/2.
It may be shown that FIR digital filters whose impulse-responses are symmetric in this way are linear phase. We can easily illustrate this for either of the two examples just given. Take the second. Its frequency-response is the DTFT of {h[n]} i.e.
= e-5jW/2 (e5jW/2 +3e3jW/2 +5ejW/2 +5e-jW/2 +3e-3jW/2 +e-5jW/2 )
= e-5jW/2 (e5jW/2 +e-5jW/2 +3e3jW/2 + 3e-3jW/2 +5ejW/2 +5e-jW/2 )
= e-5jW/2 (2cos(2.5W) + 6cos (1.5W) + 10cos(W/2) )
Since the cosine terms are purely real, when their sum is positive we can say that H(ejW) = G(W)ejf(W) with G(W) = 2cos(2.5W) + 6cos 1.5W) + 10cos(W/2) and f(W) = -5W/2. Hence f(W) / W = -5/2 which is constant and this means that H(ejW) is linear phase.
Exercise repeat this calculation for the other example where M=6 which is even.
It is possible to design FIR filters which are not linear phase, and we will consider this later. The technique described in this section is known as the “ windowing ” technique or the "Fourier series approximation technique ".
4.6 Further applications of the windowing design technique for FIR filters
The windowing technique is even more powerful than has been yet indicated; and it is not restricted to linear phase filters. Consider some further examples of its use.
4.6.1. Fractional sampling interval delay filter: An interesting and very useful example is a “fractional sampling interval delay” digital filter whose gain-response is one for all W, but whose phase-response is f(W) = - 0.5 W. Therefore:
ì e -j W / 2 : 0 £ W £ p
H(e j W ) = í
î e j W/ 2 : -p £ W £ 0
Passing any digitised signal through such a filter will not affect its magnitude spectrum. However it will be delayed by half a sampling interval. It will be just as though you had taken each sample of the original analogue signal half a sampling interval later. Applying the inverse DTFT to H(ejW) as defined above produces an impulse-response which must then be windowed and delayed in the usual way. The delay needed to make the impulse-response causal will be added to the half sample delay. So a tenth order filter will delay the input signal by 5.5 samples.