EKT 4304 DIGITAL SIGNAL PROCESSING 20072008 CHAPTER 8

EKT 430/4 DIGITAL SIGNAL PROCESSING 2007/2008 CHAPTER 8 IIR FILTER DESIGN

IIR based Digital Filters: The design uses the analog prototype of given specifications. For the given specifications, IIR filters yield much smaller order of the filter. Hence requires much less number of elements. Appropriate mapping and spectral transformation is carried out to result in a desired digital filter.

The object is to: Investigate the relationship between s- and z- domains. Workout the equivalence of analog and digital systems. To study the different mapping rules and response matching.

Equivalence of analog and digital systems: The approximate continuous time domain impulse response h(t) at the sampling interval ts (=1/fs) is: h(t) ha(t) = ts h(t) (t-kts) = ts h(kts) (t-kts); k =[- , ]. For the normalized sampled time impulse response, the sampling period ts is 1 s and replacing t by n : hs[n] = hs[k] [n-k] for k =[- , ].

Equivalence of analog and digital systems: For k [- , ], the Laplace transform of ha(t)= ts h(kts) (t-kts) is H(s) Ha(s) = ts h(kts) exp(-sts)k and z-transform H(z) of normalized hs[n] is: hs[n] = hs[k] [n-k] is Hd(z) = hs[k]z-k

Equivalence …. . Thus for zk = exp(skts), Ha(s) = ts Hd(z) the conditions can be rewritten as z exp(sts) and s ln(z)/ts. These relationships describe a mapping between variables s and z: z = exp( j )ts = exp( ts) exp(j ts) = exp( ts) exp(j ) Here = ts = 2 f/fs = 2 F, where F is the digital frequency in cycles/sample ts is the sampling rate per second is the analog frequency in radians per second. .

…. . and mapping The DTFT of hs[n] has a periodic spectrum. Therefore ensure that the analog signal h(t) is strictly band-limited and moreover, to obliterate aliasing, is sampled much above the Nyquist rate The normalized band of frequency should lie well within the principal period of Hp(f) – 0. 5 F 0. 5.

aliasing revisited: The relationships z exp(sts) s ln(z)/ts does not describe a unique mapping between the s-plane and the z-plane. Since exp(j ) has a period of 2 ; frequencies beyond , such as k fs are also mapped on the same circle after k . Thus a one to one mapping is possible only if the signal is restricted in the range f/fs= F = 0. 5. j 3 /T Im z /T Re z -1 1 - /T -3 /T s-plane z-plane

Notes … The logarithmic nature of the transformation s ln(z)/ts, being a infinite series, , does not convert a rational TF H(s) into a rational TF H(z). The unique representation is true only for band limited signals, those lie well within the principal range. The spectra of such band limited signal repeats at every sampling frequency. If the signal is not band-limited, it will cause inherent aliasing (introduction of corrupted signal in the principal region).

Notes…. We need to study the practical mappings. The methods are: Use of standard z-transform: Equivalent finite element integrators. Other methods (not discussed) w-v transform z-forms etc.

Use of standard z-transform Impulse invariant filters and Matched z-transform filters § n n These filters use standard z-transform. Therefore the intelligence has to be band limited. The Impulse invariant filter preserves time response. The Matched z-transform filters preserves frequency response. Aliasing is severe. Use high sampling rate.

Impulse Invariance filter mapping n It is designed to preserves time response of analog filters using standard z-transform. The mapping rule is: s=ln(z)/ T n For compatible frequency response, sampling frequency should be considerably high. The signal should be strictly band-limited. n

Matched z-transform n Ø Ø Ø Factored form of s-domain transfer function is used. Individual singularity is transformed using the standard z-transform. (s-a) = (1 -z-1 ea. T)= z/(z-ea. T) Characteristic polynomial (denominator) of this and Impulse invariant method are same. It is due to the use of standard z-transform. The zeros are different. The filter so obtained preserves the frequency response of the analog counterpart, though exhibits considerable distortion. Standard z-transform related filters should be band limited.

Equivalent Integrator: The basic philosophy: The simulation element in signal processing is Integrator. It is (1/s) in Laplace Transform. Different numerical methods of integration are used to create equivalence. They are grouped as : Backward difference, Forward difference, Trapezoidal integration, Simpson’s integration. [not used, increases the order of filter] “ Philips CL, + Nagle HT, “ Digital Control System Analysis and Design” 3/e, PH Ch. 11

