Filter Analysis and Design

Butterworth Filters Butterworth filters have a transfer function whose squared 1 2 magnitude has the form H a ( jω ) = 2n . 1 + (ω / ω c )

*

M. J. Roberts - All Rights Reserved

2

Butterworth Filters The poles of a lowpass Butterworth filter lie on a semicircle of radius ω c in the open left half-plane. The number of poles is n and the angular spacing between the poles is always π / n. For n odd there is one pole on the negative real axis and all the others occur in complex conjugate pairs. For even n all the poles occur in complex conjugate pairs. The form of the transfer function is n n 1 1 pk Ha ( s) = =∏ =∏ (1 − s / p1 )(1 − s / p2 )(1 − s / pn ) k=1 1 − s / pk k=1 pk − s

*

M. J. Roberts - All Rights Reserved

3

Frequency Transformation A normalized lowpass filter can be transformed into an un-normalized lowpass filter or into a highpass, bandpass or bandstop filter by frequency transformation. Lowpass to lowpass ⇒ s → s / ω c Lowpass to highpass ⇒ s → ω c / s s 2 + ω Lω H Lowpass to bandpass ⇒ s → s (ω H − ω L ) Lowpass to bandstop ⇒ s →

s (ω H − ω L ) s 2 + ω Lω H

⎛ ω c is the cutoff frequency in rad/s of a lowpass or highpass filter and ⎜ ω and ω are the lower and upper cutoff frequencies in rad/s of a H ⎜ L ⎝ bandpass or bandstop filter. *

M. J. Roberts - All Rights Reserved

⎞ ⎟ ⎟ ⎠ 4

Frequency Transformation Example! Design a first-order bandstop Butterworth filter with a cutoff frequencies of 1 kHz and 2 kHz. 1 1 s 2 + ω Lω H H norm ( s ) = ⇒ H BS ( s ) = = 2 s (ω H − ω L ) s +1 s + s (ω H − ω L ) + ω Lω H +1 2 s + ω Lω H ω L = 2π × 1000 = 2000π , ω H = 4000π s 2 + 8π 2 × 10 6 s 2 + 7.896 × 10 7 H BS ( s ) = 2 ≅ 2 2 6 s + 2000π s + 8π × 10 s + 6283s + 7.896 × 10 7 Zeros at s = ± j8886 which correspond to ω = ±8886 rad/s and f = 1414 Hz Poles at s = −3142 ± j8886.

*

M. J. Roberts - All Rights Reserved

5

Frequency Transformation Example!

*

M. J. Roberts - All Rights Reserved

6

Chebyshev, Elliptic and Bessel Filters The Butterworth filter is "maximally flat" in its passband. This means that its frequency response in the passband is monotonic and the slope approaches zero at the maximum response. The Chebyshev filter has ripple in either its passband or stopband depending on which type of Chebyshev filter it is. Ripple is a variation of the frequency response between two limits. Allowing this ripple in the passband or stopband lets the transition from pass to stop band be faster than for a Butterworth filter of the same order. The Elliptic (or Cauer) filter has ripple in both the passband and stopband and also has an even faster transition between passband and stopband than the Chebyshev filters.

*

M. J. Roberts - All Rights Reserved

7

Chebyshev, Elliptic and Bessel Filters

*

M. J. Roberts - All Rights Reserved

8

Impulse and Step Invariant Design

*

M. J. Roberts - All Rights Reserved

9

Impulse Invariant Design Let hδ ( t ) = h a ( t )δ Ts ( t ) . Then ∞

Hδ ( jω ) = fs

∑ H ( j (ω − kω )) a

s

k=−∞

L where h a ( t ) ←⎯ → Ha ( s ).

Also let h d [ n ] = h a ( nTs ) . Then, using ω /fs = Ω ⇒ ω = fs Ω

( )

H d e jΩ = fs

*

∞

∑ H ( jf (Ω − 2π k )) a

s

k=−∞

M. J. Roberts - All Rights Reserved

10

Impulse Invariant Design In the expression

( )

H d e jΩ = fs

∞

∑ H ( jf (Ω − 2π k ))

k = −∞

a

s

