University of Manchester School of Computer Science Comp30291 Digital Media Processing 2009-10 Section 4 ‘Design of FIR

1 / 114

# University of Manchester School of Computer Science Comp30291 Digital Media Processing 2009-10 Section 4 ‘Design of FIR - PowerPoint PPT Presentation

University of Manchester School of Computer Science Comp30291 Digital Media Processing 2009-10 Section 4 ‘Design of FIR digital filters’. x[n]. . z -1. z -1. z -1. z -1. a M. a M-1. a 0. a 1. y[n]. 4.1.Introduction

I am the owner, or an agent authorized to act on behalf of the owner, of the copyrighted work described.

## PowerPoint Slideshow about 'University of Manchester School of Computer Science Comp30291 Digital Media Processing 2009-10 Section 4 ‘Design of FIR' - lamond

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.

- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -
Presentation Transcript

University of Manchester

School of Computer Science

Comp30291

Digital Media Processing 2009-10

Section 4

‘Design of FIR digital filters’

Comp30291: Section 4

x[n]

...

z-1

z-1

z-1

z-1

aM

aM-1

a0

a1

y[n]

4.1.Introduction

FIR digital filter of order M 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]

Comp30291: Section 4

x[n]

z-1

z-1

z-1

z-1

z-1

z-1

z-1

z-1

z-1

z-1

a0

a10

a1

a3

a4

a5

a6

a7

a8

a9

a2

+

+

+

+

+

+

+

+

+

+

y[n]

10th order FIR digital filter

Comp30291: Section 4

Background, objective & methodology

• Its impulse-response is {h[n]} = {..., 0, ..., a0, a1, a2,..., aM, 0, ...}
• Taking DTFT of impulse-response gives the frequency-response:
• Objective is to choose a0, a1,..., aM such that H(ej) is close to some target frequency-response H’(ej) ).
• Inverse DTFT of H’(ej) gives required impulse-response :
• Methodology is to use inverse DTFT to get an impulse- response {h[n]} & then realise some approximation to it.

Comp30291: Section 4

• It is an integral
• It has complex numbers
• Range of integration is from - to  so it involves negative frequencies.

Comp30291: Section 4

x(t)

t

(Have +ve & ve areas)

a

b

Comp30291: Section 4

• Let x = a + j.b, j = [-1]
• Modulus: |x| = [a2 + b2]
• Arg(x)=tan-1(b/a) + {.sign(b) if a < 0}

= tan2(b,a) = angle(a + j.b) : range - to 

• Polar: x = Rej where R = |x| &  = Arg(x)
• De Moivre: ej = cos() + j.sin()
• e-j = cos() - j.sin()
•  ej + e-j = 2cos() & ej - e-j = 2 j.sin()
• Complex conjugate: x* = a - j.b = Re-j

Comp30291: Section 4

• Examine the DTFT formula for H(ej).
• If h[n] real then h[n]ej is complex-conj of h[n]e-j.
• Adding up terms gives H(e-j ) as complex conj of H(ej).
• G() = G(-) since G() = |H(ej)| & G() = H(e-j)|
• (Mod of a complex no. is Mod of its complex conj.)
• (-) = () since() = Arg(H(ej)) & (-) = Arg(H(e-j))
• (Arg of a complex no. is Arg of its complex conj).

Comp30291: Section 4

G()

-

0

-()

-

Gain & phase response graphs with ve frequencies

As G(-) = G(),gain-response is always symmetric about =0

As () = (),phase-response is always anti-symmetric about  = 0

Comp30291: Section 4

G()

1

-

/3

/3

0

4.2. Design of an FIR low-pass filter

Assume we require lowpass filter whose gain-response approximates the ideal 'brick-wall' gain-response:

If we take phase-response () = 0 for all ,

the required target frequency-response is:

Comp30291: Section 4

To complete the calculation:

Comp30291: Section 4

Graph of sinc(x) against x

sinc(x)

1

Main ‘lobe’

‘Ripples’

x

-1

-3

-2

1

2

3

-4

‘Zero-crossings’ at x =1, 2, 3, etc.

Comp30291: Section 4

