2.161: Lecture 20 Notes

Massachusetts Institute of Technology Department of Mechanical Engineering 2.161 Signal Processing - Continuous and Discrete Fall Term 2008 Lecture 201 Reading: 1 • Class Handout: Direct-Form Digital Filter Structures • Proakis and Manolakis: Sec. 9.1 – 9.3 • Oppenheim, Schafer, and Buck: 6.0 – 6.5 Direct-Form Filter Structures Linear shift-invariant digital filters can be represented in block diagram form in terms of the three primitive elements u n it d e la y c o e ffic ie n t ( g a in ) a fn y = a fn n z fn s u m m e r -1 y n = fn -1 f n y - g 2 n = fn - g n Transversal FIR Structure Let the FIR structure to be implemented be H(z) = N � bn z −k k=0 so that the difference equation is yn = N � bk fn−k . k=0 The following block diagram is the transversal form of this system: f z n b -1 0 f n -1 b z 1 -1 f n -2 b z -1 f n -3 b z 2 3 -1 f n -N b N y 1 c D.Rowell 2008 copyright  20–1 n n The following MATLAB code implements this structure in a point-by-point filtering function: % ------------------------------------------------------------------------­ % 2.161 Classroom Example - firdf - Demonstration FIR Direct Form % implementation. % Usage : 1) Initialization: % b = [1 2 3 4 5 4 3 2 1]; % y = iirdf1(’initial’, b); % where b are the numerator polynomial coefficients. Example: % y = iirdf1(’initial’,[1 2 5 2 1]); % Note: firdf returns y = 0 for initialization % 2) Filtering: % y_out = firdf(f); % where f is a single input value, and % y_out is the computed output value. % Example: To compute the step response: % for j=1:100 % y(j) = firdf(1); % end % ------------------------------------------------------------------------­ % function y_n = firdf(f_n,B) persistent f_register Bx N % % The following is initialization, and is executed once % if (ischar(f_n) && strcmp(f_n,’initial’)) N = length(B); Bx = B; f_register = zeros(1,N); y_n = 0; else % Filtering: y_n = 0; for J = N:-1:2 f_register(J) = f_register(J-1); y_n = y_n + Bx(J)*f_register(J); end y_n = y_n + Bx(1)*f_n; f_register(1) = f_n; end 20–2 3 IIR Direct Form Structures Let the IIR structure to be implemented be H(z) = �N 1+ k=0 � N bn z −k ak z −k k=1 where it is assumed that the orders of the numerator and denominator of H(z) are equal. The difference equation is yn = − N � N � ak yn−k + k=1 bk fn−k . k=0 Write H(z) as a pair of cascaded sub-systems, H(z) = H1 (z)H2 (z) where H1 (z) = N � bn z −k , and H2 (z) = 1+ k=0 3.1 1 �N k=1 ak z −k . Direct Form I Define an intermediate variable xn , and implement as X(z) = H1 (z)F (z) and Y (z) = H2 (z)X(z), or in difference equation form as xn = N � bk fn−k k=0 N � yn = − ak yn−k + xn k=1 as shown below: f b x 0 n n y z f -1 b a 1 -1 z 1 y n -1 z f -1 b a 2 -1 z 2 y n -2 z -1 f a 3 n -3 z f b b a N n -N 20–3 z 3 -1 z N n -1 n -2 -1 y n -3 y n -N -1 n The following MATLAB code implements the Direct Form I structure in a point-by-point filtering function. % ------------------------------------------------------------------------­ % 2.161 Classroom Example - iirdf1 - Demonstration IIR Direct Form I % implementation. % Usage : 1) Initialization: % y = iirdf1(’initial’, b, a) % where b, a are the numerator and denominator polynomial % coefficients. Example: % [b,a] = butter(7,0.4); % y = iirdf1(’initial’,b,a); % Note: iirdf1 returns y = 0 for initialization % 2) Filtering: % y_out = iirdf1(f_{in}; % where f_in is a single input value, and % y_out is the computed output value. % Example: To compute the step response: % for j=1:100 % y(j) = iirdf1(1); % end % ------------------------------------------------------------------------­ function y_n = iirdf1(f_n,B,A) persistent f_register y_register Bx Ax N % % The following is initialization, and is executed once % if (ischar(f_n) && strcmp(f_n,’initial’)) N = length(A); Ax = A; Bx = B; f_register = zeros(1,N); y_register = zeros(1,N); y_n = 0; else % Filtering: (Note that a Direct Form I filter needs two shift registers.) x = 0; y = 0; for J = N:-1:2 y_register(J) = y_register(J-1); % Move along the shift register f_register(J) = f_register(J-1); y = y - Ax(J)*y_register(J); x = x + Bx(J)*f_register(J); end x = x + Bx(1)*f_n; y_n = y + x; f_register(1) = f_n; y_register(1) = y_n; end 20–4 3.2 Direct Form II The Direct Form II structure results from reversing the order of H1 (z) and H2 (z) so that X(z) = H2 (z)F (z) and Y (z) = H1 (z)X(z), or in difference equation form as xn = − N � ak fn−k yn = k=1 N � bk xn−k . k=0 as shown below: f x n b n 0 y a a x z 3 z -1 -1 -1 z n x - b 2 b x 1 x 2 a 3 z 3 -1 b a N n -N b n -1 y n 1 z b n -2 -1 2 b x z 0 -1 z a b n -1 z a n -3 -1 x f 1 n -2 x z N b n -1 x a -1 -1 z 2 a z -1 z 1 n 3 n -3 -1 b N x N n -N From the left hand figure it can be seen that the values xn−k , k = 0, . . . N , in the two shift registers is equal, and that they can be combined to create the Direct Form II structure, as is shown on the right. The following MATLAB code implements the Direct Form II structure in a point-by-point filtering function: % ------------------------------------------------------------------------­ % 2.161 Classroom Example - iirdf2 - Demonstration IIR Direct Form II % implementation. % Usage : 1) Initialization: % y = iirdf2(’initial’, b, a) % where b, a are the numerator and denominator polynomial % coefficients. Example: % [b,a] = butter(7,0.4); % y = iirdf2(’initial’,b,a); % Note: iirdf2 returns y = 0 for initialization % 2) Filtering: % y_out = iirdf2(f_{in}; % where f_in is a single input value, and % y_out is the computed output value. % Example: To compute the step response: 20–5 % for j=1:100 % y(j) = iirdf2(1); % end % ------------------------------------------------------------------------­ % function y_n = iirdf2(f_n,B,A) persistent register Bx Ax N % % The following is initialization, and is executed once % if (ischar(f_n) && strcmp(f_n,’initial’)) N = length(A); Ax = A; Bx = B; register = zeros(1,N); y_n = 0; else % Filtering: (Note that a Direct Form II filter needs only a single % shift register.) x = 0; y = 0; for J = N:-1:2 register(J) = register(J-1); % Move along the shift register x = x - Ax(J)*register(J); y = y + Bx(J)*register(J); end x = x + f_n; y_n = y + Bx(1)*x; register(1) = x; end 4 Transposed Direct Forms The transposed forms result from the transposition theorem from signal-flow graph theory, which states that in a signal-flow graph if • The arrows on all graph branches are reversed. • Branch points become summers, and summers become branch points. • The input and output are swapped, then the input/output relationships remain unchanged. The same applies to block diagrams. 4.1 Transposed Transversal FIR Filter The transposed FIR structure is shown below: 20–6 f x n b n 0 y n f b n 0 -1 z x b n -1 y -1 z 1 b n 1 -1 z x z b n -2 2 b -1 b x x 2 3 z n -3 b z -1 z -1 3 -1 x b N n -N b z -1 N T ra n s p o s e d T ra n s v e rs a l F IR T ra n s v e rs a l F IR % ------------------------------------------------------------------------­ % 2.161 Classroom Example - firtdf - Demonstration Transposed FIR Direct % Form implementation. % Usage : 1) Initialization: % y = firtdf(’initial’, b) % where b, a are the numerator and denominator polynomial % coefficients. Example: % b = [1 2 3 4 5 4 3 2 1]; % y = firtdf(’initial’,b); % % Note: firtdf returns y = 0 for initialization % 2) Filtering: % y_out = firtdf(f_{in}); % where f_in is a single input value, and % y_out is the computed output value. % Example: To compute the step response: % for j=1:100 % y(j) = firtdf(1); % end % ------------------------------------------------------------------------­ % function y_n = firtdf(f_n,B) persistent register Bx N % % The following is initialization, and is executed once % if (ischar(f_n) && strcmp(f_n,’initial’)) N = length(B); Bx = B; register = zeros(1,N-1); 20–7 y_n = 0; else % Filtering: y_n = register(1) + Bx(1)*f_n; % Update for the next iteration for J = 1:N-2 register(J) = register(J+1) + Bx(J+1)*f_n; end register(N-1) = Bx(N)*f_n; end 4.2 Transposed Direct Form II The following diagram shows the result when the transposition theorem is applied to a Direct Form II structure. f n x - x 1 x 2 a y f n b n 0 b n -1 3 z -1 z 1 b b n -2 a 1 - 2 b 1 x 2 z n -3 b x 3 a - 3 N n -N b D ir e c t F o r m 2 -1 -1 N a - 3 b z -1 z b -1 x a y -1 z a 0 -1 z a b n II z N -1 - a N T r a n s p o s e d D ir e c t F o r m II This block diagram simply reorganizes the difference equation as yn = b0 fn + N � (bk fn−k − ak yn−k ) k=1 which is implemented in the MATLAB function iirtdf2() on the next page. 20–8 n % ------------------------------------------------------------------------­ % 2.161 Classroom Example - iirtdf2 - Demonstration Transposed IIR Direct % Form II implementation. % Usage : 1) Initialization: % y = iirtdf2(’initial’, b, a) % where b, a are the numerator and denominator polynomial % coefficients. Example: % [b,a] = butter(7,0.4); % y = iirtdf2(’initial’,b,a); % Note: iirdf2 returns y = 0 for initialization % 2) Filtering: % y_out = iirtdf2(f_{in}; % where f_in is a single input value, and % y_out is the computed output value. % Example: To compute the step response: % for j=1:100 % y(j) = iirtdf2(1); % end % ------------------------------------------------------------------------­ % function y_n = iirtdf2(f_n,B,A) persistent register Bx Ax N % % The following is initialization, and is executed once % if (ischar(f_n) && strcmp(f_n,’initial’)) N = length(A); Ax = A; Bx = B; register = zeros(1,N-1); y_n = 0; else % Filtering: (Note that a Transposed Direct Form II filter needs only a single % register.) Also note that this is not strictly a shift register. y_n = register(1) + Bx(1)*f_n; % Update for the next iteration for J = 1:N-2 register(J) = register(J+1) + Bx(J+1)*f_n - Ax(J+1)*y_n; end register(N-1) = Bx(N)*f_n - Ax(N)*y_n; end 5 Coefficient Sensitivity in Direct Form Filters In practice high-order IIR Direct Form filters are rarely used because of the sensitivity of pole and zero positions to small perturbations in the values of the coefficients ak and bk in the difference 20–9 equation. If the transfer function is �M −1 A(z) k=0 bk z , = H(z) = � −1 B(z) 1+ N k=1 ak z and the poles are clustered near the unit circle, then small perturbations in any of the ak from the desired value (perhaps because of finite precision limitations) may cause the filter to become unstable. To demonstrate this, consider a low-pass filter with A(z) = 1 + N � ak z −1 = k=1 N � � 1 − pk z −1 � k=1 where the poles pk are within the unit circle, but close to z = 1, and write pk = 1 + k , where |k |  1. Now let a single (arbitrary) coefficient ar be be perturbed by δ to ar = ar + δ so that the denominator polynomial becomes  A (z) = 1 + N � ak z −1 + δz −r . k=1 As |δ| increases, one or more of the poles may move outside the unit circle, leading to instability. It is difficult to define the general condition, but we can easily find the condition that leads to a pole migrating to z = 1, since then A (1) = A(1) + δ = 0, that is, there will be a pole at z = 1 if δ = −A(1), or alternatively, if δ= N � (−(k). k=1 Example 1 Consider a low-pass filter H(z) = 1 1 = (1 − 0.99z −1 ) 1 − 2.97z −1 + 2.9403z −2 − 0.970299z −3 with three poles at z = 0.99. Find the perturbation allowed in any coefficient that will create a marginally stable system with a pole at z = 1. Discuss some methods of decreasing the sensitivity. 20–10 Solution: For the third-order system A(1) = −10−6 , so any change of δ = −A(1) = 10−6 in any coefficient will move one of the poles from z = 0.99 to z = 1. Any perturbation larger than this will generate an unstable filter. Now consider the effect of implementing this filter as a cascade connection of two filters, a second-order filter H1 (z), and a first-order filter H2 (z), that is H(z) = H1 (z)H2 (z) = 1 1 . −1 2 (1 − 0.99z ) 1 − z −1 with a pair of difference equations xn = 1.98xn−1 − 0.9801xn−2 + fn yn = 0.99yn−1 + xn . For H1 (z), A1 (1) = −10−4 , while for H2 (z) A2 (1) = −10−2 and the sensitivity is significantly reduced. If the filter is implemented as a cascade connection of three first-order filters, 1 1 1 . . H(z) = H1 (z)H2 (z)H(3 z)(z) = 1 − z −1 1 − z −1 1 − z −1 with a set of difference equations wn = 0.99wn−1 + fn vn = 0.99vn−1 + wn yn = 0.99yn−1 + vn , for any of the first-order sections Hk (z), A1 (1) = −10−2 , and the coefficient sensitivity is significantly reduced even further. This example demonstrates that the sensitivity to coefficient precision can be often drastically reduced by implementing a filter with low order sections. 5.1 Cascade Structures If the transfer function is written in terms of its poles and zeros �M2 � M1 −1 −1 −1 K=1 (1 − ek z ) K=1 (1 − gk z )(1 − g k z ) H(z) = �N � N 1 2 −1 −1 −1 K=1 (1 − ck z ) K=1 (1 − dk z )(1 − dk z ) where the ck and ek are real poles and zeros, and dk , dk and gk , g k are complex conjugate pole and zero pairs, it is common to realize the system as a cascade chain of first- and second-order sections (usually Direct Form II): A first-order Direct Form II section, implementing a real pole ck and zero ek Hk (z) = 1 − ek z −1 , 1 − ck z −1 as yn = ck yn−1 + fn − ek fn−1 is shown below 20–11 x f n z c n y -1 x k n -e n -1 k A second-order Direct Form 2 section, implementing a conjugate pole pair dk , dk = r e±j θ has a denominator polynomial (1 − dk z −1 )(1 − dk z −1 ) = 1 − 2r cos(θ)z −1 + r2 z −2 and when paired with a pair of zeros (either real, or a complex conjugate pair) to give a transfer function 1 − b1 z −1 + b2 z −2 H(z) = 1 − 2r cos(θ)z −1 + r2 z −2 and difference equation yn = 2r cos(θ)yn−1 + r2 yn−2 + fn + b1 fn−1 + b2 fn−2 is shown below f x n 2 r c o s (G ) -r 2 n y n -1 z x z n -1 b 1 -1 b x 2 n -2 Example 2 Implement the system H(z) = 0.04756z 3 + 0.14273z 2 + 0.14273z + 0.04756 z 3 − 1.3146z 2 + 1.17043z − 0.47524 as a set of cascaded first- and second-order systems. Solution: Factor the transfer function and rewite as H(z) = 0.04756(1 + z −1 )3 . (1 − 0.6711z −1 + 0.7386z −2 )(1 − 0.6435z −1 ) 20–12 Implement the filter as a cascaded pair 0.04756(1 + 2z −1 + z −2 ) 1 − 0.6711z −1 + 0.7386z −2 1 + z −1 H2 (z) = 1 − 0.6435z −1 H1 (z) = with a pair of difference equations xn = 0.6711xn−1 − 0.7386xn−2 + 0.04756(fn + 2fn−2 + fn−2 ) yn = 0.6435yn−1 + xn + xn−1 . There is a lot of flexibility in choosing which zeros to associate with the poles of each low order section, and how to distribute the overall gain between the sections. A general (rule-of-thumb) procedure is (1) Select the poles closest to the unit circle. (2) Find the closest zeros to those poles. (3) Combine into a second-order section. (4) Repeat until all zeros are accounted for. 5.2 Parallel Structures A high order filter may also be realized as a set of parallel second- and first-order sections using partial fractions, and implemented as H(z) = H1 (z) + H2 (z) + H3 (z) . . . + HN (z) so that Y (z) = (H1 (z) + H2 (z) + H3 (z) . . . + HN (z)) F (z) f H n H (z ) 2 H H (z ) 1 N 3 (z ) (z ) 20–13 y n Example 3 Implement the system of Example 2 as a parallel realization H(z) = 0.04756z 3 + 0.14273z 2 + 0.14273z + 0.04756 z 3 − 1.3146z 2 + 1.17043z − 0.47524 Solution: Using a partial fraction expansion H(z) = 0.0476 + 0.2929 .0877z − 0.2271 − 2 . z − 0.6435 z − 0.6711z + 0.7386 Implement as three sections H1 (z) = 0.0476 0.2929z −1 H2 (z) = 1 − 0.6435z −1 .0877z −1 − 0.2271z −2 H3 (z) = − 1 − 0.6711z −1 + 0.7386z −2 with difference equations un vn wn yn = = = = 0.0476fn 0.6435vn−1 + 0.2929fn−1 0.6711wn−1 − 0.7386wn−2 + 0.0877fn−1 − 0.2271fn−2 un + vn − wn 20–14