Digital Signal Processing Moslem Amiri, Václav Přenosil Masaryk University Understanding Digital Signal Processing, Third Edition, Richard Lyons (0-13-261480-4) © Pearson Education, 2011. Infinite Impulse Response Filters 2 Introduction  Infinite impulse response (IIR) digital filters  Are fundamentally different from FIR filters  FIR filter output samples depend only on past input samples  Each IIR filter output sample depends on previous input samples and previous filter output samples  IIR filters have memory of past outputs (require feedback)  As in all feedback systems, perturbations at IIR filter input could cause filter output to become unstable and oscillate indefinitely  Infinite impulse response  IIR filters have more complicated structures (block diagrams), are harder to design and analyze, and do not have linear phase responses 3 Introduction  Why use an IIR filter?  Because they are very efficient  IIR filters require far fewer multiplications per filter output sample to achieve a given frequency magnitude response  From a hardware standpoint, IIR filters can be very fast, allowing us to build real-time IIR filters that operate over much higher sample rates than FIR filters 4 Introduction where the 19-tap FIR filter requires 19 multiplications per filter output sample, the 4th-order IIR filter requires only 9 multiplications for each filter output sample reduced passband ripple and a sharper filter rolloff, with less than half the multiplication workload 5 Introduction  IIR vs. FIR  An FIR filter’s frequency response with very steep transition regions requires a very long impulse response  The maximum number of FIR filter taps we can have (length of impulse response) depends on how fast our hardware can perform the required number of multiplications and additions to get a filter output before the next input sample arrives  IIR filters can be designed to have impulse responses longer than their number of taps  Thus, IIR filters can give much better filtering for a given number of multiplications per output sample 6 An Introduction to IIR Filters  If IIR filter’s input suddenly becomes all zeros, its output could remain nonzero forever  This is because of feedback structure of their delay units, multipliers, and adders 7 An Introduction to IIR Filters )3()3()2()2( )1()1()()0()(   nxhnxh nxhnxhny difference equation 8 An Introduction to IIR Filters )3()3()2()2()1()1( )3()3()2()2()1()1()()0()(   nyanyanya nxbnxbnxbnxbny difference equation describing this IIR filter 9 An Introduction to IIR Filters  How determine a(k) and b(k) IIR filter coefficients  Window method of lowpass FIR filter design  Define frequency response of desired FIR filter  take inverse Fourier transform  shift that transform result  we get filter’s time-domain impulse response  Due to the nature of transversal FIR filters, the desired h(k) filter coefficients turn out to be exactly equal to the impulse response sequence  Following that same procedure with IIR filters  Desired frequency response of IIR filter  inverse Fourier transform  time-domain impulse response  But there’s no direct method for computing IIR filter’s a(k) and b(k) coefficients from impulse response  FIR filter design techniques cannot be used here 10 An Introduction to IIR Filters  Standard IIR filter design techniques  Fall into three basic classes: impulse invariance, bilinear transform, and optimization methods  These design methods use a discrete sequence, mathematical transformation process known as the z-transform whose origin is Laplace transform used in analyzing of continuous systems 11 The Laplace Transform  Laplace transform  A mathematical method of solving linear differential equations  Process of using Laplace transform  A time-domain differential equation is written that describes input/output relationship of a physical system  The differential equation is Laplace transformed, converting it to an algebraic equation  Standard algebraic techniques are used to determine desired output function’s equation in Laplace domain  The desired Laplace output equation is, then, inverse Laplace transformed to yield the desired time-domain output function’s equation 12 The Laplace Transform  Laplace transform of a continuous timedomain function f(t)  Variable s is the complex number s = σ + jω  f(t) is defined only for positive time (t > 0)  Systems where system conditions for negative time (t < 0) are not needed (one-sided ) are referred to as causal systems  Causal systems may have initial conditions at t = 0 that must be taken into account, but we don’t need to know what the system was doing prior to t = 0     0 )()( dtetfsF st 13 The Laplace Transform  s = σ + jω  σ is a real number  ω is frequency in radians/second  e−st is dimensionless  s has dimension of 1/time, or frequency  Laplace variable s is often called a complex frequency 14 The Laplace Transform  Laplace transform  e−jωt is a unity-magnitude phasor rotating clockwise around origin of a complex plane at a frequency of ω radians/second  eσt is a real number whose value is one at t = 0  As t increases, eσt gets larger (when σ is positive), and complex e−st phasor’s magnitude gets smaller as phasor rotates on complex plane  The tip of that phasor traces out a curve spiraling in toward origin of complex plane ttt tj tjttjst st e t j e t e e eeee dtetfsF     )sin()cos( )()( )( 0        15 The Laplace Transform real part of F(s), for a particular value of s, is correlation of f(t) with a cosine wave of frequency ω and a damping factor of σ, and imaginary part of F(s) is correlation of f(t) with a sinewave of frequency ω and a damping factor of σ 16 The Laplace Transform  s-plane  If we associate each of different values of complex s variable with a point on a complex plane, called s-plane, we could plot real part of F(s) correlation as a surface above (or below) that s-plane and generate a second plot of imaginary part of F(s) correlation as a surface above (or below) s-plane  We can also graph magnitude |F(s)| as a function of s 17 The Laplace Transform )( )( )( )()( 01012 2 2 txb dt tdx btya dt tdy a dt tyd a  We’ll use Laplace transform toward our goal of figuring out what the y(t) output will be for any given x(t) input 18 The Laplace Transform )()()()( )( )( )( )( )(functiontransfer )()()()(or )()()()()( )( )( )( )()( )( ,..., )( , )( , )( 01 2 2 01 01 2 2 01 0101 2 2 0101 2 2 and,offunctionsbeandletweIf 01012 2 2 3 3 3 2 2 2 sHsX asasa bsb sXsY asasa bsb ex ey sX sY sH exbsbeyasasa exbesxbeyaesyaeysa txb dt tdx btya dt tdy a dt tyd a es dt ed es dt ed es dt ed se dt ed ts ts tsts tststststs )y(e)x(eey(t)x(t) tsn n tsn ts ts ts ts ts ts ststst              19 The Laplace Transform )()()()( )( )( )( )( )(functiontransfer 01 2 2 01 01 2 2 01 sHsX asasa bsb sXsY asasa bsb ex ey sX sY sH ts ts        20 The Laplace Transform  Laplace analysis technique is based on system’s x(t) input being some function of est, or x(est)  All practical x(t) input functions can be represented with complex exponentials, e.g.,  A constant, c = ce0t  Sinusoids, sin(ωt) = (ejωt − e−jωt)/2j or cos(ωt) = (ejωt + e−jωt)/2  A monotonic exponential, eat  An exponentially varying sinusoid, e−at cos(ωt) 21 The Laplace Transform  System’s transfer function H(s)  If we know H(s), we can take Laplace transform of any x(t) input to determine X(s), multiply that X(s) by H(s) to get Y(s), and then inverse Laplace transform Y(s) to yield time-domain expression for the output y(t)  Not needed because it’s H(s) in which we’re interested  Being able to express H(s) mathematically or graph the surface |H(s)| as a function of s will tell us two important properties we need to know about system:  Is system stable  And if so, what is its frequency response 22 The Laplace Transform  Poles and zeros on s-plane and stability  A system is stable if, given any bounded input, the output will always be bounded  Instability would result in a filter output that’s not at all representative of filter input—a situation we’d like to avoid if we can  Example  If s = −a0/a1, denominator equals zero and H1(s) would have an infinite magnitude  This s = −a0/a1 point on s-plane is called a pole 10 10 01 0 1 / / )( )( )(functiontransferLaplacesSystem' aas ab asa b sX sY sH     23 The Laplace Transform the pole is located exactly on negative portion of real σ axis H1(s) is stable because its y(t) output approaches zero as time passes distance of pole from σ = 0 axis, a0/a1 for our H1(s), gives the decay rate of y(t) impulse response If the system described by H1 were at rest and we disturbed it with an impulse like x(t) input at time t = 0, its continuous time-domain y(t) output would be the damped exponential curve shown in (b) 24 The Laplace Transform Laplace transform is a more general case of Fourier transform because if σ = 0, then s = jω intersection of vertical σ = 0 plane (jω axis plane) and |H1(s)| surface, gives us frequency magnitude response |H1(ω)| of system 25 The Laplace Transform  Example: 2nd-order transfer function H2(s)  Order of transfer function is the largest exponential order of s in either numerator or denominator  If s is equal to −p or −p*, one of polynomial roots in the denominator will equal zero, and H2(s) will have an infinite magnitude  two complex poles                              2 2 2 2 2 imagrealimagreal21 2 2021 2 21 01 2 2 1 2 4 4 24 4 2 )( ,,/where ))(( )( /)/( )/( )( )( )( a acb a b s a acb a b scbsassf jpppjpppabA psps As sH aasaas sab asasa sb sX sY sH 26 The Laplace Transform the two complex poles are located off the negative portion of real σ axis If H2 system were at rest and we disturbed it with an impulse-like x(t) input at t = 0, its continuous time-domain y(t) output would be the damped sinusoidal curve shown in (b) H2(s) is stable because its oscillating y(t) output approaches zero as time increases distance of poles from σ = 0 axis (−preal) gives decay rate of sinusoidal y(t) impulse response. Likewise, distance of poles from jω = 0 axis (±pimag) gives frequency of sinusoidal y(t) impulse response 27 The Laplace Transform 28 The Laplace Transform stable systems: when disturbed from their at-rest condition, they respond and, at some later time, return to that initial condition 1/s transfer function conditional stability: if an H(s) transfer function has conjugate poles located exactly on jω axis (σ = 0), the system will go into oscillation when disturbed from its initial condition instability: the poles lie to the right of jω axis. When disturbed from their initial atrest condition by an impulse input, their outputs grow without bound 29 The Laplace Transform for a system to be stable, all of its transfer-function poles must lie on the left half of s-plane 30 The z-Transform  z-transform  Discrete-time cousin of continuous Laplace transform  While Laplace transform is used to simplify analysis of continuous differential equations, z-transform facilitates analysis of discrete difference equations  z-transform is performed on a discrete h(n) sequence, converting that sequence into a continuous function H(z) of the continuous complex variable z      n n znhzH )()( 31 The z-Transform  z-transform  This Equation can be interpreted as Fourier transform of product of original sequence h(n) and exponential sequence r−n  When r = 1, it simplifies to Fourier transform  Thus on z-plane, the contour of H(z) surface for those values where |z| = 1 is Fourier transform of h(n)  If h(n) represents a filter impulse response sequence, evaluating H(z) transfer function for |z| = 1 yields frequency response of filter               n njn n njjrez n n ernhrenhreHzH znhzH j )()())(()()( )()(  32 The z-Transform 33 The z-Transform jω frequency axis on continuous Laplace s-plane is linear and ranges from − ∞ to + ∞ radians/second. The ω frequency axis on complex z-plane, however, spans only the range from −π to +π radians  z-plane frequency axis is equivalent to coiling s-plane’s jω axis about the unit circle on z-plane 34 The z-Transform  Poles, zeros, and digital filter stability  Region of filter stability is mapped to the inside of unit circle on z-plane  Given H(z) transfer function of a digital filter, we can examine that function’s pole locations to determine filter stability  If all poles are located inside unit circle, filter will be stable 35 The z-Transform  Example  If a causal filter’s H(z) transfer function has a single pole at location p on z-plane, its transfer function can be represented by  Filter’s time-domain impulse response sequence  u(n) represents a unit step (all ones) sequence beginning at time n = 0  When |p| < 1, h(n) impulse response sequence is unconditionally bounded in amplitude  |p| < 1 means that pole must lie inside z-plane’s unit circle 1 1 1 )(    pz zH )()( nupnh n  36 The z-Transform point z = −1 corresponds to π radians (or πfs radians/second = fs/2 Hz)  ωo = π/4 means that fo = fs/8 and our y(n) will have eight samples per cycle of fo 37 Using z-Transform to Analyze IIR Filters  Representing delay operation  Assume we have a sequence x(n) whose ztransform is X(z) and a sequence y(n) = x(n−1) whose z-transform is Y(z)  Thus, effect of a single unit of time delay is to multiply z-transform of undelayed sequence by z−1 )]([)( )()()( )1()()( 1)(1 1)1(1letweif zXzzkxz zzkxzkxzY znxznyzY k k k k k knk n n nn                        38 Using z-Transform to Analyze IIR Filters X(z)z−k is z-transform of x(n) delayed by k samples. So a transfer function of z−k is equivalent to a delay of kts seconds from the instant when t = 0, where ts = 1/fs Because a delay of one sample is equivalent to factor z−1, the unit time delay symbol is usually indicated by z−1 operator 39 Using z-Transform to Analyze IIR Filters 40 Using z-Transform to Analyze IIR Filters  Fig. 6-17  Is a general Mth-order IIR filter  This IIR filter structure is called Direct Form I                                    M k k N k k N k k M k k M k k N k k M N zka zkb zX zY zHzkbzXzkazY zkazYzkbzXzY zzYMazzYazzYa zzXNbzzXbzzXbzXbzY MnyManyanya NnxNbnxbnxbnxbny 1 0function transfer 01 10 21 21 )(1 )( )( )( )()()()(1)( )()()()()( )()(...)()2()()1( )()(...)()2()()1()()0()( )()(...)2()2()1()1( )()(...)2()2()1()1()()0()( order of H(z) and order of filter = the largest exponential order of z in either numerator or denominator 41 Using z-Transform to Analyze IIR Filters  Two things to know about an IIR filter: its frequency response and stability  We can evaluate denominator of H(z) to determine positions of filter’s poles on z-plane indicating filter’s stability  From H(z) we develop an expression for IIR filter’s frequency response  H(z) is a complex-valued surface above, or below, the z-plane  Intersection of H(z) surface and perimeter of a cylinder representing z = ejω unit circle is the filter’s complex frequency response 42 Using z-Transform to Analyze IIR Filters                            M k M k N k N kje M k jk N k jk jM k k N k k kkajkka kkbjkkb H eka ekb ez zHH zka zkb zX zY zH j 11 00)sin()cos( 1 0 1 0 )sin()()cos()(1 )sin()()cos()( )( )(1 )( )()( )(1 )( )( )( )(          43 Using z-Transform to Analyze IIR Filters 44 Using z-Transform to Analyze IIR Filters  Fig. 6-18  (a) is a 2nd-order lowpass IIR filter whose positive cutoff frequency is ω = π/5 (fs/10 Hz) )]2sin(436.0)1sin(194.1[)2cos(436.0)1cos(194.11 )]2sin(0605.0)1sin(121.0[)2cos(0605.0)1cos(121.00605.0 )( 436.0194.11 0605.0121.00605.0 )( 436.0194.11 0605.0121.00605.0 )( )( )( )(436.0)(194.1 )(0605.0)(121.0)(0605.0)( )2(436.0)1(194.1 )2(0605.0)1(121.0)(0605.0)( 21 210 21 210 21 21                          j j H ee eee H zz zzz zX zY zH zzYzzY zzXzzXzXzY nyny nxnxnxny jj jjj 45 Using z-Transform to Analyze IIR Filters although both filters require the same computational workload, five multiplications per filter output sample, lowpass IIR filter has the superior frequency magnitude response phase nonlinearity is inherent in IIR filters knowing that the filter’s phase response is nonlinear, we should expect the impulse response to be asymmetrical infinite impulse response: if we used infiniteprecision arithmetic in our filter implementation, h(k) impulse response would be infinite in duration 46 Using z-Transform to Analyze IIR Filters  To determine our IIR filter’s stability  So when z = p0 = 0.597 − j0.282, or when z = p1 = 0.597 + j0.282, filter’s H(z) transfer function’s denominator is zero and |H(z)| is infinite  Because those pole locations are inside the unit circle (their magnitudes are less than one), our example IIR filter is unconditionally stable )282.0597.0)(282.0597.0( )1)(1( ))(( ))(( )( 436.0194.1 0605.0121.00605.0 )( 436.0194.11 0605.0121.00605.0 )( )( )( 10 10 2 2 bymultiply 21 210 22 jzjz zz pzpz zzzz zH zz zz zH zz zzz zX zY zH /zzH(z)                47 Using z-Transform to Analyze IIR Filters if we were to unwrap the bold |H(ω)| curve and lay it on a flat surface, we would have the |H(ω)| curve in Fig. 6-19(a) 48 Using Poles and Zeros to Analyze IIR Filters  IIR filter transfer function algebra  Several ways to write H(z) = Y(z)/X(z) z-domain transfer function ))()()(( ))()()(( )( )4()3()2()1( )4()3()2()1()0( )( )4()3()2()1(1 )4()3()2()1()0( )( 3210 3210formfactored 234 234 )formpolynomial( /bygmultiplyin 4321 4321 44 pzpzpzpz zzzzzzzz zH azazazaz bzbzbzbzb zH zazazaza zbzbzbzbb zH zz              useful because we can replace z with ejω to obtain an expression for frequency response of filter necessary so we can factor (find roots of) polynomials to obtain values (locations) of numerator zeros and denominator poles 49 Using Poles and Zeros to Analyze IIR Filters  Using poles/zeros to obtain transfer functions  We can analyze an IIR filter’s frequency-domain performance based solely on poles and zeros  Given zk zeros and pk poles, we can write the factored form of filter’s transfer function as  G = G1/G2 is an arbitrary gain constant  Filter zeros are associated with decreased frequency magnitude response, and filter poles are associated with increased frequency magnitude response )...)()()()(( )...)()()()(( )...)()()()(( )...)()()()(( )( 43210 43210 432102 432101 pzpzpzpzpz zzzzzzzzzzG pzpzpzpzpzG zzzzzzzzzzG zH       50 Using Poles and Zeros to Analyze IIR Filters if a filter has no z-plane zeros, and one pole at p0 = 0.8, we can write its transfer function as H1(z)=G/(z-0.8) |H1(ω)| is normalized. P0 is closest to ω = 0 radians/sample (z = 1) on unit circle  lowpass filter. |p0| < 1  filter is unconditionally stable if a filter has a zero at z0 = 1, and a pole at p0 = −0.8, we write its transfer function as 8.0)8.0( )1( )(2       z GGz z zG zH pole is closest to ω = π radians/sample (z = −1) on unit circle  highpass filter. the zero located at z = 1 causes filter to have infinite attenuation at ω = 0 radians/sample (zero Hz) 51 Using Poles and Zeros to Analyze IIR Filters  Fig. 6-21(c)  Consider a filter having two complex conjugate zeros at −0.707 ± j0.707, as well as two complex conjugate poles at 0.283 ± j0.283  The two poles on the right side of z-plane make this a lowpass filter having a wider passband than H1(z)  Two zeros are on unit circle at angles of ω = ±3π/4 radians, causing filter to have infinite attenuation at frequencies ω = ±3π/4 radians/sample (±3fs/8 Hz) )283.0283.0()283.0283.0( )707.0707.0()707.0707.0( )]283.0283.0([)]283.0283.0([ )]707.0707.0([)]707.0707.0([ )(3 jzjz jzjzG jzjz jzjzG zH       52 Using Poles and Zeros to Analyze IIR Filters 53 Using Poles and Zeros to Analyze IIR Filters  Fig. 6-21(d)  If we add a z-plane zero at z = 1 to H3(z)  The zero at z = 1 yields infinite attenuation at ω = 0 radians/sample (zero Hz), creating a bandpass filter  Because p0 and p1 poles are oriented at angles of θ = ±π/4 radians, filter’s passbands are centered in the vicinity of frequencies ω = ±π/4 radians/sample (±fs/8 Hz) )283.0283.0()283.0283.0( )707.0707.0()707.0707.0()1( )(4 jzjz jzjzzG zH    54 Using Poles and Zeros to Analyze IIR Filters  Fig. 6-21(e)  If we increase magnitude of H4(z) filter’s poles, making them equal to 0.636 ± j0.636, we position the conjugate poles much closer to unit circle  Poles near unit circle now have a much more profound effect on filter’s magnitude response  The poles’ infinite gains cause H5(z) passbands to be very narrow (sharp)  When a pole is close to unit circle, center frequency of its associated passband can be accurately estimated to be equal to pole’s angle )636.0636.0()636.0636.0( )707.0707.0()707.0707.0()1( )(5 jzjz jzjzzG zH    55 Using Poles and Zeros to Analyze IIR Filters  Fig. 6-21(f)  Consider an FIR filter—a digital filter whose H(z) transfer function denominator is unity  For an FIR filter to have linear phase, each z-plane zero located at z = z0 = Mejα, where M ≠ 1, must be accompanied by a zero having an angle of −α and a magnitude of 1/M  z0 is accompanied by z3  If FIR filter’s transfer function polynomial has realvalued bk coefficients, then a z0 zero not on the zplane’s real axis will be accompanied by a complex conjugate zero at z = z2  Likewise, for FIR filter to have linear phase, z2 zero must be accompanied by z1 zero  z1 and z3 zeros are complex conjugates of each other 56 Using Poles and Zeros to Analyze IIR Filters  z-plane pole/zero properties  Filter poles are associated with increased frequency magnitude response (gain)  Filter zeros are associated with decreased frequency magnitude response (attenuation)  To be unconditionally stable, all filter poles must reside inside the unit circle  Filter zeros do not affect filter stability  The closer a pole (zero) is to unit circle, the stronger will be its effect on filter’s gain (attenuation) at the frequency associated with the pole’s (zero’s) angle 57 Using Poles and Zeros to Analyze IIR Filters  z-plane pole/zero properties  A pole (zero) located on unit circle produces infinite filter gain (attenuation)  If a pole is at the same z-plane location as a zero, they cancel each other  Poles or zeros located at origin of z-plane do not affect frequency response of filter  Filters whose transfer function denominator (numerator) polynomial has real-valued coefficients have poles (zeros) located on real zplane axis, or pairs of poles (zeros) that are complex conjugates of each other 58 Using Poles and Zeros to Analyze IIR Filters  z-plane pole/zero properties  For an FIR filter (transfer function denominator is unity) to have linear phase, any zero on z-plane located at z0 = Mejα, where z0 is not on unit circle and α is not zero, must be accompanied by a reciprocal zero whose location is 1/z0 = e−jα/M  If an FIR filter has real-valued coefficients, is linear phase, and has a z-plane zero not located on real z-plane axis or on unit circle, that z-plane zero is a member of a “gang of four” zeros  If we know z-plane location of one of those four zeros, then we know location of the other three zeros 59 Alternate IIR Filter Structures  Direct Form I, Direct Form II, and transposed structures  Direct Form I structure of an IIR filter can be converted to several alternate forms 60 Alternate IIR Filter Structures thinking of feedforward and feedback portions as two separate filter stages, because both stages are linear and time invariant, we can swap them with no change in y(n) output because sequence g(n) is shifted down along both delay lines in (b), we can eliminate one of delay paths and arrive at filter structure shown in (c), where only half the delay storage registers are required compared to Direct Form I structure 61 Alternate IIR Filter Structures  Transposition theorem  There is a process in DSP that allows us to change structure of an LTI digital network without changing network’s transfer function (its frequency response)  That network conversion process follows transposition theorem  A transposed version of some digital network might be easier to implement, or may exhibit more accurate processing, than the original network 62 Alternate IIR Filter Structures  Steps to transpose a digital filter (starting with Direct Form II)  1. reverse direction of all signal-flow arrows  2. convert all adders to signal nodes  3. convert all signal nodes to adders  4. swap x(n) input and y(n) output labels 63 Alternate IIR Filter Structures By convention, we flip the network in (b) from left to right so that x(n) input is on the left as shown in (c) 64 Alternate IIR Filter Structures  Fig. 6-23  Transposed filter contains the same number of delay elements, multipliers, and addition operations as the original filter, and both filters have the same transfer function given by  When implemented using infinite-precision arithmetic, Direct Form II and transposed Direct Form II filters have identical frequency response properties 21 21 )2()1(1 )2()1()0( )( )( )(      zaza zbzbb zX zY zH 65 Alternate IIR Filter Structures  Fig. 6-23  Transposed Direct Form II structure is less susceptible to errors that can occur when finiteprecision binary arithmetic is used to represent data values and filter coefficients within a filter implementation  Direct Form II filters implement (possibly high-gain) feedback pole computations before feedforward zeros computations  large intermediate data values which must be truncated  Transposed Direct Form II filters implement zeros computations first followed by pole computations  Direct Form I filter has the most resistance to coefficient quantization and stability problems 66 Impulse Invariance IIR Filter Design Method  Impulse invariance method  Is based upon the notion that we can design a discrete filter whose time-domain impulse response is a sampled version of impulse response of a continuous analog filter  If that analog filter (called prototype filter) has some desired frequency response, then our IIR filter will yield a discrete approximation of that desired response 67 Impulse Invariance IIR Filter Design Method 68 Impulse Invariance IIR Filter Design Method  Impulse invariance method  Our goal is to design a digital filter whose impulse response is a sampled version of analog filter’s continuous impulse response  We can map each pole on s-plane for analog filter’s Hc(s) transfer function to a pole on z-plane for discrete IIR filter’s H(z) transfer function  Impulse invariance method yields useful IIR filters as long as sampling rate is high relative to bandwidth of signal to be filtered  IIR filters designed using impulse invariance method are susceptible to aliasing problems because practical analog filters cannot be perfectly band-limited 69 Impulse Invariance IIR Filter Design Method we prefer to make fs as large as possible to minimize the overlap between primary frequency response curve and its replicated images spaced at multiples of ±fs Hz Due to aliasing behavior of impulse invariance design method, this filter design process should never be used to design highpass digital filters 70 Impulse Invariance IIR Filter Design Method  Two methods for designing IIR filters using impulse invariance  Method 1  Requires that an inverse Laplace transform as well as a z-transform be performed  Method 2  Uses a direct substitution process to avoid inverse Laplace and z-transformations at the expense of needing partial fraction expansion algebra necessary to handle polynomials 71 Impulse Invariance IIR Filter Design Method  Method 1  Step 1: Design a prototype analog filter with desired frequency response  In a lowpass filter design, for example, filter type (Chebyshev, Butterworth, elliptic), filter order (number of poles), and cutoff frequency are parameters to be defined in this step  Result of this step is a continuous Laplace transfer function Hc(s) expressed as ratio of two polynomials, such as  N < M, and a(k) and b(k) are constants           M k k N k k MM NN c ska skb asasMasMa bsbsNbsNb sH 0 0 1 1 )( )( )0()1(...)1()( )0()1(...)1()( )( 72 Impulse Invariance IIR Filter Design Method  Method 1  Step 2: Determine analog filter’s continuous timedomain impulse response hc(t) from Hc(s) Laplace transfer function  Can be done using Laplace tables as opposed to actually evaluating an inverse Laplace transform equation  Step 3: Determine digital filter’s fs, and calculate ts = 1/fs  fs is chosen based on absolute frequency, in Hz, of prototype analog filter  Because of aliasing problems associated with this impulse invariance design method, fs should be made as large as is practical 73 Impulse Invariance IIR Filter Design Method  Method 1  Step 4: Find z-transform of continuous hc(t) to obtain IIR filter’s z-domain transfer function H(z) in form of a ratio of polynomials in z  Step 5: Substitute the value ts for the continuous variable t in H(z) transfer function obtained in Step 4  In performing this step, we are ensuring that IIR filter’s discrete h(n) impulse response is a sampled version of continuous filter’s hc(t) impulse response so that h(n) = hc(nts), for 0 ≤ n ≤ ∞ 74 Impulse Invariance IIR Filter Design Method  Method 1  Step 6: H(z) from Step 5 will now be of the general form  Because process of sampling continuous impulse response results in a digital filter frequency response that’s scaled by a factor of 1/ts, we include ts factor in this equation  Incorporating ts makes IIR filter time-response scaling independent of sampling rate, and discrete filter will have the same gain as prototype analog filter              M k k N k k MM NN zka zkb azazMazMa bzbzNbzNb zH 1 0 1)1( 1)1( )(1 )( )0()1(...)1()( )0()1(...)1()( )(         M k k N k k s zka zkbt zX zY zH 1 0 )(1 )( )( )( )( 75 Impulse Invariance IIR Filter Design Method  Method 1  Step 7: By inspection, we can express filter’s time-domain difference equation as )()(...)2()2()1()1( )]()(...)2()2()1()1()()0([)( )(1 )( )( )()(...)2()2()1()1( )()(...)2()2()1()1()()0()( )(1 )( )( 1 0 1 0 MnyManyanya NnxNbnxbnxbnxbtny zka zkbt zH MnyManyanya NnxNbnxbnxbnxbny zka zkb zH s M k k N k k s M k k N k k                       these time-domain expressions apply only to the filter structure in Fig. 6-18. The a(k) and b(k), or ts · b(k), coefficients, however, can be applied to the improved IIR structure shown in Fig. 6-22 to complete our design 76 Impulse Invariance IIR Filter Design Method 77 Impulse Invariance IIR Filter Design Method  Method 2  Step 1: Obtain Laplace transfer function Hc(s) for prototype analog filter  Same as Method 1, Step 1  Step 2: Select an appropriate sampling frequency fs and calculate sample period as ts = 1/fs  Same as Method 1, Step 3           M k k N k k MM NN c ska skb asasMasMa bsbsNbsNb sH 0 0 1 1 )( )( )0()1(...)1()( )0()1(...)1()( )( 78 Impulse Invariance IIR Filter Design Method  Method 2  Step 3: Express analog filter’s Laplace transfer function Hc(s) as sum of single-pole filters  This requires us to use partial fraction expansion methods where M > N, individual Ak factors are constants, and the kth pole is located at −pk on s-plane  Hk(s) = kth single-pole analog filter M M M k k k MM NN c ps A ps A ps A ps A asasMasMa bsbsNbsNb sH               ... )0()1(...)1()( )0()1(...)1()( )( 2 2 1 1 1 1 1 k k k ps A sH  )( 79 Impulse Invariance IIR Filter Design Method  Method 2  Step 4: Substitute 1−e−pktsz−1 for s+pk in Hc(s) equation in Step 3  This mapping of each Hk(s) pole, located at s = −pk on s-plane, to z = e−pkts location on z-plane is how we approximate the impulse response of each single-pole analog filter by a single-pole digital filter  So, the kth analog single-pole filter Hk(s) is approximated by a single-pole digital filter whose zdomain transfer function is           M k M k tp k k tp k k ze A zHzH ze A zH sk sk 1 1 1 H(z)combinedfinal 1 1 )()( 1 )( 80 Impulse Invariance IIR Filter Design Method  Method 2  Step 5: Calculate z-domain transfer function of the sum of M single-pole digital filters in the form of a ratio of two polynomials in z  Because H(z) in Step 4 will be a series of fractions, we’ll have to combine those fractions over a common denominator to get a single ratio of polynomials in form of         M k k N k k zka zkb zX zY zH 1 0 )(1 )( )( )( )( 81 Impulse Invariance IIR Filter Design Method  Method 2  Step 6: As in Method 1, Step 6, by inspection, we express filter’s time-domain equation in general form of )()(...)2()2()1()1( )]()(...)2()2()1()1()()0([)( )(1 )( )( )()(...)2()2()1()1( )()(...)2()2()1()1()()0()( )(1 )( )( 1 0 1 0 MnyManyanya NnxNbnxbnxbnxbtny zka zkbt zH MnyManyanya NnxNbnxbnxbnxbny zka zkb zH s M k k N k k s M k k N k k                       Finally, we implement the improved IIR structure shown in Fig. 6-22 using a(k) and b(k) coefficients or a(k) and ts·b(k) coefficients from these time-domain expressions 82 Impulse Invariance IIR Filter Design Method  Impulse invariance design Method 1 example  Design an IIR filter that approximates a 2nd-order Chebyshev prototype analog lowpass filter whose passband ripple is 1 dB  fs = 100 Hz (ts = 0.01)  Filter’s 1 dB cutoff frequency = 20 Hz 83 Impulse Invariance IIR Filter Design Method  Method 1 example  Given above filter requirements, assume that analog prototype filter design effort results in  It’s this transfer function that we intend to approximate with our discrete IIR filter 145.1741094536.137 145.17410 )( 2   ss sHc 22222 22 2)( )sin( )( :)(:)(oftransformLaplace),(               ss A s A tAe s A txtxsX t 84 Impulse Invariance IIR Filter Design Method )485173.112sin(77724.154 )sin()( )485173.112()972680.68( )485173.112)(77724.154( )( 77724.154 145.17410 485173.112145.17410145.17410 972680.68 2 94536.137 2145.1741094536.137 145.17410 )( 972680.68 filteranalogprototypeof responseimpulse domain-time 22 222 2222 te tAeth s sH A ss A ss sH t t c c c                      85 Impulse Invariance IIR Filter Design Method 21 1 268972680.02168972680.0 168972680.0 201.0972680.682101.0972680.68 101.0972680.68 variable continuous for01.0 substitute 2972680.6821972680.68 1972680.68 972680.68 221 1 25171605.043278805.01 059517.70 )( )( )]12485173.1cos([21 )12485173.1sin(77724.154 )]01.0485173.112cos([21 )01.0485173.112sin(77724.154 )( )]485173.112cos([21 )485173.112sin(77724.154 )( )485173.112sin(77724.154)( )]cos([21 )sin( )sin( :)(oftransform-z),(:)(                              zz z zX zY zeze ze zeze ze zH zezte zte zH teth zezte ztCe tCe txzXtx t t tt t t c tt t t s       86 Impulse Invariance IIR Filter Design Method )2(25171605.0)1(43278805.0)1(70059517.0)( )2(25171605.0)1(43278805.0)1(059517.7001.0)( )(25171605.0)(43278805.0)(059517.70)( )059517.70()()25171605.043278805.01()( 25171605.043278805.01 059517.70 )( )( )( scalingfor010bytcoefficien)1(Multiplyfilter.IIRforexpressiondomain-Get time 211 121 21 1             nynynxny nynynxny zzYzzYzzXzY zzXzzzY zz z zX zY zH .=tn-x s these coefficients are what we use in implementing the improved IIR structure shown in Fig. 6-22 to approximate the original 2nd-order Chebyshev analog lowpass filter 87 Impulse Invariance IIR Filter Design Method  Impulse invariance design Method 2 example  Given original prototype filter’s Hc(s) as and ts = 0.01 145.1741094536.137 145.17410 )( 2   ss sHc )2/)(2/( )( 4 4 24 4 2 )( )( |4/)4|(and1 22 2 14517410 94536137 2 jRbsjRbs c sH cbb s cbb s c sH csbs c sH c cbRj c c .c = .b =                             If we substitute the values for b and c, the quantity under radical sign is negative  factors in denominator are complex 88 Impulse Invariance IIR Filter Design Method )2/( 2/ )2/( 2/ )( 222/2/ )2/( )2/)(2/( 222/2/ )2/( )2/)(2/( )2/()2/()2/)(2/( )( 2/ 2 2/ 1 21 jRbs Rjc jRbs Rjc sH R jc jR c jRbjRb c jRbs jRbsjRbs c K R jc jR c jRbjRb c jRbs jRbsjRbs c K jRbs K jRbs K jRbsjRbs c sH c jRbs jRbs c                                         89 Impulse Invariance IIR Filter Design Method 212/ 12/ 21)2/()2/( 1)2/()2/( 21)2/(1)2/( 1)2/(1)2/( 1)2/(1)2/( 1)2/(1)2/( 1)2/(1)2/( termsfor1substitute :plane-ztoplane-sfrom2and2polesmap )(1 )()2/( )( )(1 ))(2/( )( 1 )11)(2/( )( )1)(1( )1)(2/()1)(2/( )( 1 2/ 1 2/ )( )2/( 2/ )2/( 2/ )( 1 21                                     zezeee zeeeRjc zH zezee zeeRjc zH zezeze zezeRjc zH zeze zeRjczeRjc zH ze Rjc ze Rjc zH jRbs Rjc jRbs Rjc sH ssss sss sss ss sss ss ss ss ss k stkp btjRtjRtbt jRtjRtbt bttjRbtjRb tjRbtjRb bttjRbtjRb tjRbtjRb tjRbtjRb tjRbtjRb tjRbtjRb psze +jRb/=pjRb/=p c 90 Impulse Invariance IIR Filter Design Method 21 1 21 1 010and,48517112,94536137,14517410 212/ 12/ 212/ 12/ 2/)()(osand2/)()sin(:Euler 212/ 12/ 25171605.043278805.01 059517.70 25171605.0)86262058.0)(50171312.0(1 )902203655.0)(50171312.0)(77724.154( )( )]cos(2[1 )][sin()/( )]cos(2[1 )]sin(2[)2/( )( )(1 )()2/( )(                             zz z zz z zH zezRte zRteRc zezRte zRtjeRjc zH zezeee zeeeRjc zH .=t.R=.b=.c= bt s bt s bt bt s bt s bt eecjee btjRtjRtbt jRtjRtbt s ss s ss s jjjj ssss sss   91 Impulse Invariance IIR Filter Design Method )2(25171605.0)1(43278805.0)1(70059517.0)( )2(25171605.0)1(43278805.0)1(059517.7001.0)( )2(25171605.0)1(43278805.0)1(059517.70)( )(25171605.0)(43278805.0)(059517.70)( )059517.70()()25171605.043278805.01()( 25171605.043278805.01 059517.70 )( )( )( bytcoefficien)1(hemultiply t 211 121 21 1               nynynxny nynynxny nynynxny zzYzzYzzXzY zzXzzzY zz z zX zY zH stnx 92 Impulse Invariance IIR Filter Design Method )2/( 2/ )2/( 2/ )( jRbs Rjc jRbs Rjc sHc      93 Impulse Invariance IIR Filter Design Method  IIR filter’s z-plane pole locations   64.450.5017:poleconjugate 64.45-0.5017 radians )2/()2/( )1( )2/( )1( )2/( )( 1 2/ 1 2/ )( poleplane-zupper 2/2/)2/(poleplane-zlower )2/()2/( 1)2/(1)2/( 1)2/(1)2/(                         s btjRtbttjRb tjRbtjRb tjRbtjRb tjRbtjRb Rteeeez ez zRjc ez zRjc zez Rjcz zez Rjcz z z zH ze Rjc ze Rjc zH ssss ss ss ss 94 Impulse Invariance IIR Filter Design Method b(0) coefficient is zero here equivalent 95 Impulse Invariance IIR Filter Design Method  In general, Method 2 is more popular for two reasons  (1) inverse Laplace and z-transformations can be very difficult for higher-order filters  (2) unlike Method 1, Method 2 can be coded in a software routine or a computer spreadsheet 96 Impulse Invariance IIR Filter Design Method  Fig. 6-35(b)  Roll-off is not particularly steep  Attenuation slope is so gradual that it doesn’t appear to be of much use as a lowpass filter  Any frequency representation (be it a digital filter magnitude response or a signal spectrum) that has nonzero values at +fs/2, most probably has aliasing problem  We also see that passband ripple is greater than the desired value of 1 dB in Fig. 6-34  It’s not the low order of filter that contributes to its poor performance, but the sampling rate used 97 Impulse Invariance IIR Filter Design Method Increasing sampling rate to 400 Hz results in the much improved frequency response Sampling rate changes do not affect filter order or implementation structure; it only changes ts in our design equations, resulting in a different set of filter coefficients the smaller we make ts (larger fs), the better the resulting filter because the replicated spectral overlap indicated in Fig. 6-32(b) is reduced due to the larger fs impulse invariance IIR filter design techniques are most appropriate for narrowband filters, that is, lowpass filters whose cutoff frequencies are much smaller than the sampling rate 98 Bilinear Transform IIR Filter Design Method  Bilinear transform method  A popular analytical IIR filter design technique  Like impulse invariance method, this design technique approximates a prototype analog filter defined by continuous Laplace transfer function Hc(s) with a discrete filter whose transfer function is H(z) 99 Bilinear Transform IIR Filter Design Method  Bilinear transform method has great utility  It allows to substitute a function of z for s in Hc(s) to get H(z), eliminating the need for Laplace and z-transformations as well as any necessity for partial fraction expansion algebra  It maps the entire s-plane to z-plane, enabling us to completely avoid frequency-domain aliasing problems we had with impulse invariance design method  It induces a nonlinear distortion of H(z)’s frequency axis, relative to the original prototype analog filter’s frequency axis, that sharpens the final roll-off of digital lowpass filters 100 Bilinear Transform IIR Filter Design Method  Bilinear transform method  If transfer function of a prototype analog filter is Hc(s), we obtain discrete IIR filter z-domain transfer function H(z) by substituting the following for s in Hc(s) where ts is discrete filter’s sampling period (1/fs)              1 1 1 12 z z t s s 101 Bilinear Transform IIR Filter Design Method  Bilinear transform’s s- to z-plane mapping  Any pole on the left side of s-plane will map to the inside of unit circle in z-plane  If σ < 0, |z| will be less than 1  If σ > 0, |z| will be greater than 1  This confirms that when using bilinear transform, any pole located on the left side of s-plane (σ < 0) will map to a z-plane location inside unit circle 22 22 rdenominato numerator 1 1 )2/()2/1( )2/()2/1( Mag Mag 2/)2/1( 2/)2/1( 2/2/1 2/2/1 )2/(1 )2/(1 1 12 sas sas sas sas sas sasjs s s s tt tt z tjt tjt tjt tjt z st st z z z t s a                                  analog 102 Bilinear Transform IIR Filter Design Method  Bilinear transform’s s- to z-plane mapping  jωa axis of s-plane maps to perimeter of unit circle in z-plane 1 )2/()1( )2/()1( Mag Mag 2/1 2/1 2/)2/1( 2/)2/1( 22 22 rdenominato numerator 0 0              sa sa sa sa sas sas t t z tj tj z tjt tjt z         103 Bilinear Transform IIR Filter Design Method  Bilinear transform’s s- to z-plane mapping  Frequency mapping from jωa axis of s-plane to perimeter of unit circle in z-plane is not linear )2/tan( 2 )2/cos( )2/sin( 2 22 )]2/cos(2[ )]2/sin(2[2 )( )(2 1 12 1 12 2/ 2/ 2/ 2/ 2/)()cos(and 2/)()sin(:Euler 2/2/2/ 2/2/2/ 1circleuniton 1 1 d sd d j j s a d j d j s a ee jee jjj jjj s a j j s ez s t jj e e t js e je t js eee eee t js e e t s z z t s d d d djj jj ddd ddd d ddj                                                                  104 Bilinear Transform IIR Filter Design Method  Bilinear transform’s s- to z-plane mapping  Range of ωd is ±π, and dimensions of digital frequency ωd are radians/sample  Range of ωa is ±∞, and dimensions of analog frequency ωa are radians/second                      d sa d ad s a d s a t t t j js 2 tan2 )2/tan( 2 0 )2/tan( 2 1 105 Bilinear Transform IIR Filter Design Method            d sa d t 2 tan2 1 frequency warping due to bilinear transform no matter how large s-plane’s analog ωa becomes, z-plane’s ωd will never be greater than π radians/sample (fs/2 Hz) 106 Bilinear Transform IIR Filter Design Method bilinear transform maps s-plane’s entire jωa axis onto the unit circle in z-plane 107 Bilinear Transform IIR Filter Design Method  Frequency warping  ωd = π radians/sample corresponds to fd = fs/2 Hz  ωd = 2π(fd / fs) Hz )/tan( 2 tan 2 Hz )/(tan 2 tan2 2 )/(2 1 2 )/(2 1           sds a f ff d s a sas d f ff sa d fff f t fff f t aa sdd aa sdd                       108 Bilinear Transform IIR Filter Design Method Hz )/(tan 1   sas d fff f   fd frequency warping (compression) becomes more severe as fd approaches fs/2 if a bilinear-transformdesigned digital bandpass filter is desired to have an upper cutoff frequency of fd1 Hz, then the original prototype analog bandpass filter must be designed (prewarped) to have an upper cutoff frequency of fa1 Hz using: Hz )/tan(   sds a fff f  no IIR filter response aliasing can occur with bilinear transform design method 109 Bilinear Transform IIR Filter Design Method  Steps to perform an IIR filter design using bilinear transform method  Step 1: Obtain Laplace transfer function Hc(s) for the prototype analog filter  Step 2: Determine digital filter’s equivalent fs and establish ts = 1/fs  Step 3: In Laplace Hc(s) transfer function, substitute the expression for the variable s to get IIR filter’s H(z) transfer function             1 1 1 12 z z ts 110 Bilinear Transform IIR Filter Design Method  Step 4: Multiply numerator and denominator of H(z) by appropriate power of (1 + z−1) and collect terms of like powers of z in the form  Step 5: By inspection, express IIR filter’s timedomain equation in the general form of  Although this expression only applies to filter structure in Fig. 6-18, we can apply a(k) and b(k) coefficients to improved IIR structure shown in Fig. 6-22         M k k N k k zka zkb zH 1 0 )(1 )( )( )()(...)2()2()1()1( )()(...)2()2()1()1()()0()( MnyManyanya NnxNbnxbnxbnxbny   111 Bilinear Transform IIR Filter Design Method  Bilinear transform design example  Design an IIR filter that approximates 2nd-order Chebyshev prototype analog lowpass filter whose passband ripple is 1 dB  fs = 100 Hz (ts = 0.01)  Filter’s 1 dB cutoff frequency = 20 Hz  Original prototype filter’s Laplace transfer function is given as 145.1741094536.137 145.17410 )( 2   ss sHc 112 Bilinear Transform IIR Filter Design Method 22122 21 zofpowerslikecollecting 2111212 21 )1(byrdenominato andnumeratormultiply 1 1 2 1 1 2 /2 :simplifyto 1 1 2 1 1 1 12 2 14517410 94536137 2 )()22()( )21( )( )1()1)(1()1( )1( )( 1 1 1 1 )( 1 12 1 12 )( )( 145.1741094536.137 145.17410 )( 21 1 1                                                                                       zabcazaccaba zzc zH zczzabza zc zH c z z ab z z a c zH c z z t b z z t c zH csbs c sH ss sH z ta ss z z t s c .c = .b = c s s 113 Bilinear Transform IIR Filter Design Method )2(35083938.0)1(53153089.0 )2(20482712.0)1(40965424.0)(20482712.0)( 35083938.053153089.01 20482712.040965424.020482712.0 )( 35083938.053153089.01 )21(20482712.0 )( )( )( )( )22( 1 )21( )( )( filterIIRfor expression domain-time 21 21 21 21 14517410 94536137 2002 2 2 2 1 2 2 21 2 c)ab(abyrdenominatoandnumeratordivide r,denominatoinoneofermconstant tagetto 2                            nyny nxnxnxny zz zz zH zz zz zH z caba abca z caba ac zz caba c zH .c = .b = =/ta = s 114 Bilinear Transform IIR Filter Design Method bilinear-transformdesigned filter’s magnitude response approaches zero at folding frequency of fs/2 = 50 Hz 115 Bilinear Transform IIR Filter Design Method bilinear transform design method gives a much sharper roll-off for our lowpass filter for two reasons: 1) frequency warping of bilinear transform method compresses (sharpens) roll-off portion of a lowpass filter; 2) the price we pay in terms of additional complexity of implementation of our IIR filter: our new filter requires five multiplications per filter output sample where impulse invariance design filter in Fig. 6-36(a) required only three multiplications 116 Bilinear Transform IIR Filter Design Method  Prewarping  If cutoff frequency is a large percentage of fs, resultant |Hd(fd)| cutoff frequency will be below the desired value  To avoid this, we prewarp prototype analog filter’s cutoff frequency requirement before the analog Hc(s) transfer function is derived in Step 1  In that way, they compensate for the bilinear transform’s frequency warping before it happens  To determine prewarped analog filter lowpass cutoff frequency that we want mapped to the desired IIR lowpass cutoff frequency, use        2 tan 2 d s a t   we plug desired IIR cutoff frequency ωd in here to calculate ωa cutoff frequency used to derive prototype analog filter’s Hc(s) transfer function 117 Optimized IIR Filter Design Method  Optimization methods  Developed for situation when the desired IIR filter frequency response is not of standard lowpass, bandpass, or highpass form  Closed-form expressions for filter’s z-transform do not exist  no explicit equations to work with  Designer should describe, in some way, the desired IIR filter frequency response  The algorithm, then, assumes a filter transfer function H(z) as a ratio of polynomials in z and starts to calculate filter’s frequency response  Based on some error criteria, the algorithm iteratively adjusts filter’s coefficients to minimize the error between desired and actual filter frequency response 118 Optimized IIR Filter Design Method  Optimized IIR filter design routines  Are used to design the simpler lowpass, bandpass, or highpass forms even though analytical techniques exist  They only require the designer to specify a few key amplitude and frequency values, and the desired order of IIR filter (the number of feedback taps), and the software computes the final feedforward and feedback coefficients 119 Optimized IIR Filter Design Method In specifying a lowpass IIR filter, a software design routine might require us to specify the values for δp, δs, f1, and f2 120 Comparison of IIR and FIR Filters Characteristic IIR FIR (nonrecursive) Number of necessary multiplications Least Most Sensitivity to filter coefficient quantization Can be high Very low Probability of overflow errors Can be high Very low Stability Must be designed in Guaranteed Linear phase No Guaranteed Can simulate prototype analog filters Yes No Required coefficient memory Least Most Hardware filter control complexity Moderate Simple Availability of design software Good Very good Ease of design, or complexity of design software Moderately complicated Simple Difficulty of quantization noise analysis Most complicated Least complicated Supports adaptive filtering With difficulty Yes So long as FIR coefficients are symmetrical (or antisymmetrical)