Plotting this ‘sinc’ graph in MATLAB

clear all; close all; clc;

x = [-10: 0.1 : 10];

y = sinc(x);

figure(1); plot(x,y); grid on;

title('sinc function');

xlabel( 'x'); ylabel('sinc(x)');

legend('sinc(x)', 'Location','Best');

axis([-10 10 -0.3 1]);

Comp30291: Section 4

Graph of sinc(x) against x

‘Main lobe’

‘Zero-crossings’ at x =1, 2, 3, etc

‘Ripples’

Comp30291: Section 4

Plotting this ‘stem’ graph in MATLAB

n=[-20:20]; h = (1/3)*sinc(n/3);

figure(1); stem(n,h,'.:'); grid on;

title('h[n] = (1/3) sinc(n/3)');

xlabel( 'n'); ylabel('h[n]');

legend('(1/3)sinc(n/3)', 'Location','Best');

axis([-20 20 -0.1 0.35]);

for n=-20:20,

disp(sprintf('n:%2d, h[n]: %6.3f' ,n,h(n+21)));

end;

Comp30291: Section 4

‘Ideal’ impulse-response

• Reading from the graph or the MATLAB display, we get:
• {h[n]} = { ..., -0.055, -0.07, 0, 0.14, 0.28, 0.33, 0.28, 0.14, 0, -0.07, -0.055, ... }
• A digital filter with this impulse-response would have exactly the ideal frequency-response we applied to the inverse-DTFT.
• ‘Brick-wall’ low-pass gain response & phase = 0 for all .
• But {h[n]} has non-zero samples extending from n = - to ,
• Not a finite impulse-response.
• Also not causal.
• Not realisable in practice.

Comp30291: Section 4

To produce a realisable impulse-response {h[n]}:

Assume M is order of required filter & that M is even.

(2) Delay resulting sequence by M/2 samples to ensure that the first non-zero sample occurs at n = 0.

Step (1) is ‘truncation’ or ‘windowing’ applied symmetrically about n=0.

Comp30291: Section 4

n

Starting with ‘ideal’ impulse-response:

h[n] = (1/3)sinc(n/3)

Comp30291: Section 4

h[n]

M=10

Truncate to M/2

n

h[n]

Delay by M/2 samples

n

Comp30291: Section 4

Causal & finite impulse-response

• If M=10, finite impulse-response obtained is:
• { ...0, -0.055, -0.07, 0, 0.14, 0.28, 0.33, 0.28, 0.14, 0, -0.07, -0.055, 0... }
• Obtained by truncating & delaying {h[n]} for ‘ideal’ lowpass digital filter with cut-off /3 radians/sample.
• Taking M = 4, we would obtain:
• {.., 0, .., 0, 0.14, 0.28, 0.33 , 0.28 , 0.14 , 0 ,..,0,..}
• Resulting causal impulse-response now realised by setting:
• an = h[n] for n = 0,1,2,...,M.

Comp30291: Section 4

x[n]

-1

z-1

z-1

z

z-1

a4

a3

a1

a0

a2

0.14

0.28

0.33

0.28

0.14

+

y[n]

FIR filter realisation (M=4)

( Note:4th order FIR filter has 4 delays & 5 multipliers ).

• Gain & phase responses may be derived by MATLAB.

Comp30291: Section 4

-6dB

-21 dB

Gain dropsto - dB

/3 Fs/6

FC

Fs/2

Straight line

Jumpsby 180O

Plot gain & phase responses of 4th order FIR filter

Comp30291: Section 4

MATLAB prog for truncating & plotting gain & phase

n=[-20:20];

h = (1/3)*sinc(n/3);

M=4; % FIR filter order

Fs = 8000; % Sampling freq (Hz)

for n=-M/2 : M/2

a(n + M/2 + 1) = h(n+21);

end;

figure(1); freqz(a,1,200,Fs);

Comp30291: Section 4

20

-6dB

0

-21 dB

-20

Gain dropsto - dB

Magnitude (dB)

-40

Fs/2

-60

0

500

1000

1500

2000

2500

3000

3500

4000

/3 Fs/6

FC

Frequency (Hz)