it is clear that the digital filter's frequency response consists of multiple scaled aliases of the analog filter's frequency response and, to the extent that the aliases overlap, the two frequency responses must differ.

*

M. J. Roberts - All Rights Reserved

11

Impulse Invariant Design As an example of impulse invariant design let Aω c2 Ha ( s) = 2 . 2 s + 2ω c s + ω c

(

)

Then h a ( t ) = 2Aω c e−ω c t / 2 sin ω ct / 2 u ( t ) . Sample fs times

(

)

per second to form h d [ n ] = 2Aω c e−ω c nTs / 2 sin ω c nTs / 2 u [ n ]. Then H d ( z ) = 2Aω c

(

ze−ω c Ts / 2 sin ω cTs / 2 z 2 − 2e−ω c Ts /

2

(

cos ω cTs / 2 z + e−2ω c Ts /

and

( )

H d e jΩ = 2Aω c *

)

)

(

e jΩ e−ω c Ts / 2 sin ω cTs / 2 e j 2Ω − 2e−ω c Ts /

2

(

)

2

)

cos ω cTs / 2 e jΩ + e−2ω c Ts /

M. J. Roberts - All Rights Reserved

2

12

Impulse Invariant Design The digital filter's frequency response is

( )=

Hd e

jΩ

2Aω c

(

e jΩ e−ω c Ts / 2 sin ω cTs / 2 e j 2Ω − 2e−ω c Ts /

2

(

)

)

cos ω cTs / 2 e jΩ + e−2ω c Ts /

( )

and, from a previous slide, H d e jΩ = fs

2

∞

∑ H ( jf (Ω − 2π k )).

k = −∞

a

s

Therefore Aω c2

∞

( )= f ∑ ⎡ jf ( Ω − 2π k ) ⎤

Hd e

jΩ

s

k = −∞

⎣

= 2Aω c

*

s

2 + j 2 ω f Ω − 2 π k + ω ( ) c s c ⎦ 2

(

e jΩ e−ω c Ts / 2 sin ω cTs / 2 e j 2Ω − 2e−ω c Ts /

2

(

)

)

cos ω cTs / 2 e jΩ + e−2ω c Ts /

M. J. Roberts - All Rights Reserved

2

13

Impulse Invariant Design Let A = 10 and ω c =100 and fs =200. Then

( )

Hd e

jΩ

∞

1

= 2000 ∑

k = −∞ ⎡ ⎣ j2 ( Ω − 2π k ) ⎤⎦ + j2 2 ( Ω − 2π k ) + 1 2

and

( ) = 1000

Hd e

jΩ

2

(

e jΩ e−1/2 2 sin 1 / 2 2 e j 2Ω − 2e−1/2

2

(

)

)

cos 1 / 2 2 e jΩ + e−1/

2

At Ω = 0, ∞

1 = 1958.5 2 2 k = −∞ 1 − 16π k − j4 2π k

H d (1) = 2000 ∑ and H d (1) = *

343.825 = 1958.5 1 − 1.31751 + 0.49306 M. J. Roberts - All Rights Reserved

14

Impulse Invariant Design By design the digital filter's impulse response consists of samples of the analog filter's impulse response.

*

M. J. Roberts - All Rights Reserved

15

Impulse Invariant Design The frequency response of the digital filter differs from the frequency response of the analog filter due to aliasing.

*

M. J. Roberts - All Rights Reserved

16

Impulse Invariant Design Direct Form II Realization of the Digital Filter Designed by the Impulse Invariant Technique

*

M. J. Roberts - All Rights Reserved

17

Impulse Invariant Design Let H a ( s ) =

Aω c Aω c Aω c ⇒ H a ( j2π f ) = ⇒ H a ( jω ) = . s + ωc j2π f + ω c jω + ω c

Then h a ( t ) = Aω c e−ω c t u ( t ) . Sample at rate fs to form h d [ n ] = Aω c e−ω c nTs u [ n ] z e jΩ jΩ and H d ( z ) = Aω c ⇒ H d e = Aω c jΩ − ω c Ts z−e e − e−ω c Ts jΩ ∞ A ω e c H d e jΩ = fs ∑ = Aω c jΩ e − e−ω c Ts k = −∞ jf s ( Ω − 2π k ) + ω c

