660 likes | 922 Views
University of Manchester School of Computer Science Comp30291: Digital Media Processing 2009-10 Section 3 : Discrete-time LTI systems. 3.1 Introduction Consider a DSP system as implemented by a digital processor. Takes discrete time input signal { x[ n ] },
E N D
University of Manchester School of Computer Science Comp30291: Digital Media Processing 2009-10 Section 3 : Discrete-time LTI systems Comp30291 Section 3
3.1 Introduction • Consider a DSP system as implemented by a digital processor. • Takes discrete time input signal { x[ n ] }, • & produces an output signal { y[ n ] }. Comp30291 Section 3
Sequences representing digital signals • {x[n]} is sequence whose value at t=nT is x[n]. • Similarly for {y[n]}. • T is sampling interval in seconds (T = 1/FS). • {x[n-N]} is sequence whose value at t=nT is x[n-N]. • {x[n-N]} is {x[n]} with every sample delayed by N sampling intervals. Comp30291 Section 3
A x[n] y[n] Examples of discrete time systems • (i) Discrete time ‘amplifier’: • y[n] = A . x[n]. • Described by‘difference equation’: y[n] = A x[n]. • Represented in diagram form by a ‘signal flow graph’: Comp30291 Section 3
(ii) Non-recursive digital filter • Output at t=nT obtained by weighting & summing present & previous input samples: e.g. • y[n] = A0 x[n] + A1 x[n-1] + A2 x[n-2] + A3 x[n-3] + A4x[n-4] • This is a ‘non-recursive difference equation’ • Represented by signal flow graph below. • Boxes marked ‘ z -1 ‘ produce a delay of one sampling interval. Comp30291 Section 3
A1 y[n] x[n] z-1 B2 • (iii) Recursive digital filter • Output at t= nT calculated from a recursive difference equation: • e.g. y[n] = A 1 x[n] - B 2 y[n-1] • Represented by signal flow graph below. • Recursive means that previous values of y[n] as well as present & previous values of x[n] are used to calculate y[n]. Comp30291 Section 3
y[n] x[n] • (iv) A non-linear system • Output at t=nT calculated from some non-linear equation, e.g. • y[n] = (x[n]) 2 • Represented below: Comp30291 Section 3
3.2. Linearity & time-invariance • (I) A DSP system is linear if: • Given any two discrete time signals {x 1 [n]} & {x 2 [n]}, • if {x 1 [n]}{y 1 [n]} & {x 2 [n]}{y 2 [n]} • then response to k 1{x 1[n]} + k 2{x 2[n]} • must bek 1{y1[n]} + k 2{y 2 [n]} • for any values of k 1 and k 2 , • To multiply a sequence by k, multiply each element by k, • k{x[n]} = {k x[n]}. • To add two sequences together, add corresponding samples, • {x[n]} + {y[n]} = {x[n] + y[n]}.) Comp30291 Section 3
(II) A DSP system is ‘time-invariant’ if: • Given any discrete time signal {x[n]}, • if response to {x[n]} is {y[n]}, • response to {x[n-N]} must be {y[n-N]} for any N. • Delaying input by N samples only delays output by N samples. • An LTI system is both linear & time-invariant • Examples (i), (ii) and (iii) are LTI whereas (iv) is not LTI . Comp30291 Section 3
3.3. Discrete time unit impulse • Useful to consider response of LTI systems to a discrete time unit • impulse, or in short an impulse denoted by {d [n] } with: {d[n-N]} is delayed impulse where the only non-zero sample occurs at n=N rather than at n=0. Comp30291 Section 3
Impulse-response • When input is {d[n]}, output is impulse-response {h[n]}. • If impulse-response of an LTI system is known, its response to any other input signal may be obtained. Comp30291 Section 3
X3 X1 X4 X2 X5 - 1 - 1 - 1 - 1 z z z z x[n] A3 A4 A5 A2 A1 Y y[n] • 3.4. Implementing signal-flow-graphs • Consider the non-recursive signal flow graph below with A1, A2, A3, A4, A5 set to specific constants. • Notice the labels X1, X2, etc. • Realised by MATLAB program & flow-diagram on next slides. Comp30291 Section 3
Flow-diagram for non-recursive signal-flow-graph Set values of A1,A2, A3, A4, A5 Set X1, X2, X3, X4, X5 to zero INPUT X1 Y = A1*X1 + A2*X2 + A3*X3 +A3*X4 + A5*X5 OUTPUT Y X5 = X4; X4 = X3; X3 = X2; X2 =X1 Comp30291 Section 3
MATLAB program for non-recursive signal flow graph clear all; A1=1; A2=2; A3=3; A4=-4; A5=5; X1=0; X2=0; X3=0; X4=0; X5=0; while 1 X1 = input( 'X1 = '); Y= A1*X1 + A2*X2 + A3*X3 + A4*X4 + A5*X5 ; disp([' Y = ' num2str(Y)]); X5 = X4 ; X4 = X3; X3 = X2 ; X2 =X1; end; Comp30291 Section 3
More efficient version A = [1 2 3 -4 5 ]’ ; x = [0 0 0 0 0 ]’ ; while 1 x(1) = input( ‘x(1) = '); Y=0; for k = 1 : 5 Y = Y + A(k)*x(k); end; disp([' Y = ' num2str(Y)]); for k=5:-1:2 x(k) = x(k-1); end; end; Comp30291 Section 3
Even more efficient version A = [1 2 3 -4 5 ]' ; x = [0 0 0 0 0 ]' ; while 1 x(1) = input( 'x(1) = '); Y = A(1)*x(1); for k = 5 : -1: 2 Y = Y + A(k)*x(k); x(k) = x(k-1); end; disp(['Y = ' num2str(Y)]); end; Comp30291 Section 3
Comments • Ready to be converted to DSP assembler. Y = Y + A(k)*x(k); x(k) = x(k-1); One DSP instruction:Mult-acc-shift - MACS • The ‘while 1’ statement initiates an infinite loop. • Program runs for ever or until interrupted by ‘CONTROL+C’ • Either of the following prints out value of Y: • disp(['Y = ' num2str(Y)]); • disp(sprintf(‘Y=%d’,Y)); • A = [1 2 3 -4 5]’ makes A a column vector (not a row). Note ’. Comp30291 Section 3
Use of ‘filter’ to implement non-rec signal-flow-graph - & thus to produce its impulse-response clear all; x=[0 1 0 0 0 0 0 0 0 0]'; a = [1 2 3 -4 5]'; y=filter(a,1,x); y • Just writing ‘y’ without ‘;’ prints out the array. Comp30291 Section 3
Impulse-response for non-rec signal-flow-graph • Use any of previous versions & enter values for X1 :0, 0, 0, 1, 0, 0, 0, 0, .... • Sequence of output samples printed out will be : • 0, 0, 0, A1, A2, A3, A4, A5, 0, 0, .... • Impulse-response can also be obtained by tabulation (later). • Output must be zero until input becomes equal to 1 at n=0 • Impulse response is: • {..., 0, ..., 0, A1, A2, A3, A4, A5, 0, 0, ... ,0, ...} • where the sample at n=0 is underlined. • Only five non-zero output samples are observed. • This is a ‘ finite impulse-response ‘ (FIR). Comp30291 Section 3
Exercise 3.1 • Calculate impulse-responses, by tabulation, for: (i) y[n] = x[n] + 2 x[n-1] + 3 x[n-2] - 4 x[n-3] + 5 x[n-4] (ii) y[n] = 4 x[n] - 0.5 y[n-1] • Difference eqn (i) will produce a finite impulse-response. • Difference eqn (ii) produces infinite response whose samples gradually reduce in amplitude but never quite become zero. Comp30291 Section 3
Impulse-response for example (i) by tabulation y[n] = x[n] + 2 x[n-1] + 3 x[n-2] - 4 x[n-3] + 5 x[n-4] n x[n] x[n-1] x[n-2] x[n-3] x[n-4] y[n] -1 0 0 0 0 0 0 0 1 0 0 0 0 1 1 0 1 0 0 0 2 2 0 0 1 0 0 3 3 0 0 0 1 0 -4 4 0 0 0 0 1 5 5 0 0 0 0 0 0 :::::: : Impulse response is: {.. 0, .., 0, 1, 2, 3, -4, 5, 0, .., 0, ...} Comp30291 Section 3
n x[n] y[n-1] y[n] 0 1 0 4 1 0 4 -2 2 0 -2 1 3 0 1 -0.5 4 0 -0.5 0.25 5 0 0.25 -0.125 :: : : Impulse-response by tabulation for Example (ii) y[n] = 4 x[n] - 0.5 y[n-1] Impulse response is: {.., 0, .., 0, 4, -2, 1, -0.5, 0.25, -0.125, ...) Comp30291 Section 3
n x[n] y[n-1] y[n] 0 1 0 4 1 0 4 -8 2 0 -8 16 3 0 16 -32 4 0 -32 64 5 0 64 128 :: : : Further example: Impulse-response by tabulation y[n] = 4 x[n] - 2 y[n-1] • Impulse response is: {.., 0, .., 0, 8, 16, -32, 64, 128, ...) • This IIR filter is ‘unstable’ (see later) Comp30291 Section 3
3.5. FIR & IIR digital filters • A ‘digital filter’ is a digitally implemented LTI discrete time system governed by a difference equation of finite order; e.g. : (i) y[n] = x[n] + 2 x[n-1] + 3 x[n-2] - 4 x[n-3] + 5 x[n-4] (ii) y[n] = 4 x[n] - 0.5 y[n-1] • Difference equation (i) is "non-recursive" • & produces a finite impulse response (FIR). • Difference equation (ii) is " recursive " . • Impulse-response of a recursive difference equation can have an infinite number of non-zero terms. • In this case it is an infinite impulse-response (IIR). Comp30291 Section 3
MATLAB program for recursive signal flow graph - example (ii) clear all; Y2=0; while 1 X1 = input( 'X1 = '); Y1= 4*X1 - 0.5*Y2 ; Y2 = Y1; % for next time round disp([' Y1 = ' num2str(Y)]); end; Comp30291 Section 3
Use of ‘filter’ for FIR & IIR digital filters y = filter(A, B, x) filters signal in array x to create array y. For FIR example (i), A = [ 1 2 3 -4 5 ] & B = [1]. For IIR example (ii), A = [4], B = [1 0.5] Consider a third IIR example: y[n] = 2x[n] + 3x[n-1] + 4x[n-2] -0.5 y[n-1] - 0.25 y[n-2] In this case set A = [2 3 4] and B = [1 0.5 0.25]. Why are A & B defined in this way? Comp30291 Section 3
Definition of A & B in ‘y = filter(A,B,x);’ • A digital filter has a ‘system function’ which is with b0 = 1 for the difference equation: • A contains [a0 a1 ... aN] & B contains [b0, b1, ..., bM] . • Reasons for this & more details will be given later in course. • ‘filter’ can be used without knowing why H(z) is defined in this way. Comp30291 Section 3
3.6. Discrete time convolution If impulse-response of an LTI system is {h[n]} its response to any input {x[n]} is an output {y[n] } whose samples are given by the following ‘convolution’ formulae: • Formulae are equivalent. • Clearly, if we know impulse-response {h[n]} we can produce the response to any other input sequence, from either of these formulae. • Proof of convolution in Appendix 3A. Comp30291 Section 3
Example Calculate response of a system with impulse response: {h[n]} = { ..., 0,..., 0, 1, 2, 3, -4, 5, 0, .....0, .... } to {x[n]} = { ... 0, ... , 0 , 1, 2, 3, 0, ..., 0, ....} Solution: By first discrete time convolution formula, This is difference equation for an LTI system with impulse response {.., 0, .., 0, 1, 2, 3, -4, 5, 0, .., 0, ..} Comp30291 Section 3
Completing the Example • Program discussed earlier implements this difference equation, • We could run it, and enter {0 1 2 3 0 0 0 0...}, • Output sequence produced is what we want. • Alternatively, we could use tabulation as follows: Comp30291 Section 3
Response of: y[n] = x[n] +2x[n-1]+3x[n-2]-4x[n-3] +5x[n-4] to input sequence {...,0,1,2,3,0,...} by tabulation n x[n] x[n-1] x[n-2] x[n-3] x[n-4] y[n] : : : : : : : -1 0 0 0 0 0 0 0 1 0 0 0 0 1 1 2 1 0 0 0 4 2 3 2 1 0 0 10 3 0 3 2 1 0 8 4 0 0 3 2 1 6 5 0 0 0 3 2 -2 6 0 0 0 0 3 15 7 0 0 0 0 0 0 : : : : : : : {y[n]} = { .... 0, ....., 0, 1, 4, 10, 8, 6, -2, 15, 0, ...., 0, ....} Comp30291 Section 3
3.7. Stability An LTI system is stable if its impulse-response {h[n]} satisfies: This means that {h[n]} must be either an FIR or an IIR whose samples decay towards zero as n . Comp30291 Section 3
3.8. Causality An LTI system operating in real time must be ‘causal’ which means that its impulse-response {h[n]} must satisfy: h[n] = 0 for n < 0. Non-causal system would need “crystal ball ” to predict future. Comp30291 Section 3
h[n] Causal, but Looks stable h[n] not stable but is not causal. n n h[n] Causal & looks stable. n Illustration of stability & causality Comp30291 Section 3
3.9. Relative Frequency • Study effect of digital filters on sinusoids of different frequencies. • Discrete time sinusoid obtained by sampling A cos(t + ). • If sampling frequency is Fs Hz, and T=1/Fs, we obtain: • x[n] = A cos(nT + ) • = A cos(n + ) • = T is ‘relative frequency’ of sampled sinusoid. • Units of are 'radians per sample'. Comp30291 Section 3
Radians per sample • To convert back to true frequency (radians/s ) multiply by Fs. • radians / sample samples / second = radians / second • Analogue signals in range 0 to FS/2 Hz =1/(2T) • Restricts to the range 0 to . Comp30291 Section 3
cos( t ) =2 / (8T) t -4T -3T 3T 4T -T T cos( T n ) = 2 / 8 n -4 -3 3 4 -1 1 Comp30291 Section 3
Values of & corresponding true frequencies Relative frequency True frequency (radians/sample) (radians/s) (Hz) 0 0 0 /6 fS/6 fS/12 /4 fS/4 fs/8 /3 fS/3 fs/6 /2 fS/2 fs/4 2/3 2fS/3 fs/3 fSfS/2 ‘radians / sample’ ‘samples / second’ = ‘radians / second’ Comp30291 Section 3
3.10. Relative frequency response It us useful to analyse the response to a sampled sinusoid: x[n] = Acos(n + ) To begin with, set A=1, =0 and remember de Moivre’s Theorem: Easier to calculate response to the complex signal x[n] = ejnthan to cos(n) directly. Comp30291 Section 3
Relative frequency response (cont) If x[n] = ejn is applied to a system with impulse-response {h[n]}, output would be, by convolution : Comp30291 Section 3
Discrete time Fourier Transform (DTFT) • H( e j ) is the DTFT of { h[n] }. • It is called the ‘relative frequency-response’ • Complex number for any value of . • Note similarity to the analogue Fourier transform. Comp30291 Section 3
e j nH(e j ) e j n LTI {h[n]} Recap • If input is x[n]=e j n, output is same sequence with each element multiplied by H( e j ). Comp30291 Section 3
3.11. Gain & phase responses • G() = |H( e j )| is ‘gain’ • () = arg ( H( e j ) ) is ‘phase lead’. • Both vary with . • Can express: H( e j ) = G() e j () • If input is {A cos(n)}, output is: { G()A cos(n + ()) } When input is sampled sinusoid of relative frequency , output is sinusoid of same frequency , but with amplitude scaled by G() & phase increased by (). Comp30291 Section 3
Gain & phase response graphs again 20log10[G()] dB -() -() G() in dB 0 /4 /2 3/4 Comp30291 Section 3
Graphs of G() & () against . • G() often converted to dBs by calculating 20 log10( G() ). • Restrict to lie in range 0 to • Adopt a linear horizontal frequency scale. Comp30291 Section 3
Example • Derive frequency-response of FIR digital filter below. • Impulse-response is: {.., 0, .., 0, 1, 2, 3, -4, 5 ,0, .., 0,..}. • By the formula established above, • H( e j ) = 1 + 2 e - j + 3 e -2 j - 4 e -3 j + 5 e - 4 j Comp30291 Section 3
Example: Plot gain & phase responses of this FIR filter • To do this ‘by hand’, we would need to take modulus & phase of expression for H( e j ). • Fortunately we almost never need to do this ‘by hand’ • (except for very simple expressions). • Best done by MATLAB. • To calculate gain & phase responses for • H( e j ) = 1 + 2 e - j + 3 e -2 j - 4 e -3 j + 5 e - 4 j • start MATLAB and type: • freqz( [1 2 3 -4 5] ); Comp30291 Section 3
‘freqz’ graphs of gain & phase responses • Frequency scale normalised to fS/2 & labelled 0 to 1 • instead of 0 to . • To plot freq-response for an IIR filter equally straightforward using ‘freqz’ if you know how to define the system function. • See earlier. • More about this in a later section. Comp30291 Section 3
FIR & IIR digital filter design • In next section we see how FIR filter coefficients can be chosen to achieve a particular type of gain & phase response. • Tells us how the MATLAB function ‘fir1’ works in principle. • In later sections, principles of IIR digital filter design will be considered also Comp30291 Section 3