Gain-response of 4th order FIR filter (zoomed)

• Gain starts at 0 dB in pass-band & falls to -6 dB at cut-off frequency.
• There are two ‘stop-band ripples’ & gain is  -21 dB at peak of first.

Comp30291: Section 4

Phase-lag

2

300O

200

100

Fs/2

/2

F Hz

0

1000

FC

2000

3000

4000

Phase-lag response of 4th order FIR filter

• Linear phase in pass-band
• Don’t worry about 180 degree jumps in stop-band for now.

Comp30291: Section 4

()

-

-()

-

Why do we not get a zero phase-response?
• We started by specifying that phase = 0 for all .
• We ended up with a phase-response as follows:

Only pass-band part is shown here

• Phase response affected by the delay of M/2 samples.

Comp30291: Section 4

Estimation of phase-delay from slope
• For a linear phase LTI system, -() /  is ‘phase-delay’ in sampling intervals.
• Look at phase-response graph in pass-band.
• Phase decreases by 180o in 2000Hz (Fs = 8000 Hz)
• -() /     /2 = 2 sampling intervals
• So the phase-delay is 2.
• Not a surprise as we delayed the impulse-response by 2 samples to make it causal (M=4).

Comp30291: Section 4

Using MATLAB functions ‘fir1’ & ‘freqz’
• Same result obtained using ‘fir1’ as follows:

c = fir1(4, 0.33, rectwin(5), 'noscale');

• Reason for rectwin(5) & ‘noscale’ will be clear later.
• To plot gain & phase response:

freqz(c, 1, 500, Fs);

• Plots 500 points & frequency-axis to range 0-Fs/2 !!
• Plots () against  rather than phase lag -().
• Phase ‘unwrapped’ to avoid 360o jumps.

Comp30291: Section 4

• Gain & phase responses different from those originally specified.
• Gain-response: cut-off rate not ‘brick-wall’,
• 0 dB at 0 Hz & drops to -6 dB at cut-off frequency,
• ‘ripples’ and ‘zeros’ appear in stop-band,
• peak of the first ripple at about -21dB.
• Phase-response: not zero for all  as originally specified,
• linear phase in pass-band
• -(  )/ = M/2 for |  | /3;
• phase-delay = M/2 sampling intervals.
• Jumps of 180o occur in stop-band.
• Jumps of 360o avoided by unwrapping.

Comp30291: Section 4

Can we improve low-pass filter by increasing order to ten?

• Taking 11 terms of { (1/3)sinc(n/3)} we get,
• after delaying by 5 samples:
• {...0,-0.055,-.069, 0,.138,.276,.333,.276,.138,0,-.069,-.055,0,...}.
• To obtain same result from MATLAB7:
• c = fir1(10, 0.33, rectwin(11), 'noscale');
• freqz(c);

Comp30291: Section 4

x[n]

z-1

z-1

z-1

z-1

z-1

z-1

z-1

z-1

z-1

z-1

0

-.07

-.07

.28

.28

0.14

0

-.055

0.14

.33

-.055

+

+

+

+

+

+

+

+

+

+

y[n]

10th order FIR digital lowpass filter(/3 rect)

When coeffs are symmetric, the filter is linear phase.

Comp30291: Section 4

Comp30291: Section 4

Comp30291: Section 4

Comp30291: Section 4

Comp30291: Section 4

Comp30291: Section 4

Effect of increasing order on gain & phase responses

• Cut-off rate increases as order (M) increases.
• Number of stop-band ripples increases with order.
• Gain at peak of first ripple after cut-off remains at -21 dB.
• Remains exactly ‘linear phase’ in pass-band but slope of phase-response (=phase-delay) increases since -() /  = M/2.
• NB: freqz ‘unwraps’ phase-response avoiding 360o jumps!
• Filter not really improving with increasing order.
• To improve matters we need to discuss ‘windowing’.

Comp30291: Section 4

rM+1[n]

n

-M/2

M/2

4.3. Windowing

By truncating {h[n]} for ideal filter, we effectively multiplied {h[n]} by a rectangular window sequence {rM+1[n]} to produce {h[n]} where

Comp30291: Section 4

