Discrete-Time Signals And Systems - University Of Colorado Colorado Springs

CONTENTS x[k] h2[n-k] n k -1 0 1 2 3 4 5 6 7 8 k Signal Creation and Convolution Using Python Creating xŒn , h1Œn , and h2Œn at the IPython qtconsole or in the Jupyter notebook with module sigsys.py imported with alias ss The function ss.dstep(n) is useful in this problem The starting point is to create an Python ndarray to hold the axes values of n over some support interval Here I choose n 2 Œ 1; 11/ In In In In In In In In In In In In 2-20 [22]: [23]: [24]: [25]: [26]: [29]: [30]: [31]: [32]: [33]: [34]: [35]: import sk dsp comm.sigsys as ss n arange(-1,11) x ss.dstep(n) - ss.dstep(n-8) h1 ss.dstep(n) - ss.dstep(n-4) h2 0.8**n * ss.dstep(n) subplot(311) stem(n,x) subplot(312) stem(n,h1) subplot(313) stem(n,h2) tight layout() ECE 5650/4650 Modern DSP

2.4. DISCRETE-TIME SYSTEMS u[ n] - u[ n - 8] u[ n] - u[ n - 4] 0.8n u[ n] The input signal and system impulse responses in Python. Note I have made use of the helper function dstep(n) in creating the sequences The function description can be found to by typing ss.dstep? The Python code module sigsys.py (aliased as ss) contains many useful functions for doing signals and systems modeling and simulation; see scikit-dsp-comm The functions in ss include: In [82]: dir(ss) Out[82]: [BPSK tx, CIC, NRZ bits, NRZ bits2, OA filter, OS filter, PN gen, am rx, am rx BPF, am tx, biquad2, bit errors, cascade filters, conv integral,conv sum, cpx AWGN, cruise control, deci24, delta eps, dimpulse, downsample, drect, dstep, env det, ex6 2, eye plot, fft, fir iir notch, from wav, fs approx, ECE 5650/4650 Modern DSP 2-21

CONTENTS fs coeff, interp24, line spectra, lms ic, lp samp, lp tri, m seq, my psd, peaking, plot na, position CD, prin alias, rc imp, rect, rect conv, scatter, signal, simpleQuant, simple SA, sinusoidAWGN, soi snoi gen, splane, sqrt rc imp, step, ten band eq filt, ten band eq resp, to wav, tri, upsample, wavfile, zplane] Obtaining System Outputs In block diagram form we have the following: h1[n] y1[n] h2[n] y2[n] x[n] At the IPython prompt (assuming pylab is loaded and you have imported sk dsp comm.sigsys as ss) simply enter In [117]: n arange(-1,11) In [118]: x ss.dstep(n)-ss.dstep(n-8) In [119]: h1 ss.dstep(n) - ss.dstep(n-4) In [120]: y1,ny1 ss.conv sum(x,n,h1,n) Output support: (-2, 20) In [121]: subplot(211) In [122]: stem(ny1,y1) In [123]: n arange(-1,31) In [124]: x ss.dstep(n)-ss.dstep(n-8) In [125]: h2 0.8**n * ss.dstep(n) In [126]: y2,ny2 ss.conv sum(x,n,h2,n,extent (’f’,’r’)) Output support: (-2, 29) In [127]: subplot(212) In [128]: stem(ny2,y2) The function ss.conv sum performs discrete-time convolution using y signal.conv(x1,x2) from the scipy.signal package 2-22 ECE 5650/4650 Modern DSP

2.4. DISCRETE-TIME SYSTEMS The function ss.conv sum() wraps signal.conv() and adds sequence index manipulation for accurate plotting and properly truncates the convolution if the sequences have infinite extent to the left or right In the above convolution operation it is assumed that both sequences start at n D 1, thus the output should start at n D 2 and end at n D 11 C 31 D 42 For y1Œn the length of the output sequence will be the sum of the two input lengths minus one (why?), e.g., here we have 8 C 4 1 D 11 The theoretical length of y2Œn is infinite because h2Œn starts at n D 0 and extends to right out to C1 11 samples Finite Duration output since duration of input and system impulse response are finite Semi-infinite duration output since the system impulse response is semi-infinite ECE 5650/4650 Modern DSP 2-23

