Digital Signal Processing - University of Cambridge

1 Digital Signal ProcessingMarkus KuhnComputer 2004 Part IITextbooks Lyons:Understanding Digital Signal , 2004. ( 45) Oppenheim, Schafer:Discrete-time Signal ed., Prentice-Hall, 1999. ( 47) J. Stein: Digital Signal Processing a computer science , 2000. ( 74) Smith: Digital Signal Processing a practical guide forengineers and , 2003. ( 40) K. Steiglitz:A Digital Signal Processing primer with appli-cations to Digital audio and computer ,1996. ( 40) Sanjit K. Mitra: Digital Signal Processing a , 2002. ( 38)2 Digital Signal processingConcerned with algorithms to interpret, transform, and model wave-forms and the information they contain. Some typical applications: geophysicsseismology, oil exploration astronomyspeckle interferometry, VLBI experimental physicssensor-data evaluation communication systemsmodulation/demodulation, channelequalization, echo cancellation aviationradar, radio navigation musicsynthetic instruments, audio effects,noise reduction medical diagnosticsECG, EEG, MEG, audiology,computer tomography, magnetic-resonance imaging consumer electronicsperceptual coding of audio and videoon DVDs, speech synthesis, speechrecognition securitysteganography, Digital watermarking,biometric identification, SIGINT engineeringcontrol systems, feature extractionfor pattern recognition3 Sequences and systemsAdiscrete sequence{xn}is a sequence of numbers.

2 , x 2, x 1, x0, x1, x2, ..wherexndenotes then-th number in the sequence (n Z). A discretesequence maps integer numbers onto real (or complex) notation is not well standardized. Some authors writex[n]instead ofxn, othersx(n).Where a discrete sequence{xn}samples a continuous functionx(t)asxn=x(ts n) =x(n/fs),we calltsthesampling periodandfs= 1/tsthesampling systemTreceives as input a sequence{xn}and transformsit into an output sequence{yn}=T{xn}:.. , x2, x1, x0, x 1, .. , y2, y1, y0, y 1, ..discretesystemT4 Properties of sequencesA sequence{xn}isabsolutely summable Xn= |xn|< square summable Xn= |xn|2< periodic k >0 : n Z:xn=xn+kA square-summable sequence is also called anenergy Signal , and Xn= |xn|2its energy. This terminology reflects that ifUis a voltage supplied to a load resistorR, thenP=U I=U2/Ris the power even where we drop physical units ( , volts) for simplicity in sequence calcu-lations, it is still customary to refer to the squared values of a sequenceaspowerand to its sum or integral over time non-square-summable sequence is apower signalif its average powerlimk 11 + 2kkXn= k|xn| sequencesUnit-step sequence:un= 0, n <01, n 0 Impulse sequence: n= 1, n= 00, n6= 0=un un 16 Types of discrete systemsAcausal systemcannot look into the future:yn=f(xn, xn 1, xn 2.)

3 Amemory-less systemdepends only on the current input value:yn=f(xn)Adelay systemshifts a sequence in time:yn=xn dTis atime-invariant systemif for anyd{yn}=T{xn} {yn d}=T{xn d}.Tis alinear systemif for any pair of sequences{xn}and{x n}T{a xn+b x n}=a T{xn}+b T{x n}.7 Examples:Theaccumulator systemyn=nXk= xkis a causal, linear, time-invariant system with memory, as are theback-ward difference systemyn=xn xn 1,theM-point moving average systemyn=1MM 1Xk=0xn k=xn M+1+ +xn 1+xnMand theexponential averaging systemyn= xn+ (1 ) yn 1= Xk=0(1 )k xn for time-invariant non-linear memory-less systems areyn=x2n, yn= log2xn, yn= max{min{ 256xn ,255},0},examples for linear but not time-invariant systems areyn= xn, n 00,n <0=xn unyn=x n/4 yn=xn (e jn)and examples for linear time-invariant non-causal system areyn=12(xn 1+xn+1)yn=9Xk= 9xn+k sin( k ) k [ + cos( k/10)]9 Constant-coefficient difference equationsOf particular practical interest are causal linear time-invariantsystemsof the formyn=b0 xn NXk=1ak yn kz 1z 1z 1ynxnb0yn 1yn 2yn 3 a1 a2 a3 Block diagram representationof sequence operations:z 1xnxnxnx nxn 1axnaxn+x nDelay:Addition:Multiplicationby constant.