( )

( )

Let A = 10, ω c = 50 and fs = 100 and check equality at Ω = 0. ∞

fs

∑

k = −∞

∞ Aω c 50000 = ∑ = 1020.7 jfs ( Ω − 2π k ) + ω c k = −∞ − j200π k + 50

e jΩ 1 Aω c jΩ = 500 = 1270.7 − ω c Ts −1/2 e −e 1− e *

M. J. Roberts - All Rights Reserved

18

Impulse Invariant Design Why are the two results ∞

fs

∑

k = −∞

∞ Aω c 50000 = ∑ = 1020.7 jfs ( Ω − 2π k ) + ω c k = −∞ − j200π k + 50

and e jΩ 1 Aω c jΩ = 500 = 1270.7 − ω c Ts −1/2 e −e 1− e different by almost 25%? The difference comes from sampling h a ( t ) to form

h d [ n ] = Aω c e−ω c nTs u [ n ]. h a ( t ) has a discontinuity at t = 0, a sample time. If we

( )

let the first sample be Aω c / 2 instead of Aω c the two expressions for H d e jΩ

are now exactly equal. Aω c / 2 is the average of the two limits at the point of discontinuity. This is consistent with Fourier transform theory which says that the Fourier representation at a discontinuity is the average of the two limits. *

M. J. Roberts - All Rights Reserved

19

Impulse Invariant Design 9.87 × 10 4 s 2 Let H a ( s ) = 4 s + 444.3s 3 + 2.467 × 10 6 s 2 + 5.262 × 10 8 s + 1.403 × 1012 and let fs = 1000. h a ( t ) = ⎡⎣ 246.07e−122.41t cos(1199.4t − 1.48) + 200.5e−99.74t cos(977.27t + 1.683) ⎤⎦ u ( t )

h d [ n ] = ⎡⎣ 246.07e−0.12241n cos(1.1994n − 1.48) + 200.5e−0.09974 n cos(0.97727n + 1.683) ⎤⎦ u [ n ] 48.4z 3 − 107.7z 2 + 51.46z Hd ( z) = 4 z − 1.655z 3 + 2.252z 2 − 1.319z + 0.6413

*

M. J. Roberts - All Rights Reserved

20

Impulse Invariant Design

*

M. J. Roberts - All Rights Reserved

21

Impulse Invariant Design

*

M. J. Roberts - All Rights Reserved

22

Impulse Invariant Design

*

M. J. Roberts - All Rights Reserved

23

Step Invariant Design In step invariant design the analog unit step response is sampled to form the digital unit sequence response. ⎛ Ha ( s) ⎞ h −1a ( t ) = L ⎜ ⇒ h −1d [ n ] = h −1,a ( nTs ) ⎟ ⎝ s ⎠ −1

z z −1 Z ( h −1d [ n ]) = Hd ( z) ⇒ Hd ( z) = Z ( h −1d [ n ]) z −1 z 9.87 × 10 4 s 2 Let H a ( s ) = 4 s + 444.3s 3 + 2.467 × 10 6 s 2 + 5.262 × 10 8 s + 1.403 × 1012 and let fs = 1000.

*

M. J. Roberts - All Rights Reserved

24

Step Invariant Design ⎡ 0.2041e−122.408t cos(1199.4t + 3.1312) ⎤ h −1a ( t ) = ⎢ ⎥ u (t ) −99.74t cos(977.27t + 0.01042) ⎦ ⎣ +0.2041e ⎡ 0.2041( 0.8847 )n cos(1.1994n + 3.1312) ⎤ h −1d [ n ] = ⎢ ⎥ u[n] n ⎢⎣ +0.2041( 0.9051) cos(0.97727n + 0.0102) ⎥⎦ 0.03443z 3 − 0.03905z 2 − 0.02527z + 0.02988 Hd ( z) = z 4 − 1.655z 3 + 2.252z 2 − 1.319z + 0.6413 y [ n ] = 0.03443x [ n − 1] − 0.03905 x [ n − 2 ] − 0.02527 x [ n − 2 ] + 0.02988 x [ n − 4 ] +1.655y [ n − 1] − 2.252 y [ n − 2 ] + 1.319 y [ n − 3] − 0.6413y [ n − 4 ]