CONTENTS The tabulated sequence values are: In [129]: ny1 Out[129]: array([-2, -1, 0, 1, 2, 3, 4, 15, 16, 17, 18, 19, 20]) In [130]: y1 Out[130]: array([ 0., 0., 0., 0., 1., 0., 2., 0., 3., 0., 5, 4., 0., 6, 4., 0., 7, 8, 4., 0., 4., 0., 9, 10, 11, 12, 13, 14, 4., 3., 0.]) 2., 1., In [131]: ny2 Out[131]: array([-2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29]) In [132]: y2 Out[132]: array([ 0. , 2.952 , 3.32891136, 1.09081767, 0.35743914, 0.11712566, 0.03837973, 0. , 3.3616 , 2.66312909, 0.87265414, 0.28595131, 0.09370052, 0.03070379]) 1. , 3.68928 , 2.13050327, 0.69812331, 0.22876105, 0.07496042, 1.8 , 3.951424 , 1.70440262, 0.55849865, 0.18300884, 0.05996834, 2.44 , 4.1611392 , 1.36352209, 0.44679892, 0.14640707, 0.04797467, The fact that y1Œn is trapezoidal in shape should come as no surprise since convolving two rectangles in the continuoustime domain is known to produce a trapezoid Convolving Other Types of Sequences The Python function conv sum, found in the module sigsys.py, will actually properly convolve any two finite duration sequences, as well semi-infinite left- and right-sided sequences The help listing for the function, shown below, explains how this is done 2-24 ECE 5650/4650 Modern DSP

2.4. DISCRETE-TIME SYSTEMS def conv sum(x1,nx1,x2,nx2,extent (’f’,’f’)): """ Discrete convolution of x1 and x2 with proper tracking of the output time axis. Convolve two discrete-time signals using the SciPy function signal.convolution. The time (sequence axis) are managed from input to output. y[n] x1[n]*x2[n]. Parameters ---------x1 : ndarray of signal x1 corresponding to nx1 nx1 : ndarray time axis for x1 x2 : ndarray of signal x2 corresponding to nx2 nx2 : ndarray time axis for x2 extent : (’e1’,’e2’) where ’e1’, ’e2’ may be ’f’ finite, ’r’ right-sided, or ’l’ left-sided Returns ------y : ndarray of output values y ny : ndarray of the corresponding sequence index n Notes ----The output time axis starts at the sum of the starting values in x1 and x2 and ends at the sum of the two ending values in x1 and x2. The default extents of (’f’,’f’) are used for signals that are active (have support) on or within n1 and n2 respectively. A right-sided signal such as a n*u[n] is semi-infinite, so it has extent ’r’ and the convolution output will be truncated to display only the valid results. Examples ------- nx arange(-5,10) x drect(nx,4) y,ny conv sum(x,nx,x,nx) stem(ny,y) # Consider a pulse convolved with an exponential (’r’ type extent) h 0.5**nx*dstep(nx) y,ny conv sum(x,nx,h,nx,(’f’,’r’)) # note extents set stem(ny,y) # expect a pulse charge and discharge sequence """ ECE 5650/4650 Modern DSP 2-25

CONTENTS 2.5 Properties of LTI Systems The LTI system properties considered here are for the most part based on relationships involving the convolution sum and the impulse response Consider a cascade connection of systems with impulse responses h1Œn and h2Œn yŒn D fxŒn h1Œn g h2Œn Since convolution is commutative, it follows that yŒn D xŒn hŒn ; where hŒn D h1Œn h2Œn x[n ] h1 [ n ] h2 [ n ] y[n] x[n ] h2 [ n ] h1 [ n ] y [ n ] Equivalent Systems x[n ] h [ n ] h 1 [ n ] h 2 [ n ] y[n] Consider a parallel connection of systems with impulse responses h1Œn and h2Œn yŒn D xŒn h1Œn C xŒn h2Œn Since convolution is distributive, it follows that yŒn D xŒn hŒn ; where hŒn D h1Œn C h2Œn 2-26 ECE 5650/4650 Modern DSP