4 Theakandbmareconstant xn mz 1z 1z 1xnynb0b1b2b3xn 1xn 2xn 3or the combination of both:NXk=0ak yn k=MXm=0bm xn mz 1z 1z 1z 1z 1z 1b0yn 1yn 2yn 3xna 10b1b2b3xn 1xn 2xn 3 a1 a2 a3ynThe MATLAB functionfilteris an efficient implementation of the last linear time-invariant (LTI) systems can be represented in the formyn= Xk= ak xn kwhere{ak}is a suitably chosen sequence of operation over sequences is calledconvolutionand defined as{pn} {qn}={rn} n Z:rn= Xk= pk qn {yn}={an} {xn}is a representation of an LTI systemT, with{yn}=T{xn}, then we call the sequence{an}theimpulse responseofT, because{an}=T{ n}.12 Convolution examplesABCDEFA BA CC AA ED EA F13 Properties of convolutionFor arbitrary sequences{pn},{qn},{rn}and scalarsa,b: Convolution is associative({pn} {qn}) {rn}={pn} ({qn} {rn}) Convolution is commutative{pn} {qn}={qn} {pn} Convolution is linear{pn} {a qn+b rn}=a ({pn} {qn}) +b ({pn} {rn}) The impulse sequence (slide 6) is neutral under convolution{pn} { n}={ n} {pn}={pn} Sequence shifting is equivalent to convolving with a shiftedimpulse{pn d}={pn} { n d}14 Can all LTI systems be represented by convolution?

5 Any sequence{xn}can be decomposed into a weighted sum of shiftedimpulse sequences:{xn}= Xk= xk { n k}Let s see what happens if we apply a linear( )time-invariant( )systemTto such a decomposed sequence:T{xn}=T Xk= xk { n k}!( )= Xk= xk T{ n k}( )= Xk= xk { n k} T{ n}= Xk= xk { n k}! T{ n}={xn} T{ n} The impulse responseT{ n}fully characterizes an LTI 1 Afinite-length sequenceis non-zero only at a finite number ofpositions. Ifmandnare the first and last non-zero positions, respectively,then we calln m+1thelengthof that sequence. What maximum lengthcan the result of convolving two sequences of lengthkandlhave?Exercise 2 The length-3 sequencea0= 3,a1= 2,a2= 1is convolvedwith a second sequence{bn}of length 5.(a) Write down this linear operation as a matrix multiplication involving amatrixA, a vector~b R5, and a result vector~c.(b) Use MATLAB to multiply your matrix by the vector~b= (1,0,0,2,2)and compare the result with that of using theconvfunction.(c) Use the MATLAB facilities for solving systems of linear equations toundo the above convolution 3(a) Find a pair of sequences{an}and{bn}, where eithercontains at least three different values and where the convolution{an} {bn}results in an all-zero sequence.

6 (b) Does every LTI systemThave an inverse LTI systemT 1such that{xn}=T 1T{xn}for all sequences{xn}? Why?16 Direct form I and II implementationsz 1z 1z 1z 1z 1z 1b0b1b2b3a 10 a1 a2 a3xn 1xn 2xn 3xnyn 3yn 2yn 1yn=z 1z 1z 1a 10 a1 a2 a3xnb3b0b1b2ynThe block diagram representation of the constant-coefficient differenceequation on slide 11 is called thedirect form I implementation. Thenumber of delay elements can be halved by using the commutativityof convolution to swap the two feedback loops, leading to thedirectform II implementationof the same LTI are sine waves useful?Adding together two sine waves of equal frequency, but arbitrary am-plitude and phase, results in another sine wave of the same frequency:A1 sin( t+ 1) +A2 sin( t+ 2) =A sin( t+ )withA=qA21+A22+ 2A1A2cos( 2 1)tan =A1sin 1+A2sin 2A1cos 1+A2cos 2 tA2AA1 2 1A1 sin( 1)A2 sin( 2)A2 cos( 2)A1 cos( 1)AlsoA sin( t+ ) =a sin( t) +b cos( t)withA= a2+b2andtan = of a discrete sequence{xn}with another sequence{yn}is nothing but adding together scaled and delayed copies of{xn},according to a decomposition of{yn}into a sum of impulses.

7 If{xn}was a sampled sine wave of frequencyf, so will{xn} {yn} Sine-wave sequences form a family of discrete sequences that isclosed under convolution with arbitrary waves are orthogonal to each otherZ sin( 1t+ 1) sin( 2t+ 2) dt= 0 16= 2 1 2= (2k+ 1) (k Z)and therefore can be used to form an orthogonal function basis phasorsAmplitude and phase are two distinct characteristics of a sine functionthat are inconvenient to keep separate functions (and discrete sequences) of the formA ej t+ =A [cos( t+ ) + j sin( t+ )](wherej2= 1) are able to represent both amplitude and phase inone single algebraic to complex multiplication, we can also incorporate in onesinglefactor both a multiplicative change of amplitude and an additivechangeof phase of such a function. This makes discrete sequences of the formxn= ej neigensequenceswith respect to an LTI systemT, because for each ,there is a complex number (eigenvalue)H( )such thatT{xn}=H( ) {xn}20 Recall: Fourier transformThe Fourier integral transform and its inverse are defined asF{g(t)}( ) =G( ) = Z g(t) e j tdtF 1{G( )}(t) =g(t) = Z G( ) e j td where and are constants chosen such that = 1/(2 ).

