Digital Filter Design

Chapters 9 and 10 Appendix A and Appendix B Dr. Iyad Jafar

Outline y Introduction I t d ti y Common C Id l Filter Ideal Filt T Types y Selecting Filter T Type pe y Digital IIR Filter Design y Digital FIR Filter Design 2

Introduction y What is filtering? y Remove unwanted frequency content in the signal y For LTI systems, systems this can be done using convolution in

the with the impulse response of the filter h[n] in the time domain or multiplication in the frequency domain with the frequency response Filter h[n]

3

x[n]

y[n]

Introduction y Digital filter design is the process of identifying the

impulse h[n] response or the transfer function H(z) of a system that is capable of passing certain frequency content of an input signal x[n] and blocking others with the following considerations y H(z) is realizable y H(z) is causal y H(z) H( ) meets the h frequency f response specifications f y H(z) meets the phase response specifications y Complexity of H(z) y Is H(z) an IIR or FIR system? y 4

Common Ideal Filter Types y The simplest way to design filters is to consider the

frequency response in the ideal case y Frequency response of an ideal digital lowpass filter Passband

H(ejω) H(

Stopband

1

-π

-ω ωc

ωc

ω π Cutoff frequency

5

Common Ideal Filter Types y Frequency response of an ideal digital highpass filter |H(ejω)| 1

-π

-ωc

ωc

ω π

y Frequency response of an ideal digital band band-pass pass filter |H(ejω)| 1

6

-π

-ωc1 -ωc2

ωc1 ωc2

ω π

Common Ideal Filter Types y Frequency response of an ideal digital band-reject filter |H(e | ( jω)| 1

-π

-ω ωc2

-ωc1

ωc1

ωc2

ω π

y These filter are easy to specify ! y They have unity magnitude in the passband and zero in

the h stopband, b d andd they h have h zero phase h (real ( l functions!) f i !)

7

y Problems?!!!

Common Ideal Filter Types y Example 1 – the impulse response of an ideal lowpass filter 0 05 0.05

1

0.6

h[n]

H(ejZ)

0.8

0

0.4 0.2 0 -4

-0.05 -2

0

Z

2

4

-200

-100

0

n

100

y The impulse response is y Non causal, causal not absolutely summable, summable and doubly infinite??

8

Î Not realizable or stable!? y Solution y Magnitude Approximation and ignore the phase!

200

Magnitude Approximation of Digital Filters y Magnitude g approximation pp

includes y The magnitude g response p specifications in the stop and pass bands are given with acceptable tolerance l

y The transition between the pass

and stop bands spans a range of frequencies 9

Magnitude Approximation of Digital Filters

10

Magnitude Approximation of Digital Filters y Alternatively, y, magnitude g specifications p mayy byy ggiven in

normalized form where the maximum magnitude in the ppassband is assumed to be unityy

11

Magnitude Approximation of Digital Filters

12

Magnitude Approximation of Digital Filters Example 2 2

13

Selection of Filter Type y The desired filter could be IIR or FIR! y For IIR filters,, the transfer function should be a real rational

function of z-1

In this case, H(z) should be stable with lowest order N to reduce computational complexity y For FIR filters,, the transfer function is a p polynomial in z-1 with y real coefficients

p y N should be the smallest For reduced complexity, 14

Selection of Filter Type y Comparison between FIR and IIR

15

Digital IIR Filter Design y Given the digital filter specifications (δp, δs, ωc, N, FT….), we

want to find

that approximates the magnitude response, stable, with lowest order N. y The most common approach is to convert the digital filter

specifications f into analog l lowpass l specifications f to determine d the h transfer function H(s) of the analog lowpass filter. Then, the obtained H(s) is converted into the digital filter transfer function H(z). y Why this approach? y Analog approximation techniques are highly advanced y Theyy result in closed-form solution 16

y Many applications require the digital simulation of analog filters

Digital IIR Filter Design Digital g ta Filter te Specifications

Convert to Analog Filter S Specifications ifi i Using i Bilinear Transformation

Spectral Transformation into to aan Analog a og Lowpass owpass Filter Prototype (If needed)

Convert H(s) into G(z) U i Bili Using Bilinear Transformation

Spectral Transformation to the Required q Analog g Filter (If needed)

Design g the analog g Lowpass Filter H(s)

• We need to know how to • Design lowpass analog filters • Perform spectral transformation between different types of 17

analogg filters • Perform bilinear transformation

Analog Filter Specifications y Typical magnitude response of an analog lowpass filter

18

Analog Filter Specifications y Magnitude specifications can be given in normalized form

19

Analog Lowpass Filter Prototypes y Butterworth Approximation

20

Analog Lowpass Filter Prototypes y Butterworth Approximation

21

Analog Lowpass Filter Prototypes y Butterworth Approximation

22

Analog Lowpass Filter Prototypes y Butterworth Approximation

E ample 3 Example

23

Analog Lowpass Filter Prototypes y Butterworth Approximation y Example 3 – continued

y In Matlab, Matlab this can be obtained using

24

[N,Wc] = buttord(1000,5000,1,40,'s') N=4 Wc = 1581 Hz

Analog Lowpass Filter Prototypes y Butterworth Approximation

25

Analog Lowpass Filter Prototypes y Butterworth Approximation y Example 4 - Determine the transfer function H(s) of a

analog butterworth lowpass filter of order 2 and :c 2 —1000 rad / sec y Solution N 2 ⇒ Two roots U1 and U2

U1

2 —1000 — e jS (2 21)/4

1000  j1000

U2

2 —1000 — e jS (2 41)/4

1000  j1000

H a (s) ( )

26

( 2 —1000) 2 (s  1000  j1000)(s  1000  j1000) H a (s)

2 —106 s 2  2000s  2 —106

Analog Lowpass Filter Prototypes y Butterworth Approximation y Example 4 – continued

In Matlab (N = 2, :c

27

2 —1000 rad / sec )

% compute the numerator and denumerator of HLP(s) >> [num,den] = butter(2,sqrt(2)*1000,'s') ; % determine the transfer function >> tf(num,den) 2e006 -------------------s^2 + 2000 s + 2e006 % plot the magnitude of the frequency response >> W = 0: 0.1 : 6000 ; >> H = freqs(num, eqs( u , den, de , W); ); >> plot(W,20*log10(abs(H));

Analog Lowpass Filter Prototypes y Butterworth Approximation y Example 4 – continued

In Matlab 0 -5

|H(: )|| (dB)

-10 -15 15 -20 -25 -30 0 28

1000

2000

3000 4000 : (rad/sec)

5000

6000

Analog Lowpass Filter Prototypes y Other Approximations

29

Design of Other Analog Filters y How about highpass, bandpass, band-stop filter? Do

we need to specify magnitude approximations for these filters? y Given the analog lowpass filter prototypes, can we design other h ffilters? l ? y We consider frequency or spectral transformations! (Appendi B) (Appendix

30

y Basically, y The filter specifications are mapped into a lowpass filter specifications. y The lowpass p filter is designed g y The inverse spectral transformation is applied to obtain the desired filter y We consider here the design of analog highpass filter as an example. (refer to Appendix B for other filters).

Design of Other Analog Filters

31

Design of Analog Highpass Filter

32

Design of Analog Highpass Filter y Example 5. Design an analog highpass filter such that y Stopband edge frequency of 2 Krad/sec with minimum attenuation of 5 dB y Passband edge frequency 4 Krad/sec with maximum attenuation of 1 dB y Solution

* Let : p

:

ˆ :p: p  ˆ :

ˆ maps to : 1 rad/sec, then : s s

2 rad/sec for

lowpass filter * Now, we compute N and :c for the lowpass filter as we did in Example 3. This gives N

2 and :c

1.6493 rad / sec

* The transfer function is obtained as in Example 4 to be

H LP (s) 33

2.72 .7 s ^ 2  2.332 s  2.72

Design of Analog Highpass Filter y Example 5. Continued

ˆ by substituting * Now, we transform H LP (s) into H D (s) ˆ :p: 4000 p s sˆ sˆ * This gives the desired transfer function of the analog highpass filter

H D (sˆ )

34

sˆ 2  1.108 —1012 sˆ  2.97 — 10025 sˆ 2  3430 sˆ  5.882 —106

Design of Analog Highpass Filter y Example 5. Continued

In Matlab % Step-by-step approach wp_HP = 4000 ; ws_HP = 2000 ; wp_LP = 1 ; ws_LP LP = 1 * wp_HP HP / ws_HP HP ; [N,Wc] = buttord(wp_LP,ws_LP,1,5,'s'); [numLP,denLP] = butter(N,Wc,'s'); [numHP denHP] = lp2hp(numLP [numHP,denHP] lp2hp(numLP,denLP,4000); denLP 4000); tf(numHP,denHP) W = 0:0.1:6000 ; H HP = freqs(numHP,denHP,W); H_HP freqs(numHP denHP W); plot(W,20*log10(abs(H_HP)),'linewidth',2) axis([0 max(W) -20 0]) xlabel('\Omega','fontweight','b','fontsize',14) xlabel( \Omega , fontweight , b , fontsize ,14) ylabel('|H_D(\Omega)| (dB)','fontweight','b','fontsize',14) (g , , , g , ) set(gca,'fontsize',14,'fontweight','b') grid on 35

% Direct approach [N,Wc] = buttord(4000,2000,1,5,'s') [numHP denHP] = butter(N [numHP,denHP] butter(N,Wc,'high','s'); Wc 'high' 's'); tf(numHP,denHP) W = 0:0.1:6000 ; hHP = freqs(numHP,denHP,W); freqs(numHP denHP W); plot(W,20*log10(abs(hHP)),'linewidth',2) axis([0 max(W) -20 0]) xlabel('\Omega' xlabel( \Omega ,'fontweight' fontweight ,'b' b ,'fontsize' fontsize ,14) 14) ylabel('|H_D(\Omega)| (dB)','fontweight','b','fontsize',14) set(gca,'fontsize',14,'fontweight','b') set(gca, fontsize ,14, fontweight , b ) grid on

Design of Analog Highpass Filter y Example 5. Continued 0

-10

D

|H H (: )| (dB)

-5

-15

-20 0 36

1000

2000

3000 4000 : (rad/sec)

5000

6000

Bilinear Transformation

37

Bilinear Transformation

38

Bilinear Transformation

Z

3.1416

0

-3.1416

39

-10

-5

0

:

5

10

Bilinear Transformation

40

Bilinear Transformation y Example 6. Design a digital lowpass IIR filter with the

following specifications y ωp = 0.2π rad/sample with maximum passband ripple of 0.5

dB y ωs = 0.6π 0 6 rad/sample d/ l with ith minimum i i stopband t b d attenuation tt ti off 15 dB

y Solution

(1) Prewrap frequencies to obtain the analog lowpass filter specs : p tan(Zp / 2) tan(0.2S / 2) 0.3249 rad/sec :s

41

tan(Zs / 2)

tan(0 tan(0.6 6S / 2) 1.3764 1 3764 rad/sec

Bilinear Transformation y Example 6. Continued

(2) Compute N and :c of the lowpass filter ⎛ 1 ⎞ 2 20 log10 ⎜  0.5 ⇒ H 0.1220 ⎟ 2 ⎝ 1 H ⎠ ⎛1⎞ 20 log10 ⎜ ⎟ 15 ⇒ A 2 31.6228 ⎝A⎠ ⎡ 1 log l 10 ((A 2  1) / H 2 ) ⎤ N= ⎢ ⎥ ⎢ 2 log10 (:s / : p ) ⎥ 1 2 H a ( j: s ) 1  (:s / :c ) 2N 42

⎡ ⎛ 1/ k1 ⎞ ⎤ ⎢log10 ⎜ 1/ k ⎟ ⎥ ⎝ ⎠⎥ ⎢ 1 ⇒ :c 2 A

2

0.5804 rad/sec

Bilinear Transformation y Example 6. Continued

(3) Find H a (s) H a (s)

: cN N

¸ (s p )

0.58042 (s  0.4104 0 4104 + 0.4104i)(s 0 4104i)(s 0.4104 0 4104 - 00.4104i) 4104i)

k

k 1

H a (s)

0.3396 0 3396 s 2  0.8202s  0.3369

(4) Find H(z) H(z)

43

H a (s) s

1 z

1

1 z 1

0.1561  0.3122 z 1  0.1561z 2 1  0.6147z 0 6147z 1  0.2393z 0 2393z 2

(5) The system as a difference equation y[n]  00.6147y[n 6147y[n  1]  0.2393y[n 0 2393y[n  2] 0.561x[n]  0.3122 x[n  1]  0.1561x[n  2]

X(z) Y(z)

Bilinear Transformation y Example 6. Continued

((6)) Test the stability y of the system y Æ Examine the locations of the poles. Since all poles are inside the unit circle, then we can define a stable and causal system 1

I Imaginary Part

0.5

2

0

-0.5

-1 44

-1

-0.5

0 Real Part

0.5

1

Design of Analog Highpass Filter y Example 6. Continued

In Matlab

45

% Step-by-step approach Wp = tan(0.2*pi/2); Ws = tan(0.6 tan(0 6*pi/2); pi/2); [N,Wc] = buttord(Wp,Ws,0.5,15,'s') [numA,denA] = butter(N,Wc,'s') [numD denD] = bilinear(numA,denA,0.5) [numD,denD] bilinear(numA denA 0 5) w = 0:0.001:pi; h = freqz(numD,denD,w); figure Zplane(numD,denD) figure pplot(w,20*log10(abs(h)),'linewidth',2) ( , g ( ( )), , ) axis([0 max(w)+1 -20 .1]) xlabel('\omega','fontweight','b','fontsize',14) ylabel('|H(z)| (dB)','fontweight','b','fontsize',14) set(gca,'fontsize',14,'fontweight','b') grid on

% Direct approach [[N,wc] = buttord(0.2,0.6,0.5,15,'z') [numD denD] = butter(N,wc,'z') [numD,denD] butter(N wc 'z') w = 0:0.001:pi; h = freqz(numD,denD,w); figure zplane(numD,denD) figure plot(w 20*log10(abs(h)) plot(w,20 log10(abs(h)),'linewidth' linewidth ,2) 2) axis([0 max(w)+1 -20 .1]) xlabel('\omega','fontweight','b','fontsize',14) ylabel('|H(z)| ylabel( |H(z)| (dB)','fontweight','b','fontsize',14) (dB) , fontweight , b , fontsize ,14) set(gca,'fontsize',14,'fontweight','b') grid on

Design of Analog Highpass Filter y Example 6. Continued 0

Z

|H H(ej )| (dB)

-5

-10

-15

-20 0 46

02 0.2

04 0.4

0.6 0 6 Z (/S )

08 0.8

1

Digital FIR Filter Design y For FIR filters, the FIR transfer function is a

polynomial in z-1 with real coefficients

For reduced complexity, complexity N should be the smallest y There are many techniques for the design of digital FIR filters y Windowed Fourier series y Frequency sampling y Computer-based p 47

Digital FIR Filter Design

48

Digital FIR Filter Design

49

Digital FIR Filter Design

50

Digital FIR Filter Design

51

Digital FIR Filter Design

52

1

1

0.8

0.8

0.6

0.6

hd[n]

Hd(ejZ)

Digital FIR Filter Design

0.4

0.4

0.2

0.2

0 -4

0 -2

0

Z

2

4

-100 100

1

Delay by M samples ht[n]]

ht[n n]

0.6 0.4

n

50

100

5

10

Truncation 2M+1 samples

0.8 0.6 0.4

0.2

0.2

0

53

0

1

0.8

0

-50 50

5

10

n

15

20

0 -10

-5

0

n

Digital FIR Filter Design y Impulse response of typical filters 1

Hd((ejZ)

08 0.8 0.6 0.4 0.2 0

54

-3 3.14 14

-wc wc

0

Z

wc

3 14 3.14

Digital FIR Filter Design y Impulse response of typical filters 1

Hd(ejZ)

08 0.8 0.6 0.4 02 0.2 0

55

-3.14

-wc

0

Z

wc

3.14

Digital FIR Filter Design y Impulse response of typical filters 1

Hd(ejZ)

0.8 0.6 0.4 0.2 0

56

-3.14-wc2

-wc1

0

Z

wc1

wc2 3.14

Digital FIR Filter Design y Impulse response of typical filters 1

Hd((ejZ)

08 0.8 0.6 0.4 0.2 0

57

-3 3.14 14-wc2 wc2

-wc1 wc1

0

Z

wc1

wc2 3.14 3 14

Digital FIR Filter Design y How to perform truncation? y Truncation of the target impulse response ht[n] is represented

by multiplying the desired impulse response hd[n] with a window function w[n]

y The simplest window function is the rectangular window 1

wR[n]

0.8 0.6 0.4 0.2 0

58

-M

0

n

y Any problems?

M

Digital FIR Filter Design y Recall that multiplication in time is equivalent to

q y convolution in frequency

y Thus, Ht(ejω) is obtained by periodic convolution between

Hd(ejω) and ψ(ejω) 1

1

0.8

WR(ejZ)

Hd(ejZ)

0.8 0.6 0.4

59

0.4 0.2

0.2 0 -4

0.6

0

-2

0

Z

2

4

-0.2 02 -50

0

Z

50

Digital FIR Filter Design 0.8

Ht(ejZ)

0.6 0.4 02 0.2 0 -0.2 -15

60

-10

-5

0

Z

5

10

15

Digital FIR Filter Design

61

Digital FIR Filter Design

62

Digital FIR Filter Design

63

Digital FIR Filter Design

64

Digital FIR Filter Design y Estimating the filter order y Kaiser’s Formula

N#

20 log10 ( Gp Gs )  13 14 6(Z s Z 14.6( Zp ) / 2S

y Bellanger’s g Formula

N#

65

2 log10 (10Gp Gs ) 3(Z s Zp ) / 2S

1

Digital FIR Filter Design y Example 7. Design a digital lowpass FIR filter with

the following specifications y δp = 0.1 , δs = 0.1, Fp = 1 KHz , Fs = 4KHz y FT = 10KHz 10KH Solution ωp = 2π * 1000 / 10000 = 0.2 π ωs = 2π * 4000 / 10000 = 0.8 π 20 log10 ( 0.01)  13 N# ⎡⎢1.59 ⎤⎥ 2 14 6(0 8S  00.2 14.6(0.8 2S) / 2S Thus, the minimum length of the filter is L = N+1 = 3 Let ωc = (ωs + ωp)/2 = 0.5 05π 66

Digital FIR Filter Design y Example 7. Continued

For L = 3, F 3 2M+1 = 3, 3 so M = 1. 1 The Th truncated t t d filter filt starts from n = -1. E l i hLP[n] Evaluating [ ] , for f n = -1, 1 0, 0 1

gives hLP[n] = {0.3183 0.5000 0.3183} Delaying by M = 1 hLP[n] = {0.3183 {0 3183 0.5000 0 5000 0.3183} 0 3183} 67

Digital FIR Filter Design 1.4

y Example 7. Continued

N=2 1.2

t

Z

|H (ej )|

1

0.6 0.4 0.2 0

0

0.2

0.4

Z

0.6

0.8

1

0.8

1

1 N=2

05 0.5 0 -0.5

jZ

I (e )

wc = 0.5 0 5*pi pi ; N=2; n = -floor((N+1)/2):floor((N+1)/2); y = sin(wc*n) i ( * ) ./pi./n; / i/ y(find(isnan(y))) = wc/pi; [h,w] = freqz(y); figure plot(w/pi,abs(h)); figure g plot(w/pi,unwrap(angle(h)))

0.8

-1 -1.5 -2 2

68

-2.5 0

0.2

0.4

Z

0.6

Digital FIR Filter Design y Example 7. Continued

Eff t off filter Effect filt order d 1.4 N=2 N=8 N=20 N=50

1.2

0.8

t

Z

|H ((ej )|

1

0.6 0.4 0.2 0

69

0

0.2

0.4

Z

0.6

0.8

1

Digital FIR Filter Design y Example 7. Continued

Eff t off truncation Effect t ti with ith Hamming H i window i d 1.4 N=2 N=8 N=20 N=50

1.2

0.8

t

Z

|H ((ej )|

1

0.6 0.4 0.2 0

70

0

0.2

0.4

Z

0.6

0.8

1

Digital FIR Filter Design y Matlab functions y kaiserord k i d y fir1 y fir2 y hamming y hann y balckman y kaiser y fdatool

71