Effect of rectangular window on freq-response

• Sudden transitions to zero at rectangular window edges cause the stop-band ripples we saw in the gain-response graphs.
• Why is this?
• It may be shown that DTFT of {rM+1[n]} is:

Comp30291: Section 4

DTFT of rect-window {rM+1[n]} with M=20

• Purely real
• Looks a bit like sinc.
• Height of main lobe = M+1
• Area of main lobe  2 Width = 4/(M+1)

Comp30291: Section 4

MATLAB program to plot RM+1(ej)

M = 20;

W = -3.2 : 0.001 : 3.2;

if W==0, R=M+1;

else R=sin((M+1)*W/2)./sin(W/2);

end;

figure(1); plot(W,R); grid on;

title(sprintf('Dirichlet K of order %d ', M));

axis([-3.2 3.2 -M*0.3 M+2]);

Comp30291: Section 4

Frequency-domain convolution
• Multiplying two sequences {x[n]} & {[y[n]} to produce {x[n].y[n]} is called ‘time-domain multiplication’.
• It may be shown that if {x[n]} has DTFT X(ej) & [y[n]} has DTFT Y(ej) then the DTFT of {x[n].y[n]} is:
• This is frequency domain convolution.
• ‘Multiply, find the area then scale by 1/ 2’
• Time-domain multiplication  freq-domain convolution.

Comp30291: Section 4

Apply this formula to rectangular window
• Multiplying an ideal impulse-resp {h[n]} by the rect window {rM+1[n]} causes H(ej) to be ‘convolved’ with RM+1(ej).
• If H(ej) has ideal brick-wall gain-response, & RM+1(ej) is as shown in previous graph, convolving H(ej) with RM+1(ej) reduces cut-off rate & introduces stop-band ripples.
• The sharper the main-lobe of RM+1(ej) and the lower the ripples, the better.

Comp30291: Section 4

Illustration for ideal /3 lowpass filter
• What happens to area under curve as  increases from 0 to  ?

Comp30291: Section 4

=0

Area under curve /3 to +/3 when =0

RM+1(ej)

 2



/3

/3

Comp30291: Section 4

RM+1(ej)



/6

/2

=/6

Area under curve /3 to +/3 when =/6

 2

Comp30291: Section 4

RM+1(ej)



2/3

0

=/3

Area under curve /3 to +/3 when =/3

 

Comp30291: Section 4

RM+1(ej)



/6

5/6

Area under curve /3 to +/3 when =/2

=/2

 0

Comp30291: Section 4



/3

=2/3

Area under curve /3 to +/3 when =2/3

RM+1(ej)

 0

Comp30291: Section 4

Conclusions from these graphs
• Ripples arise from freq-domain convolution between ideal freq-response & RM+1(ej).
• At =0, area is 2 so gain is 1 which is 0 dB.
• At  = C (= /3), area is half of 2, so gain  0.5 ( 6 dB) as half of main lobe lies between   /3 &  + /3.
• As  is increased from /3, main lobe no longer included,so area comes from the ripples.
• Ripples are small & alternately +ve & negative, so area & therefore the gain will be small.
• But area will increase & decrease (ripple) as  increases.

Comp30291: Section 4

wM+1[n]

n

M/2

M/2

4.4. Non-rectangular windows

• Levels of ripples reduced if {rM+1[n]} replaced by non-rectangular window sequence { wM+1[n] }.
• Produces a more gradual transition at the window edges.
• Simple non-rectangular window sequence is Hann window
• It’s a ‘raised cosine’ with M+1 non-zero samples centred on n=0.

Comp30291: Section 4

Hann window w21[n]

Comp30291: Section 4

Hann window w41[n]

Comp30291: Section 4

‘Hamming’ & other ‘Hann’ windows
• Slightly different formulae exist for the Hann window.
• In practice, the difference is usually unimportant.
• MATLAB has 2 functions: ‘hann’ & ‘hanning’.
• ‘hanning(M+1)’ gives our Hann window {wM+1[n]}
• But it’s shifted by 1+M/2 samples to start at n=1.
• Best known non-rectangular window is ‘Hamming’.
• Its name & its formula are fairly similar to ‘Hann’:

Comp30291: Section 4

Hamming window w21[n]

Comp30291: Section 4

Hamming window w41[n]

Comp30291: Section 4

Effect of non-rectangular windows

• Multiplying {h[n]} by {wM+1[n]} instead of {rM+1[n]} gradually tapers impulse-response towards zero at window edges.
• To understand why this reduces stop-band ripples, compare DTFT of {wM+1[n]} with DTFT of {rM+1[n]} with M=20.

Comp30291: Section 4

DTFT of Hann window {wM+1[n]}
• Height of main lobe = 1 + M/2 ( half that of rect window)
• Area of main lobe  2 (same as rect window)
• Width of main lobe 8/M ( twice that of rect window)
• It may be shown that:

WM+1(ej) = 0.5RM+1(ej) + 0.25RM+1(ej(-/(1+M/2))) + 0.25RM+1(ej(+/(1+M/2)))

• Purely real only for a window that is symmetric about =0.

Comp30291: Section 4

DTFT of Hamming window {hwM+1[n]}
• Height of main lobe = 1 + 0.54M ( half that of rect window)
• Area of main lobe  2 (same as rect window)
• Width of main lobe  8/M ( twice that of rect window)
• It may be shown that:

HWM+1(ej) = 0.54RM+1(ej) + 0.23RM+1(ej(-/(1+M/2))) + 0.23RM+1(ej(+/(1+M/2)))

• Let’s stay with Hann for now as the differences are small.

Comp30291: Section 4

To plot DTFT of Hann & Rect windows

clear all; close all; clc

M=20;

W = -3.2:0.001: 3.2; W=W+0.0000001; % Sorry

R=sin((M+1)*W/2)./sin(W/2);

W1=W-pi/(1+M/2); W2=W+pi/(1+M/2);

R1=sin((M+1)*W1/2)./sin(W1/2); R2=sin((M+1)*W2/2)./sin(W2/2);

Hann=0.5*R + 0.25*R1+0.25*R2;

figure(1); plot(W,Hann,W,R); grid on;

title(sprintf('DTFT of Hann & Rect windows, order %d ', M));

axis([-3.2 3.2 -M*0.3 (M+2)]);

legend('Hann','Rect','Location','Best');

Comp30291: Section 4

Plot on dB scale

20.log10(abs(DTFT))

Comp30291: Section 4

Comparing DTFT of Rect & Hann windows

• Ripples in W21(ej) greatly reduced in comparison to R21(ej).
• This is good!
• Main lobe W21(ej) less sharp & lower in comparison to R21(ej).
• Reduces sharpness of the cut-off
• Price paid for reducing stop-band ripples.

Comp30291: Section 4

Ideal /3 lowpass filter with Hann window
• As before, consider what happens to area under curve as  increases from 0 to , but now with Hann window ?

Comp30291: Section 4

=0

Area under curve /3 to +/3 when =0

WM+1(ej)

DTFT of Hann window

 2



/3

/3

Comp30291: Section 4

=/3

Area under curve /3 to +/3 when =/3

WM+1(ej)

DTFT of Hann window

 



0

2/3

Comp30291: Section 4

=/2

Area under curve /3 to +/3 when =/2

WM+1(ej)

DTFT of Hann window

 0.1



/6

5/6

Comp30291: Section 4

DTFT of Hann window

 0



=2/3

Area under curve /3 to +/3 when =2/3

WM+1(ej)

/3

Comp30291: Section 4

Conclusions from these graphs
• Have plotted freq-domain convolution between an ideal freq-response & DTFT of Hann window.
• At =0, area is 2 so gain is 1 as for rect window.
• At  = /3, area is  (half of main lobe), so gain  0.5 (-6dB) as for rect window.
• As  is increased from /3, area of main lobe disappears more gradually than for rect window.
• Ripples of DTFT of Hann window are much smaller than for rect window, so ripples in the area & therefore the gain will be much reduced.
• Slower cut-off (bad) but reduced stopband ripples (good)

Comp30291: Section 4

Applying 4th order Hann window

• Consider again low-pass filter with cut-off /3 radians/sample
• Ideal impulse-response was found to be:
• {h[n]} = { …..... , 0.14, 0.28, 0.33, 0.28, 0.14, ………}
• When M = 4, Hann window is:
• {w5[n]} = {..,0,..,0, 0.25, 0.75, 1, 0.75, 0.25, 0,..,0,..}
• Multiplying term by term & delaying by M/2 = 2 samples we get:
• {..,0,..,0, 0.04, 0.21, 0.33, 0.21, 0.04, 0,..,0,..}

Comp30291: Section 4

x[n]

-1

z-1

z-1

z

z-1

a4

a3

a1

a0

a2

0.04

0.21

0.33

0.21

0.04

+

y[n]

Resulting ‘Hann-windowed’ FIR filter of order 4:

Its gain-response is shown on next slide.

Comp30291: Section 4

Effect of increasing order to 10

• To design 10th order FIR lowpass filter with Hann window & cut-off frequency /3 ( Fs/6):

clear all;

M=10

for n= -M/2 : M/2

w = 0.5 *(1+cos(n*pi/(1+M/2)) );

c(1+n+M/2) = (1/3)*sinc(n/3)*w;

end;

Fs = 8000; freqz(c,1,500,Fs);

axis([0 4000 -60 0]);

Comp30291: Section 4

x[n]

z-1

z-1

z-1

z-1

z-1

z-1

z-1

z-1

z-1

z-1

0

-.017

-.017

.257

.257

0.1

0

-.004

0.1

.333

-.004

+

+

+

+

+

+

+

+

+

+

y[n]

10th order FIR lowpass (/3, Hann)

Since coeffs are symmetric, filter is exactly linear phase.

Comp30291: Section 4

Use of MATLAB function ‘fir1’

• ‘fir1’ uses ‘windowing method’ we have just seen.
• By default it uses a Hamming window.
• Also scales coeffs to make gain exactly 0 dB at 0 Hz
• To design 10th order FIR lowpass filter with Hamming window & cut-off frequency /3 ( Fs/6):
• c=fir1(10, 0.33);

Comp30291: Section 4

x[n]

z-1

z-1

z-1

z-1

z-1

z-1

z-1

z-1

z-1

z-1

0

-.01

-.01

.25

.25

0.1

0

-.005

0.1

.33

-.005

+

+

+

+

+

+

+

+

+

+

y[n]

10th order FIR lowpass (/3, Hamming)

Again, coeffs are symmetric so filter is exactly linear phase.

Comp30291: Section 4

• Effect is to gradually reduce amplitude of ideal impulse-response towards zero at edges of window rather than to abruptly truncate.
• Effect on gain-response of FIR filter obtained is:
• i) to greatly reduce stop-band ripples ( good ).
• ii) to reduce the cut-off rate ( bad ).
• Phase-response is not affected in the pass-band.
• We can improve the cut-off rate by going to higher orders.
• Graphs next are for Hann window with M=10, 20 & 60:

