1 / 69

University of Manchester School of Computer Science C omp 3 0 291 : Digital Media Processing

University of Manchester School of Computer Science C omp 3 0 291 : Digital Media Processing Section 5 z-transforms & IIR-type digital filters. Introduction. General causal digital filter has difference equation:. Order is maximum of N and M.

emlyn
Download Presentation

University of Manchester School of Computer Science C omp 3 0 291 : Digital Media 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 School of Computer Science Comp 30291 : Digital Media Processing Section 5 z-transforms & IIR-type digital filters Comp30291 DMP Section 5

  2. Introduction • General causal digital filter has difference equation: • Order is maximum of N and M. • Recursive if any b j is non-zero. • A 2nd order recursive filter has the difference-equation: • 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] Comp30291 DMP Section 5

  3. y[n] x[n] + + a0 z-1 z-1 z-1 z-1 + + a1 -b1 a2 -b2 Signal-flow-graph for 2nd order recursive 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] Comp30291 DMP Section 5

  4. Derivation of ‘z-transform’ If {x[n]} with x[n] = zn is applied to an LTI system with impulse-response {h[n]}, output would be, by convolution : Comp30291 DMP Section 5

  5. ‘z-transform’ of {h[n]} • If input is {zn}, output is {H(z).zn} • z may be any real or complex number for which series converges. • For causal & stable IIR, the series converges when |z|  1. • Replacing z by ej gives frequency-response: Comp30291 DMP Section 5

  6. Imaginary part of z z = ej  1 Real part of z Visualising H(z) On Argand diagram (‘z-plane’), z = ej lies on unit circle. Comp30291 DMP Section 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 Comp30291 DMP Section 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:If x[n] = z n then y[n] = H(z) z n , • y[n-1] = H(z) z n - 1 • Substitute to obtain: • H(z) z n = a 0 z n + a 1 z n - 1 - b 1 H(z) z n - 1 • H(z) = a 0 + a 1 z - 1 - b 1 H(z) z - 1 • (1 + b1z-1) H(z) = a 0 + a 1 z - 1 • When z = -b1, H(z) =  Comp30291 DMP Section 5

  9. Find H(z) for the 2nd order difference-equation y[n] = a0 x[n] + a1x[n-1] + a2x[n-2] - b1 y[n-1] - b2 y[n-2] The same method gives: H(z) z n = a0zn + a1zn - 1 +a2zn - 2 - b1H(z) zn - 1 - b2H(z)zn-2 a 0 + a 1 z - 1 + a 2 z - 2  H(z) =  b 0 + b 1 z - 1 + b 2 z - 2 with b0 = 1. Comp30291 DMP Section 5

  10. Consider difference-equation for general digital filter • N M • y[n] =  a i x[n  i]  b j y[n  j] • i=0 j=1 • The same method gives: • This is the ‘system function’ of the digital filter. • Also referred to as ‘transfer function’ Comp30291 DMP Section 5

  11. System function • Ratio of 2 polynomials in z-1 • Equal to z-transform of impulse-response when this converges. • Easily derived from difference-equation & signal-flow graph. • Replacing z by ej gives frequency-response Comp30291 DMP Section 5

  12. z-1 z-1 z-1 z-1 y[n] x[n] + + a0 + + a1 -b1 -b2 a2 Example 5.3 Give a signal-flow graph for the system function: Solution: The difference-equation is: y[n] = a0 x[n] + a1 x[n-1] + a2 x[n-2] - b1 y[n-1] - b2 y[n-2] Represented by the signal-flow graph below: 2nd order or ‘bi-quadratic’ IIR section in ‘direct form 1’. Comp30291 DMP Section 5

  13. Y + + X a0 z-1 z-1 z-1 z-1 X1 + + Y1 a1 -b1 X2 Y2 a2 -b2 To implement ‘direct form 1’ biquad as a program Label ‘z-1’ box inputs & outputs as shown: Comp30291 DMP Section 5

  14. In MATLAB using floating point 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; Comp30291 DMP Section 5

  15. In MATLAB using Signal Processing toolbox: y = filter([a0 a1 a2], [b0 b1 b2], x ); Comp30291 DMP Section 5

  16. In MATLAB using fixed point 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; Comp30291 DMP Section 5

  17. y[n] x[n] + + a0 + + a1 -b1 z-1 z-1 z-1 z-1 -b2 a2 • Look again at ‘Direct Form 1’ signal-flow-graph • It may be thought of as two signal-flow-graphs: Non-recursive part Recursive part x[n] y[n] Comp30291 DMP Section 5

  18. L1 L2 x[n] y[n] Re-ordering LTI systems • It may be shown than if we have two LTI systems as shown: then re-ordering L1 & L2 does not change the behaviour of the overall system. L2 L1 x[n] y[n] • Only guaranteed to work for LTI systems Comp30291 DMP Section 5

  19. + + a0 + + a1 -b1 z-1 z-1 z-1 z-1 z-1 z-1 z-1 z-1 z-1 z-1 -b2 a2 + + + + x[n] x[n] y[n] y[n] a0 a0 + + + + a1 -b1 -b1 a1 a2 a2 -b2 -b2 Alternative signal-flow-graph Look again at ‘Direct Form 1’: y[n] x[n] Re-order the two ‘halves’ & then simplify to ‘direct form 2’: Comp30291 DMP Section 5

  20. W W1 z-1 z-1 W2 + + x[n] y[n] a0 + + -b1 a1 a2 -b2 ‘Direct Form II’ signal-flow-graph • Direct form II economises on ‘delay boxes’. • Notice labels: W,W1 & W2 It is a 2nd order (bi-quad) section whose system function is: 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 Comp30291 DMP Section 5

  21. Program to implement Direct Form II using normal arithmetic 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 Comp30291 DMP Section 5

  22. Direct Form II in fixed point arithmetic with 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 right-shift Y = W*A0+W1*A1+W2*A2; % Non-rec. part W2 = W1; W1 = W; %For next time Y = round(Y/K); %By arith right-shift disp(sprintf( ‘Y=%f’, Y)); %Output Y using disp end;% Back for next sample Comp30291 DMP Section 5

  23. Poles & zeros of H(z) • For a discrete time filter: • Re-express as: • Now factorise numerator & denominator: Comp30291 DMP Section 5

  24. Poles & zeros of H(z) continued • z 1 , z 2 ,..., z N, are ‘zeros’. p 1 ,p 2 ,..., pN, 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 all poles must satisfy  p i < 1. • i.e. on Argand diagram: poles must lie inside unit circle. • No restriction on the positions of zeros. Comp30291 DMP Section 5

  25. Design of IIR ‘notch’ filter • Design a 4th order 'notch' filter to eliminate an unwanted sinusoid at 800 Hz without severely affecting rest of signal. The sampling rate is FS = 10 kHz. • One way is to use the MATLAB function ‘butter’ as follows: • FL = 800 – 25 ; FU = 800+25; • [a b] = butter(2, [FL FU]/(FS/2),’stop’); • a = [0.98 -3.43 4.96 -3.43 0.98] • b= [ 1 -3.47 4.96 -3.39 0.96] • freqz(a, b); • freqz(a, b, 512, FS); % Better graph • axis([0 FS/2 -50 5]); % Scales axes • Notch has -3 dB frequency band: 25 + 25 = 50 Hz. MATLAB response Comp30291 DMP Section 5

  26. Gain/phase response of notch filter Comp30291 DMP Section 5

  27. 0 Magnitude (dB) -20 -40 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 Frequency (Hz) 0 -100 Phase (degrees) -200 -300 -400 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 Frequency (Hz) Gain/phase responses of notch filter Comp30291 DMP Section 5

  28. Details • How sharp is the notch? • Can answer this question by specifying notch’s -3 dB bandwidth. • Have just designed what Barry calls a 4th order band-stop filter. • MATLAB calls it a 2nd order band-stop filter. • For a sharper notch, decrease -3 dB bandwidth • But this will decrease its ‘depth’, i.e. the attenuation at the ‘notch’ frequency. • If necessary, increase the order to 6 (3) or 8 (4). Comp30291 DMP Section 5

  29. Gain 0 dB 1 0.5 -3 dB F 0 FS/2 800 -  800 800 +  • Sketch the same gain-response • The -3dB frequencies are at ( 800 +  ) and ( 800 -  ) Hz. • Given  (= 25 Hz say) can sketch gain-response: Comp30291 DMP Section 5

  30. Implement the 4th order ‘notch’ filter MATLAB function gave us: a = [0.98 -3.43 4.96 -3.43 0.98] b= [ 1 -3.47 4.96 -3.39 0.96] Transfer (System) Function is: Comp30291 DMP Section 5

  31. + x[n] y[n] 0.98 3.47 z-1 z-1 z-1 z-1 + + + + + + + -3.43 -4.96 4.96 3.39 -3.43 -0.96 0.0.98 A direct Form 2 implementation of the 4th order IIR notch filter Comp30291 DMP Section 5

  32. Problems with ‘direct form’ IIR implementations • Implementation on previous slide works fine in MATLAB. • But ‘direct form’ IIR implementations of order >2 are rarely used. • Sensitivity to round-off error in coeff values will be high. • Also range of ‘intermediate’ signals in z-1 boxes is high. • High wordlength floating point arithmetic hides this problem • But in fixed point arithmetic, great difficulty occurs. • Instead we use ‘cascaded biquad sections’ Comp30291 DMP Section 5

  33. x[n] y[n] H1(z) H2(z) H(z) G x[n] y[n] Using biquad (2nd order) IIR sections Given 4th order H(z), instead of: we prefer to arrange two biquad sectns as follows: Comp30291 DMP Section 5

  34. Converting to the new implementation • Get a & b for 4th order H(z) as before: [a b] = butter(2, [FL FU]/(FS/2),’stop’); • Then execute: [SOS G] = tf2sos(a,b) • MATLAB responds with: SOS = 1 -1.753 1 1 -1.722 0.9776 1 -1.753 1 1 -1.744 0.9785 G = 0.978 Transfer function to 2nd order sectns First sectn 2nd sectn Comp30291 DMP Section 5

  35. x[n] y[n] 0.978 -1.753 1.722 -1.753 1.744 -0.978 -0.979 Fourth order IIR notch filter realised as two biquad (SOS) sections H(z) may now be realised as: Comp30291 DMP Section 5

  36. Example A digital filter with a sampling rate of 200 Hz is required to eliminate an unwanted 50 Hz sinusoidal component of an input signal without affecting the magnitudes of other components too severely. Design a 4th order "notch" filter for this purpose whose 3dB bandwidth is not greater than 3.2 Hz. Solution method: FS=200; FL=50-1.6; FU=50+1.6; [a b]=butter(2,[FL,FU]/(FS/2), ‘stop’); [SOS G] = tf2sos(a,b) Comp30291 DMP Section 5

  37. IIR digital filter design by bilinear transformation • Many design techniques for IIR digital filters have adopted ideas of analogue filters. • Can transform analogue ‘prototype’ transfer function Ha(s) into H(z) for an IIRdigital filter. • Analogue filters have infinite impulse-responses. • Many gain-response approximations exist which are realisable by analogue filters • e.g. Butterworth low-pass approximation which can be transformed to high-pass, band-pass & band-stop. Comp30291 DMP Section 5

  38. Butterworth low-pass gain approximation of order n At C , Gain  0.71 i.e -3 dB Comp30291 DMP Section 5

  39. Transformations from Ha(s) to H(z) for IIR digital filter • Can transform Ha(s), with gain-response Ga(), to H(z) for an IIR digital filter with similar gain-response G(). • Many ways exist. • Most famous is ‘bilinear transformation’. • Replaces s by 2(z-1)/(z+1) to transform Ha(s) to H(z). • Fortunately MATLAB does all this for us. Comp30291 DMP Section 5

  40. Properties of bilinear transformation • (i) Order of H(z) = order of Ha(s) • (ii) If Ha(s) is causal & stable, so is H(z). • (iii) G() = Ga() where  = 2 tan(/2) • So gain of analog filter at  radians/second becomes gain of digital filter at  radians/sample where  = 2tan( /2). Comp30291 DMP Section 5

  41. Frequency warping By (iii),  from - to  mapped to  in range - to .   Comp30291 DMP Section 5

  42. Frequency warping (cont) • Shape of G() will change under the transformation. • If Ga() is Butterworth, G() will not have exactly the same shape, but we still call it Butterworth. • Mapping approx linear for  in the range -2 to 2. • As  increases above 2 , a given increase in  produces smaller and smaller increases in . Comp30291 DMP Section 5

  43. G() Ga() W w p/2 p p (a): Analogue gain response (b): Effect of bilinear transformation Comparing Ga() with G() • G() becomes more and more compressed as . • Illustrate for an analog gain-response with ripples: Comp30291 DMP Section 5

  44. ‘Prototype’ analogue transfer function • Although the shape changes, we would like G() at its cut off C to the same as Ga() at its cut-off frequency. • If Ga() is Butterworth, it is -3dB at its cut-off freq • So we would like G() to be -3 dB at its cut-off C. • Achieved if analogue prototype is designed to have its cut-off frequency at C = 2 tan(C/2). • C is the ‘pre-warped’ cut-off frequency. • Designing analog prototype with cut-off freq 2 tan(C/2) guarantees that the digital filter will have its cut-off at C. Comp30291 DMP Section 5

  45. Design 2nd order IIR lowpass digital filter by bilinear transfm • Let required cut-off frequency C = /4 radians/sample. • Need prototype transfer fn Ha(s) for 2nd order Butt low-pass filter with 3 dB cut-off at 2tan(C/2) = 2 tan(/8) radians/second. • C = 2 tan(/8) = 0.828 • I happen to remember that the transfer fn for a 2nd order Butt low-pass filter with cut-off C is: • If you don’t believe me, check that replacing s by j and taking the modulus gives G() = 1/[1+(/C)2n] with n=2. • Set C = 0.828 in this formula, then replace s by 2(z-1)/(z+1). • Gives us H(z) for an IIR digital filter. That’s it! Comp30291 DMP Section 5

  46. x[n] y[n] 0.098 2 0.94 -0.33 - Resulting IIR digital filter Comp30291 DMP Section 5

  47. Design of 2nd order IIR low-pass digital filter by bilinear transform using MATLAB • If required cut-off freq is /4 radians/sample, type: • [a b] = butter(2, 0.25) • MATLAB gives us: • a = [0.098 0.196 0.098] • b = [1 -0.94 0.33] • The required expression for H(z) is therefore: To save multipliers it is a good idea to re-express this as: Comp30291 DMP Section 5

  48. x[n] y[n] 0.09 8 2 0.94 - 0.33 Realise by ‘direct form 2’ signal-flow graph Comp30291 DMP Section 5

  49. Higher order IIR digital filters Example: Design 4th order Butterwth-type IIR low-pass digital filter with 3 dB c/o at fS / 16. . Solution: Relative cut-off frequency is /8. Typing: [a b] = butter(4, 0.125) gives the response: a = 0.0009 0.0037 0.0056 0.0037 0.0009 b = 1 -2.9768 3.4223 -1.7861 0.3556 Comp30291 DMP Section 5

  50. + x[n] y[n] 0.00093 2.977 z-1 z-1 z-1 z-1 + + + + + + + 0.0037 -3.422 0.0056 1.79 0.0037 -0.356 0.00093 A direct Form 2 implementation of 4th order IIR filter Comp30291 DMP Section 5

More Related