8 Many equivalent forms of the Fourier transform are used in the literature, and there is no strongconsensus of whether the forward transform usese j tand the backwards transformej t, orvice versa. Some authors set = 1and = 1/(2 ), to keep the convolution theorem free of aconstant prefactor, others use = = 1/ 2 , in the interest of substitution = 2 fleads to a form without prefactors:F{h(t)}(f) =H(f) =Z h(t) e 2 jf tdtF 1{H(f)}(t) =h(t) =Z H(f) e 2 jf tdf21 Another notation is in the continuous caseF{h(t)}( ) =H(ej ) =Z h(t) e j tdtF 1{H(ej )}(t) =h(t) =12 Z H(ej ) ej td and in the discrete-sequence caseF{hn}( ) =H(ej ) = Xn= hn e j nF 1{H(ej )}(t) =hn=12 Z H(ej ) ej nd which treats the Fourier transform as a special case of thez-transformto be introduced theoremContinuous form:F{(f g)(t)}=F{f(t)} F{g(t)}F{f(t) g(t)}=F{f(t)} F{g(t)}Discrete form:{xn} {yn}={zn} X(ej ) Y(ej ) =Z(ej )Convolution in the time domain is equivalent to (complex) scalar mul-tiplication in the frequency domain, and convolution in the frequencydomain corresponds to scalar multiplication in the time s delta functionThe continuous equivalent of the impulse sequence{ n}is known asDirac s delta function (x).

9 It is a generalized function, defined suchthat (x) = 0,x6= 0 , x= 0Z (x) dx= 1and can be thought of as the limit of function sequences such as (x) = limn 0,|x| 1/nn/2,|x|<1/nor (x) = limn n e n2x2 The delta function is mathematically speaking not a function, but adistribution, that is anexpression that is only defined when properties of Dirac s delta function:Z f(x) (x a) dx=f(a)Z e 2 jf tdf= (t)12 Z e j td = (t)Fourier transform:F{ (t)}( ) =Z (t) e j tdt= e0= 1F 1{1}(t) =12 Z 1 ej td = (t)25 Sine and cosine in the frequency domaincos(2 f t) =12e2 jf t+12e 2 jf tsin(2 f t) =12je2 jf t 12je 2 jf t ff ffAs any real-valued signalx(t)can be represented as a combinationof sine and cosine functions, the spectrum of any real-valued signalwill show the symmetryX(ej ) = [X(e j )] , where denotes thecomplex conjugate ( , negated imaginary part).26 Sampling using a Dirac combThe loss of information in the sampling process that converts a con-tinuous functionx(t)into a discrete sequence{xn}defined byxn=x(ts n) =x(n/fs)can be modelled through multiplyingx(t)by a comb of Dirac impulsess(t) = Xn= (t ts n)to obtain the sampled function x(t) =x(t) s(t)The function x(t)now contains exactly the same information as thediscrete sequence{xn}, but is still in a form that can be analysed usingthe Fourier transform on continuous Fourier transform of a Dirac combs(t) = Xn= (t ts n)is another Dirac combS(f) =F( Xn= (t tsn))(f) = Z Xn= (t tsn) e2 jf tdt=1ts Xn= f nts.

10 Tss(t)S(f)fs 2ts ts2ts 2fs fs2fs00ft28 Sampling, aliasing and Nyquist limit0i fs ffA wavecos(2 tf)sampled at frequencyfscannot be distinguishedfromcos(2 t(kfs f))for anyk Z, therefore ensure|f|< view of samplingffs 2fs fs02fsffs 2fs fs02fsf fs0f0fsWithout anti-aliasing filterWith anti-aliasing filterX(f) X(f)X(f) X(f)30 Exercise 4 Generate a one second long Gaussian noise sequence{rn}(usingMATLAB functionrandn) with a sampling rate of 300 Hz. Use thefir1(50, 45/150)function to design a finite impulse re-sponse low-pass filter with a cut-off frequency of 45 Hz. Use thefiltfiltfunction in order to apply that filter to the generated noisesignal, resulting in the filtered noise Signal {xn}. Then sample{xn}at 100 Hz by setting all but every third samplevalue to zero, resulting in sequence{yn}. Generate another lowpass filter with a cut-off frequency of 50 Hz andapply it to{yn}, in order to interpolate the reconstructed filterednoise Signal {zn}. Multiply the result by three, to compensate theenergy lost during sampling.