2.5. PROPERTIES OF LTI SYSTEMS h1 [ n ] x[n] y[ n] Equivalent Systems h2 [ n ] x[n] y[ n] h [ n ] h1 [ n ] h2 [ n ] Theorem: LTI systems are BIBO stable if and only if SD 1 X jhŒk j 1 kD 1 proof ( ˇ 1 ˇX ˇ jyŒn j D ˇ hŒk xŒn ˇ kD 1 ˇ 1 ˇ X ˇ k ˇ jhŒk jjxŒn ˇ k j kD 1 Now since jxŒn j Bx it follows that jyŒn j Bx 1 X jhŒk j kD 1 QED proof ) First assume that S D 1 and show that a (just one) bounded input can be found that will result in an unbounded output. One such sequence is ( h Œ n ; hŒn 0 xŒn D jhŒ n j 0; hŒn D 0 ECE 5650/4650 Modern DSP 2-27

CONTENTS Clearly Bx D 1, but yŒ0 is yŒ0 D 1 X kD 1 1 X jhŒk j2 xŒ k hŒk D DS !1 jhŒk j kD 1 QED If an LTI system is causal it then follows that hŒn D 0; n 0 A causal sequence is one which is zero for n 0 An impulse response with only a finite number of nonzero samples is called a finite-duration impulse response (FIR) Clearly FIR systems are always stable so long as the impulse response values have finite magnitude An impulse response with infinite duration is said to have an infinite-duration impulse response (IIR) 2-28 ECE 5650/4650 Modern DSP

2.5. PROPERTIES OF LTI SYSTEMS Example 2.7: Impulse, Step, and Exponential An ideal delay has impulse response hŒn D ıŒn nd ; nd 0 P Since n jhŒn j D 1 we know that the system is stable, and since hŒn D 0 for n 0 provided nd 0, we also know that the system is causal. If hŒn D uŒn then the system is an accumulator. Clearly an accumulator is IIR and causal, but since 1 X juŒn j ! 1 nD 1 we see that the accumulator is also unstable. Consider a system with exponential impulse response hŒn D anuŒn , a real. The system is stable provided SD 1 X jajn 1 nD0 – Recall the infinite geometric series formula 1 X rn D nD0 Thus SD 1 1 1 1 jaj and the system is stable provided r ; jrj 1 1 jaj 1 ECE 5650/4650 Modern DSP 2-29

CONTENTS 2.6 Linear Constant-Coefficient Difference Equations A subclass of LTI systems which is of considerable interest are those systems which have input output relationship satisfying an N th-order constant coefficient difference equation (LCCDE) of the form N X ak yŒn k D kD0 M X bk xŒn k kD0 Example 2.8: Simple LCCDE Examples As a specific example consider an accumulator yŒn D n X xŒk kD 1 Note that the first backward difference formed at the output returns the original input yŒn yŒn 1 D xŒn ; or yŒn D yŒn 1 C xŒn The form yŒn D yŒn 1 C xŒn shows that the accumulator can be written as an LCCDE with N D 1, a0 D 1, a1 D 1, M D 0, and b0 D 1 2-30 ECE 5650/4650 Modern DSP

2.6. LINEAR CONSTANT-COEFFICIENT DIFFERENCE EQUATIONS x[n] y[n] One-Sample Delay y[n – 1] Accumulator LCCDE block diagram As a second example consider the moving average system with M1 D 0 so that the system is causal M 2 X 1 yŒn D xŒn M2 C 1 k kD0 This system fits the LCCDE form given above if we let N D 0, M D M2, a0 D 1, and bk D 1 .M2 C 1/, for 0 k M2 The impulse response for the causal moving average system, as obtained by setting xŒn D ıŒn , is given by hŒn D 1 .uŒn M2 C 1 uŒn M2 1 / Note that since uŒn M D ıŒn M uŒn , we can also write hŒn D 1 .ıŒn M2 C 1 ıŒn M2 1 / uŒn Recalling that uŒn is the impulse response of an accumulator, we observe that the above form of the moving average system impulse response is in the form of a cascade system ECE 5650/4650 Modern DSP 2-31

CONTENTS x[n ] Attenuator 1 ---------------M2 1 ( M2 1 ) Sample Delay – x1 [ n ] Accumulator System y[ n] Causal moving average, recursive LCCDE block diagram The first stage produces output x1Œn D 1 .xŒn M2 C 1 xŒn M2 1 / The second stage is the accumulator, thus the total input/output relationship is yŒn D yŒn 1 C 1 .xŒn M2 C 1 xŒn M2 1 / which corresponds to an LCCDE with coefficients N D 1, a0 D 1, a1 D 1, b0 D bM2C1 D 1 .M2 C 1/, and bk D 0 otherwise From the above we conclude that the simple causal moving average system does not have a unique LCCDE representation 2-32 ECE 5650/4650 Modern DSP

2.6. LINEAR CONSTANT-COEFFICIENT DIFFERENCE EQUATIONS Example 2.9: A Nonlinear System Example An example of a nonlinear system with difference equation representation (not an LCCDE however) is the square root algorithm 1 xŒn yŒn D yŒn 1 C 2 yŒn 1 If we let the input be an amplitude step (i.e. xŒn D AuŒn ), then p as the initial estimate pthe initial condition yŒn 1 serves of A, and the response yŒn tends to A as n increases – Let A D 2, and yŒ 1 D 1, then yŒ0 D 3 2, yŒ1 D 1:416667, and yŒ2 D 1:4142157 – Let A D 2, and yŒ 1 D 1:5, then yŒ0 D 1:41667, and yŒ1 D 1:4142157 p – Note that 2 D 1:4142136 2.6.1 Classical Solution of LCCDEs The solution of a LCCDE can be written as the sum yŒn D yhŒn C yp Œn where yhŒn is the natural response which is a solution to the homogeneous equation N X ak yŒn k D 0 kD0 ECE 5650/4650 Modern DSP 2-33

CONTENTS and yp Œn is the forced response which is a particular solution of the nonhomogeneous equation N X ak yŒn k D xŒn kD0 To solve the homogeneous equation we assume a solution of the form yhŒn D z n If we substitute this solution into the homogeneous equation we obtain N X ak yhŒn k D 0 kD0 or zn N .a0z N C a1z N 1 C C aN / D 0 ƒ‚ „ characteristic polynomial The system characteristic polynomial, as indicated above, will have N roots, zi ; i D 1; : : : ; N , which may be real/complex and distinct/repeated. Assuming the coefficients a0; a1; : : : ; aN are real, then complex roots will always occur in complexconjugate pairs. If the roots are distinct, then the general homogeneous solution is of the form yhŒn D A1z1n C A2z2n C C AN zNn where the coefficients A1; A2; : : : ; AN are determined from the system initial conditions. 2-34 ECE 5650/4650 Modern DSP

2.6. LINEAR CONSTANT-COEFFICIENT DIFFERENCE EQUATIONS If a repeated root of multiplicity m N should occur, then the homogeneous solution is of the form yhŒn D A1z1n C A2nz1n C C Amnm 1z1n n C AmC1zmC1 C C AN zNn To solve the nonhomogeneous equation and obtain the forced response, we can use the method of undetermined coefficents which involves guessing a particular solution, typically having the same form as xŒn Particular solutions for several types of inputs xŒn (uŒn implied) yp Œn (uŒn implied) A (const.) B AM n BM n AnM B0nM C K1nM 1 C C KM An n M An.B0nM C K1nM 1 C C KM / A cos !on B1 cos !on C B2 sin !on A sin !on B1 cos !on C B2 sin !on Example 2.10: Classical solution: Consider the first order system yŒn D ayŒn 1 C xŒn jaj 1 with auxiliary (initial) condition yŒ 1 D c. Let xŒn D KuŒn . The homogeneous equation is yŒn ayŒn 1 D 0 and the assumed solution is yhŒn D z n, thus we must solve zn ECE 5650/4650 Modern DSP az n 1 D 0 ) z1 D a 2-35

