ECE 445S Real-Time Digital Signal Processing Lab Fall 2023 Infinite Impulse Response Filters Prof. Brian L. Evans Dept. of Electrical and Computer Engineering The University of Texas at Austin Lecture 6 http://www.ece.utexas.edu/~bevans/courses/realtime

Outline IIR biquad filter structures Stability Bounded-Input Bounded Output Equalization Example X X O O IIR filter design IIR filter implementation Filter design demos Conclusion 6 - 2

IIR Biquad Filter Structures Discrete-Time IIR Filters Infinite Impulse Response (IIR) filter has impulse response of infinite duration, e.g. Z H(z)= 2 n=0 n n n =1+1 1 1 1 1 n h 2z-1+...= z-n= 2z-1 = [ ] u n 1-1 2 2z-1 n=0 How to implement the IIR filter by computer? Let x[n] be the input signal and y[n] the output signal, y[n]-1 1 H(z)=Y(z) 2y[n-1]= x[n] Y(z)= X(z) 1-1 X(z) 2z-1 y[n]=1 Y(z)-1 2y[n-1]+x[n] 2 z-1Y(z)= X(z) Y(z)= H(z)X(z) Recursively compute output y[n], n 0, given y[-1] and x[n] 6 - 3

IIR Biquad Filter Structures Many Equivalent Representations Difference equation [ 2 Recursive computation needs y[-1] and y[-2] For the filter to be LTI, y[-1] = 0 and y[-2] = 0 Transfer function Assumes LTI system Block diagram representation 1 1 = ] 1 + + [ ] [ ] 2 [ ] y n y n y n x n 8 x[n] y[n] Unit Delay 1/2 y[n-1] Unit Delay y[n-2] 1/8 Second-order filter section (a.k.a. biquad) with 2 poles and 0 non-trivial zeros 1 1 = + + 1 2 ( ) ( ) ( ) ( ) Y z z Y z z Y z X z 2 8 ( ) 1 Y z = = ( ) H z 1 1 ( ) X z 1 2 1 z z Poles at 0.183 and +0.683 6 - 4 2 8

IIR Biquad Filter Structures IIR Biquad Block Diagrams Biquad: short for biquadratic Transfer function is ratio of two quadratic polynomials Two poles, and zero, one or two non-trivial zeros v[n] x[n] y[n] b0 Unit Delay Feedforward paths a1 b1 Feedback paths all-pole section from previous slide v[n-1] Unit Delay v[n-2] a2 b2 Express above as cascade of two familiar filters v[n] All-pole IIR Filter FIR Filter x[n] y[n] 6 - 5

IIR Biquad Filter Structures IIR Biquad Transfer Function Assume LTI system (zero initial conditions) V(z) FIR Filter All-pole IIR Filter X(z) Y(z) H1(z) H2(z) V(z) = H1(z) X(z) and Y(z) = H2(z) V(z) Y(z) = H1(z) H2(z) X(z) and hence H(z) = H1(z) H2(z) Overall transfer function + + 1 2 ( ) ( ) ( ) b b z b z Y z V z Y z = = = 0 1 z 2 ( ) H z 1 2 ( ) ( ) ( ) 1 X z X z V z a a z 1 2 Real a1, a2: poles are conjugate symmetric ( j ) or real Real b0, b1, b2: zeros are conjugate symmetric or real 6 - 6

IIR Biquad Filter Structures IIR Biquad Difference Equation Derive difference equation from transfer function Take inverse z-transform of both sides: Collect terms to isolate y[n] on the left-hand side: How many multiplications and additions per output sample? How to convert biquad to a first-order IIR section?

IIR Biquad Filter Structures First-Order Frequency Response Pole at p0= 0.9 and zero at z0= 0.0 Angle of pole near unit circle indicates frequency at which peak occurs in magnitude response Angle of zero on or near unit circle indicates frequency at which valley occurs in magnitude response 6 - 8

Impact of Zeros on Magnitude Resp. ? ??? = ??? ?0 ??? ?1 ??? ?2 6 - 9

