Digital Signal Processing lecture notes

lecture notes on advanced digital signal processing and advanced digital communications systems and signal processing techniques pdf free
JuliyaMadenta Profile Pic
Published Date:15-07-2017
Your Website URL(Optional)
Signal processing Signals may have to be transformed in order to Digital Signal Processing → amplify or filter out embedded information → detect patterns Markus Kuhn → prepare the signal to survive a transmission channel → prevent interference with other signals sharing a medium → undo distortions contributed by a transmission channel → compensate for sensor deficiencies → find information encoded in a different domain Computer Laboratory To do so, we also need → methods to measure, characterise, model and simulate trans- mission channels → mathematical tools that split common channels and transfor- Lent 2009 – Part II mations into easily manipulated building blocks 3 Signals Analog electronics → flow of information Passive networks (resistors, capacitors, R inductances, crystals, SAW filters), → measured quantity that varies with time (or position) non-linear elements (diodes, ...), U L C U in out (roughly) linear operational amplifiers → electrical signal received from a transducer Advantages: (microphone, thermometer, accelerometer, antenna, etc.) • passive networks are highly linear U in → electrical signal that controls a process over a very large dynamic range U in and large bandwidths Continuous-time signals: voltage, current, temperature, speed, ... U out • analog signal-processing circuits √ 0 1/ LC ω (= 2πf) t require little or no power Discrete-time signals: daily minimum/maximum temperature, Z t lap intervals in races, sampled continuous signals, ... • analog circuits cause little addi- U −U 1 dU in out out = U dτ +C out tional interference R L dt −∞ Electronics (unlike optics) can only deal easily with time-dependent signals, therefore spatial signals, such as images, are typically first converted into a time signal with a scanning process (TV, fax, etc.). 2 4 UoutDigital signal processing Syllabus Analog/digital and digital/analog converter, CPU, DSP, ASIC, FPGA. Signals and systems. Discrete sequences and systems, their types and proper- Advantages: ties. Linear time-invariant systems, convolution. Harmonic phasors are the eigen functions of linear time-invariant systems. Review of complex arithmetic. Some → noise is easy to control after initial quantization examples from electronics, optics and acoustics. MATLAB. Use of MATLAB on PWF machines to perform numerical experiments → highly linear (within limited dynamic range) and visualise the results in homework exercises. → complex algorithms fit into a single chip Fourier transform. Harmonic phasors as orthogonal base functions. Forms of the Fourier transform, convolution theorem, Dirac’s delta function, impulse combs in → flexibility, parameters can easily be varied in software the time and frequency domain. → digital processing is insensitive to component tolerances, aging, Discrete sequences and spectra. Periodic sampling of continuous signals, pe- environmental conditions, electromagnetic interference riodic signals, aliasing, sampling and reconstruction of low-pass and band-pass signals, spectral inversion. But: Discrete Fourier transform. Continuous versus discrete Fourier transform, sym- metry, linearity, review of the FFT, real-valued FFT. → discrete-time processing artifacts (aliasing) Spectralestimation. Leakageandscallopingphenomena,windowing,zeropadding. → can require significantly more power (battery, cooling) → digital clock and switching cause interference 5 7 Finite and infinite impulse-response filters. Properties of filters, implementa- Typical DSP applications tion forms, window-based FIR design, use of frequency-inversion to obtain high- pass filters, use of modulation to obtain band-pass filters, FFT-based convolution, → communication systems → astronomy polynomial representation, z-transform, zeros and poles, use of analog IIR design modulation/demodulation, channel VLBI, speckle interferometry techniques (Butterworth, Chebyshev I/II, elliptic filters). equalization, echo cancellation Random sequences and noise. Random variables, stationary processes, autocor- → experimental physics relation, crosscorrelation, deterministic crosscorrelation sequences, filtered random → consumer electronics sensor-data evaluation sequences, white noise, exponential averaging. perceptual coding of audio and video on DVDs, speech synthesis, speech Correlation coding. Random vectors, dependence versus correlation, covariance, recognition → aviation decorrelation, matrix diagonalisation, eigen decomposition, Karhunen-Lo`eve trans- form, principal/independent component analysis. Relation to orthogonal transform radar, radio navigation → music coding using fixed basis vectors, such as DCT. synthetic instruments, audio effects, Lossy versus lossless compression. What information is discarded by human noise reduction → security senses and can be eliminated by encoders? Perceptual scales, masking, spatial steganography,digitalwatermarking, biometric identification, surveillance resolution, colour coordinates, some demonstration experiments. → medical diagnostics systems, signals intelligence, elec- magnetic-resonance and ultrasonic Quantization, image and audio coding standards. A/μ-law coding, delta cod- tronic warfare imaging, computer tomography, ing, JPEG photographic still-image compression, motion compensation, MPEG ECG, EEG, MEG, AED, audiology video encoding, MPEG audio encoding. → engineering Note: The last three lectures on audio-visual coding were previously part of the course “Informa- → geophysics control systems, feature extraction tion Theory and Coding”. A brief introduction to MATLAB was given in “Unix Tools”. seismology, oil exploration for pattern recognition 6 8Objectives Sequences and systems ∞ By the end of the course, you should be able to A discrete sequence x is a sequence of numbers n n=−∞ → apply basic properties of time-invariant linear systems ...,x ,x ,x ,x ,x ,... −2 −1 0 1 2 → understand sampling, aliasing, convolution, filtering, the pitfalls of spectral estimation wherex denotesthen-thnumberinthesequence(n∈Z). Adiscrete n → explain the above in time and frequency domain representations sequence maps integer numbers onto real (or complex) numbers. ∞ We normally abbreviate x to x , or to x if the running index is not obvious. n n n n n=−∞ → use filter-design software The notation is not well standardized. Some authors write xn instead of x , others x(n). n → visualise and discuss digital filters in the z-domain Where a discrete sequencex samples a continuous functionx(t) as n → use the FFT for convolution, deconvolution, filtering x =x(t ·n) =x(n/f ), n s s → implement,applyandevaluatesimpleDSPapplicationsinMATLAB → applytransformsthatreducecorrelationbetweenseveralsignalsources we call t the sampling period and f = 1/t the sampling frequency. s s s → understand and explain limits in human perception that are ex- A discrete system T receives as input a sequencex and transforms n ploited by lossy compression techniques it into an output sequencey =Tx : n n → provide a good overview of the principles and characteristics of sev- discrete eral widely-used compression techniques and standards for audio- ...,x ,x ,x ,x ,... ...,y ,y ,y ,y ,... 2 1 0 −1 2 1 0 −1 system T visual signals 9 11 Textbooks Properties of sequences A sequencex is n → R.G. Lyons: Understanding digital signal processing. Prentice- ∞ Hall, 2004. (£45) X absolutely summable ⇔ x ∞ n → A.V. Oppenheim, R.W. Schafer: Discrete-time signal process- n=−∞ ∞ ing. 2nd ed., Prentice-Hall, 1999. (£47) X 2 square summable ⇔ x ∞ n → J. Stein: Digital signal processing – a computer science per- n=−∞ spective. Wiley, 2000. (£74) periodic ⇔ ∃k 0 :∀n∈Z :x =x n n+k → S.W. Smith: Digital signal processing – a practical guide for A square-summable sequence is also called an energy signal, and engineers and scientists. Newness, 2003. (£40) ∞ X 2 x n → K. Steiglitz: A digital signal processing primer – with appli- n=−∞ cations to digital audio and computer music. Addison-Wesley, is its energy. This terminology reflects that if U is a voltage supplied to a load 1996. (£40) R 2 resistorR, thenP =UI =U /Risthepowerconsumed, and P(t)dttheenergy. So even where we drop physical units (e.g., volts) for simplicity in calculations, it → Sanjit K. Mitra: Digital signal processing – a computer-based is still customary to refer to the squared values of a sequence as power and to its approach. McGraw-Hill, 2002. (£38) sum or integral over time as energy. 10 126 A non-square-summable sequence is a power signal if its average power Examples: The accumulator system k X 1 2 lim x n n X k→∞ 1+2k n=−k y = x n k k=−∞ exists. is a causal, linear, time-invariant system with memory, as are the back- ward difference system Special sequences y =x −x , n n n−1 Unit-step sequence:  0, n 0 the M-point moving average system u = n 1, n≥ 0 M−1 X 1 x +···+x +x n−M+1 n−1 n y = x = n n−k Impulse sequence: M M k=0  1, n = 0 and the exponential averaging system δ = n 0, n = 0 ∞ X k = u −u y =α·x +(1−α)·y =α (1−α) ·x . n n−1 n n n−1 n−k k=0 13 15 Examples for time-invariant non-linear memory-less systems: Types of discrete systems A causal system cannot look into the future: 2 y =x , y = log x , y = maxmin⌊256x ⌋,255,0 n n n n n 2 n y =f(x ,x ,x ,...) n n n−1 n−2 Examples for linear but not time-invariant systems: A memory-less system depends only on the current input value:  x , n≥ 0 n y = =x ·u n n n y =f(x ) n n 0, n 0 y = x n ⌊n/4⌋ A delay system shifts a sequence in time: ωjn y = x ·ℜ(e ) n n y =x n n−d Examples for linear time-invariant non-causal systems: T is a time-invariant system if for any d 1 y = (x +x ) n n−1 n+1 y =Tx ⇐⇒ y =Tx . n n n−d n−d 2 9 X ′ sin(πkω) T is a linear system if for any pair of sequencesx andx n n y = x · ·0.5+0.5·cos(πk/10) n n+k πkω ′ ′ k=−9 Ta·x +b·x =a·Tx +b·Tx . n n n n 14 16Constant-coefficient difference equations Convolution Of particular practical interest are causal linear time-invariant systems All linear time-invariant (LTI) systems can be represented in the form of the form x b y n 0 n ∞ X N y = a ·x X n k n−k −1 z y =b ·x − a ·y n 0 n k n−k k=−∞ −a 1 k=1 y n−1 wherea is a suitably chosen sequence of coefficients. k −1 z This operation over sequences is called convolution and defined as −a 2 Block diagram representation y n−2 ∞ of sequence operations: X −1 ′ z p ∗q =r ⇐⇒ ∀n∈Z :r = p ·q . n n n n k n−k x n −a 3 k=−∞ ′ y n−3 x x +x n n n Addition: If y = a ∗x is a representation of an LTI system T, with n n n Multiplication x a ax n n The a and b are k m y =Tx , then we call the sequence a the impulse response n n n by constant: constant coefficients. of T, becausea =Tδ . n n x x n n−1 −1 Delay: z 17 19 or Convolution examples x x x x n n−1 n−2 n−3 −1 −1 −1 z z z M A B C D X y = b ·x b b b b 0 1 2 3 n m n−m m=0 y n or the combination of both: −1 A∗B A∗C x b a y E F n 0 n 0 −1 −1 z z b −a 1 1 N M X X x y n−1 n−1 a ·y = b ·x k n−k m n−m −1 −1 z z k=0 m=0 C∗A A∗E D∗E A∗F b −a 2 2 x y n−2 n−2 −1 −1 z z b −a 3 3 x y n−3 n−3 The MATLAB function filter is an efficient implementation of the last variant. 18 20Properties of convolution Exercise1 Whattypeofdiscretesystem(linear/non-linear,time-invariant/ non-time-invariant, causal/non-causal, causal, memory-less, etc.) is: For arbitrary sequencesp ,q ,r and scalars a, b: n n n → Convolution is associative 3x +x n−1 n−2 (a) y =x n n (e) y = n x n−3 (p ∗q )∗r =p ∗(q ∗r ) n n n n n n n/14 (f) y =x ·e (b) y =−x +2x −x n n n n−1 n n+1 → Convolution is commutative 8 Y (g) y =x ·u p ∗q =q ∗p n n n n n n n (c) y = x n n−i ∞ X i=0 → Convolution is linear (h) y = x ·δ n i i−n+2 1 (d) y = (x +x ) n 2n 2n+1 2 i=−∞ p ∗a·q +b·r =a·(p ∗q )+b·(p ∗r ) n n n n n n n Exercise 2 → The impulse sequence (slide 13) is neutral under convolution Prove that convolution is (a) commutative and (b) associative. p ∗δ =δ ∗p =p n n n n n → Sequence shifting is equivalent to convolving with a shifted impulse p =p ∗δ n−d n n−d 21 23 Can all LTI systems be represented by convolution? Exercise 3 A finite-length sequence is non-zero only at a finite number of positions. If m and n are the first and last non-zero positions, respectively, Any sequencex can be decomposed into a weighted sum of shifted n then we call n−m+1 the length of that sequence. What maximum length impulse sequences: can the result of convolving two sequences of length k and l have? ∞ X Exercise 4 The length-3 sequence a =−3, a = 2, a = 1 is convolved x = x ·δ n k n−k 0 1 2 with a second sequence b of length 5. k=−∞ n (a) Write down this linear operation as a matrix multiplication involving a (∗) (∗∗) 5 Let’s see what happens if we apply a linear time-invariant system matrix A, a vector b∈R , and a result vectorc. T to such a decomposed sequence: (b) Use MATLAB to multiply your matrix by the vector b = (1,0,0,2,2) and compare the result with that of using the conv function. ∞ ∞ X X (∗) Tx = T x ·δ = x ·Tδ (c) Use the MATLAB facilities for solving systems of linear equations to n k n−k k n−k undo the above convolution step. k=−∞ k=−∞ ∞ ∞ X X (∗∗) Exercise 5 (a) Find a pair of sequences a and b , where each one = x ·δ ∗Tδ = x ·δ ∗Tδ n n k n−k n k n−k n containsatleastthreedifferentvaluesandwheretheconvolutiona ∗b k=−∞ k=−∞ n n results in an all-zero sequence. = x ∗Tδ q.e.d. n n −1 (b) Does every LTI system T have an inverse LTI system T such that −1 ⇒ The impulse response Tδ fully characterizes an LTI system. n x =T Tx for all sequences x ? Why? n n n 22 24Direct form I and II implementations Convolution: electronics example −1 −1 x b a y x a b y n 0 n n 0 n 0 0 U in R U U in in √ 2 −1 −1 −1 z z z U C U in out b −a −a b 1 1 1 1 x y n−1 n−1 U out = −1 −1 −1 0 z z z 1/RC ω (= 2πf) t 0 b −a −a b 2 2 2 2 x y n−2 n−2 Any passive network (R,L,C) convolves its input voltage U with an in −1 −1 −1 impulse response function h, leading to U =U ∗h, that is z z z out in b −a −a b 3 3 3 3 Z ∞ x y n−3 n−3 U (t) = U (t−τ)·h(τ)·dτ out in Theblockdiagramrepresentationoftheconstant-coefficientdifference −∞ equation on slide 18 is called the direct form I implementation. In this example: The number of delay elements can be halved by using the commuta-  −t 1 U −U dU RC in out out ·e , t≥ 0 tivity of convolution to swap the two feedback loops, leading to the RC =C· , h(t) = R dt 0, t 0 direct form II implementation of the same LTI system. These two forms are only equivalent with ideal arithmetic (no rounding errors and range limits). 25 27 Convolution: optics example Why are sine waves useful? If a projective lens is out of focus, the blurred image is equal to the 1) Adding together sine waves of equal frequency, but arbitrary ampli- original image convolved with the aperture shape (e.g., a filled circle): tude and phase, results in another sine wave of the same frequency: A ·sin(ωt+ϕ )+A ·sin(ωt+ϕ ) = A·sin(ωt+ϕ) 1 1 2 2 with q ∗ = 2 2 A = A +A +2A A cos(ϕ −ϕ ) 1 2 2 1 1 2 A sinϕ +A sinϕ 1 1 2 2 tanϕ = A A ·sin(ϕ ) 2 2 2 A cosϕ +A cosϕ 1 1 2 2 A as Point-spread function h (disk, r = ): 2f image plane focal plane ϕ 2 Sine waves of any phase can be  A A ·cos(ϕ ) 2 2 1 1 2 2 2 , x +y ≤r formed from sin and cos alone: 2 r π h(x,y) = 2 2 2 ωt A ·sin(ϕ ) 1 1 0, x +y r ϕ ϕ 1 a A·sin(ωt+ϕ) = A ·cos(ϕ ) 1 1 Original image I, blurred image B =I∗h, i.e. a·sin(ωt)+b·cos(ωt) ZZ ′ ′ ′ ′ ′ ′ s √ B(x,y) = I(x−x ,y−y )·h(x ,y )·dx dy b 2 2 with a =A·cos(ϕ), b =A·sin(ϕ) and A = a +b , tanϕ = . a f 26 28 U out6 Note: Convolution of a discrete sequencex with another sequence Why are complex numbers so useful? n y is nothing but adding together scaled and delayed copies ofx . n n 1) They give us all n solutions (“roots”) of equations involving poly- √ (Think ofy decomposed into a sum of impulses.) n nomials up to degree n (the “ −1 = j ” story). Ifx is a sampled sine wave of frequency f, so isx ∗y n n n 2) They give us the “great unifying theory” that combines sine and =⇒ Sine-wave sequences form a family of discrete sequences exponential functions: that is closed under convolution with arbitrary sequences.  1 jωt −jωt cos(ωt) = e +e The same applies for continuous sine waves and convolution. 2  1 jωt −jωt sin(ωt) = e −e 2) Sine waves are orthogonal to each other: 2j Z ∞ or  1 jωt+ϕ −jωt−ϕ sin(ω t+ϕ )·sin(ω t+ϕ )dt “=” 0 1 1 2 2 cos(ωt+ϕ) = e +e −∞ 2 or ⇐⇒ ω =ω ∨ ϕ −ϕ = (2k+1)π/2 (k∈Z) 1 2 1 2 jωn+ϕ jω n jϕ cos(ωn+ϕ) = ℜ(e ) = ℜ(e ) ·e They can be used to form an orthogonal function basis for a transform. jωn+ϕ jω n jϕ sin(ωn+ϕ) = ℑ(e ) = ℑ(e ) ·e The term “orthogonal” is used here in the context of an (infinitely dimensional) vector space, where the “vectors” are functions of the form f :R→R (or f :R→C) and the scalar product R ∞ 2 is defined as f ·g = f(t)·g(t)dt. Notation: ℜ(a+ jb) :=a and ℑ(a+ jb) :=b where j =−1 and a,b∈R. −∞ 29 31 We can now represent sine waves as projections of a rotating complex Why are exponential functions useful? vector. This allows us to represent sine-wave sequences as exponential Adding together two exponential functions with the same base z, but jω sequences with basis e . different scale factor and offset, results in another exponential function Aphaseshiftinsuchasequencecorrespondstoarotationofacomplex with the same base: vector. t+ϕ t+ϕ t ϕ t ϕ 1 2 1 2 A ·z +A ·z = A ·z ·z +A ·z ·z 1 2 1 2 3)Complexmultiplicationallowsustomodifytheamplitude andphase ϕ ϕ t t 1 2 = (A ·z +A ·z )·z =A·z 1 2 of a complex rotating vector using a single operation and value. Likewise, if we convolve a sequencex of values Rotation of a 2D vector in (x,y)-form is notationally slightly messy, n 2 but fortunately j =−1 does exactly what is required here: −3 −2 −1 2 3 ...,z ,z ,z ,1,z,z ,z ,...       n n x x −y x 3 2 2 1 x =z with an arbitrary sequenceh , we gety =z ∗h , n n n n = · (x ,y ) 3 3 y y x y 3 2 2 1 ∞ ∞ ∞   X X X x x −y y (−y ,x ) n−k n −k n 1 2 1 2 2 2 y = x ·h = z ·h =z · z ·h =z ·H(z) = n n−k k k k x y +x y 1 2 2 1 (x ,y ) 2 2 k=−∞ k=−∞ k=−∞ (x ,y ) 1 1 z =x + jy , z =x + jy 1 1 1 2 2 2 where H(z) is independent of n. z ·z =x x −y y + j(x y +x y ) 1 2 1 2 1 2 1 2 2 1 Exponential sequences are closed under convolution with arbitrary sequences. The same applies in the continuous case. 30 32Complex phasors Another notation is in the continuous case Z ∞ Amplitude and phase are two distinct characteristics of a sine function jω −jωt Fh(t)(ω) = H(e ) = h(t) ·e dt that are inconvenient to keep separate notationally. −∞ Complex functions (and discrete sequences) of the form Z ∞ 1 −1 jω jω jωt j(ωt+ϕ) F H(e )(t) = h(t) = H(e )·e dω A·e =A·cos(ωt+ϕ)+ j·sin(ωt+ϕ) 2π −∞ 2 (where j = −1) are able to represent both amplitude and phase in and in the discrete-sequence case one single algebraic object. Thankstocomplexmultiplication, wecanalsoincorporateinonesingle ∞ X jω −jωn factorbothamultiplicativechangeofamplitudeandanadditivechange Fh (ω) = H(e ) = h ·e n n of phase of such a function. This makes discrete sequences of the form n=−∞ Z π jωn 1 x = e −1 jω jω jωn n F H(e )(t) = h = H(e )·e dω n 2π −π eigensequences with respect to an LTI system T, because for each ω, there is a complex number (eigenvalue) H(ω) such that which treats the Fourier transform as a special case of thez-transform (to be introduced shortly). Tx =H(ω)·x n n jω In the notation of slide 30, where the argument of H is the base, we would write H(e ). 33 35 Recall: Fourier transform Properties of the Fourier transform The Fourier integral transform and its inverse are defined as If Z ∞ ∓jωt x(t) •−◦ X(f) and y(t) •−◦ Y(f) Fg(t)(ω) = G(ω) = α g(t) ·e dt −∞ are pairs of functions that are mapped onto each other by the Fourier Z ∞ transform, then so are the following pairs. −1 ±jωt F G(ω)(t) = g(t) = β G(ω)·e dω −∞ Linearity: ax(t)+by(t) •−◦ aX(f)+bY(f) where α and β are constants chosen such that αβ = 1/(2π). Many equivalent forms of the Fourier transform are used in the literature. There is no strong −jωt jωt consensus on whether the forward transform uses e and the backwards transform e , or Time scaling:   vice versa. We follow here those authors who set α = 1 and β = 1/(2π), to keep the convolution √ 1 f theorem free of a constant prefactor; others prefer α =β = 1/ 2π, in the interest of symmetry. x(at) •−◦ X a a The substitution ω = 2πf leads to a form without prefactors: Z ∞ Frequency scaling: ∓2πjft Fh(t)(f) = H(f) = h(t) ·e dt   −∞ 1 t Z ∞ x •−◦ X(af) a a −1 ±2πjft F H(f)(t) = h(t) = H(f)·e df −∞ 34 366 Time shifting: Convolution theorem Continuous form: −2πjfΔt x(t−Δt) •−◦ X(f)·e F(f∗g)(t) = Ff(t)·Fg(t) Frequency shifting: Ff(t)·g(t) = Ff(t)∗Fg(t) 2πjΔft x(t)·e •−◦ X(f−Δf) Discrete form: jω jω jω x ∗y = z ⇐⇒ X(e )·Y(e ) = Z(e ) n n n Parseval’s theorem (total power): Z Z ∞ ∞ Convolution in the time domain is equivalent to (complex) scalar mul- 2 2 x(t) dt = X(f) df tiplication in the frequency domain. −∞ −∞ Convolution in the frequency domain corresponds to scalar multiplica- tion in the time domain. R R R R −jωr −jωr Proof: z(r) = x(s)y(r−s)ds ⇐⇒ z(r)e dr = x(s)y(r−s)e dsdr = s r r s R R R R t:=r−s −jωr −jωs −jω(r−s) x(s) y(r−s)e drds = x(s)e y(r−s)e drds = s r s r R R R R P R −jωs −jωt −jωs −jωt x(s)e y(t)e dtds = x(s)e ds· y(t)e dt. (Same for instead of .) s t s t 37 39 Fourier transform example: rect and sinc Dirac’s delta function The continuous equivalent of the impulse sequence δ is known as n The Fourier transform of the “rectangular function” Dirac’s delta function δ(x). It is a generalized function, defined such  1  1 ift that  2 1 1 rect(t) = ift =  2 2  1  0, x = 0 0 otherwise δ(x) = ∞, x = 0 Z ∞ is the “(normalized) sinc function” δ(x)dx = 1 1 Z −∞ 2 sinπf 0 x −2πjft Frect(t)(f) = e dt = = sinc(f) 1 πf and can be thought of as the limit of function sequences such as − 2  0, x≥ 1/n and vice versa δ(x) = lim n→∞ n/2, x 1/n Fsinc(t)(f) = rect(f). or Some noteworthy properties of these functions: n 2 2 −n x R R δ(x) = lim √ e ∞ ∞ • sinc(t)dt = 1 = rect(t)dt n→∞ π −∞ −∞ • sinc(0) = 1 = rect(0) The delta function is mathematically speaking not a function, but a distribution, that is an • ∀n∈Z\0 : sinc(k) = 0 expression that is only defined when integrated. 38 40Some properties of Dirac’s delta function: Fourier transform symmetries Z ∞ We call a function x(t) f(x)δ(x−a)dx = f(a) −∞ odd if x(−t) = −x(t) Z ∞ even if x(−t) = x(t) ±2πjft e df = δ(t) −∞ ∗ ∗ and· is the complex conjugate, such that (a+ jb) = (a− jb). Z ∞ 1 ±jωt Then e dω = δ(t) 2π −∞ ∗ x(t) is real ⇔ X(−f) = X(f) ∗ Fourier transform: x(t) is imaginary ⇔ X(−f) =−X(f) Z x(t) is even ⇔ X(f) is even ∞ −jωt 0 Fδ(t)(ω) = δ(t) ·e dt = e = 1 x(t) is odd ⇔ X(f) is odd −∞ x(t) is real and even ⇔ X(f) is real and even Z ∞ x(t) is real and odd ⇔ X(f) is imaginary and odd 1 −1 jωt F 1(t) = 1 ·e dω = δ(t) x(t) is imaginary and even ⇔ X(f) is imaginary and even 2π −∞ x(t) is imaginary and odd ⇔ X(f) is real and odd 41 43 Sine and cosine in the frequency domain Example: amplitude modulation 1 1 1 1 Communication channels usually permit only the use of a given fre- 2πjf t −2πjf t 2πjf t −2πjf t 0 0 0 0 cos(2πf t) = e + e sin(2πf t) = e − e 0 0 quency interval, such as 300–3400 Hz for the analog phone network or 2 2 2j 2j 590–598 MHz for TV channel 36. Modulation with a carrier frequency 1 1 Fcos(2πf t)(f) = δ(f−f )+ δ(f +f ) 0 0 0 f shifts the spectrum of a signal x(t) into the desired band. c 2 2 Amplitude modulation (AM): j j Fsin(2πf t)(f) = − δ(f−f )+ δ(f +f ) 0 0 0 2 2 ℜ ℜ y(t) =A·cos(2πtf )·x(t) c 1 1 2 2 X(f) Y(f) ℑ ℑ 1 1 j j 2 2 ∗ = −f f f −f f f −f 0 f f −f f f −f 0 f f 0 0 0 0 l l c c c c The spectrum of the baseband signal in the interval −f f f is l l As any real-valued signal x(t) can be represented as a combination of sine and cosine functions, jω −jω ∗ ∗ shifted by the modulation to the intervals±f −f f ±f +f. c l c l the spectrum of any real-valued signal will show the symmetry X(e ) = X(e ) , where denotes the complex conjugate (i.e., negated imaginary part). How can such a signal be demodulated? 42 44Sampling using a Dirac comb Sampling and aliasing The loss of information in the sampling process that converts a con- sample tinuous function x(t) into a discrete sequencex defined by n cos(2π tf) cos(2π t(k⋅ f ± f)) s x =x(t ·n) =x(n/f ) n s s can be modelled through multiplyingx(t) by a comb of Dirac impulses ∞ X 0 s(t) =t · δ(t−t ·n) s s n=−∞ to obtain the sampled function xˆ(t) =x(t)·s(t) The function xˆ(t) now contains exactly the same information as the discrete sequencex , but is still in a form that can be analysed using Sampled at frequency f , the function cos(2πtf) cannot be distin- s n guished from cos2πt(kf ±f) for any k∈Z. the Fourier transform on continuous functions. s 45 47 The Fourier transform of a Dirac comb Frequency-domain view of sampling ∞ ∞ X X x(t) s(t) xˆ(t) 2πjnt/t s s(t) =t · δ(t−t ·n) = e s s n=−∞ n=−∞ · = is another Dirac comb ... ... ... ... ( ) ∞ X 0 t −1/f 0 1/f t −1/f 0 1/f t s s s s S(f) =F t · δ(t−tn) (f) = s s ˆ X(f) S(f) X(f) n=−∞ ∞ Z   ∞ ∞ X X n 2πjft t · δ(t−tn)e dt = δ f− . ∗ = s s t s n=−∞ n=−∞ −∞ ... ...... ... 0 f −f f f −f 0 f f s s s s s(t) S(f) Sampling a signal in the time domain corresponds in the frequency domain to convolving its spectrum with a Dirac comb. The resulting copies of the original signal spectrum in the spectrum of the sampled −2t −t 0 t 2t t −2f −f 0 f 2f f s s s s s s s s signal are called “images”. 46 48Nyquist limit and anti-aliasing filters Reconstruction of a continuous band-limited waveform If the (double-sided) bandwidth of a signal to be sampled is larger than Theidealanti-aliasingfilterforeliminatinganyfrequencycontentabove the sampling frequencyf , the images of the signal that emerge during s f/2 before sampling with a frequency of f has the Fourier transform s s sampling may overlap with the original spectrum. ( f s 1 iff 2 Such an overlap will hinder reconstruction of the original continuous H(f) = = rect(tf). s f s 0 iff signal by removing the aliasing frequencies with a reconstruction filter. 2 Therefore, it is advisable to limit the bandwidth of the input signal to This leads, after an inverse Fourier transform, to the impulse response the sampling frequencyf before sampling, using an anti-aliasing filter.   s sinπtf 1 t s In the common case of a real-valued base-band signal (with frequency h(t) =f · = ·sinc . s πtf t t s s s content down to 0 Hz), all frequencies f that occur in the signal with The original band-limited signal can be reconstructed by convolving non-zero power should be limited to the interval−f/2f f/2. s s this with the sampled signal xˆ(t), which eliminates the periodicity of The upper limit f/2 for the single-sided bandwidth of a baseband s the frequency domain introduced by the sampling process: signal is known as the “Nyquist limit”. x(t) =h(t)∗xˆ(t) Note that sampling h(t) gives the impulse function: h(t)·s(t) =δ(t). 49 51 Impulse response of ideal low-pass filter with cut-off frequency f/2: Nyquist limit and anti-aliasing filters s Without anti-aliasing filter With anti-aliasing filter single-sided Nyquist X(f) bandwidth X(f) limit = f /2 s anti-aliasing filter 0 f −f 0 f f s s double-sided bandwidth reconstruction filter ˆ ˆ X(f) X(f) 0 −2f −f 0 f 2f f −2f −f 0 f 2f f s s s s s s s s −3 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5 3 Anti-aliasing and reconstruction filters both suppress frequencies outside ff /2. s t⋅ f s 50 52Reconstruction filter example Exercise 6 Digital-to-analog converters cannot output Dirac pulses. In- stead, for each sample, they hold the output voltage (approximately) con- stant, until the next sample arrives. How can this behaviour be modeled sampled signal mathematically as a linear time-invariant system, and how does it affect the interpolation result spectrum of the output signal? scaled/shifted sin(x)/x pulses Exercise 7 Many DSP systems use “oversampling” to lessen the require- ments on the design of an analog reconstruction filter. They use (a finite approximation of) the sinc-interpolation formula to multiply the sampling frequency f of the initial sampled signal by a factor N before passing it to s the digital-to-analog converter. While this requires more CPU operations and a faster D/A converter, the requirements on the subsequently applied analog reconstruction filter are much less stringent. Explain why, and draw schematic representations of the signal spectrum before and after all the relevant signal-processing steps. Exercise 8 Similarly, explain how oversampling can be applied to lessen the requirements on the design of an analog anti-aliasing filter. 1 2 3 4 5 53 55 Reconstruction filters Band-pass signal sampling The mathematically ideal form of a reconstruction filter for suppressing Sampled signals can also be reconstructed if their spectral components aliasing frequencies interpolates the sampled signalx =x(t ·n) back n s remain entirely within the interval n·f/2 f (n+1)·f/2 for s s into the continuous waveform some n∈N. (The baseband case discussed so far is just n = 0.) ∞ X In this case, the aliasing copies of the positive and the negative frequencies will interleave instead sinπ(t−t ·n) s of overlap, and can therefore be removed again with a reconstruction filter with the impulse x(t) = x · . n π(t−t ·n) response s n=−∞ „ « sinπtf /2 2n+1 sinπt(n+1)f sinπtnf s s s h(t) =f ·cos 2πtf = (n+1)f −nf . s s s s πtf /2 4 πt(n+1)f πtnf s s s Choice of sampling frequency Due to causality and economic constraints, practical analog filters can only approx- ˆ X(f) anti-aliasing filter X(f) reconstruction filter imate such an ideal low-pass filter. Instead of a sharp transition between the “pass band” ( f /2) and the “stop band” ( f /2), they feature a “transition band” s s in which their signal attenuation gradually increases. The sampling frequency is therefore usually chosen somewhat higher than twice the highest frequency of interest in the continuous signal (e.g., 4×). On the other 5 5 −f f s s − f 0 f f −f 0 f f s s s s hand, the higher the sampling frequency, the higher are CPU, power and memory 4 4 2 2 requirements. Therefore, the choice of sampling frequency is a tradeoff between n = 2 signal quality, analog filter cost and digital subsystem expenses. 54 56Exercise 9 Reconstructing a sampled baseband signal: Spectrum of a periodic signal A signalx(t) that is periodic with frequencyf can be factored into a • Generate a one second long Gaussian noise sequence r (using p n single period x˙(t) convolved with an impulse comb p(t). This corre- MATLAB function randn) with a sampling rate of 300 Hz. sponds in the frequency domain to the multiplication of the spectrum • Use the fir1(50, 45/150) function to design a finite impulse re- of the single period with a comb of impulses spaced f apart. p sponse low-pass filter with a cut-off frequency of 45 Hz. Use the filtfilt function in order to apply that filter to the generated noise x(t) x˙(t) p(t) signal, resulting in the filtered noise signalx . n • Then sample x at 100 Hz by setting all but every third sample n = ∗ value to zero, resulting in sequence y . n ... ... ... ... • Generate another low-pass filter with a cut-off frequency of 50 Hz −1/f 0 1/f t t −1/f 0 1/f t p p p p and apply it toy , in order to interpolate the reconstructed filtered n ˙ X(f) X(f) P(f) noise signal z . Multiply the result by three, to compensate the n energy lost during sampling. = · • Plot x , y , andz . Finally compare x andz . n n n n n ... ... Why should the first filter have a lower cut-off frequency than the second? −f 0 f f f −f 0 f f p p p p 57 59 Exercise 10 Reconstructing a sampled band-pass signal: Spectrum of a sampled signal • Generate a 1 s noise sequence r , as in exercise 9, but this time n A signalx(t) that is sampled with frequencyf has a spectrum that is s use a sampling frequency of 3 kHz. periodic with a period of f . s • Apply to that a band-pass filter that attenuates frequencies outside theinterval31–44Hz,whichtheMATLABSignalProcessingToolbox x(t) s(t) xˆ(t) function cheby2(3, 30, 31 44/1500) will design for you. • Then sample the resulting signal at 30 Hz by setting all but every · = 100-th sample value to zero. ... ... ... ... • Generate with cheby2(3, 20, 30 45/1500) another band-pass 0 t −1/f 0 1/f t −1/f 0 1/f t s s s s filter for the interval 30–45 Hz and apply it to the above 30-Hz- ˆ X(f) S(f) X(f) sampled signal, to reconstruct the original signal. (You’ll have to multiply it by 100, to compensate the energy lost during sampling.) ∗ = • Plot all the produced sequences and compare the original band-pass ... ... ... ... signal and that reconstructed after being sampled at 30 Hz. 0 f −f f f −f 0 f f s s s s Why does the reconstructed waveform differ much more from the original if you reduce the cut-off frequencies of both band-pass filters by 5 Hz? 58 60Continuous vs discrete Fourier transform Discrete Fourier Transform visualized       • Sampling a continuous signal makes its spectrum periodic x X 0 0       x      X  • A periodic signal has a sampled spectrum 1 1             x X 2 2       We sample a signalx(t) withf , gettingxˆ(t). We taken consecutive s       x X 3 3       · = samplesofxˆ(t)andrepeattheseperiodically, gettinganewsignalx¨(t)       x X 4 4       ¨ with period n/f . Its spectrum X(f) is sampled (i.e., has non-zero s       x X 5 5       value) at frequency intervalsf/n and repeats itself with a periodf . s s       x X    6   6  ¨ Now both x¨(t) and its spectrum X(f) are finite vectors of length n. x X 7 7 ¨ x¨(t) X(f) Then-point DFT of a signalx sampled at frequencyf contains in i s the elementsX toX of the resulting frequency-domain vector the 0 n/2 frequency components 0, f/n, 2f/n, 3f/n, ..., f/2, and contains s s s s in X downto X the corresponding negative frequencies. Note n−1 n/2 ... ... ... ... that for a real-valued input vector, bothX andX will be real, too. 0 n/2 −1 −1 −n/f f 0f n/f t −f −f /n 0 f /n f f s s s s s s s s Why is there no phase information recovered at f /2? s 61 63 Discrete Fourier Transform (DFT) Inverse DFT visualized n−1 n−1 X X ik 1 ik −2πj 2πj n n X = x ·e x = · X ·e k i k i       n i=0 i=0 X x 0 0       The n-point DFT multiplies a vector with an n×n matrix   X   x   1 1         1 1 1 1 ··· 1    X   x  2 2 1 2 3 n−1       −2πj −2πj −2πj −2πj   n n n n 1 e e e ··· e       1 X x   3 3 2(n−1) 2 4 6         −2πj −2πj −2πj −2πj · · = n n n n 1 e e e ··· e         X x 8 4 4  3(n−1)  3 6 9       F = n −2πj −2πj −2πj −2πj  n n n n  1 e e e ··· e         X x    5   5   . . . . . .  . . . . . .        .  . . . . . X x    6   6  n−1 2(n−1) 3(n−1) (n−1)(n−1) −2πj −2πj −2πj −2πj n n n n 1 e e e ··· e X x 7 7         x X X x 0 0 0 0         x X X x 1 1 1 1         1         ∗ x X X x 2 2 2 2 F ·  = , ·F ·  =  n n  .   .  n  .   .  . . . .         . . . . x X X x n−1 n−1 n−1 n−1 62 64Fast Fourier Transform (FFT) Fast complex multiplication Calculating the product of two complex numbers as n−1 X  ik n−1 −2πj n F x = x ·e n i i i=0 k (a+ jb)·(c+ jd) = (ac−bd)+ j(ad+bc) i=0 n n −1 −1 2 2 X X ik ik k involves four (real-valued) multiplications and two additions. −2πj −2πj −2πj n/2 n/2 n = x ·e + e x ·e 2i 2i+1 The alternative calculation i=0 i=0      n n k −1 −1  −2πj n 2 2  n n n F x + e · F x , k α = a(c+d)  2i 2i+1 i=0 i=0 2 2 2 k k = (a+ jb)·(c+ jd) = (α−β)+ j(α+γ) with β = d(a+b)     n n k  −1 −1 2 −2πj 2 n  n n n  F x + e · F x , k≥ 2i 2i+1 γ = c(b−a) i=0 i=0 2 n n 2 2 k− k− 2 2 provides the same result with three multiplications and five additions. The DFT over n-element vectors can be reduced to two DFTs over ThelattermayperformfasteronCPUswheremultiplicationstakethree n/2-element vectors plusn multiplications andn additions, leading to or more times longer than additions. log n rounds and nlog n additions and multiplications overall, com- 2 2 This trick is most helpful on simpler microcontrollers. Specialized signal-processing CPUs (DSPs) 2 pared to n for the equivalent matrix multiplication. feature 1-clock-cycle multipliers. High-end desktop processors use pipelined multipliers that stall A high-performance FFT implementation in C with many processor-specific optimizations and where operations depend on each other. support for non-power-of-2 sizes is available at 65 67 Efficient real-valued FFT FFT-based convolution m−1 n−1 The symmetry properties of the Fourier transform applied to the discrete Calculatingtheconvolutionoftwofinitesequencesx andy i i i=0 i=0 n−1 n−1 Fourier transformX =F x have the form i n i i=0 i=0 of lengths m and n via ∗ ∀i :x = ℜ(x ) ⇐⇒ ∀i :X = X i i n−i i minm−1,i X ∗ ∀i :x = j·ℑ(x ) ⇐⇒ ∀i :X =−X i i n−i z = x ·y , 0≤im+n−1 i i j i−j j=max0,i−(n−1) These two symmetries, combined with the linearity of the DFT, allows us to calculate two real-valued n-point DFTs takes mn multiplications. ′ n−1 ′ n−1 ′′ n−1 ′′ n−1 Can we apply the FFT and the convolution theorem to calculate the X =F x X =F x i n i i n i i=0 i=0 i=0 i=0 convolution faster, in just O(mlogm+nlogn) multiplications? simultaneously in a single complex-valued n-point DFT, by composing its −1 input as z =F (Fx·Fy) i i i ′ ′′ x =x + j·x i i i There is obviously no problem if this condition is fulfilled: and decomposing its output as x andy are periodic, with equal period lengths 1 1 i i ′ ∗ ′′ ∗ X = (X +X ) X = (X −X ) i i i n−i i n−i 2 2 Inthiscase, thefactthattheDFTinterpretsitsinputasasingleperiod of a periodic signal will do exactly what is needed, and the FFT and Tooptimizethecalculationofasinglereal-valuedFFT,usethistricktocalculatethetwohalf-size inverse FFT can be applied directly as above. real-value FFTs that occur in the first round. 66 68In the general case, measures have to be taken to prevent a wrap-over: Each block is zero-padded at both ends and then convolved as before: −1 F F(A)⋅F(B) A B ∗ ∗ ∗ −1 A’ B’ F F(A’)⋅F(B’) = = = The regions originally added as zero padding are, after convolution, aligned Bothsequencesarepaddedwithzerovaluestoalengthofatleastm+n−1. to overlap with the unpadded ends of their respective neighbour blocks. Thisensuresthatthestartandendoftheresultingsequencedonotoverlap. The overlapping parts of the blocks are then added together. 69 71 Zero padding is usually applied to extend both sequence lengths to the Deconvolution ⌈log (m+n−1)⌉ 2 next higher power of two (2 ), which facilitates the FFT. A signal u(t) was distorted by convolution with a known impulse re- With a causal sequence, simply append the padding zeros at the end. sponseh(t)(e.g., throughatransmissionchannelorasensorproblem). The “smeared” result s(t) was recorded. With a non-causal sequence, values with a negative index number are wrapped around the DFT block boundaries and appear at the right Can we undo the damage and restore (or at least estimate) u(t)? end. In this case, zero-padding is applied in the center of the block, between the last and first element of the sequence. Thanks to the periodic nature of the DFT, zero padding at both ends ∗ = has the same effect as padding only at one end. IfbothsequencescanbeloadedentirelyintoRAM,theFFTcanbeap- plied to them in one step. However, one of the sequences might be too large for that. It could also be a realtime waveform (e.g., a telephone signal) that cannot be delayed until the end of the transmission. In such cases, the sequence has to be split into shorter blocks that are ∗ = separately convolved and then added together with a suitable overlap. 70 72The convolution theorem turns the problem into one of multiplication: Spectral estimation Z s(t) = u(t−τ)·h(τ)·dτ Sine wave 4×f /32 s Discrete Fourier Transform 1 s = u∗h 15 Fs = Fu·Fh 10 0 5 Fu = Fs/Fh −1 −1 0 u = F Fs/Fh 0 10 20 30 0 10 20 30 In practice, we also record some noise n(t) (quantization, etc.): Sine wave 4.61×f /32 s Discrete Fourier Transform Z 1 15 c(t) =s(t)+n(t) = u(t−τ)·h(τ)·dτ +n(t) 10 0 Problem – At frequencies f where Fh(f) approaches zero, the 5 noise will be amplified (potentially enormously) during deconvolution: −1 0 −1 −1 u˜ =F Fc/Fh =u+F Fn/Fh 0 10 20 30 0 10 20 30 73 75 Typical workarounds: We introduced the DFT as a special case of the continuous Fourier transform, where the input is sampled and periodic. → ModifytheFouriertransformoftheimpulseresponse,suchthat If the input is sampled, but not periodic, the DFT can still be used Fh(f)ǫ for some experimentally chosen threshold ǫ. to calculate an approximation of the Fourier transform of the original → If estimates of the signal spectrum Fs(f) and the noise continuous signal. However, there are two effects to consider. They spectrum Fn(f) can be obtained, then we can apply the are particularly visible when analysing pure sine waves. “Wiener filter” (“optimal filter”) Sine waves whose frequency is a multiple of the base frequency (f/n) s 2 Fs(f) of the DFT are identical to their periodic extension beyond the size W(f) = 2 2 Fs(f) +Fn(f) of the DFT. They are, therefore, represented exactly by a single sharp before deconvolution: peak in the DFT. All their energy falls into one single frequency “bin” −1 in the DFT result. u˜ =F W ·Fc/Fh Sine waves with other frequencies, which do not match exactly one of Exercise 11 Use MATLAB to deconvolve the blurred stars from slide 26. the output frequency bins of the DFT, are still represented by a peak The files stars-blurred.png with the blurred-stars image and stars-psf.png with the impulse at the output bin that represents the nearest integer multiple of the response (point-spread function) are available on the course-material web page. You may find the MATLAB functions imread, double, imagesc, circshift, fft2, ifft2 of use. DFT’s base frequency. However, such a peak is distorted in two ways: Try different ways to control the noise (see above) and distortions near the margins (window- ing). The MATLAB image processing toolbox provides ready-made “professional” functions → Its amplitude is lower (down to 63.7%). deconvwnr, deconvreg, deconvlucy, edgetaper, for such tasks. Do not use these, except per- haps to compare their outputs with the results of your own attempts. → Much signal energy has “leaked” to other frequencies. 74 76The reason for the leakage and scalloping losses is easy to visualize with the help of the convolution theorem: 35 The operation of cutting a sequence of the size of the DFT input vector out 30 of a longer original signal (the one whose continuous Fourier spectrum we 25 try to estimate) is equivalent to multiplying this signal with a rectangular 20 function. This destroys all information and continuity outside the “window” that is fed into the DFT. 15 Multiplication with a rectangular window of length T in the time domain is 10 equivalent to convolution with sin(πfT)/(πfT) in the frequency domain. 5 The subsequent interpretation of this window as a periodic sequence by 0 16 0 the DFT leads to sampling of this convolution result (sampling meaning 5 10 15.5 15 20 multiplication with a Dirac comb whose impulses are spaced f /n apart). s 25 30 15 input freq. Wherethewindowlengthwasanexactmultipleoftheoriginalsignalperiod, DFT index sampling of the sin(πfT)/(πfT) curve leads to a single Dirac pulse, and The leakage of energy to other frequency bins not only blurs the estimated spec- thewindowingcausesnodistortion. Inallothercases, theeffectsofthecon- trum. The peak amplitude also changes significantly as the frequency of a tone volution become visible in the frequency domain as leakage and scalloping changes from that associated with one output bin to the next, a phenomenon losses. known as scalloping. In the above graphic, an input sine wave gradually changes from the frequency of bin 15 to that of bin 16 (only positive frequencies shown). 77 79 Windowing Some better window functions Sine wave Discrete Fourier Transform 1 300 1 200 0.8 0 100 0.6 −1 0 0 200 400 0 200 400 0.4 Sine wave multiplied with window function Discrete Fourier Transform 100 0.2 1 Rectangular window Triangular window 0 0 50 Hann window Hamming window 0 0.2 0.4 0.6 0.8 1 −1 0 All these functions are 0 outside the interval 0,1. 0 200 400 0 200 400 78 80

Advise: Why You Wasting Money in Costly SEO Tools, Use World's Best Free SEO Tool Ubersuggest.