Comp30291: Section 4

Effect of Hann window on first stop-band ripple

• At frequency of first stop-band ripple:
• since 20log10(1/10) = 20 dB,
• amplitude reduced by factor 10 with rect window.
• since 20 log10(1/100) = 40 dB,
• amplitude reduced by factor >100 with Hann window.
• i.e. 5 sin(t) becomes 0.5 sin(t) or 0.05 sin(t) .
• What do we do if this is not good enough?
• Answer, use a different window (e.g. Hamming, Kaiser)
• Hamming window similar to Hann but slightly better.

Comp30291: Section 4

4.5. Kaiser window

• Offers a range of options from rectangular, through Hamming towards even lower stop-band ripples.
• Ripple reduction at expense of less sharp cut-off rate.
• MATLAB command:
• KW = kaiser(M+1,beta)
• gives a Kaiser window array of length M+1 for any value of beta > 0.
• When beta () = 0, this is a rectangular window & when beta = 5.4414 we get a Hamming window.
• Increasing beta further gives further reduced stop-band ripples with a reduced cut-off sharpness.

Comp30291: Section 4

Comp30291: Section 4

Comp30291: Section 4

Comp30291: Section 4

Comp30291: Section 4

My MATLAB test program

clear all;

beta = 5; N=60;

kw = kaiser(N+1,beta);

