1 / 72

University of Manchester Department of Computer Science CS 3291 : Digital Signal Processing

University of Manchester Department of Computer Science CS 3291 : Digital Signal Processing Section 5: z-transforms & IIR-type discrete time filters Remember: A discrete time filter is described by a linear difference equation. For example:

raglandj
Download Presentation

University of Manchester Department of Computer Science CS 3291 : Digital Signal Processing

An Image/Link below is provided (as is) to download presentation 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. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. University of Manchester Department of Computer Science CS 3291 : Digital Signal Processing Section 5: z-transforms & IIR-type discrete time filters Remember: A discrete time filter is described by a linear difference equation. For example: 1) y[n] = x[n] + x[n-1] 2) y[n] = x[n] + x[n-1] - b y[n-2] Example 2 is recursive & can have infinite impulse-response. CS3291 DSP Sectn 5

  2. Introduction: • General causal discrete time filter has difference equation: • N M • y[n] =  a i x[n-i] -  b j y[n-j] • i=0 j=1 • which is of order max{N,M}. • Recursive if any b j is non-zero. • A 2nd order recursive filter has the diff equn: • y[n] = a 0 x[n] + a 1 x[n-1] + a 2 x[n-2] - b 1 y[n-1] - b 2 y[n-2] • Frequency response of a discrete time filter:- •  • H(e j  ) =  h[n]e - j  n • n=- CS3291 DSP Sectn 5

  3. Consider response of causal stable LTI system to {z n } • i.e. { …, z-2, z-1, 1, z, z2, z3, ……… } • By discrete time convolution, output is {y[n]} where •  • y[n] =  h[k] z n - k = z n h[k] z - k • k=- k=- •  • = z n H(z) with H(z) =  h[k] z - k • k=- • H(z) is "z-transform" of impulse-response. • Equivalent of Laplace transform. • H(z) is a complex number for any z. CS3291 DSP Sectn 5

  4. For causal stable system, H(z) must be finite for z 1 •  • H(z) =  h[k] z - kh[k] z - k by causality • k=0 k=0 •  • =  h[k] z - k • k=0 •  • h[k] when  z  1 • k=0 • finite by defn of stability CS3291 DSP Sectn 5

  5. Some properties of H(z): (1) H(z) is complex no. for any z with z1. (2) Replacing z by e j  , gives H(e j  ) rel frequency-resp. On Argand diagram ("z-plane"), z = e j  lies on unit circle. (3) If input is {z n }, the output is {H(z) z n }. CS3291 DSP Sectn 5

  6. System function is z-transform of {h[n]} : •  • H(z) =  h[k] z - n • n=- • Frequency-response is DTFT of {h[n]} : •  • H(e j  ) =  h[n]e - j  n • n=- • Replace z by e j  CS3291 DSP Sectn 5

  7. Example 5.1: Find H(z) for the non-recursive difference equation: y[n] = x[n] + x[n-1] Solution: {h[n]} = { ... , 0, 1, 1, 0, ... }, therefore H(z) = 1 + z - 1 CS3291 DSP Sectn 5

  8. Example 5.2: Find H(z) for the recursive difference equation: • y[n] = a 0 x[n] + a 1 x[n-1] - b 1 y[n-1] • Solution: Method used above becomes hard. • Remember that if x[n] = z n then y[n] = H(z) z n , • y[n-1] = H(z) z n - 1 etc. • Substitute to obtain: • H(z) z n = a 0 z n + a 1 z n - 1 - b 1 H(z) z n - 1 • a 0 + a 1 z - 1 • H(z) =  • 1 + b 1 z - 1 except when z = - b 1 . • When z = -b 1 , H(z) = . CS3291 DSP Sectn 5

  9. For the general discrete time filter: • N M • y[n] =  a i x[n-i] -  b j y[n-j] • i=0 j=1 • the same method gives: • a 0 + a 1 z - 1 + a 2 z - 2 + ... + a N z - N • H(z) =  • b 0 + b 1 z - 1 + b 2 z - 2 + ... + b M z - M • with b0 = 1. • Given H(z) in this form, we can easily get its difference equn • & hence its signal flow graph. CS3291 DSP Sectn 5

  10. Example 5.3: Give a signal flow graph for: • a 0 + a 1 z -1 + a 2 z - 2 • H(z) =  • 1 + b 1 z - 1 + b 2 z - 2 • Solution: The difference equation is: • y[n] = a 0 x[n] + a 1 x[n-1] + a 2 x[n-2] - b 1 y[n-1] - b 2 y[n-2] • Signal flow graph readily deduced. • Referred to as 2nd order or "bi-quadratic" IIR section • in "direct form 1". CS3291 DSP Sectn 5

  11. Bi-quad IIR section in ‘direct form 1’ y[n] = a 0 x[n] + a 1 x[n-1] + a 2 x[n-2] - b 1 y[n-1] - b 2 y[n-2] CS3291 DSP Sectn 5

  12. To implement this signal flow-graph as a program: a0 X Y z - 1 z -1 a1 - b1 X1 Y1 z - 1 z - 1 a2 - b2 X2 Y2 CS3291 DSP Sectn 5

  13. In MATLAB using floating pt arithmetic: X1=0; X2=0; Y1=0; Y2=0; while 1 X = input(‘X=’); Y = a0*X + a1*X1 + a2*X2 - b1*Y1 - b2*Y2; disp(sprintf(‘Y = %f’ , Y) ) ; % output Y X2 = X1; X1 = X ; Y2 = Y1; Y1 = Y ; end; CS3291 DSP Sectn 5

  14. In MATLAB using fixed pt arith & shifting: K=1024; A0=round(a0*K); A1=round(a1*K); A2=round(a2*K); B1=round(b1*K); B2=round(b2*K); X1=0; X2=0; Y1=0; Y2=0; while 1 X = input(‘X=’) ; Y = A0*X + A1*X1 + A2*X2 - A1*Y1 - A2*Y2 ; Y = round(Y/K); % Divide by arith right shift disp(sprint(‘Y=%f’,Y)); % Output Y X2 = X1; X1 = X ; % Prepare for next time Y2 = Y1; Y1 = Y ; end; CS3291 DSP Sectn 5

  15. Alternative signal flow graphs: Re-order the two halves of Direct Form I Impulse response not changed (Problem 5.10). Simplify to : Direct Form II CS3291 DSP Sectn 5

  16. "Direct Form II" signal flow graph:- It is a 2nd order (bi-quad) section whose system function is: a 0 + a 1 z -1 + a 2 z - 2 H(z) =  1 + b 1 z - 1 + b 2 z - 2 Its difference equation is:- y[n] = a 0 x[n] + a 1 x[n-1] + a 2 x[n-2] - b 1 y[n-1] - b 2 y[n-2] i.e. exactly the same as Direct Form 1 CS3291 DSP Sectn 5

  17. Direct Form II has minimum no. of delay boxes. Example: Given values for a 1 , a 2 , a 0 , b 1 & b 2 , write program to implement Direct Form II using normal arithmetic. Solution: W1 = 0; W2 = 0; %For delay boxes while 1 X = input(‘X=’) ; % Input to X W =X - b1*W1 - b2*W2; % Recursive part Y = W*a0 + W1*a1 + W2*a2; % Non-rec. part W2 = W1; W1 = W; % For next time disp(sprintf(‘Y=%f’,Y)); % Output Y using disp end;% Back for next sample CS3291 DSP Sectn 5

  18. % Direct Form II in fixed point arithmetic & shifting. K=1024; A0=round(a0*K); A1=round(a1*K); A2=round(a2*K); B1=round(b1*K); B2=round(b2*K); W1 = 0; W2 = 0; %For delay boxes while 1 X = input(‘X=’) ; % Assign X to input W =K*X - B1*W1 - B2*W2; % Recursive part W =round( W / K); % By arith shift Y = W*A0+W1*A1+W2*A2; % Non-rec. part W2 = W1; W1 = W; %For next time Y = round(Y/K); %By arith shift disp(sprintf( ‘Y=%f’, Y)); %Output Y using disp end;%Back for next sample CS3291 DSP Sectn 5

  19. System function: • Expression for H(z) is ratio of polynomials in z - 1 . • Remove restriction z 1, & call H(z) "system function". • When z < 1 , H(z) need not be z-transform of {h[n]} • & need not be finite. • Nevertheless it is still useful to us. CS3291 DSP Sectn 5

  20. Poles and zeros of H(z) • For a general discrete time filter:- • a 0 + a 1 z - 1 + a 2 z - 2 + ... + a N z - N • H(z) =  • b 0 + b 1 z - 1 + b 2 z - 2 + ... + b M z - Mwith b0 = 1. • Re-express as: • (z N + a 1 z N -1 + ... + a N ) • H(z) = K z M - N • (z M + b 1 z M-1 + ... + b M )where K = a 0/b 0 • Factorise denom & numerator polynomials: • (z - z 1 )(z - z 2 )...( z - z N ) • H(z) = K z M - N • (z - p 1 )( z - p 2 )...(z - p M ) CS3291 DSP Sectn 5

  21. So z 1 , z 2 ,..., z N, are "zeros". p 1 ,p 2 ,..., p M , are "poles”. • H(z) generally infinite with z equal to a pole. • H(z) generally zerowith z equal to a zero. • For a causal stable system: H(z) must be finite for  z  1. •  Cannot be a pole with modulus  1. •  All poles must satisfy  p i < 1. • On Argand diagram: poles must lie inside unit circle. • no restriction on the positions of zeros. CS3291 DSP Sectn 5

  22. Gain response in terms of poles & zeros: • e j - z1e j - z2...ej - zN • H(e j ) = K • e j - p1e j - p2...e j - pM • Phase response: • Arg(H(e j )) = (M - N)  + arg(e j - z 1 )+ ... + arg(ej - z N ) • - arg(e j - p1 ) - ... - arg(ej - p M) • On Argand diagram: e j - p1 is distance from p1 to e j  • & Arg(e j - p1) is "pole angle" shown. • Similarly for e j - z1, e j - p2, etc. CS3291 DSP Sectn 5

  23. Imag part of z exp(j) |ej- p| Arg(ej- p) P Real part of z Fig 5.5: z - plane CS3291 DSP Sectn 5

  24. “Distance rule” Product distances from zeros to e j  H(e j ) = K  Product distances from poles to e j  “Phase rule”: Arg(H(e j )) = (M-N) + sum of zero angles - sum of pole angles. CS3291 DSP Sectn 5

  25. CS3291 DSP Sectn 5

  26. Example: Calculate gain-response of a filter with 1 + z - 1 H(z) =  1 - 1.273z - 1 + 0.81z - 2 Then plot poles & zeros of H(z) & comment on their effects: Before we start, please give (1) its difference equation & (2) its signal flow graph CS3291 DSP Sectn 5

  27. x[n] y[n] z-1 z - 1 1.273 z - 1 -0.81 • Difference equation is: • y[n] = x[n] + x[n-1] + 1.273 y[n-1] - 0.81 y[n-2] • Signal flow graph in (Direct Form 1) is: CS3291 DSP Sectn 5

  28. x[n] y[n] z-1 1.273 z-1 -0.81 • In Direct Form II, signal flow graph for:: • y[n] = x[n] + x[n-1] + 1.273 y[n-1] - 0.81 y[n-2] CS3291 DSP Sectn 5

  29. Gain response: z (z + 1) This is H(z) = hard work! z 2 - 1.273 z + 0.81 [ e j W ] [1 + e j W ]  H(e j  ) =  e 2 j  - 1.273 e j W + 0.81 H(e j )2 = H(e j ) H * (e j ) = H(e j W) H(e - j W) [1 + e j  ] [1 + e - j  ] =  [ e 2 j  - 1.27 e j  + 0.8] [ e - 2 j  - 1.27 e - j  + 0.81] 2 + 2 cos() =  3.28 - 4.61 cos() + 1.62 cos(2) CS3291 DSP Sectn 5

  30. Plot the gain response against  CS3291 DSP Sectn 5

  31. Second part of problem: • z (z + 1) • H(z) =  • z 2 - 1.27 z + 0.81 • Zeros of H(z): z = 0 & z = -1. • Poles by solving: z 2 - 1.27 z + 0.81 = 0 ( in polar form). • Assume poles are R e j  and R e - j  . Then: • z 2 - 1.27 z + 0.81 = (z - R e j  ) ( z - R e- j  ) • = z 2 - R( e j  + e - j  ) z + R 2 • = z 2 - 2 R cos () z + R 2 •  R 2 = 0.81 & 2R cos() = 1.2 Hence R = 0.9 &  = /4. • Poles of H(z): z = 0.9 e j /4 & z = 0.9 e - j/4 CS3291 DSP Sectn 5

  32. Poles & zeros plotted on z-plane Imag(z) P1 Z2 Z1 Real(z) P2 CS3291 DSP Sectn 5

  33. Pole P1 causes gain response to peak at =/4. • Zero Z 2 causes gain to be zero at  = . • Z 1 has no effect on gain (only phase) • P 2 will have little dramatic effect. • Instead of calculating gain response exactly, • we can estimate its shape by "distance rule". CS3291 DSP Sectn 5

  34. Estimating gain-response from poles & zeros:  Estimated distances Estimated gain 0 2 * 1 / [0.7 * 0.7] 4 /4 1.8 * 1 / [1.3 * 0.1] 12 /2 1 * 1.5 / [0.7 * 1.7] 1  1 * 0 / [1.5 * 1.5] 0 CS3291 DSP Sectn 5

  35. 3 dB points • Gain resp. most interesting around  = /4. • Peak caused by P 1 close to unit circle. • If  = /4 & increases slightly, • distance P1 to e j  changes dramatically. • distances from P2, Z1 & Z2 do not. • Consider effect of increasing  by 0.1 CS3291 DSP Sectn 5

  36. Increasing  by 0.1 increases pole distance from 0.1 to 0.12. • Effect of multiplying pole distance by 2: • Assuming other distances change little, gain will decrease by 2. • This is a 3 dB reduction in gain. • Same if  is decreased from /4 to /4 - 0.1. •  “3 dB pts” around peak at  = /4 are: • /4 + 0.1 & /4 - 0.1 radians/sample • Can now draw reasonably accurate sketch of G(W) against W: CS3291 DSP Sectn 5

  37. CS3291 DSP Sectn 5

  38. Comments: • Sketch not as accurate as Excel graph, but useful for design. • Assume we wanted slightly different gain response • Re-position poles and/or zeros. • Effect on gain of a zero close to the unit circle:- • Assume zero is at z = (1 - ) e j  with  small. • 3 dB pts around spectral trough at  = : at  + & -  • Gain increased to 3 dB above the minimum. CS3291 DSP Sectn 5

  39. Example: Calculate phase response of H(z). Then estimate it from pole/zero diagram. Solution: e j  / 2 [ 2 cos (/2) ] H(ej) =  e j  - 1.27 + 0.81 e - j  e j  / 2 [ 2 cos (/2) ] =  [1.81 cos( ) - 1.27 ] - j [0.19 sin( ) ] As cos (/2)  0 for 0  :  0.19 sin() Arg[H(e j )] =  - arctan  2 1.81 cos() - 1.27 Hard work as well! CS3291 DSP Sectn 5

  40. CS3291 DSP Sectn 5

  41. Estimating phase response from pole/zero plot: CS3291 DSP Sectn 5

  42. Example: An FIR filter has the difference eqn: y[n] = x[n] - 0.95 x[n-1] + 0.9 x[n-2] Plot its poles & zeros & estimate its gain (& phase) resp. Solution: H(z) = 1 - 0.95 z - 1 + 0.9 z - 2 z 2 - 0.95 z + 0.9 =  z 2 (z - 0.95 e j  / 3)(z - 0.95 e - j  / 3) =  ( z - 0 ) 2 CS3291 DSP Sectn 5

  43. Zeros at z = 0.95 e j  / 3 & 0.95 e - j  / 3 . Two poles, both at z = 0. Poles or zeros at z = 0 do not affect gain response. They do affect the phase response. Complete as an exercise. CS3291 DSP Sectn 5

  44. Design of IIR notch by pole/zero placement • By strategically locating poles & zeros in the z-plane. • 2nd order 'notch' filter to eliminate an unwanted sinusoidal component without severely affecting rest of signal. • Assume frequency of sinusoid: =/4, • Zero on unit circle at z = exp(j/4) = z1. • Must have complex conjugate at z2 = exp(-j/4). CS3291 DSP Sectn 5

  45. CS3291 DSP Sectn 5

  46. System function with these zeros is: • H(z) = (z - z 1)( z - z 2 ) = z 2 - 1.414z + 1. • Introduce double pole at z = 0 to obtain: • z 2 - 1.414z + 1. • H(z) =  • z 2 • = 1 - 1.414 z -1 + z - 2 • This is the system function of an FIR filter. • Use “distance rule” to estimate its gain-response: CS3291 DSP Sectn 5

  47.  Dist. from Z1 Dist. from Z2 Product 0 0.75 0.75 0.6 /4 0 1.4 0 /2 0.75 1.85 1.4 3/4 1.4 2 2.8  1.85 1.85 3.4 CS3291 DSP Sectn 5

  48. This is not a very good notch filter. CS3291 DSP Sectn 5

  49. Eliminates tone at  = /4, but severely affects spectrum. • Improved by placing poles p1 & p2 close to z1 & z2 CS3291 DSP Sectn 5

  50. Take p1 = 0.9 e j  / 4 , p2 = 0.9e - j  / 4 • (z - z1)(z - z2) • H(z) =  • (z - p1)(z - p2) • (z - e j  / 4 )(z - e - j  / 4 ) • =  • (z - 0.9 e j  / 4 )(z - 0.9 e - j  / 4 ) CS3291 DSP Sectn 5

More Related