CONTENTS The homogeneous solution is yhŒn D A1an The particular solution of the nonhomogeneous equation is obtained by assuming yp Œn D BuŒn where B is a scale factor Substituting the assumed solution into the nonhomogeneous equation we obtain BuŒn 1 D KuŒn aBuŒn To solve for B we evaluate the above equation where none of the terms vanish, that is for n 1 B aB D K or B D so yp Œn D K 1 a K 1 a uŒn The total solution can now be written as yŒn D yhŒn C yp Œn K D A1 a n C ; n 0 1 a To solve for A1 use the fact that yŒ 1 D c and write yŒ0 D A1 C 2-36 K 1 a ECE 5650/4650 Modern DSP

2.6. LINEAR CONSTANT-COEFFICIENT DIFFERENCE EQUATIONS also yŒ0 D ayŒ 1 C K thus K A1 D ayŒ 1 C K 1 a Rewriting, the total solution becomes yŒn D ca nC1 D ca Ka 1 a anC1 CK ; n 0 1 a 1 Note that the zero state response (i.e. c D 0) is anC1 uŒn yŒn D K 1 a 1 Recursive solution: In general a recursive solution for yŒn can be obtained using the recurrence formula yŒn D N X ak kD1 a0 yŒn k C M X br rD0 a0 xŒn r along with the auxiliary conditions yŒ 1 ; : : : ; yŒ N . We start by finding yŒ0 . A similar procedure can be used to find yŒn for n 0. For the example we are considering we can write yŒ0 yŒ1 yŒ2 ::: yŒn ECE 5650/4650 Modern DSP D D D D D ac C K a2c C aK C K a3c C a2K C aK C K ::: anC1c C anK C C aK C K 1 anC1 nC1 D ca CK 1 a 2-37