hw=hamming(N+1);

a=fir1(N, 0.33, rectwin(N+1), 'noscale');

for n=1:N+1

akw(n)=a(n)*kw(n); ahw(n) = a(n)*hw(n);

end;

figure(1); freqz(akw); grid on;

figure(2); freqz(ahw); grid on;

Comp30291: Section 4

4.6. Highpass, bandpass & bandstop linear phase FIR filters

• Can be designed almost as easily as low-pass.
• Remember to define required gain-response G() from - to +
• Make G(-) = G().
• Band-pass filter with pass-band from FS/8 to FS/4 has following gain response ideally:-

Comp30291: Section 4

FIR band-pass filter (/4 to /2)

• Taking () = 0 for all  initially as before, we obtain:
• Applying inverse DTFT (not forgetting negative ):
• Can evaluate this, apply window & delay as before.
• Use MATLAB to get {h[n]}

Comp30291: Section 4

• c = fir1(10,0.33,'high') designs a 10th order high-pass filter.
• c = fir1(10,[0.2 0.4],'bandpass') designs a 10th order band-pass filter with cut-off frequencies 0.2 and 0.4.
• c = fir1(20,[0.2 0.4],'stop') designs a 20th order band-stop filter with cut-off frequencies 0.2 and 0.4.
• By default ‘fir1’ uses a ‘Hamming’ window (very similar to a Hann) and scales the pass-band gain to 0 dB.
• Linear phase response obtained.