IIR Biquad Filter Structures IIR Biquad Frequency Responses Magnitude response Zeros z0& z1and poles p0& p1 |a b| is distance between complex numbers a and b |ej p0| is distance from point on unit circle ej and pole p0 ( )( )( ) z z z z = 0 1 ( ) H z C ( ) z p z p )( )( 0 1 ( e ) j j e z e z = j 0 1 ( ) H e C ( ) j j p e p 0 1 j j e z e z 0 1 = j ( ) H e C j j e p e p 0 1 6 - 10

Frequency Response? Poles have radius r Zeros have radius 1/r handout Im(z) Im(z) Im(z) O Zeros are on the unit circle O X X X X O O Re(z) Re(z) Re(z) X X O O Lowpass, highpass, bandpass, bandstop, allpass, notch? zeroAngle = pi/4; rz = 1/0.8; z0 = rz*exp(j*zeroAngle); z1 = rz*exp(-j*zeroAngle); numer = [1 -(z0+z1) z0*z1]; zeroAngle = pi/4; z0 = exp(j*zeroAngle); z1 = exp(-j*zeroAngle); numer = [1 -(z0+z1) z0*z1]; zeroAngle = 15*pi/16; z0 = exp(j*zeroAngle); z1 = exp(-j*zeroAngle); numer = [1 -(z0+z1) z0*z1]; r = 1/rz; poleAngle = pi/4; p0 = r * exp(j*poleAngle); p1 = r * exp(-j*poleAngle); denom = [1 -(p0+p1) p0*p1]; r = 0.9; poleAngle = pi/4; p0 = r * exp(j*poleAngle); p1 = r * exp(-j*poleAngle); denom = [1 -(p0+p1) p0*p1]; r = 0.9; poleAngle = pi/16; p0 = r * exp(j*poleAngle); p1 = r * exp(-j*poleAngle); denom = [1 -(p0+p1) p0*p1]; %%% Normalize the DC response %%% to 1 in linear units by %%% setting H(z) evaluated %%% at z = 1 to be 1 C = (denom * [1 1 1]') / (numer * [1 1 1]'); %%% Normalize the DC response %%% to 1 in linear units by %%% setting H(z) evaluated %%% at z = 1 to be 1 C = (denom * [1 1 1]') / (numer * [1 1 1]'); %%% Normalize the DC response %%% to 1 in linear units by %%% setting H(z) evaluated %%% at z = 1 to be 1 C = (denom * [1 1 1] ) / (numer * [1 1 1]'); freqz(C*numer, denom); freqz(C*numer, denom); freqz(C*numer, denom);

IIR Biquad Filter Structures Demonstrations DSP First, 2nded., ch. 10 IIR filtering tutorial (Link) Connection between the Z and frequency domains (Link) Time/frequency/z domain movies for IIR Filters (Link) link IIR filter with one pole and a zero at the origin Movement of Poles 6 - 12

Stability Stability A discrete-time LTI system is bounded-input bounded-output (BIBO) stable if for any bounded input x[n] such that | x[n] | B1< , then the filter response y[n] is also bounded | y[n] | B2< Proposition: A discrete-time filter with an impulse response of h[n] is BIBO stable if and only if = n | [ | ] n h handout Every finite impulse response LTI system (even after implementation) is BIBO stable A causal infinite impulse response LTI system is BIBO stable if and only if its poles lie inside the unit circle 6 - 13

Stability BIBO Stability Rule #1: For a causal sequence, poles are inside the unit circle (applies to z-transform functions that are ratios of two polynomials) OR Rule #2: Unit circle is in the region of convergence. (In continuous-time, imaginary axis would be in region of convergence of Laplace transform.) 1 Z n Example: n u for a z a pole at z=a 1 1 a z Stable if |a| < 1 by rule #1 or equivalently Stable if |a| < 1 by rule #2 because |z|>|a| and |a|<1 6 - 14

Stability Equalization Design equalizer to compensate frequency distortion LTI LTI System Model (Filter) h[m] x[m] y[m] Equalizer g[m] Applications: Audio, image and communication systems Goal: Make cascade all-pass: Example: H(z) is an FIR filter Im(z) notch O X cancellation O X Re(z) Stability issues? O X O X cancellation all-pass Zeros of H(z) become poles of G(z) X O notch

Cascade of Biquads IIR Filter Design Classical IIR filter design algorithms Design Method Butterworth Chebyshev type I Chebyshev type II Elliptic Passband Monotonic Monotonic Ripples Ripples Stopband Monotonic Ripples Monotonic Ripples Manual placement of poles and zeros Automated adaptation of pole and zero locations Center frequency of bandpass resonator HW 1.1(d) & 2.1(d) Notch filter with parameterized notch frequency HW 3.1(c)

Cascade of Biquads IIR Filter Implementation Cascade of biquads for N poles N even: N/2 conjugate pairs of poles. Cascade of N/2 biquads. N odd: One real pole and (N-1)/2 conjugate pairs of poles. Cascade of first-order filter (real pole) and (N-1)/2 biquads For M biquads, M! orderings of conjugate pairs of poles and for each ordering, M! orderings of conjugate pairs of zeros Exhaustive search: simulate all combinations to pick best one Quick pairing: Order biquads in ascending order of quality factors (see next slide) and for each each biquad, choose the conjugate zero pair that is closest in Euclidean distance Single section (single difference equation) Could make a BIBO stable filter become BIBO unstable

Cascade of Biquads Quality Factors for Digital Biquads Measures sensitivity in pole locations to perturbations = + 2 2 j , where For poles at a is pole radius (r < 1 for stability), with y = 2 a: + = r ) 1 ( 2 Real poles: b = 0 and 1 < a < 1, so r = |a| and y = 2 a and Q = (impulse response is C0anu[n] + C1n anu[n]) Poles on unit circle: r = 1 so Q = (oscillatory response) Imaginary poles: a = 0 so r = |b| and y = 0, and Processors with 16-bit x 16-bit multipliers and 40-bit accumulators: Qmax 40 j b = r e r a b 2 2 2 1 ( ) r y 1 where Q Q 2 2 + + 2 2 1 1 1 1 r b = = Q 2 2 2 1 2 Filter design programs use r as approximation of quality factor 1 r b

Cascade of Biquads Quality Factors for Analog Biquads Second-order filter section Conjugate symmetric poles a j b, where a < 0 With no zeros, impulse response is which is pure decay when b = 0 and pure sinusoid when a = 0 Quality factor: measures sensitivity of pole locations to perturbations Real poles: b = 0 so Q = (exponential decay response) Imaginary poles: a = 0 so Q = (oscillatory response) Maximum Q values: 25 for board-level RC circuits, 40 for switched capacitor designs, and 80 for analog IC designs Classical designs have biquads with high Q factors = + a t ( ) cos( ) h t C e b t a2+b2 -2a Q= 6 - 19

Common IIR Filter Structures IIR Filter as Single Section IIR filters have rational transfer functions 1 ) M N + + + 1 N ... ( ) ( ) b 1 b z b z Y z B z = m = k = m k ( ( ) = = = Y z a z X z b z 0 1 ( ) N H z m k 1 M ( ) ( ) X z A z a z a z 1 0 1 M Direct form realization Dot product of vector of N +1 coefficients and vector of current input and previous N inputs (FIR section) Dot product of vector of M coefficients and vector of previous M outputs ( FIR filtering of previous output values) Computation: M + N + 1 multiply-accumulates (MACs) Memory: M + N words for previous inputs/outputs and M + N + 1 words for coefficients M N = = = + [ ] [ ] [ ] y n a y n m b x n k m k 1 0 m k 6 - 20

Common IIR Filter Structures IIR Filter as Single Section y[n] x[n] b0 N = k = + [ ] [ ] y n b x n k k Unit Delay Unit Delay 0 M = m [ ] a y n m b1 a1 x[n-1] y[n-1] m 1 Unit Delay Unit Delay Full Precision Wordlength of y[0] is 2 words. Wordlength of y[n] increases with n for n > 0. b2 a2 x[n-2] y[n-2] Feed- forward Feedback Unit Delay Unit Delay M and N may be different bN aM x[n-N] y[n-M] 6 - 21

Stability As Single Section IIR filter design to meet a specification Manual or automatic placement of poles and zeros Poles are inside unit circle for stability Factored transfer function in poles, zeros and gain Expanded to obtain coefficients for difference equation Expansion for feedback coefficients can lead to instability Second-order example: zeros z0& z1, poles p0& p1 Each addition and multiplication in 64-bit floating point worst-case loss of one bit 6 - 22 worst-case loss of 53 bits

Stability As Single Section 6 - 23

Cascade of Biquads MATLAB Filter Designer Demo #1 Filter design/analysis Lowpass filter design specification (all demos) fpass = 9600 Hz fstop = 12000 Hz fsampling = 48000 Hz Apass = 1 dB Astop = 80 dB Under analysis menu Show magnitude, phase and group delay responses FIR filter equiripple Also called Remez Exchange or Parks-McClellan design Minimum order is 50 Change Wstop to 80 Order 100 gives Astop 100 dB Order 200 gives Astop 175 dB Order 300 does not converge how to get higher order filter? FIR filter Kaiser window Minimum order 101 meets spec

Cascade of Biquads MATLAB Filter Designer Demo #2 IIR filter elliptic Use second-order sections Filter order of 8 meets spec Achieved Astop of ~80 dB Poles/zeros separated in angle Zeros on or near unit circle indicate stopband Poles near unit circle indicate passband Two poles very close to unit circle but still BIBO stable Show group delay response IIR filter elliptic Use second-order sections Increase filter order to 9 Eight complex symmetric poles and one real pole: Same observations on left 6 - 25

Cascade of Biquads MATLAB Filter Designer Demo #3 Single section (Edit menu) Oscillation frequency ~9 kHz appears in passband BIBO unstable: two pairs of poles outside unit circle IIR filter elliptic Use second-order sections Increase filter order to 20 Two poles very close to unit circle but BIBO stable IIR filter design algorithms return poles-zeroes-gain (PZK format): Impact on response when expanding polynomials in transfer function from factored to unfactored form (see slide 6-25)

Cascade of Biquads MATLAB Filter Designer Demo #4 IIR filter - constrained least pth-norm design Use second-order sections Limit pole radii 0.92 Order 8 does not meet both passband/stopband specs for different weights Order 9 does indeed meet specs using Wpass = 1.5 and Wstop = 800 Filter order might increase but worth higher complexity for being more robust perturbations in pole locations 6 - 27

Conclusion FIR Filters Higher IIR Filters Implementation complexity (1) Lower (sometimes by factor of four) Elliptic design algorithm Minimum order design algorithm Parks-McClellan (Remez exchange) algorithm (2) BIBO stable? Always May become unstable when implemented (3) Linear phase If impulse response is symmetric or anti- symmetric about midpoint No, but phase may made approximately linear over passband (or other band) (1) For same piecewise constant magnitude specification (2) Algorithm to estimate minimum order for Parks-McClellan algorithm by Jim Kaiser may be off by 10%. Search for minimum order is often needed. (3) Choice of IIR filter structure matters -- use cascade of biquads instead of a single section. Some IIR filter design algorithms can constrain pole locations.

Conclusion Keep poles computed by filter design algorithms Polynomial deflation (rooting) reliable in floating-point Polynomial inflation (expansion) may degrade roots Direct form IIR structures expand zeros and poles May become unstable for large order filters (order > 12) due to degradation in pole locations from polynomial expansion Cascade of biquads (second-order sections) Only poles and zeros of second-order sections expanded Biquads placed in order of ascending quality factors Optimal ordering of biquads requires exhaustive search and each ordering can be simulated in parallel for speed 6 - 29