CONTENTS For systems with input/output satisfying an LCCDE we can make the following observations: – The system will in general not be linear. Recall that in the previous example when K D 0 (i.e. xŒn D 0) we found that yŒn D canC1 0 unless c D 0 – The system will in general not be time-invariant. Again with reference to the previous example we note that with x1Œn D xŒn no D KuŒn no the output is y1Œn D ca nC1 CK 1 anC1 no uŒn .1 a/ no – If the system is initially at rest, then the system will be linear, time-invariant, and causal Example 2.11: Find the impulse response, hŒn , of the second order system yŒn 3yŒn 1 4yŒn 2 D xŒn C 2xŒn 1 Since hŒn D yŒn when xŒn D ıŒn we must solve yŒn 3yŒn 4yŒn 2 D ıŒn C 2ıŒn 1 1 C 4yŒn 2 C ıŒn C 2ıŒn 1 1 or yŒn D 3yŒn For n 2 we have the homogeneous equation yŒn 2-38 3yŒn 1 4yŒn 2 D 0 ECE 5650/4650 Modern DSP

2.6. LINEAR CONSTANT-COEFFICIENT DIFFERENCE EQUATIONS and characteristic equation z 2 3z 4 D 0, which has solution yŒn D A1. 1/n C A24n; n 2 A1 and A2 are determined by satisfying the nonhomogeneous equation for n D 0 and 1 Assuming zero initial conditions we must have yŒ0 D 1 D A1 C A2 yŒ1 D 3 C 2 D 5 D A1 C 4A2 which implies that A2 D 1 A1 D 1 5 and A2 D 6 5 A1 and 5 D 5A1 C 4, thus Substituting the constants we obtain 6 1 . 1/n C 4n uŒn hŒn D 5 5 Note, if M D N in the general case, then in order to satisfy the nonhomogeneous equation we would need to modify the solution by including the term A0ıŒn – For the more general case of M N we add to the solution a term of the form M XN Ak ıŒn k kD0 A case in point being the LCCDE yŒn ECE 5650/4650 Modern DSP yŒn 1 D xŒn C xŒn 1 2-39

CONTENTS Since M D 1 D N and yŒn yŒn 1 D 0 for n 2 we have as general solution for hŒn hŒn D A1.1/nuŒn C A0ıŒn D A1 C A0ıŒn ; n 0 To find A1 and A0 we solve yŒ0 D 1 D A1 C A0 yŒ1 D 2 D A1 Thus A1 D 2 and A0 D 1 giving hŒn D 2uŒn ıŒn As a check note that yŒn yŒ0 yŒ1 yŒ2 ::: yŒn D D D D yŒn 1 C ıŒn C ıŒn 0C1C0D1 1C0C1D2 2C0C0D2 ::: D 2; n 1 1 Example 2.12: The Python Filter Function for LCCDEs Given an LCCDE description of a system, the Python scipy.signal function lfilter() can be used to process arbitrary signals through the system. We begin with a description of the function lfilter() 2-40 ECE 5650/4650 Modern DSP

2.6. LINEAR CONSTANT-COEFFICIENT DIFFERENCE EQUATIONS Definition: signal.lfilter(b, a, x, axis -1, zi None) Docstring: Filter

2.4. DISCRETE-TIME SYSTEMS 2.4 Discrete-Time Systems Definition: A discrete-time system is an operator that maps an input sequence into an output sequence y„n"DTfx„n"g x[n] y[n] T{ ! } Example 2.3:Moving Average (MA) Operator Define Tfgsuch that y„n"D 1 M 1CM 2C1 XM2 kDM1 x„n k" The averaging is causal if we set M 1 D0so that .