Comp30291: Section 4

4.7. Summary of ‘windowing method’

• To design FIR filter of even order M, with gain-response approximating G() & linear phase,
• 1) Set H(ej) = G(). This assumes () = 0.
• 2) I-DTFT to produce ideal impulse-response {h[n]}.
• 3) Window to M/2 using chosen window.
• 4) Delay windowed impulse-response by M/2 samples.
• Realise by setting multipliers of FIR filter.
• MATLAB routine fir1 does all this for LP, HP, BP, BS

Comp30291: Section 4

Properties of resulting FIR digital filter

• Instead H(ej) = ej0G(  ), we get H(ej) = e-jM/2G()
• G() is distorted version of G() due to windowing.
• Phase-response is () = -M/2 in pass-band
• Linear phase with phase-delay: -() /  = M/2 samples.
• Filter coeffs are symmetric about M/2.
• e.g. {…2, -3, 5, 7, 5, -3, 2, …} M =6 (even)
• {…, 1, 3, 5, 5, 3, 1, …} M=5 (odd)
• This is because the FIR filter is linear phase.

Comp30291: Section 4

Demonstrate thatan FIR filter whose impulse-response is symmetric is linear phase.

• Let {h[n]} = {…, 1, 3, 5, 5, 3, 1, …} (symmetric)
• = e-5j/2 (e5j/2 +3e3j/2 +5ej/2 +5e-j/2 +3e-3j/2 +e-5j/2 )
• = e-5j/2 (2cos(2.5) + 6cos (1.5) + 10cos(/2) )
• = G()ej() with () = -5/2.
• Hence () /  = -5/2 = constant, so H(ej) is linear phase.

Comp30291: Section 4

4.8 Using MATLAB SP Toolbox
• See notes

Comp30291: Section 4

4.9. Remez Exchange Algorithm method

• Better than windowing technique, but more complicated.
• Available in MATLAB.
• Design 40th order FIR lowpass filter whose gain is unity (0 dB) in range 0 to 0.3 radians/sample & zero in range 0.4 to .
• a = remez (40, [0, 0.3, 0.4,1],[1, 1, 0, 0] );
• The 41 coefficients will be found in array ‘a’.
• Produces 'equi-ripple' gain-responses where peaks of stop-band ripples are equal rather than decreasing with increasing frequency.
• Highest peak in stop-band lower than for FIR filter of same order designed by windowing technique to have same cut-off rate.
• There are 'equi-ripple' pass-band ripples.

Comp30291: Section 4

Comp30291: Section 4