*

M. J. Roberts - All Rights Reserved

25

Step Invariant Design

*

M. J. Roberts - All Rights Reserved

26

Step Invariant Design

*

M. J. Roberts - All Rights Reserved

27

Step Invariant Design

*

M. J. Roberts - All Rights Reserved

28

Finite-Difference Design Finite-difference design is based on approximating continuous-time 1 derivatives by finite differences. The transfer function H a ( s ) = s+a d implies the differential equation ( y a ( t )) + a y a ( t ) = x a ( t ) . The dt y d [ n + 1] − y d [ n ] d derivative can be approximated by ( y a ( t )) ≅ . Then dt Ts

y d [ n + 1] − y d [ n ] + a y d [ n ] = x d [ n ] and the differential equation in Ts continuous time has become a difference equation in discrete time. Then Yd ( z ) Ts we can form the discrete-time transfer function H d ( z ) = = . X d ( z ) z − (1 − aTs ) *

M. J. Roberts - All Rights Reserved

29

Finite-Difference Design A derivative can be approximated by a forward difference y d [ n + 1] − y d [ n ] Z z − 1 d y a ( t )) ≅ ←⎯→ Yd ( z ) ( dt Ts Ts a backward difference y d [ n ] − y d [ n − 1] Z 1 − z −1 d z −1 y ( t )) ≅ ←⎯→ Yd ( z ) = Yd ( z ) ( dt Ts Ts zTs or a central difference y d [ n + 1] − y d [ n − 1] Z z − z −1 d z2 − 1 y a ( t )) ≅ ←⎯→ Yd ( z ) = Yd ( z ) . ( dt 2Ts 2Ts 2zTs Since every multiplication by s represents a differentiation in time, finite difference design can be done by simply replacing s with a z-transform expression for a difference that approximates a derivative. *

M. J. Roberts - All Rights Reserved

30

Finite-Difference Design Using a forward difference, 1 1 ⎡ 1 ⎤ Hd ( z) = Ha ( s) = ⎢ = = Ts ⎥ z − 1 z −1 z − (1 − aTs ) ⎣ s + a ⎦ s→ + a Ts Ts Using a backward difference, 1 Ts z ⎡ 1 ⎤ Hd ( z) = Ha ( s) = ⎢ = = ⎥ z − 1 1 z −1 1 + aT ⎣ s + a ⎦ s→ s +a z− zTs zTs 1 + aTs Using a central difference, 1 z ⎡ 1 ⎤ Hd ( z) = Ha ( s) = ⎢ = 2 = 2Ts 2 2 ⎥ z −1 z −1 z + 2aTs z − 1 ⎣ s + a ⎦ s→ + a 2 zTs 2zTs *

M. J. Roberts - All Rights Reserved

31

Finite-Difference Design Using the finite-difference technique it is possible to approximate a stable continuous-time filter by an unstable discrete-time filter. 1 In H a ( s ) = the pole is at s = −a. In the digital filter s+a 1 H d ( z ) = Ts the pole is at z = (1 − aTs ) . If H a ( s ) is z − (1 − aTs ) stable, a > 0 and 1 − aTs is on the real axis of the z plane. If aTs ≥ 2 the digital filter is unstable. One way to stabilize the design is to use a smaller value for Ts , which implies a higher sampling rate. This design was based on using a forward difference.

*

M. J. Roberts - All Rights Reserved

32

Finite-Difference Design If we use a backward difference instead of a forward difference we get the digital filter ⎡ 1 ⎤ Hd ( z) = Ha ( s) = ⎢ ⎣ s + a ⎥⎦ s→ z −1 zTs

=

Ts z . 1 + aTs z − 1 / (1 + aTs )

The pole is now at z = 1 / (1 + aTs ) and is guaranteed to be inside the unit circle in the z plane, making this filter stable for any a > 0 and any Ts . Using a backward difference in any digital filter design guarantees stability, but that guarantee comes at a price. The pole locations are restricted to a small region inside the unit circle. *

M. J. Roberts - All Rights Reserved

33

Finite-Difference Design 9.87 × 10 4 s 2 Let H a ( s ) = 4 s + 444.3s 3 + 2.467 × 10 6 s 2 + 5.262 × 10 8 s + 1.403 × 1012 and design a digital filter using the finite-difference technique with fs = 1000. 0.169z 2 ( z − 1) Hd ( z) = 4 z − 1.848z 3 + 1.678z 2 − 0.7609z + 0.1712 2

*

M. J. Roberts - All Rights Reserved

34

Finite-Difference Design

*

M. J. Roberts - All Rights Reserved

35

Finite-Difference Design

*

M. J. Roberts - All Rights Reserved

36

Matched z and Direct Substitution The matched-z and direct substitution digital filter design techniques use the correspondence z = e sTs to map the poles and zeros of the analog filter to corresponding locations in the z plane. If the analog filter has a pole at s = a, the corresonding digital filter pole would be at z = eaTs . The matched-z method aTs z − e replaces every factor of the form s − a with the factor 1 − eaTs z −1 or and z the direct substitution method replaces it with the factor z − eaTs . Zeros are mapped in the same way.

s=a *

M. J. Roberts - All Rights Reserved

z = eaTs 37

Matched z and Direct Substitution 9.87 × 10 4 s 2 Let H a ( s ) = 4 . 3 6 2 8 12 s + 444.3s + 2.467 × 10 s + 5.262 × 10 s + 1.403 × 10 Then, using the matched-z technique z ( z − 1) z 4 − 1.655z 3 + 2.252z 2 − 1.319z + 0.6413 y [ n ] = 98700 x [ n ] − 197400 x [ n − 1] + 98700 x [ n − 2 ] H d ( z ) = 98700

2

2

+1.655 y [ n − 1] − 2.252 y [ n − 2 ] + 1.319 y [ n − 3] − 0.6413y [ n − 4 ]

*

M. J. Roberts - All Rights Reserved

38

Matched z and Direct Substitution

*

M. J. Roberts - All Rights Reserved

39

Matched z and Direct Substitution

*

M. J. Roberts - All Rights Reserved

40

Matched z and Direct Substitution

*

M. J. Roberts - All Rights Reserved

41

The Bilinear z Method One approach to digital filter design would be to use the mapping z = e sTs ⇒ s → (1 / Ts ) ln ( z ) to convert the s-domain expression to the z-domain expression H d ( z ) = H a ( s ) s→ (1/T ) ln ( z ) . s

Unfortunately this creates a transcendental function of z with infinitely many poles. We can simplify the result by making an approximation to the exponential function. The exponential function is defined by 2 3 k ∞ x x x ex = 1 + x + + + = ∑ 2! 3! k = 0 k!

*

M. J. Roberts - All Rights Reserved

42

The Bilinear z Transform We can approximate the exponential function by truncating its series definition to two terms e x ≅ 1 + x. Then z = 1 + sTs or z −1 s→ . This is exactly the same as the finite difference technique Ts using forward differences and has the same drawback, the possibility of converting a stable analog filter into an unstable digital filter. A simple but very significant variation on this idea is to express e sTs in e sTs /2 the form − sTs /2 and then approximate both exponentials with the first e 1 + sTs / 2 2 z −1 two terms of the series. Then z = or s → . This is 1 − sTs / 2 Ts z + 1 the basis of the bilinear digital filter design technique. *

M. J. Roberts - All Rights Reserved

43

The Bilinear z Transform 2 z −1 let s = jω . Then Ts z + 1 2 + sTs 2 + jω Ts j 2 tan −1 (ω Ts /2 ) −1 z= = = 1 tan (ω Ts / 2 ) = e 2 − sTs 2 − jω Ts

In the mapping s →

The ω axis of the s plane maps into the unit circle of the z plane. If s = σ + jω , σ < 0, the mapping is from the left half of the s plane to the interior of the unit circle in the z plane. If s = σ + jω , σ > 0, the mapping is from the right half of the s plane to the exterior of the unit circle in the z plane.

*

M. J. Roberts - All Rights Reserved

44

The Bilinear z Transform For real discrete-time frequencies Ω the correspondence between 2 e jΩ − 1 2 ⎛ Ω⎞ the s and z planes is s = = j tan ⎜ ⎟ . Then, since jΩ ⎝ 2⎠ Ts e + 1 Ts

s = σ + jω , σ = 0 and ω = ( 2 / Ts ) tan ( Ω / 2 ) or Ω = 2 tan −1 (ω Ts / 2 ) .

*

M. J. Roberts - All Rights Reserved

45

The Bilinear z Transform

9.87 × 10 4 s 2 Let H a ( s ) = 4 . 3 6 2 8 12 s + 444.3s + 2.467 × 10 s + 5.262 × 10 s + 1.403 × 10 Then H d ( z ) = 12.38

*

( z + 1)2 ( z − 1)2

z 4 − 1.989z 3 + 2.656z 2 − 1.675z + 0.711

M. J. Roberts - All Rights Reserved

46

The Bilinear z Transform

*

M. J. Roberts - All Rights Reserved

47

The Bilinear z Transform

*

M. J. Roberts - All Rights Reserved

48

The Bilinear z Transform

*

M. J. Roberts - All Rights Reserved

49

FIR Filter Design An analog filter can be approximated by an FIR digital filter. The impulse response of the analog filter has infinite duration in time but beyond some time it has dropped to a very low value and the remainder beyond that time can be truncated with little change in the filter characteristics. We can create a digital FIR filter by sampling over that time. *

M. J. Roberts - All Rights Reserved

50

FIR Filter Design

We can also approximate ideal filters, which are non-causal, by truncating that part of the impulse response occurring before time n = 0 and later when the impulse response has fallen to a very low value.

*

M. J. Roberts - All Rights Reserved

51

FIR Filter Design The truncated impulse response has the form hN [n] =

N −1

∑ a δ [n − m]

m=0

m

where N is the number of samples retained. This type of FIR digital filter can be realized in Direct Form II. Its transfer function is Hd ( z) =

*

N −1

−m a z ∑ m .

m=0

M. J. Roberts - All Rights Reserved

52

FIR Filter Design ⎧h d [ n ] , 0 ≤ n < N ⎫ Let h N [ n ] = ⎨ ⎬ = h d [ n ] w [ n ] be the , otherwise ⎭ ⎩0 impulse response of the digital filter where h d [ n ] is the

sampled ideal analog filter impulse response and w [ n ] is a window function. Then the frequency response is

( )

( ) ( ) The periodic convolution with W ( e ) changes the H N e jΩ = H d e jΩ  W e jΩ . jΩ

frequency response from ideal and introduces ripple in the frequency response. *

M. J. Roberts - All Rights Reserved

53

FIR Filter Design Impulse response and frequency response of an FIR lowpass digital filter designed by truncating the ideal impulse response.

*

M. J. Roberts - All Rights Reserved

54

FIR Filter Design Impulse response and frequency response of an FIR lowpass digital filter designed by truncating the ideal impulse response.

*

M. J. Roberts - All Rights Reserved

55

FIR Filter Design Impulse response and frequency response of an FIR lowpass digital filter designed by truncating the ideal impulse response.

*

M. J. Roberts - All Rights Reserved

56

FIR Filter Design The digital filter impulse and frequency responses on the previous three slides were all formed by using a rectangular ⎧1 , 0 ≤ n < N ⎫ window of the form, w [ n ] = ⎨ ⎬ . Other window ⎩0 , otherwise ⎭ functions could be used to reduce the ripple in the frequency response.

*

M. J. Roberts - All Rights Reserved

57

FIR Filter Design Frequently-Used Window Functions von Hann

Bartlett

Hamming

w[n] =

1⎡ ⎛ 2π n ⎞ ⎤ 1 − cos ⎜⎝ ⎟⎠ ⎥ , 0 ≤ n < N ⎢ 2⎣ N −1 ⎦

N −1 ⎧ 2n , 0 ≤ n ≤ ⎪⎪ N − 1 2 w[n] = ⎨ ⎪2 − 2n , N − 1 ≤ n < N ⎪⎩ N −1 2 ⎛ 2π n ⎞ w [ n ] = 0.54 − 0.46 cos ⎜ , 0≤n