Digital Signal Processing Moslem Amiri, Václav Přenosil Masaryk University Understanding Digital Signal Processing, Third Edition, Richard Lyons (0-13-261480-4) © Pearson Education, 2011. The Discrete Fourier Transform (1) 2 The Discrete Fourier Transform  Discrete Fourier transform (DFT)  DFT is a mathematical procedure used to determine harmonic, or frequency, content of a discrete signal sequence  DFT’s origin is continuous Fourier transform X(f) where x(t) is some continuous time-domain signal  DFT equation (exponential form) x(n) is sequence of time-domain sampled values dtetxfX ftj 2 )()(     Nmnj N n enxmX /2 1 0 )()(     3 Understanding the DFT Equation  DFT equation (rectangular form)  From Euler’s relationship, e−jø = cos(ø) − jsin(ø)  m = index of DFT output in frequency domain m = 0, 1, 2, 3, . . ., N−1  n = time-domain index of input samples n = 0, 1, 2, 3, . . ., N−1  N = number of samples of input sequence and number of frequency points in DFT output )]/2sin()/2[cos()( )()( 1 0 /2 1 0 NmnjNmnnx enxmX N n Nmnj N n            4 Understanding the DFT Equation  X(m) DFT output  Each X(m) DFT output term is sum of point-forpoint product between an input sequence of signal values and a complex sinusoid of the form cos(ø) − jsin(ø)  Exact frequencies of different sinusoids depend on both sampling rate fs at which the original signal was sampled, and number of samples N  The N separate DFT analysis frequencies are N mf mf s analysis )( 5 Understanding the DFT Equation 6 Understanding the DFT Equation  Magnitude and power contained in each X(m) term 222 1 22 )()()()(:spectrumpower )( )( tan)( )()()()( )(ofangleanat)()()()( mXmXmXmX mX mX mX mXmXmXmX mXmXmjXmXmX imagrealmagPS real imag imagrealmag magimagreal                7 Understanding the DFT Equation  Example  We want to sample and perform an 8-point DFT on a continuous input signal containing components at 1 kHz and 2 kHz, expressed as  The 2 kHz term is shifted in phase by 3π/4 radians relative to the 1 kHz sinewave )4/320002sin(5.0)10002sin()(   tttxin 8 Understanding the DFT Equation  Example (cont.)  N = 8  we need 8 input sample values on which to perform DFT  Sample rate = fs  sampling every 1/fs = ts sec.  If we sample xin(t) at fs = 8000 samples/second, DFT results will indicate what signal amplitude exists in x(n) at analysis frequencies of mfs/N, or 0 kHz, 1 kHz, 2 kHz, . . ., 7 kHz  We use DFT equation for m = 0, …, 7 )4/320002sin(5.0)10002sin()()(   sssin ntntntxnx )]/2sin()/2[cos()()( 1 0 NmnjNmnnxmX N n      9 Understanding the DFT Equation 10 Understanding the DFT Equation 11 Understanding the DFT Equation  Example (cont.)        9040.40.0)7(:kHz7 8 kHz87 452414.1414.1)6(:kHz6 8 kHz86 000.00.0)5(:kHz5 8 kHz85 000.00.0)4(:kHz4 8 kHz84 000.00.0)3(:kHz3 8 kHz83 452414.1414.1)2(:kHz2 8 kHz82 9040.40.0)1(:kHz1 8 kHz81                      jmX N mf jmX N mf jmX N mf jmX N mf jmX N mf jmX N mf jmX N mf s s s s s s s 12 Understanding the DFT Equation  Example (cont.)  When m = 0  X(0) is sum of x(n) samples  X(0) is proportional to average of x(n)  = N times x(n)’s average value  X(0) frequency term is the non-time-varying (DC) component of x(n)  Our x(n) has no DC component       1 0 1 0 )()]0sin()0[cos()()0( N n N n nxjnxX  000.00.0)0(:kHz0 8 kHz80    jmX N mfs 13 Understanding the DFT Equation 14 Understanding the DFT Equation  Fig. 3-4  Indicates that xin(t) has signal components at 1 kHz (m = 1) and 2 kHz (m = 2)  1 kHz tone has a magnitude twice that of 2 kHz tone  DFT phase at frequency mfs/N is relative to a cosine wave at that same frequency of mfs/N Hz where m = 1, 2, 3, ..., N−1  E.g., phase of X(1) is −90 degrees, so input sinusoid whose frequency is 1 · fs/N = 1000 Hz was a cosine wave having an initial phase shift of −90 degrees (or a sinewave having an initial phase of zero) 15 Understanding the DFT Equation  Fig. 3-4  When DFT input signals are real-valued, DFT phase at 0 Hz (m = 0, DC) is always zero because X(0) is always real-only     1 0 )()0( N n nxX 16 DFT Symmetry  Fig. 3-4  There is a symmetry in DFT results  When input sequence x(n) is real, complex DFT outputs for m = 1 to m = (N/2) − 1 are redundant with frequency output values for m > (N/2) for m = 1, 2, 3, . . . , N−1  To obtain DFT of x(n), we need only compute the first N/2+1 values of X(m) where 0 ≤ m ≤ (N/2); X(N/2+1) to X(N−1) DFT output terms provide no additional information about spectrum of real sequence x(n) )()( degrees)(at)( degrees)(at)()( mNXmX mNXmNX mXmXmX       17 DFT Symmetry  Proving symmetry of DFT of real input sequences  X(N−m) is merely X(m) with the sign reversed on X(m)’s exponent                        1 0 /21)2sin()2cos( /22 1 0 /)(2 1 0 onsubstituti /2 1 0 )( )( )()( )()( 2 N n Nmnjnjne Nmnjnj N n NmNnj N n Nmnj N n enx eenx enxmNX enxmX nj      18 DFT Symmetry  An additional symmetry property of DFT  Determine DFT of real input functions where input index n is defined over both positive and negative values  If real input function is even, x(n) = x(−n), then X(m) is always real and even  Xreal(m) is in general nonzero and Ximag(m) is zero  If real input function is odd, x(n) = −x(−n), then Xreal(m) is always zero and Ximag(m) is, in general, nonzero 19 DFT Linearity  DFT has linearity property  Proof )()()()()()( )()( )()( 2121 22 11 mXmXmXnxnxnx mXnx mXnx sum DFT sum DFT DFT         )()( )()( )()()( )()( 21 /2 1 0 2 /2 1 0 1 /2 1 0 21 /2 1 0 mXmX enxenx enxnxmX enxmX Nmnj N n Nmnj N n Nmnj N n sum Nmnj N n                       20 DFT Magnitudes  Output magnitude of DFT  When a real input signal contains a sinewave component of peak amplitude Ao with an integral number of cycles over N input samples, output magnitude of DFT for that sinewave is  For real inputs, hardware memory registers must be able to hold values as large as N/2 times the maximum amplitude of input sample values  If DFT input is a complex sinusoid of magnitude Ao (i.e., Aoej2πfnts) with an integer number of cycles over N samples, output magnitude of DFT for that sinewave is 2/NAM or  NAM oc  21 DFT Magnitudes  DFT is occasionally defined as  Some commercial software packages use  Forward and inverse DFTs  Scale factors are used so that there’s no scale change when transforming in either direction Nmnj N n enx N mX /2 1 0 )( 1 )('     Nmnj N m Nmnj N n emX N nx enx N mX /2 1 0 /2 1 0 )('' 1 )( )( 1 )(''            22 Summary  To recap what we’ve learned so far  Each DFT output term is sum of term-by-term products of an input time-domain sequence with a sine and a cosine wave sequences  For real inputs, an N-point DFT’s output provides only N/2+1 independent terms  DFT is a linear operation  Magnitude of DFT results is directly proportional to N  DFT’s frequency resolution (spacing) is fs/N  X(m = N/2) corresponds to mfs/N = fs/2 = folding (Nyquist) frequency 23 DFT Shifting Theorem  Shifting theorem property of DFT  A shift in time of a periodic x(n) input sequence manifests itself as a constant phase shift in angles associated with DFT results  If we sample x(n) starting at n = k (integer), DFT of those time-shifted sample values is Xshifted(m) where  If the point where we start sampling x(n) is shifted to right by k samples, DFT output spectrum of Xshifted(m) is X(m) with each of X(m)’s complex terms multiplied by linear phase shift ej2πkm/N, which is merely a phase shift of 2πkm/N radians )()( /2 mXemX Nmkj shifted   24 DFT Shifting Theorem  Example  Suppose we sampled preceding DFT example input sequence later in time by k = 3 samples xin(t) = sin(2π1000t) + 0.5sin(2π2000t+3π/4) 25 DFT Shifting Theorem 26 DFT Shifting Theorem 27 DFT Shifting Theorem  Fig. 3-6  Magnitude of Xshifted(m) should be unchanged from that of X(m)  Fig. 3-6(a) is identical to Fig. 3-4(a)  Phase of DFT result does change depending on the instant at which we started to sample xin(t)  E.g., X(1) from preceding DFT Example had a magnitude of 4 at a phase angle of −π/2  k = 3 and N = 8  Xshifted(1): magnitude of 4, phase angle of π/4 4/)8/48/6( 2/8/132/2 44 4)1()1(   jj jjNmkj shifted ee eeXeX     28 Inverse DFT  Inverse discrete Fourier transform (IDFT)  We can obtain the original time-domain signal by performing IDFT on X(m) frequency-domain values  A discrete time-domain signal can be considered the sum of various sinusoidal analytical frequencies and that the X(m) outputs of the DFT are a set of N complex values indicating the magnitude and phase of each analysis frequency comprising that sum )]/2sin()/2[cos()( 1 )( 1 )( 1 0 /2 1 0 NnmjNnmmX N emX N nx N m Nnmj N m           29 DFT Leakage  Leakage  Causes DFT results to be only an approximation of the true spectra of the original input signals prior to digital sampling  Although there are ways to minimize leakage, we can’t eliminate it entirely 30 DFT Leakage  Leakage  DFT produces correct results only when input data sequence contains energy precisely at integral multiples of fundamental frequency fs/N  If input has a signal component at some intermediate frequency between analytical frequencies of mfs/N, say 1.5fs/N, this input signal will show up to some degree in all of N output analysis frequencies of DFT  We say that input signal energy shows up in all of DFT’s output bins  DFT samples = “bins” 1,...,2,1,0where,)(  Nm N mf mf s analysis 31 DFT Leakage 32 DFT Leakage an input sequence having 3.4 cycles over N = 64 samples because input sequence does not have an integral number of cycles over 64-sample interval, input energy has leaked into all the other DFT output bins m = 4 bin, for example, is not zero because sum of products of input sequence and m = 4 analysis frequency is no longer zero 33 DFT Leakage  Cause of leakage  For a real cosine input having k cycles (k need not be an integer) in N-point input time sequence, amplitude response of an N-point DFT bin in terms of bin index m is approximated by sinc function  Ao is peak value of DFT’s input sinusoid  For our examples here, Ao is unity )( )](sin[ 2 )( mk mkNA mX o      34 DFT Leakage the curve comprises a main lobe and periodic peaks and valleys known as sidelobes DFT output will be a sampled version of the continuous spectrum when DFT’s input sequence has exactly an integral k number of cycles (centered exactly in the m = k bin), no leakage occurs 35 DFT Leakage  Example (Fig. 3-10(a))  A real 8 kHz sinusoid, having unity amplitude, is sampled at a rate of fs = 32000 samples/second  If we take a 32-point DFT of samples, DFT’s frequency resolution, or bin spacing, is fs/N = 32000/32 Hz = 1.0 kHz  We can predict DFT’s magnitude response by centering input sinusoid’s spectral curve at positive frequency of 8 kHz  DFT output is a sampled version of continuous spectral curve  DFT outputs reside on continuous spectrum at its peak and exactly at curve’s zero crossing points 36 DFT Leakage frequency-domain sampling results in nonzero magnitudes for all DFT output bins results in the leaky DFT output shown 37 DFT Leakage If the continuous spectra that we’re sampling are symmetrical, why does DFT output in this Figure look so asymmetrical? Bins to the right of third bin are decreasing in amplitude faster than bins to the left of third bin With k = 3.4 and m = 3, from sinc function the X(3) bin’s magnitude should be 24.2— but Fig. (b) shows X(3) bin magnitude to be greater than 25. Why? 38 DFT Leakage 39 DFT Leakage  Fig. 3-11  When examining a DFT output, we’re normally interested only in m = 0 to m = (N/2−1) bins  Thus, only the first 32 bins are shown in Fig. 3-8(b)  DFT is periodic in frequency domain  Upon examining DFT’s output for higher frequencies, we end up going in circles, and spectrum repeats itself forever  Fig. 3-11 shows cyclic representation of 64-point DFT shown in Fig. 3-8(b)  The more conventional way to view a DFT output is to unwrap the spectrum in Fig. 3-11 to get the spectrum in Fig. 3-12 40 DFT Leakage 41 DFT Leakage  Fig. 3-12  As some of the input 3.4-cycle signal amplitude leaks into 2nd bin, 1st bin, and 0th bin, leakage continues into −1st bin, −2nd bin, −3rd bin, etc  63rd bin = −1st bin, 62nd bin = −2nd bin, and so on  These bin equivalencies allow us to view DFT output bins as if they extend into negative-frequency range, as shown in Fig. 3-13(a)  Result is that the leakage wraps around m = 0, as well as around m = N  m = 0 frequency is m = N frequency  Leakage wraparound at m = 0 accounts for the asymmetry around DFT’s m = 3 bin in Fig. 3-8(b) 42 DFT Leakage 43 DFT Leakage  Fig. 3-13(b)  when a DFT input sequence x(n) is real, DFT outputs from m = 0 to m = (N/2−1) are redundant with frequency bin values for m > (N/2)  |X(m)| = |X(N−m)|  This means that leakage wraparound also occurs around m = N/2 bin  This can be illustrated using an input of 28.6 cycles per sample interval (32 − 3.4) whose spectrum is shown in Fig. 3-13(b)  Figs. 3-13(a) and (b) are similar 44 DFT Leakage 45 DFT Leakage  Fig. 3-14  DFT exhibits leakage wraparound about m = 0 and m = N/2 bins  Minimum leakage asymmetry will occur near N/4th bin 46 Windows  Windowing  Windowing reduces DFT leakage by minimizing magnitude of sinc function’s sidelobes  Done by forcing amplitude of input time sequence at both beginning and end of sample interval to go smoothly toward a single common amplitude value 47 Windows 48 Windows  Fig. 3-15  Considering infinite-duration time signal shown in (a), a DFT can only be performed over a finitetime sample interval like that shown in (c)  We can think of DFT input signal in (c) as product of (a), and rectangular window whose magnitude is 1 over sample interval shown in (b)  Anytime we take DFT of a finite-extent input sequence, we are, by default, multiplying that sequence by a window of all ones  Sinc function’s sin(x)/x shape is caused by this rectangular window because CFT of rectangular window in (b) is sinc function 49 Windows  Fig. 3-15  Rectangular window’s abrupt changes between one and zero cause sidelobes in sinc function  To minimize spectral leakage caused by those sidelobes, we have to reduce sidelobe amplitudes by using window functions other than rectangular window  If we multiply DFT input, (c), by triangular window function, (d), to obtain windowed input signal, (e), values of final input signal appear to be the same at beginning and end of sample interval in (e)  Reduced discontinuity decreases level of relatively high-frequency components in overall DFT output; that is, DFT bin sidelobe levels are reduced in magnitude using a triangular window 50 Windows  Fig. 3-15  There are window functions that reduce leakage even more than triangular window  Hanning window in (f)  Product of window in (f) and input sequence provides signal shown in (g) as input to DFT  Hamming window in (h)  It’s much like Hanning window, but it’s raised on a pedestal 51 Windows 1,...,2,1,0for 2 cos46.054.0)(:windowHamming 1,...,2,1,0for 2 cos5.05.0)(:windowHanning 1,...,22/,12/for 2/ 2 2/,...,2,1,0for 2/)(:windowTriangular 1,...,2,1,0for,1)(:r windowRectangula )()()( 1 0 /2                             Nn N n nw Nn N n nw NNNn N n Nn N n nw Nnnw enxnwmX N n Nmnj w    52 Windows 53 Windows  Fig. 3-16  Hamming, Hanning, and triangular give reduced sidelobe levels relative to rectangular  Because Hamming, Hanning, and triangular reduce time-domain signal levels, their main lobe peak values are reduced relative to rectangular  Log magnitude response (normalized)          )0( )( log20)( 10 W mW mWdB 54 Windows  Fig. 3-16(b)  Main lobe of rectangular window’s magnitude response is the most narrow, fs/N  However, its first sidelobe level is only −13 dB below main lobe peak, which is not good  Various nonrectangular windows’ wide main lobes degrade the windowed DFT’s frequency resolution by almost a factor of two  However, important benefits of leakage reduction usually outweigh the loss in DFT frequency resolution  Hanning has rapid sidelobe roll-off  Hamming has lower first sidelobe levels, but its sidelobes roll off slowly relative to Hanning 55 Windows shape of Hanning window’s response looks broader and has a lower peak amplitude, but its sidelobe leakage is noticeably reduced from that of rectangular window 56 Windows 57 Windows  Windows  Overall frequency resolution and signal sensitivity are affected much more by size and shape of window function than mere size of DFTs  There are many different window functions described in literature of DSP  Window selection is a trade-off between main lobe widening, first sidelobe levels, and how fast the sidelobes decrease with increased frequency  Use of any particular window depends on application 58 DFT Scalloping Loss  Scalloping  Fluctuations in overall magnitude response of an N-point DFT  When no input windowing function is used, sin(x)/x shape of sinc function’s magnitude response applies to each DFT output bin 59 DFT Scalloping Loss 60 DFT Scalloping Loss  Fig. 3-19  (a) shows a DFT’s aggregate magnitude response by superimposing several sin(x)/x bin magnitude responses  Sinc function’s sidelobes are not key here  In (b), overall DFT frequency-domain response is indicated by bold envelope curve  This rippled curve, also called picket fence effect, illustrates processing loss for input frequencies between bin centers  Magnitude of DFT response fluctuates from 1.0, at bin center, to 0.637 halfway between bin centers  This envelope ripple exhibits a scalloping loss of −4 dB halfway between bin centers 61 DFT Scalloping Loss  Fig. 3-19  Illustrates a DFT output when no window (i.e., a rectangular window) is used  Because nonrectangular window functions broaden DFT’s main lobe, their use results in a scalloping loss that will not be as severe as with rectangular window  Their wider main lobes overlap more and fill in valleys of envelope curve in (b)  Scalloping loss isn’t a severe problem in practice  Real-world signals normally have bandwidths that span many frequency bins so that DFT magnitude response ripples can go almost unnoticed