4.10. Fixed point implementatn of FIR digital filters

• FIR digital filters often implemented in mobile equipment.
• Low power fixed point DSP processors are used.
• Typically with a basic 16-bit word-length.
• Must be programmed using only integer arithmetic.
• Take 4th order FIR filter with impulse response:
• {….. 0.04, 0.21, 0.33, 0.21, 0.04, …...}.
• Rounding each coeff to nearest integer clearly a mistake.
• Multiply each coeff by a large constant, e.g. 100, then round: .
• A0 = 4, A1 = 21, A2 = 33, A3 = 21 , A4 = 4.
• We must divide the output by same constant.
• Instead of 100, choose a power of two for the constant.

Comp30291: Section 4

Fixed point implementation (cont)

• Dividing by a power of two (e.g. 1024) is very simple
• It’s just an ‘arithmetic right-shift’ operation.
• Available instruction on DSPs.
• The larger the constant, the more accurate the coefficients.
• Careful not to choose too large a constant
• If integers produced get too large, we risk overflow
• Difficult balancing act between inaccuracy & overflow.
• If constant is 2^10 (=1024), rounded 4th order coeffs become:
• A0 = 35 A1 = 212 A2 = 341 A3 = 212 A4 = 35
• MATLAB prog on next slide uses integer arithmetic
• Ready to be ported to a DSP

Comp30291: Section 4

4th order low-pass filter using integer arithmetic only

A = [35 212 341 212 35] ;

x = [0 0 0 0 0 ] ;

while 1

x(1) = input( 'X = ');

Y = A(1)*x(1);

for k = 5 : -1: 2

Y = Y + A(k)*x(k);

x(k) = x(k-1);

end;

Y = round( Y/1024); %Arith right-shift 10 places

disp([' Y = ' num2str(Y)]);

end;

Comp30291: Section 4

4.11. Advantages of FIR filters compared with IIR

• FIR filters easy to program in fixed point arithmetic.
• Never become unstable as there is no feedback.
• Can be exactly linear phase
• In some cases, overflows can be allowed to occur since if gain is never greater than 1, you know that a +ve overflow will eventually be cancelled out by a ve overflow or vice versa.
• This works if you do not use ‘saturation mode’ arithmetic which avoids ‘wrap-round’.
• Can risk overflow more readily with FIR digital filters than with IIR digital filters, & thus have greater coefficient accuracy.
• Disadvantage: FIR need higher orders than IIR (later).

Comp30291: Section 4

DSP language: ‘Q-format’

• Scaling by 1024, is adopting a 'Q-format' of ten.
• Programmer assumes a binary point to exist ten bit positions
• from the right within the 16-bit word.

Comp30291: Section 4

PROBLEMS

1. Design a 100th order FIR low-pass digital filter with cut-off atfS/4 with & without a Hann window. Use MATLAB to compare gain responses obtained.

2. Design 10th order FIR bandpass filter with cut-off frequencies at /4 &/2.

3. Write MATLAB program for one of these filters using integer arithmetic only.

4. Design a 4th order FIR high-pass filter with cut-off at /3.

5. Do all FIR filters have exactly linear phase responses?

6. Explain why IIR digital filters cannot have exactly linear phase responses.

7. Explain why it is impossible to implement an ideal ‘brick-wall’ LPF.

8. Show that (  ) = -k corresponds to a delay of k samples.

9. Rearrange even order linear phase FIR filter to reduce no. of multipliers.

10. Do FIR filters have passband ripples as well as stopband ripples?

11. Explain why FIR filters have stopband ripples where IIR filters do not.

12. Explain why the gain of an FIR low pass filter with 0 dB gain at zero frequency reduces to -6 dB at its cut-off frequency.

Comp30291: Section 4

h[n]

If h[n] symmetric & IIR it must be non-causal

n

M

Questions 10 & 6

10. Yes but you can hardly see them.

6. For linear phase, impulse-response must be symmetric about some value of n, say n=M. If it is an IIR it goes on for ever as n  .

So it must go on for ever backwards as n  -.

Would have to be non-zero for values on n<0; i.e. non-causal.

Comp30291: Section 4