Z TRANSFORMS DR. VETON KEPUSKA Z-TRANSFORM In mathematics and signal processing, the Z-transform converts a discrete-time signal, which is a sequence of real or complex numbers, into a complex frequency domain representation. It can be considered as a discrete-time equivalent of the Laplace transform. This similarity is explored in the theory of time scale calculus. The z-transform is useful for the manipulation of discrete data sequences and has acquired a new significance in the formulation and analysis of discrete-time systems. It is used extensively today in the areas of applied mathematics, digital signal processing, control theory, population science, economics. These discrete models are solved with difference equations in a manner that is analogous to solving continuous models with differential equations. The role played by the ztransform in the solution of difference equations corresponds to that played by the Laplace transforms in the solution of differential equations. 02/27/2020
2 Z-TRANSFORM Definition: Consider a function x(t) defined for t0. As we should know we could sample this continuous function at times t = T, 2T, 3T, where T is the sampling period (or sampling rate). We can write the sample as a sequence using the notation Given a finite length sequence x[n], defended in the interval [0,N] and z-any complex number, the z-transform is defined as: 02/27/2020 3 Z-TRANSFORM
This transformation produces a new representation of x[n] denoted by X(z). Returning to the original sequence (via inverse z-transform) x[n] requires finding the coefficients associated with nth power of z-1. 02/27/2020 4 Z-TRANSFORM Formally transforming from the time/sequence/n-domain to z domain is represented as: 02/27/2020 5
Z-TRANSFORM A sequence and its z-transform are said to form a z-transform pair: n is operator in the domain of sequence indicating the sample is is considered an independent variable. In the z domain the independent variable is z. 02/27/2020 6 EXAMPLE 1 Appling the definition of the z transform: 02/27/2020
7 EXAMPLE 2 Appling the definition of the z transform: 02/27/2020 8 EXAMPLE 3 Appling the definition of the z transform: 02/27/2020 9
THE Z TRANSFORM AND LINEAR SYSTEMS THE Z TRANSFORM AND LINEAR SYSTEMS The z-transform is specifically useful in the analysis and design of LTI systems 02/27/2020 11 FIR FILTER 1. A system that has only for k=1,,N is said to be a Finite Impulse Response (FIR) filter.
The name reflects the fact that FIR filters have finite impulse (e.g., unit sample) response. FIR filters are also called moving average (MA) filters considering the fact their output is simply a weighted average of the input values. 02/27/2020 12 THE Z-TRANSFORM OF AN FIR FILTER We should recall that for any LTI system with input x[n] and impulse response h[n], the output y[n] is:
We are interested in the z-transform of h[n], where h[n] is an FIR filter 02/27/2020 13 Consider the input: The output y[n] is: 02/27/2020 14 The term under the parenthesis is the z-transform of h[n] This is also termed the system function This function is defined as:
02/27/2020 15 The z-transform pair that was just established is 02/27/2020 16 System Function is an Mth degree polynomial in complex variable z. As with any polynomial, it will have M roots or zeros, that is there are M values, z0, such that H(z0) = 0;
These M zeros completely define the polynomial, that is: Where zk, k=1,,M denote zeros of polynomial. 02/27/2020 17 EXAMPLE 4 Find the zeros of: The z-transform is 02/27/2020 18 The zeros of H(z) are -1/2 and +1/3
The difference equation has the same zeros, but a different scale factor; 02/27/2020 19 PROPERTIES OF THE ZTRANSFORM PROPERTIES OF THE Z-TRANSFORM The z-Transform has a few very useful properties, and is defintio extends to infinite signals/impulse responses (IIR). 02/27/2020 21 THE SUPERPOSITION (LINEARITY)
PROPERTY The Superposition (Linearity) Property 02/27/2020 22 Proof: 02/27/2020 23 THE TIME-DELAY PROPERTY Proof:
02/27/2020 24 THE TIME-DELAY PROPERTY Let Hence 02/27/2020 25 Similarly 02/27/2020 26
A GENERAL Z-TRANSFORM FORMULA We have seen that for a sequence x[n] defined over the interval 0 n N the z-transform is This definition extends for sequences having interval from - n 02/27/2020 27 PROPERTIES OF Z-TRANSFORM z-Transform Convolution Theorem: Cascading Systems x[ n] X( z)
LT 1 H1(z), h1[n] w[n ] W(z ) LT 2 H2(z), h2[n] y[ n] Y( z) 02/27/2020
28 THE Z-TRANSFORM AS AN OPERATOR THE Z-TRANSFORM AS AN OPERATOR The z-transform can be considered as an operator. Unit-Delay Operator x[ n] x[ n] Unit Delay
z-1 y[n]= x[n1] y[n]= x[n1] 02/27/2020 30 UNIT DELAY OPERATOR In the case of the unit delay, we observe that: Which is derived from the fact that Y(z) = z-1X(z) 02/27/2020 31 THE Z-TRANSFORM AS AN OPERATOR
The filter described with the following equation: This expression can be viewed as the operator: This is so because: 02/27/2020 32 EXAMPLE: TWO-TAP FIR FILTER x[ n] X(z) b0 z-1 z1
X(z) b1 x[n1] y[n ] 02/27/2020 33 EXAMPLE: TWO-TAP FIR FILTER Using the operator convention, we can write by inspection: 02/27/2020 34
FIR FILTERS FIR FILTERS General constant-coefficients equation: M y n k x n k k 0 Impulse (unit sample) response, h[n], of the filter is obtained when x[n]=[n]. y n 0 n 1 n 1 M n M
1, n n0 0 n n0 n n0 0 otherwise y 0 0 y1 1 y M M y k 0, k M 02/27/2020 36
LINEAR PHASE OF FIR FILTERS If FIR filter has coefficients that are symmetric, as depicted in the following relationship: 0 M 1 M 1 2 M 2 Then, it can be shown that the resulting filter has linear phase constant delay for all frequency components of the input signal. This property is very important in many communications data streams (speech, data, ) and image processing applications. 02/27/2020 37 FIR FILTER STRUCTURES
For a given difference equation there are different ways to implement a digital filter. Selection of a particular filter structure to be implemented is dependent on many factors: Programming considerations Hardware Sensitivity of Quantizing Coefficients Quantization noise of the input signal. 02/27/2020 38
DIRECT STRUCTURE OF AN FIR FILTER M y n k x n k k 0 x[n] 0 D D x[n-1] 1 D
D x[n-2] 2 D D x[n-M] M y[n] 02/27/2020 39 SYMMETRIC LINEAR PHASE FIR FILTER
0 Better Implementation: y[n] x[n] D D D D D D
D D 1 x[n-1] x[n-2] x[n-M/2+1] D D L/2 x[n-M/2]
02/27/2020 40 MATLAB >> filterDesigner 02/27/2020 41 FILTER DESIGNER Equiripple Least-Squares Window Constrained Least-Squares Complex Equiripple
Maximally Flat Least Pth Norm Constrained Equiripple Generalized Equiripple Const. Band Equiripple Interpolated FIR 02/27/2020 42 WINDOWING METHOD Choose Chebyshev Filter 02/27/2020 43 CHEBYSHEV WINDOW METHOD
02/27/2020 44 CONVOLUTION AND THE Z-TRANSFORM CONVOLUTION AND THE Z-TRANSFORM The impulse response of the unity delay system is: The system output written in terms of a convolution is: The system function (z-transform of h[n]) is: 02/27/2020 46
CONVOLUTION AND THE Z-TRANSFORM Hence in general, when the system is described by H(z) we can optain the output Y(z) knowing what the input is X(z) To show this we start with: 02/27/2020 47 CONVOLUTION AND THE Z-TRANSFORM
( ) ( )= h [ ] ( ( ) = h [ ] ( )= ( ) () =0 =0 02/27/2020 48 EXAMPLE: CONVOLVING A FINITE DURATION SEQUENCES
Suppose that: Find first Y(z) by applying z transform to both x[n] and h[n] Then we find the Y(z) by applying the multiplication of z 02/27/2020 transform rule: 49 EXAMPLE: CONVOLVING A FINITE DURATION SEQUENCES Then we find y[n] by inverse z transformation: 02/27/2020
50 FACTORIZATION USING MATLAB H(z) is a 3rd degree polynomial. We can use the MALTAB function roots() to obtain the roots. >> p = [1 3 -2 1]; >> r = roots(p) r= -3.6274 + 0.0000i 0.3137 + 0.4211i 0.3137 - 0.4211i 02/27/2020 51 >> conv([1 -r(2)], [1 -r(3)])
ans = 1.0000 -0.6274 0.2757 02/27/2020 52 FACTORING Cascade of the system is thus: x[n ] X(z )
1+3.6274z-1 H1(z ) w[n] W(z) 1 10.6274 +0.2757 H2(z ) 2 02/27/2020 y[n]
Y(z) 53 >> conv([1 -r(1)],conv([1 -r(2)],[1 -r(3)])) ans = 1.0000 3.0000 -2.0000 1.0000 The difference equations for each subsystem are: 02/27/2020 54 DECONVOLUTION/INVERSE FILTERING
In a two subsystems cascade can the second system undo the action of the first subsystem? For the output to equal the input we need H(z) = 1 We thus desire 02/27/2020 55 EXAMPLE: The inverse filter is: This expression is no longer an FIR filter, it is an infinite impulse response (IIR) filter. We can approximate
as n FIR filter via long devision: 02/27/2020 56 02/27/2020 57 An M+1 term approximation is: 02/27/2020 58 RELATIONSHIP BETWEEN THE Z-DOMAIN
AND THE FREQUENCY DOMAIN ( )= ( )= =0 =0
02/27/2020 59 THE Z-PLANE AND THE UNIT CIRCLE Considering z-plane by evaluating H(z) which becomes H(). 02/27/2020 60 From this interpretation we also can see why H() is periodic with period 2: As increases it continues to sweep around the unit circle over and over again 02/27/2020
61 THE ZEROS AND POLES OF H(Z) Conisder: Factoring H(z) results in The zeros are the locations where H(z) = 0; i.e., , , 02/27/2020 62 POLE-ZERO PLOT The poles are where H(z) I m ZPlane
z1 z2 Triple Pole 3 R e z3 02/27/2020 63 EXAMPLE . MATLAB has a function that supports the creation of a polezero plot given the system function coefficients
>> zplane([1 2 2 1], 1) 02/27/2020 64 1 0.8 0.6 0.4 0.2 3
0 -0.2 -0.4 -0.6 -0.8 -1 02/27/2020 -1 -0.5
0 Real Part 0.5 1 65 The Significance of the Zeros of H(z) The difference equation specifies the actual time domain of the filter. It can be used for calculating the filter output for a given filter input The difference equation coefficients are the polynomial coefficients in H(z) For . We know that
so in particular if is one of the zeros of , If a zero lies on the unit circle then the output will be zero for a sinusoidal input of the form where 02/27/2020 is the angle of the zero relative to the real axis, which is also the frequency of the corresponding complex sinusoid; why? 66 EXAMPLE: Solution: Factoring H(z): Expanding x[n] we see that:
02/27/2020 67 The action of H(z) at will remove the signal from the filter output. Setting up a simple MATLAB script will validate this: >> n = 0:100; >> w0 = pi/4; >> x = cos(w0*n); >> y = filter([1 -2*cos(w0) 1],1,x); >> stem(n,x,'filled') >> hold Current plot held >> stem(n,y,'filled','r') >> axis([0 50 -1.1 1.1]); grid 02/27/2020
68 Input (blue ) 1 0.8 0= 4 0.6
0.4 0.2 0 Outpu t (red) -0.2 -0.4 -0.6 -0.8
-1 02/27/2020 0 5 10 15 20 25 30 35
40 45 69 50 Since the input is applied at n=0, we see a small transient while the filter settles to the final output, which in this case is zero. >> zplane([1 -2*cos(w0) 1],1)% check the pole-zero plot 02/27/2020
70 1 0.8 0.6 0.4 0= 0.2 4 2
0 -0.2 -0.4 -0.6 -0.8 -1 02/27/2020 -1 -0.5
0 Real Part 0.5 1 71 Graphical Relation Between z and When we make the substitution z=ej in H(z) we know that we are evaluating the z-transform on the unit circle and thus obtain the frequency response If we plot say |H(z)| over the entire z-pane we could visualize how we could get the frequency response magnitude based on evaluation of H(z) on unit circle.
02/27/2020 72 EXAMPLE L=9, Moving Average Filter (9 taps/8th-order) 02/27/2020 73 >> ZPLANE([ONES(1,9)]/9,1) 11 0.8 0.8 0.6 0.6
0.4 0.4 0.2 0.2 88 00 -0.2 -0.2 -0.4 -0.4 -0.6 -0.6 -0.8 -0.8 -1 -1 -1
-1 -0.5 -0.5 00 0.5 0.5 11 Real Real Part Part 02/27/2020
74 >> W = -PI:(PI/500):PI; >> H = FREQZ([ONES(1,9)/9],1,W); Magnitude Plot 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 100
200 300 400 500 600 700 800 900 1000
w 02/27/2020 75 FAZE RESPONSE Phaze 2.5 2 1.5 1 0.5 0 -0.5 -1 -1.5 -2
-2.5 100 200 300 400 500 600 700 800 900
1000 w 02/27/2020 76 USEFUL FILTERS THE L-POINT MOVING AVERAGE FILTER The L-point moving average (running sum) filter has the following definition: and system function (z-transform of the impulse response): 02/27/2020
78 THE L-POINT MOVING AVERAGE FILTER The sum can be simplified using the geometric series sum formula: Note that zeros of H(z) are determined by the roots of the equation: A geometric series is a series for which the ratio of each two consecutive terms / is a constant function of the summation index k. http://mathworld.wolfram.com/GeometricSeries.html VETON KEPUSKA 02/27/2020 79 THE L-POINT MOVING AVERAGE FILTER The roots of this equation can be found by noting that ej2 = 1
for any integer k. Those roots are referred to as the L roots of unity. 02/27/2020 80 THE L-POINT MOVING AVERAGE FILTER One of the zeros sits at z=1, but there is also a pole at z=1, so there is a pole-zero cancellation, meaning that the polezero plot of H(z) corresponds to the L-roots of unity, less the root at z = 0; 02/27/2020 81 1
Z-Plane 0.8 2 0.6 0.4 0.2 8 0 -0.2
-0.4 -0.6 -0.8 L=8 -1 -2.5 -2 -1.5 -1
-0.5 0 Real Part 0.5 1 1.5 2 02/27/2020 2.5
82 Magnitude Plot 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 100 200
300 400 500 w 0 600 700 2 4
800 900 1000 02/27/2020 83 PRACTICAL FILTER DESIGN USING MATLAB >> fdatool, that will be replaced in future release with >> filterDesigner Requires signal processing toolbox Will be using it to design FIT filters
02/27/2020 84 PROPERTIES OF LINEAR-PHASE FILTERS FIR filter can be designed to have a linear-phase It requires that the FIR filters to have symmetrical coefficients: 02/27/2020 85 EXAMPLE Symmetrical FIR filter Moving to frequency by
02/27/2020 86 EXAMPLE Since we have M = 4 the linear phase term is Real function R() is of the form: 02/27/2020 87 LOCATIONS OF THE ZEROS OF FIR LINEARPHASE SYSTEMS Further investigation of for the case of symmetric coefficients reveals that: Consequently, for
having a zero at it will also have a zero at 02/27/2020 88 LOCATIONS OF THE ZEROS OF FIR LINEARPHASE SYSTEMS Assuming the filter has real coefficients, complex zeros occur in conjugate pairs, so the even symmetry condition further implies that the zeros occur as quadruplets: 02/27/2020 89 zPlane
1 0 0 0 1 0 02/27/2020 90 Quadruplet Zeros for Linear Phase 1
2 3 ( )=12 +4 2 + 4 >> zplane([1 -2 4 -2 1],1) 1.5 1.5 1 1 0.5 0.5 4 0
4 0 -0.5 -0.5 -1 -1 02/27/2020 -1.5 -1.5 -2 -2 -1.5 -1.5
-1 -1 -0.5 0 0.5 -0.5 0 0.5 Real Part 1 1 1.5 1.5 2
2 91 IIR Filters The General IIR Difference Equation After rearranging this expression: This expression has N+M+1 coefficients, and that is number of multiplies and adds operation are required to produce each output. 02/27/2020 93 Direct-Form I of IIR Filter
y[n] 0x[n] x[n] D D D D 0 1x[n1] 1 2x[n-2] 1 2y[n-2]
2 D D Mx[n-M] M D D 1y[n-1] D D 2
D D Ny[n-N] N 02/27/2020 94 Time-Domain Response Consider a first-order filter (N=1) with M=0. 02/27/2020 95
Impulse Response of a First-Order IIR System The impulse response is obtained by setting x[n]=[n]. Initial Rest Conditions of IIR Filter: 1) The input is zero prior to the start time n0, that is x[n] = 0; for n < n0 2) The output is zero prior to the start time, that is y[n] = 0 for n < n 0 Find the impulse response via direct recursion from the difference equation 02/27/2020 96
0 0 0 0 1-st order FIR filter response to the impulse is: 02/27/2020 97 EXAMPLE: First-Order IIR with
The impulse response of: is Or 02/27/2020 98 First-Order with beta=1, alpha=0.8 1 0.9 0.8 0.7
0.6 0.5 0.4 0.3 0.2 0.1 0 5 10
15 20 Time 25 02/27/2020 30 99 Linearity and Time Invariance of IIR Filters Linearity of FIR filters and the conditions under which it is exhibited was covered in our previous discussions.
Using linearity and time invariance we can find the output of the first-order system to a linear combination of time shifter impulses 02/27/2020 100 2 [ ] = [ ] h [ ] = 1 02/27/2020
101 Example Using result presented in previous slides, it follows that: Plotting the y[n] results in the following graph: 02/27/2020 102 2 1.5 1
0.5 0 -0.5 5 10 15 20 Time 25
30 02/27/2020 103 Linearity & Time-Invariance Linearity and Time-Invariance are used to find the impulse response of this IIR filter: Superposition of un-delayed and delayed input: 02/27/2020 104 Based on this observation, the impulse response is:
02/27/2020 105 Step Response of a Frist-Order Recursive System The step response allows us to see how a filter (system) responds to an infinitely long input Considering the step response of the system: 02/27/2020 106 The for in the previous slide is a geometric series, which has a solution:
Using this result and assuming that the step response of the fist-order filter is: 02/27/2020 107 Three conditions for exist: 1. When the term grows without bounds as n becomes large, resulting in an unstable condition. 2. When the term decays to zero as n becomes large, resulting in an stable condition.
3. When , we have the special case output where it is of the form which also grows without bound; with the output alienates in sign, hence this case defines so called a marginally stable condition 02/27/2020 108 Example: The step response of the this filter: The parameters of the filter are: Hence: 02/27/2020
109 2 1.5 1 0.5 0 5 10 15
20 Time 25 02/27/2020 30 110 Direct Evaluation The step response can also be obtained by direct evaluationof the convolution sum: For the problem at hand:
02/27/2020 111 To evaluate this requires careful attention to details The product tells us how to set the sum limits u[nk] n>0 n<0 k 0 u[k] 0
k 02/27/2020 112 The result is: 02/27/2020 113 SYSTEM FUNCTION OF AN IIR FILTER IIR Filter System Function From our study of the z-transform we know that convolution in the time (sequence)-domain corresponds to multiplication in the z-domain
For the case of IIR filters H(z) will be a fully rational function, meaning in general both poles and zeros (more than at z=0) Applying z-transform on both sides of the general IIR difference equation: 02/27/2020 115 Forming the ration Y(z)/X(z) gives us transfer function H(z) The coefficients of the numerator polynomial, denoted B(z), correspond to the feed-forward terms of the difference equation The coefficients of the denominator polynomial, denoted A(z), for z-l, l>0 correspond to the feedback terms of the difference
equation 02/27/2020 116 MATLAB We have used various MATLAB functions that take as input b and a coefficient vectors, e.g., filter(b,a,...), freqz(b,a,...), and zplane(b,a) In terms of the general IIR system we now identify those vectors as: 02/27/2020 117
The General First-Order Case As a special case consider N=M=1: For a1=0.5, b0=-3, and b1=2: 02/27/2020 118 >> n = 0:20; >> x = [1 zeros(1,20)]; % impulse sequence input >> y = filter([-3,2],[1 -0.5],x); >> stem(n,y,'filled') >> axis([0 10 -3.1 .6]) >> grid >> ylabel('Impulse Response h[n]') >> xlabel('Time Index (n)') 02/27/2020
119 Impulse Response 0.5 0.5 00 -0.5 -0.5 -1 -1 -1.5 -1.5 -2 -2 -2.5 -2.5 -3
-3 00 11 22 33 44 55 66 77 88
99 10 10 Time Time Index Index (n) (n) 02/27/2020 120 Example: y = filter([1 1],[1 -0.8],x)
Find the System function, Impulse response, Difference equation That corresponds to the given filter() expression. By inspection: 02/27/2020 121 By inspection: The impulse response using is obtained by applying inverse
z-transform The difference equation is: 02/27/2020 122 SYSTEM FUNCTIONS AND BLOCK-DIAGRAM STRUCTURES THE DIRECT-FORM I 0x[n] x[n] D
D D D y[n] 0 1x[n1] 1 2x[n-2] 1 2y[n-2] 2 D D
Mx[n-M] D D 1y[n-1] D D 2 D D Ny[n-N]
02/27/2020 M N 124 Direct-Form II w[n] x[n] 0 y[n] D
D 1d[n-1] 1d[n-1] D D 2d[n-2] 2d[n-2] D D Md[n-M] Md[n-M] D
D Nd[n-N] 02/27/2020 125 THE TRANSPOSED STRUCTURES A property of filter block diagrams: When all the arrows are reversed All branch points become summing nodes All summing nodes become branch points
The input and output are interchanged. The system function is unchanged: 02/27/2020 126 Transposed DirectForm II x[n ] b0 b1 b2
bm-1 y[n ] z-1 z-1 z-1 z-1 a1 a2 an-1
02/27/2020 bm an 127 Relation to the Impulse Response Consider the input-output relationship: Impulse response of which is: Z transform of h[n] is: 02/27/2020 128
The previous expression summation defines a geometric series: Applying the sum formula to the previously obtained result with give us: The condition, specifies the range of values for which this transformation exists. The z-pane region for which transformation exits specifies the region of convergence. 02/27/2020 129 Z-Domain Time-Domain Relationship
1 1 1 [] Utilizing this result we can find the z-transform of: Using linearity and delay property of z transform we get: 02/27/2020 130 z-Transform 0+1 1
1 1 1 0 1 +1 = 1 1 11 11 11 02/27/2020 131 Poles and Zeros This transfer function has: The single pole , and
Singe zero 02/27/2020 132 Pole-Zero Plot First Order Filter 1 0.8 0.6 0.4 0.2
0 -0.2 -0.4 -0.6 -0.8 -1 -1 -0.5 0 0.5
1 Real Part 02/27/2020 133 POLES OR ZEROS AT THE ORIGIN OR INFINITY For the general IIR filter/system the number of poles always equals the number of zeros For FIR systems we saw that all of the poles were at z=0. It is also possible to have poles or zeros at z=. 02/27/2020
134 Example: Zero at z= Consider: This system has a pole at z = 0.8 and a zero at z = . 02/27/2020 135 Example: Zero at z= Consider: This system has a pole at z = , and a zero at z = -0.5. 02/27/2020
136 POLE LOCATIONS AND STABILITY We know that these expressions give a time z transform pair. Note that the system has a pole at z = a and a zero at z = 0; The impulse response decays to zero so long as |a|<1,which is equivalent to requiring that the pole lies inside the unit circle. System Stability: Causal LTI IIR systems, initially at rest, are stable if all of the poles of the system function lie inside the 02/27/2020 unit circle 137
EXAMPLE Converting to positive powers of z: 02/27/2020 138 EXAMPLE: SECOND-ORDER H(Z) Suppose that H(z) is 02/27/2020 139 In polar form the poles are hence the poles are inside the
unit circle and th system is stable. Stability can be checked using MATLAB tool zplane() to plo tht epolse and zeros: >> zplane([1 0.2],[1 -1.4 0.81]) 1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 02/27/2020
-1 -1 -0.5 0 Real Part 0.5 1 140 FREQUENCY RESPONSE OF AN IIR FILTER The frequency response is found by letting
in the system function (provided the system is stable) 02/27/2020 141 EXAMPLE Making the substitution we get: Use MATLAB freqz() to plot the magnitude and phase 02/27/2020 142
9 8 7 6 5 4 3 2 1 -1 Passband -0.8 -0.6 -0.4
-0.2 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6
0.8 1 0.2 0.4 0.6 0.8 1 1.5 1 0.5
0 -0.5 -1 -1.5 -1 0 Normalized Frequency ( rad/sample) 02/27/2020 143 This particular filter is a bandpass filter because it has a
relative large magnitude response over a narrow band of frequencies and small response otherwise From the earlier pole-zero analysis, the peak gain is near the angle the poles make to the real axis, 02/27/2020 144 EXAMPLE The frequency response is obtained by first substituting Plotting using MATLAB freqz() function we get: 02/27/2020
145 >> w = -pi:(pi/500):pi; >> H = freqz(1,[1 -0.8],w); 5 4.5 4 3.5 3 2.5 2 1.5 1 -1 Passband -0.8
-0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1
0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 02/27/2020 -0.8 -1 -0.8 -0.6
-0.4 -0.2 0 Normalized Frequency ( 0.2 rad/sample) 0.4 0.6 0.8
146 1 EXAMPLE As we did in previous example we substitute Plotting using MATLAB freqz() function we get: 02/27/2020 147 >> w = -pi:(pi/500):pi; >> H = freqz(1,[1 0.8],w); 5 4.5
4 3.5 3 2.5 2 1.5 1 -1 Hipass Filter -0.8 -0.6 -0.4 -0.2
0 0.2 0.4 0.6 0.8 1 0.8 0.6 0.4 0.2 0
-0.2 -0.4 -0.6 -0.8 -1 02/27/2020 Hipass Filter Phase -0.8 -0.6 -0.4 -0.2
0 0.2 0.4 0.6 0.8 148 1 THE INVERSE Z-TRANSFORM AND APPLICATIONS 02/27/2020
149 THE INVERSE Z-TRANSFORM AND APPLICATIONS Finding the impulse response of a first-order IIR system was not too difficult using difference equation recursion, but for N>1 system order, this process becomes too difficult We need an inverse z-transform approach that will allow us to work from any rational system function, H(z), backwards to h[n] Useful z-transform properties and pairs available to help this cause are listed below 02/27/2020 150 Z-TRANSFORM RELATIONS
x[n] X(z) 1. 2. 3. 4. 1 5. 6. 02/27/2020 151
A GENERAL PROCEDURE FOR INVERSE Z-TRANSFORMATION The technique we develop here uses an algebraic decomposition known as partial fraction expansion The function we wish to inverse transform, H(z), is assumed for now to be proper rational, meaning that M
Step 3: The inverse z transform is: 02/27/2020 153 The limitation of this approach is that the are distinct In general there may be repeated poles, in which case the partial fraction expansion takes a slightly different form from step 2 Hence, at present we will only consider non-repeated poles 02/27/2020 154
EXAMPLE M=1, N=2 Finding the roots of the H(z): 02/27/2020 155 Solving for Solving for 02/27/2020 156 Hence the overall solution is
Inverse z-transform is: MALTAB single processing toolbox has a function that can perform partial fraction expansion: 02/27/2020 157 >> help esidues --- help for residue --- residue Partial-fraction expansion (residues). [R,P,K] = residue(B,A) finds the residues, poles and direct term of a partial fraction expansion of the ratio of two polynomials B(s)/A(s). If there are no multiple roots, B(s)
R(1) R(2) R(n) ---- = -------- + -------- + ... + -------- + K(s) A(s) s - P(1) s - P(2) s - P(n) Vectors B and A specify the coefficients of the numerator and denominator polynomials in descending powers of s. The residues are returned in the column vector R, the pole locations in column vector P, and the direct terms in row vector K. The number of poles is n = length(A)-1 = length(R) = length(P). The direct term
coefficient vector is empty if length(B) < length(A), otherwise length(K) = length(B)-length(A)+1. 02/27/2020 158 If P(j) = ... = P(j+m-1) is a pole of multplicity m, then the expansion includes terms of the form R(j) R(j+1) R(j+m-1) -------- + ------------ + ... + -----------s - P(j) (s - P(j))^2 (s - P(j))^m
[B,A] = residue(R,P,K), with 3 input arguments and 2 output arguments, converts the partial fraction expansion back to the polynomials with coefficients in B and A. Warning: Numerically, the partial fraction expansion of a ratio of polynomials represents an ill-posed problem. If the denominator polynomial, A(s), is near a polynomial with multiple roots, then small changes in the data, including roundoff errors, can make arbitrarily large changes in the resulting poles and residues. 02/27/2020 159 Problem formulations making use of state-space or zero-pole representations are preferable. Class support for inputs B,A,R: float: double, single
See also poly, roots, deconv. >> 02/27/2020 160 Using residuez() we find: >> [A,p,K] = residuez([1 2],[1 -3/4 1/8]) A= % The partial fraction coefficients agree 10 %
-9 % p= % The pole factoring agree 0.5000 % 0.2500 % K= here)
[] % Results from long division to make proper ration (NA % 02/27/2020 161 As a further check we can plot h[n] directly and compare it to the results obtained by direct evaluation of the difference equation via filter() 3 2.5 2 1.5
1 0.5 0 0 2 4 6 8 10
12 14 16 18 20 3 2.5 2 1.5
1 0.5 0 02/27/2020 0 2 4 6 8
10 12 14 16 18 20 162 EXAMPLE Find
for IIR system having input and system function: Find Y(z): 02/27/2020 163 Use a partial fraction expansion over the three real poles: 1,1/3 and : 02/27/2020 164 Solving for coefficients:
02/27/2020 165 Combining the results: And using inverse transform: Checking this derivation using MATLAB residuez() function 02/27/2020 166 >> [A, p, K] = residuez([2 2], conv([1 -1], [1 -1/6 -1/6])) A= 6.0000 -3.6000
-0.4000 p= 1.0000 0.5000 -0.3333 K= [] >> The results agree. 02/27/2020 167 The partial fraction expansion technique requires that M
Performing long devision: 02/27/2020 168 Now we have reduced Y(z) in to the form: And 02/27/2020 169 Finally:
And time domain solution is: Checking this result with MATLABs resudiex(); function which automaticaly perform long division 02/27/2020 170 >> [A,p,K] = residuez([2 -2.4 -0.4],[1 -0.3 -0.4]) A= -1.0000 2.0000 p= 0.8000 -0.5000
K= 1 02/27/2020 171 STEADY-STATE RESPONSE AND STABILITY 02/27/2020 172 STEADY-STATE RESPONSE AND STABILITY It can be shown that for: And H(z) an IIR system with N > M, the output will be
Where are the poles of H(z) and are corresponding partial fraction coefficients. 02/27/2020 173 The first term represents the transient response, which provided all the poles of H(z) lie inside the unit circle, willdecay to zero, leaving just the sinusoidal term. For the output to reach sinusoidal steady-state we must have |pk| < 1 to insure that the transient term (first term) decays to zero, this then insures that the system is stable 02/27/2020 174
EXAMPLE Input is a cosine with starting a n=0 to the system: This is a second-order IIR filter, with the transient term that consists of two exponentials Second-order system will be studied in more detail in the next section. 02/27/2020 175 MATLAB
02/27/2020 176 >> n = -5:50; >> x = cos(2*pi/20*n).*ustep(n,0); >> y = filter(5,[1 -1.5 0.8],x); >> subplot(211) >> stem(n,x,'filled') >> axis([-5 50 -1 1]); grid >> ylabel('x[n]') >> subplot(212) >> stem(n,y,'filled') >> axis([-5 50 -30 30]); grid >> ylabel('y[n]') >> xlabel('Time Index (n)') 02/27/2020
177 What is the filter gain at From the plot above we see that in steady-state output peak amplitude over the input peak amplitude is about This result can be checked by finding the filter frequency response magnitude: >> w = [0 0.1*pi]; >> H = abs(freqz(5,[1 -1.5 0.8],w)) H= 16.6667 22.6520 02/27/2020 178
SECOND-ORDER FILTERS 02/27/2020 179 SECOND-ORDER FILTERS Second-order IIR filters allow for the possibility of complex conjugate pole and zero pairs, yet still have real coefficients The general second-order system function is Corresponding Difference Equation is 02/27/2020 180
The direct forms I and II can be used to implement this filter. Poles and Zeros To identify the poles and zeros of H(z) we can first convertto positive powers of z and then factor into poles and zeros: The coefficients are related to the roots: 02/27/2020 181 Given real coefficients, numerator and denominator, the roots occur either as two real values or as a complexconjugate pair Poles: From the quadratic formula: Real poles occur when
Complex-conjugate poles occur otherwise, and are given by 02/27/2020 182 Where And 02/27/2020 183 Zeros: Similar results hold for the zeros. If we factor out and
then replace with and replace with / and replace with / Example of Complex Poles and Zeros: First we need to apply the quadratic formula to the numerator and denominator to find the zeroes and poles: 02/27/2020 184 We can use MATLAB function tf2zp() to convert the system function form to a zero pole form, plus a gain term: >> [z, p, k] = tf2zp([3 2 2.5], [1 -1.5 0.8]) z=
-0.3333 + 0.8498i -0.3333 - 0.8498i p= 0.7500 + 0.4873i 0.7500 - 0.4873i k= 3 02/27/2020 185 IMPULSE RESPONSE Using a partial fraction expansion we can inverse transform any rational H(z) back to the impulse response h[n]. Example: We first perform long division to reduce the numerator order
02/27/2020 186 It can be shown that the partial fraction expansion will be of the form: Where: The impulse response is: This is a very general result, because the poles may be either real or complex conjugates: 02/27/2020 187 Note that if (no long division is required)
then the delta function term in is not needed. For complex conjugates poles further simplification is possible because and will also be complex conjugates. Complex conjugate Poles: In this case we first write and in polar form: And hence we can now write: h(n) becomes: 02/27/2020
188 EXAMPLE Find the impulse response corresponding to: MATLAB residuez() 02/27/2020 189 >> [A,p,K] = residuez([3 1], [1 -2*3/4*cos(pi/4) 9/16]) A= 1.5000 - 2.4428i 1.5000 + 2.4428i p= 0.5303 + 0.5303i 0.5303 - 0.5303i
K= [] 02/27/2020 190 >> [abs(A(1)) angle(A(1))] % Get mag and angle of A1 ans = 2.8666 -1.0201 02/27/2020 191 Plot from the above analysis and compare it with a MATLAB simulation using filter() n = 0:20;
h_anal = 2*2.867*(3/4).^n.*cos(pi/4*n - 1.02); subplot(211); stem(n,h_anal,'filled'); grid; ylabel('h_{analysis}[n]'); subplot(212); x = [1 zeros(1,20)]; h_sim = filter([3 1], conv([1 -3/4*exp(j*pi/4)],... [1 -3/4*exp(-j*pi/4)]),x); stem(n,h_sim,'filled'); grid; ylabel('h_{simulation}[n]'); xlabel('Time Index (n)'); 02/27/2020 192 5
5 4 4 3 3 2 2 1 1 0 0 -1 -1 -2 0 -2
2 0 4 2 6 4 8 6 10 12
12 . ( . ) . [] 8 ( 14 ) 10 16 14
18 16 20 18 20 5 5 4 4 3 3 2 2 1
1 0 0 -1 -1 -2 0 -2 2 0 4 2 6
4 8 6 10 8 12 Time Index (n) 10 14 12 16 14
18 16 20 18 20 Time Index (n) 02/27/2020 193 FREQUENCY RESPONSE With just a second-order filter we can realize a frequency
response that is lowpass, highpass, bandpass, or bandstop! Idealized filter responses are shown below: 02/27/2020 194 02/27/2020 195 Frequency response of second order filters is somewhat limited. Lowpass/Highpass: Place zeros on the unit circle as a conjugate pair and poles inside the: unit circle as a conjugate pair:
02/27/2020 196 Bandpass: {lace zeros at z1=1 and z2=-1 and poles inside the unit circle as a conjugate pair: 02/27/2020 197 Bandstop Filter: 02/27/2020 198
EXAMPLE OF AN IIR LOWPASS FILTER 02/27/2020 199 IIR FILTER DESIGN Lowpass Filter Amilitude Response Requriements 20 10| ( )|[ ] 1 -
0 50 3 2 02/27/2020
200 LOWPASS FILTER SPECIFICATION In the above the filter amplitude response is to remain within a 1 dB tolerance band of 0 dB or unity gain, from to , and have a gain less than -50 dB for This requirement specification defines a lowpass filter. A variety of IIR filters can be designed to meet these requirements, one such filter type is an elliptic filter, which is available in the MATLAB signal processing toolbox. 02/27/2020 201 MATLAB >>
ellipord(2*(pi/3/(2*pi)),2*(pi/2/(2*pi)),1,50) ans = 5 % The required filter order, N, to meet the design requirements. >>b = 0.0194 0.0211 0.0377 0.0377 0.0211
0.0194 % M = 5 >> a = 1.0000 -2.7580 4.0110 -3.3711 1.6541 -0.3796 % N = 5
02/27/2020 202 >> ZPLANE(B,A) 1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 -1
-0.5 0 0.5 1 Real Part 02/27/2020 203 MALTAB w = -pi:(pi/500):pi; H = freqz(b,a,w); subplot(211); plot(w,20*log10(abs(H)));
xlim([-pi +pi]); ylim([-60 1]); xlabel('\fontsize{20}-\pi \leq \omega \leq \pi'); ylabel('\fontsize{20}|H(e^{j\omega})|'); grid on; subplot(212); plot(w,angle(H)); xlim([-pi +pi]); ylim([-3.5 3.5]); xlabel('\fontsize{20}-\pi \leq \omega \leq \pi'); ylabel('\fontsize{20}\angleH(e^{j\omega})'); grid on; 02/27/2020 204 1dB 00
-10 -10 -20 -20 -30 -30 -40 -40 -50 -50 -60 -60
-3-3 -2-2 -1-1 00 3 00 11 -- 11
2 22 33 33 22 11 00 -1-1 -2-2 -3-3 -3-3 -2-2
-1-1 -- 22 02/27/2020 33 205 TRANSFER FUNCTION The fifth-order elliptic design has more than met the amplitude response requirements as it achieves the -50 dB gain level before
The complete H(z) is: 02/27/2020 206