Backward Difference: (right side rectangular rule) Difference of present and the past sample. y[t] = y(t) – y(t-T) n It’s Laplace transform: s. Y(s) = [Y(s) – e-s. TY(s)]/T + Y(0+) letting Y(0+)=0, we get: s = (1 - e-s. T)/T It is not the z transform of 1/s = (1 -z-1)/T Or, 1/s = T/(1 -z-1) n z = 1/(1 -Ts) = ½ [(1 -s. T)+(1+s. T)]/(1 -s. T) = ½ [(1 - (1+s. T)/(1 -s. T)];

Backward Difference: (right side rectangular rule) z = ½ [(1 - (1+s. T)/(1 -s. T)]; replacing s = j returns: = ½ [1+ e j 2 arctan( T)] = ½ [1+ej ]. Mapping left half of j axis of s-plane into a circle of radius ½ centered at (½ , 0) in z plane. The region enclosed is very small. If our frequency extends little beyond s=0 or, z=1, there will be severe distortion. Sampling Frequency should be increased to improve the performance.

Forward Difference: left side rectangular rule It is the difference between present pulse and forth coming pulse. y[t] =[y(t+T)- y(t)]/T n It’s Laplace transform: s = (z-1)/T z = 1+Ts = 1+ j T This mapping only shifts j axis to right by 1. It can make a stable filter into unstable one. n

In literature: Same mapping rule for Impulse Invariance Integrator, backward difference Integration and right side rectangular integration. Step invariance Integrator, forward difference integration and Left side rectangular integration.

Bilinear z-transform Mapping n n It does not require band limitation. Pre warping of analog frequency is essential before bilinear z-transformation for suitable matching of frequency response. Most suitable for polynomial based flat response filters. Uses Trapezoidal Integration rule.

Algorithms for Bilinear Transformation This transformation follows the trapezoidal interpolation rule, or, Tustin’s rule, and is represented by the numerical integration algorithm: y[n] = y[n-1] + (ts /2){x[n] + x[n-1]}. So that: H(z) = 0. 5 ts [z+1]/[z-1] And equating with integration rule in s-domain. H(s) 1/s. Thus s = (2/ts){(z– 1) / (z+1)}

The time-response matching The idea behind response matching is to match the time-domain analog and digital response at the sampling instants t = ts for a given input. One way to achieve the response is as under. t-domain input sampled @ ts Input X[z] s-domain input analog filter in s-domain output in t-domain output sampled @ ts output Y[z] hence: H[z] = Y[z]/X[z]. The example below will illustrate that the H[z] is excitation specific.

Time response matching…

Example 10. 01 Convert analog filter of transfer function H(s) = 1/(s+1) into a digital filter H[z] for (a) impulse input (b) unit step input. Soln: (a) Impulse input : x(t) = (t) X(s) = 1; H(s) = 1/(s+1), output Y(s) = 1/(s+1) or, y(t) = e-t U(t). the sampled version is y(n. T) = e-n. TU[n]. And at n=0, y(n. T) = 1.

Example 7. 01 n Letting sampling time ts = 1, for impulse input, we get H[z] =Z e-n. TU[n] =1/[1 -e-1 z-1] =1/[1 - 0. 3679 z-1]. We thus note: § The DC gain H[z=1] = 1. 582 while it at H(s=0) =1. § This increase is solely due to aliasing. § Even if the DC gain is normalized, the frequency characteristics does not match. It is shown in Magnitude plot of Fig. 7. 02 (a) and (b). § This can be remedied by higher rate of sampling. But it makes it unsuitable for HP and BP filters.

Normalized amplitude response:

Phase response

Time response of Impulse Invariance digital filters and equivalent analog filter: note similarity in time response.

Response with step input: Step response of the filter H[z] = z / [z-0. 3679] is Y[z] = z 2 / [z-1][z – 0. 3679] Inverse Z-transform is y[n] = 1. 5820 U[n]– 0. 214 (0. 3679)-n U[n] =1. 5820[1 -0. 1353(0. 3679)-n]U[n].

Magnitude-frequency response for step input

Remarks on impulse invariant method: Comparison of response in analog domain filters and impulse invariant method in digital domain while the frequency and phase plots are not alike, the time response at the sampling intervals is same in analog and impulse invariant digital filters the time response is remarkably different from its analog counterpart for step input. This clearly illustrates that this scheme is excitation specific.

Why Bilinear Transformation Due to high rate of sampling, Impulse invariant filters (IIF) are not suitable for high pass/band reject/all-pass modes while matched-z filters can not realize LPF. Frequency response does not match with those of analog filters in IIF while there is a distortion in matched-z. Aliasing problem in both cases.

Bilinear Transformation Alternatively a transformation is attempted which converts 0 to 2 change in [0 to ]. Thus the circle of z-domain is converted onto entire j -axis with one to one mapping. This eliminates the aliasing problem. It is BILINEAR TRANSFORMATION

s- and z- planes mapping- we want

The relationship between analog and digital frequency : s = (2/ts)[(z-1)/(z+1)] Putting s = j , and z = ej j =(2/ts)[(ej -1) / (ej +1)] =(2/ts)[{ej /2 (ej /2 – e-j /2)} / {ej /2(ej /2+ e-j /2)}] =(2/ts)[(ej /2 – e-j /2)/ (ej /2+ e-j /2)] =j(2/ts)sin( /2)/cos( /2) Or, = (2/ts)tan( /2) = 2 tan-1( ts/2)

Frequency Warpage in BLT: The graph between and for ts=1, red color, is non-linear implying that the warpage has been introduced in the frequency translation. Blue color graph for is drawn for understanding.

Working out frequency in BLT: = (2/ts)tan( /2) and = 2 tan-1( ts/2) The above is a non linear or warped characteristic. It has one to one mapping between s and z domains. For ; /2. To undo this frequency warping, we pre-warp the analog frequency p in s-plan so that: pts/2 = tan( ts/2).

Linearization of relationship After frequency prewarping pts/2 = tan( ts/2) The Digital frequency becomes: = 2 tan-1( p ts/2) = 2 tan-1 tan( ts/2) = (2/ts)( ts/2); we get a linear relation: = .

More over due to carrying out the prewarping, the term 2/ts turns out to be redundant as it cancels out in calculations. Therefore we can individually use: n s = (z-1)/(z+1) n p = (2/ts)tan( ts/2). n Or, combining the two: s = cot ( ts/2) (z-1)/(z+1). [Ifeacher et al, ”Digital signal Processing’ 2/e, Prentice Hall, pp 473 -500.

Example: 10. 02 Using un-prewarped bilinear transformation and impulse invariance technique synthesize: H(s) = c /(s+ c ), c = 105. 1. Bilinear Transformation Method: Putting the transformation s =(2/ts)[(z-1) / (z+1)]; we get; H[z} = c/[(2/ts)[(z-1) / (z+1)] )+ c] = c. T(z+1)/(2+ c. T)z – (2 - c. T). Selecting the value of T such that T << / c say = 0. 1 / c, H[z] = 0. 1357(z +1)/[z- 0. 7284] Lathi B. P. , ”Signal Processing and Linear Systems” 1998, Oxford

Ex. 9. 02 contd. Thus |H(ej T)| = 0. 024(1+cos ( 10 -6)/ (1 -0. 9518 cos( 10 -6) and H(e j T) = tan-1[sin( 10 -6)/ [1+ cos( 10 -6)] - tan-1[sin( 10 -6)/[cos( 10 -6)-0. 7285].

Ex. 9. 02: Impulse Invariance Technique x(t) = (t), X(s) = 1; H(s) = c /(s+ c), output Y(s) = c /(s+ c) y(t) = c exp(- ct) U(t). Sampled version is y (n. T) = exp(- cn. T)U[n]. Therefore H[z] = Y[z]/X[z] Y[z] = c. Tz /[z-exp(- c. T)] For T = 0. 1 / c, as above, H[z] = 0. 3142 z/[z- 0. 7304] Thus |H(ej T)| = 0. 1 [1. 533 – 1. 4608 cos(( 10 -6)] and H(ej T) = -tan-1[0. 7304 sin( 10 -6)/ [10. 7304 cos( 10 -6)]

Ex. 09. 02: Magnitude function:

Ex. 9. 02: : Phase function

Ex. 10. 03: Design an IIR digital filter for T = 2 sec and the transfer function of analog filter given as: H(s) = (s+a)/{(s+a)2 + b 2}. where a = 0. 1, b = 3. Impulse invariance method: Inverse Laplace transform is: h(t) = e-at cos(bt) U(t). Letting T = ts, the discrete time output is: h[n. T] = e-an. Tcos[bn. T] U[n. T]; Z-transform of above is H[z] = (z 2 -e-0. 2 cos[6] z) / {z 2 - 2 e-0. 2 cos[6] z +e-0. 4}.

Ex. 9. 03: Bilinear transformation method without (pre-warping H(s) = (s+0. 1)/{(s+0. 1)2 +9} replace s by 2/T) {(z-1)/(z+1)}, H(z) = 0. 1077{(z 2 + 0. 18 z – 0. 82) / (z 2 + 1. 569 z + 0. 96)}

Ex. 9. 03: Magnitude plots

Ex. 10. 03: Design an IIR digital filter for the TF of analog filter given as: H(s) = (s+a)/{(s+a)2 + b 2}. where [a, b] =[0. 1, 3] and the resonant frequency r= /4. with pre-warping. n n n The sampling time, given that b = c= 3; Using the relation: =2 tan-1( ts/2) and replacing c and c; ts = (2/ c) tan( c/2) = (2/3)tan( /8)=0. 276. H(z) = N(z)/D(z) = H(s) for s = (2/T)(z-1)/(z+1) N(z) = (2/0. 276){(z-1)/(z+1)} + 0. 1 D(z) = N 2(z) + 9. H(z) = [1+0. 027 z-1 -0. 973 z-2]/[8. 572 – 11. 84 z-1+8. 177 z-2] [Salivahanan, “DSP” Mcgraw Hills, p. 429 -430]

BZT Bilinear transformation method are useful when the gains are constant over certain frequency bands. This condition exists in lowpass, band-pass, high-pass and band-stop filters. This method requires lower sampling rate compared to Impulse Invariance Method because of absence of aliasing. This absence of aliasing is the result of one-one mapping from a s-plane to z-plane.

Warping The digital frequency is limited in the range of -0. 5 F 0. 5, and sampling time period is decided by the highest frequency of the signal Fh is to be preserved. Warping of Digital Frequency scale is inherent in Bilinear Transformation. This can be reduced by pre warping the frequencies in analog domain.

Steps to Design IIR BZT filters n n n From specifications, workout the desired discretetime filter requirements and attributes. Pre-warp the critical discrete time frequencies into corresponding analog frequencies and estimate the order of filter and -3 d. B critical frequency. Design the analog prototype filter from the continuous-time parameters obtained from above. Convert the analog prototype into an analog filter using frequency-frequency transform. Design a digital filter using a BZTof H(s) that warps the frequency axis which was previously prewarped. [1. Taylor, ”Hands-on digital signal processing”, Mc. Graw Hills, pp 414 -420 ISBN: 0 -07 -912965 x 2. Chen, ”Digital Signal Processing” Oxford, PP. 386 -394: ISBN 0 -19 -513638 -1]

Example 10. 04: Using the bilinear transformation and a low pass analog Butterworth prototype, design a low pass filter operating at a rate of 20 k. Hz and having a pass band extending to 4 k. Hz with maximum attenuation of 0. 5 d. B and stop band starting at 5 k. Hz with a minimum stop band attenuation of 10 d. B. Soln : Given that sampling rate is fs=20 k. S/s, Pass band =fpass= 4 k. Hz. Fpass= 4/20 Hz/s at -0. 5 d. B. Stop band =fstop = 5 k. Hz. Fstop= 5/20 Hz/s at +10 d. B. [Orfanidis: Introduction to signal Processing pp 612 -617]

Soln. . : n n Prewarped Digital frequencies are: pass= tan (2 Fpass/2) = 0. 7265 rad/sample stop= tan (2 Fstop/2) = 1. 00 rad/sample. ln( )= ln ( stop/ pass) = ln(1. 0/0. 7265)

Determination of order n of the filter: Given that H(s)H(-s) = [1+(s/j)2 n]-1, Gx = 20 log |H(j x)| = -10 log[1 + ( x/ c)2 n] Where s = x, c = cut-off frequency, n is order of filter. Letting Pass band gain Gp d. B occur at p rad/sec and stop band gain Gs d. B occur at s rad/sec. We get expressions as: Gp = -10 log[1 + ( p/ c)2 n] ( p/ c)2 n = 10 -Gp/10 -1 Gs = -10 log[1 + ( s/ c)2 n] ( s/ c)2 n = 10 -Gs/10 – 1 Hence: ( s/ p)2 n = (10 -Gs/10 – 1)/ (10 -Gp/10 – 1)

Conclusion From here we conclude: n, the order of the filter = log {(10 -Gs/10 – 1)/ (10 -Gp/10 – 1)} / 2 log( s/ p). c, the cut off frequency = p/{10 -Gp/10 – 1}1/2 n ALTERNATIVELLY AS = s/{10 -Gs/10 – 1}1/2 n

Hence: ( s/ p)2 n = (10 -Gs/10 – 1)/ (10 -Gp/10 – 1) 2 nlog(( s/ p)=log[(10 -Gs/10 – 1)/ (10 -Gp/10 – 1)] Or, n = 0. 5 log[(10 -Gs/10 – 1)/ (10 -Gp/10 -1)]/log( s/ p) In this case, n = 0. 5 log[(10 -Gs/10 – 1)/ (10 -Gp/10 -1)]/log( s/ p) = 6. 73 say 7. The -3 d. B cutoff frequency calculated as c, = p/{10 -Gp/10 – 1}1/2 n = 0. 8443 f = (f / ) tan-1( )= 4. 464 k. Hz. o s c

The filter is a 7 th order filter One real pole and three complex poles. The composite TF is: = 0. 4578 (1+z-1)/1 -0. 0844 z-1 + 0. 3413(1+z-1)2/[1 -0. 2749 z-1 + 0. 6402 z-2] + 0. 2578(1+z-1)2/[1 -0. 2076 z-1+0. 2386 z-2] + 0. 2204(1+z-1)2/[1 -0. 1775 z-1+0. 0592 z-2]

Example_10. 05 [ extention of example 9. 4] A band stop filter to satisfy Butterworth and Tchebyshev Criteria is required to meet following specifications. to meet the following specifications. 1. Stop band 100 to 600 Hz (assumed at 3 d. B) 2. At and between 200 and 400 Hz, the magnitude should be at least -20 d. B. 3. Maximum Gain (at the zero) =1. 4. Pass band ripple 1. 1 d. B. 5. Sampling frequency 2000 Hz.

Conversion to normalized LPF specs. n Calculation of prewarp frequencies: pr_x= tan ( x. T/2) = tan ( f/fs) n fx Hz fx/fs rad/s since fs. T=1 pr_x rad/s fx 1 100 0. 1571 0. 1584 fx 2 200 0. 3141 0. 3249 fx 3 300 0. 6283 0. 7265 fx 4 400 0. 9425 1. 3764 Normalized Bandwidth = B = fx 4 -fx 1 = 1. 3764 -0. 1584 = 1. 218 Normalized center frequency = 02= fx 4 fx 1 = 0. 2180

Calculation of equivalent LPF n B x /( 02 - x 2) n Basic equation is: n Calculating n for various values of x in the equation =1. 218 x / (0. 2180 - x 2) n n n x 0. 1584 0. 3249 0. 7265 1. 3764 soln 1. 218 x x 0. 1584/[0. 218 – (0. 1584) 2] 0. 3249 /[0. 218 – (0. 3249) 2] 0. 7265 /[0. 218 – (0. 7265) 2] 1. 3764 /[0. 218 – (1. 3764 ) 2] n =1 = 3. 519 =-2. 856 = -1 s/ p=min( 2/ 1 3/ 4) = 2. 856 3. and 1/cosh-1(2. 856) = 0. 5846

Calculation of number of Tchebyshev poles n The control equations is where Given: Gs =-20 d. B, R=-1. 1 d. B and Calculated s/ p= 2. 856 1/cosh-1(2. 856) =0. 5846 cosh-1( )=cosh-1(18. 53) = 3. 612 n= 2. 11 3 for better transit characteristic

Chebyshev polynomial for n=3 and R = 1. 1 d. B In the polynomial calculated in 9. 4 Hlp(s) =0. 4656/(s+0. 4766)(s 2+0. 4766 s+0. 977) n To convert this normalized LP to desired BS, substitute s 1. 218 s / (s 2+ 0. 2180). n n Use BZT in the so obtained polynomial in s.

1. Butterworth filter. Number of poles of LP filter = log 10 {(10 -Gs/10 – 1)/ (10 -Gp/10 – 1)} / 2 log 10 ( s/ p) = log 10{(102 -1)/(100. 3 -1)/2 log 10(2. 856) = 2. 199 3 (for better transition characteristic). The normalized TF = 1/ (s+1)(s 2+s+1) = 1/ 1+2 s+2 s 2 +s 3 Transform normalized LP to required BS by s 1. 218 s / (s 2+ 0. 2180) USE BZT on the HBS(s).

Butterworth filter