1 / 48

UNIVERSITY of MANCHESTER Department of Computer Science CS3291 Digital Signal Processing ‘05

UNIVERSITY of MANCHESTER Department of Computer Science CS3291 Digital Signal Processing ‘05 Section 8 Introduction to the DFT. DTFT of {x[n]} is:. If {x[n]} obtained from x a (t), correctly bandlimited then:. for -  <  <  :  =  T

haru
Download Presentation

UNIVERSITY of MANCHESTER Department of Computer Science CS3291 Digital Signal Processing ‘05

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. UNIVERSITY of MANCHESTER Department of Computer Science CS3291 Digital Signal Processing ‘05 Section 8 Introduction to the DFT CS3291 : Section 8

  2. DTFT of {x[n]} is: If {x[n]} obtained from xa(t), correctly bandlimited then: for - <  < :  = T DTFT is convenient way of calculating Xa(j) by computer. CS3291 : Section 8

  3. Two difficulties: • (i)Infinite range of summation • (ii) X(e j ) is continuous function of  • Solutions: • (i) Time-domain windowing: • Restrict {x[n]} to { x[0], x[1], … x[N-1]}  {x[n]} 0, N-1 • (ii) Frequency-domain sampling: • Store values of X(ej) for - <  <  • For real signals we only need 0  <  but generalise to complex signals • Instead of - <  <  take 0  < 2 CS3291 : Section 8

  4. Why take 0  < 2 ? • X(e j ) = X( e j ( + 2 ) ) for any  • So for X(e - j  / 3 ) look up X(e j 5  / 3 ) . • Same information, & it is convenient for  to start off at 0. • In many cases, not interested in  >  CS3291 : Section 8

  5. Taking M equally spaced samples over 0  < 2 we get : { X[k] } 0, M-1  { X[0], X[1],…, X[M-1] } where X[k] = X(exp(jk)) with k = 2k/M CS3291 : Section 8

  6. For spectral analysis, the larger M, the better for drawing accurate graphs etc. • But, if we need minimum M for storage of unambiguous spectrum, take M=N. • DFT: {x[n]} 0, N-1 {X[k]} 0, N-1 •  • (complex) (complex) k = 2k/N CS3291 : Section 8

  7. DFT transforms a finite sequence to another finite sequence. • DTFT transforms infinite sequence to continuous functn of  Inverse DFT: {X[k]}0, N-1 {x[n]}0, N-1 Note Similarity with DFT: CS3291 : Section 8

  8. Programming the DFT & its inverse: k = 2k/N • Similarity exploited by programs able to perform DFT or its inverse using same code. • Programs to implement these equations in a ‘direct’ manner given in MATLAB (using complex arith) & C (using real arith only). • These ‘direct’ programs are very slow & FFT is much faster. CS3291 : Section 8

  9. Direct DFT using complex arithmetic in MATLAB % Given N complex time-domain samples in array x[1:N] E = 2*pi/N ; for k=0 : N-1 X(1+k) = 0 + j*0 ; Wk =k*E ; for L = 0 : N-1 C = cos(L*Wk) + j *sin(L*Wk); X(1+k) = X(1+k) + x(1+L) * C; end; end; % Now have N complex freq-dom samples in array X[1:N] CS3291 : Section 8

  10. Direct inverse DFT using complex arithmetic in MATLAB % Given N complex freq-domain samples in array x[1:N] E = -2*pi/N ; for k=0 : N-1 X(1+k) = 0 + j*0 ; Wk =k*E ; for L = 0 : N-1 C = cos(L*Wk) + j *sin(L*Wk); X(1+k) = X(1+k) + x(1+L) * C; end; X(1+k) = X(1+k)/N ; end; % Now have N complex time-dom samples in array X[1:N] CS3291 : Section 8

  11. Direct forward/inverse DFT using complex arith in MATLAB % Given N complex samples in array x[1:N] if (Invers == 1) E = -2*pi/N else E = 2*pi/N ; for k=0 : N-1 X(1+k) = 0 + j*0 ; Wk =k*E ; for L = 0 : N-1 C = cos(L*Wk) + j *sin(L*Wk); X(1+k) = X(1+k) + x(1+L) * C; end; if (Inverse == 1) X(1+k) = X(1+k)/N ; end; % Now have N complex samples in array X[1:N] CS3291 : Section 8

  12. // Direct fwd/inverse DFT using real arith only in ‘C’ void directdft(void) // DFT or Inverse DFT by direct method. { // Order=N, Real & imag parts of input in arrays xr & xi // Output:- Real part in array X, Imag part in Y // Invers is 0 for DFT, 1 for IDFT int k, L; float c,e,s,wk; if(Invers==1) e = -2.0*PI/(float)N; else e = 2.0*PI/(float)N; for(k=0;k<N;k++) { x[k]=0.0; y[k]=0.0 ; wk=(float)k*e ; for(L=0; L<N; L++) { c=cos(L*wk); s=sin(L*wk); x[k] = x[k] + xr[L]*c + xi[L]*s; y[k] = y[k] + xi[L]*c - xr[L]*s;} if(Invers==1) { x[k]=x[k]/(float)N; y[k]=y[k]/(float)N;} } } CS3291 : Section 8

  13. Fast Fourier Transform (FFT) • An FFT algorithm, programmed in "C ", is presented in Table 2. • Gives ‘exactly’ same results as DFT only much faster. • Its detail & how speed is achieved is outside our syllabus. • We are interested in how to use DFT & interpret its results. • MATLAB has efficient fft procedure in its ‘SP tool-box'. • Don’t need to know how it’s programmed, only how to use it! • Direct DFT programs of academic interest only. • Table 4 is MATLAB program which reads music from a 'wav' file, splits it up into 512 sample sections & • performs a DFT (by FFT) analysis on each section. CS3291 : Section 8

  14. Effect of windowing (old) • DFT of {x[0], x[1], ..., x[N-1] } is frequency-sampled version of DTFT of infinite sequence {x[n]} with all samples outside range n= 0 to N-1 set to zero. • {x[n]} effectively multiplied by "rectangular window" {r[n]} When N is even, the DTFT, R(ej) of {r[n]} is: Dirichlet Kernel CS3291 : Section 8

  15. Effect of windowing (newer) • Given {x[n]}0,N-1 {x[0], x[1], ..., x[N-1] } • Assumed obtained by multiplying infinite sequence {x[n]} by "rectangular window" {r[n]} When N is even, the DTFT, R(ej) of {r[n]} is: Dirichlet Kernel CS3291 : Section 8

  16. n 0 N • Effect of windowing • Given {x[n]}0,N-1 {x[0], x[1], ..., x[N-1] } • Assumed obtained by multiplying infinite sequence {x[n]} by "rectangular window" {r[n]} r[n] CS3291 : Section 8

  17. When N is even, the DTFT, R(ej) of {r[n]} is: Dirichlet Kernel CS3291 : Section 8

  18. sincs10(/(2)) plotted against : -  CS3291 : Section 8

  19. -  CS3291 : Section 8

  20. Notation “sincsM(/(2))” is not in textbooks. • This function is also encountered when designing FIR filters. • Causes the stop-band ripples which appear when FIR low-pass filters are designed with rectangular windows. • Magnitude of R(ej  ) shown above when M=10. • Note relatively narrow main lobe & side-lobes. • Zero-crossings occur at  = 2 / M , 4 / M, etc. • Like a "sinc" function in many ways. CS3291 : Section 8

  21. Frequency-domain convolution: • DTFT of product of {x[n]} & {r[n]} is complex ( freq.-domain)convolution between X(ej) & R(ej) denoted X(ej)  R(ej): • Observe the form of this expression & its similarity with time-domain convolution formulae. • Also note the (1/2) factor & limits - to  which may remind us of the inverse DTFT. CS3291 : Section 8

  22. Duality of time- & frequency-domain convolution (i) Fourier transform of h(t)  x(t) is H(j).X(j). (ii) DTFT of {h[n]} {x[n]} is H(e j  ).X(e j  ) Time-domain Frequency-domain Convolution (filtering) Multiplication Multiplication (windowing) Complex convolution CS3291 : Section 8

  23. Duality of time- & frequency-domain convolution Time-domain Frequency-domain Convolution (filtering) Multiplication Multiplication (windowing) Complex convolution CS3291 : Section 8

  24. Spectral analysis of sampled sine-waves • DTFT of {cos(0n)} cannot be found directly. • But inverse DTFT of • X(e j  ) = ( - 0) + ( + 0) is: •  • (1/2)((-0) + (+0) ) ejn d •  • = (1/2)[ e j o n + e -j o n ] = cos ( 0 n) for - < n <  • It may be inferred that DTFT of {cos(0n) } is: • ( - 0) + ( + 0) CS3291 : Section 8

  25. X(ej) DTFT of {cos(n)}    -0 - 0  CS3291 : Section 8

  26. X(e j  )    -0 - 0  R(e j( - ) ) • DTFT of rectangularly windowed sampled sine-wave • DTFT of {cos0n}0,N-1 = DTFT of {x[n].r[n]} • = X(ej)  R(ej) = P(ej) say. • By frequency-domain convolution formula,  CS3291 : Section 8

  27.  X(e j )  0 -0 -  R(e j(  ) )  CS3291 : Section 8

  28. X(e j ) = ( - 0) + (  + 0) Now draw this & its modulus: CS3291 : Section 8

  29. P(e j  ) - 0  0 CS3291 : Section 8

  30. | P(e j  ) |  0  0 CS3291 : Section 8

  31. When rectangular window is of width N samples: • Ampl of main peak = ampl of sine-wave*N/2 check! • N-1 zero-crossings between 0 and . • 1 main peak & N-2 ripples in magnitude spectrum. • Increasing N gives sharper peak & more ripples. CS3291 : Section 8

  32. Effect of windowing on DFT : • Effect on DTFT is frequency spreading & side lobes. • How does windowing with frequency sampling affect the DFT? • Effect is rather confusing • Consider 64 pt DFT of { cos(0.7363n)} in Fig 1. • 0 = 0.7363 lies between 7 = (2/64)*7 & 8 = (2/64)*8 • Samples of rectangular window seen. • X[7] and X[8] strongly affected by sinusoid. CS3291 : Section 8

  33. Amplitude  20 Fig 1 : Magnitude of 64 pt DFT spectrum of cos(0.7363n) CS3291 : Section 8

  34. Now consider 64 pt DFT of { cos (n/4) } in Fig 2 • 0 = /4 coincides exactly with 8. • Only X[8] affected. CS3291 : Section 8

  35. Ampl. = 32 Approx. 36% difference in ampl. For same ampl. of sinusoid. Fig2 : Magnitude of 64 point DFT of cos(pn/4). CS3291 : Section 8

  36. What has happened to the ‘sincs’ function in Figure 2 ? • We don’t see it because all the samples apart from one • occur at the ‘zero-crossings of the ‘sincs function’. • Effect of the rectangular window is no longer seen because the sampling points happen to coincide exactly with the zero-crossings of the sincs function. • This always occurs when the frequency of the cosine wave coincides exactly with a frequency sampling point. CS3291 : Section 8

  37. Difference between Figs 1 & 2 undesirable. Use non-rectangular windows e.g. Hann {w[n]}0,N-1 with (Slightly different definition from the one we had for FIR filters). CS3291 : Section 8

  38. CS3291 : Section 8

  39. In frequency-domain : • (i) broader main-lobe for Hann window whose width is approx doubled as compared with rectr window • (ii) reduced side-lobe levels. • For sine-wave, at least 3 spectral 'bins' strongly affected even when its frequency coincides with a sampling point. • Often take highest of the 3 peaks as amplitude of sine-wave. • With rectr : just one bin affected when frequency of sine-wave coincides with a bin. If frequency then shifted towards next bin,  35% variation amplitude seen. • With Hann: Variation of only about 15% seen. • Still undesirable, but not as bad. CS3291 : Section 8

  40. Reduction in amplitude estimation error for sine-waves is at expense of some loss of spectral resolution • With Hann window, 3 bins strongly affected by one sine-wave • We will only know that the frequency of the sine-wave lies within the range of these 3 frequencies. CS3291 : Section 8

  41. DFT & FFT have many applications in DSP esp in comms. • e.g. filtering by (1) performing an FFT, • (2) zeroing unwanted spectral components • (3) performing an inverse FFT. • FFT is ‘Swiss army knife' of signal processing. • Most ‘spectrum analysers’ use an FFT algorithm. • Some applications of spectral estimation are: • determining frequency & loudness of a musical note, • deciding whether a signal is periodic or non-periodic, • analysing speech to recognise spoken words, • looking for sine-waves buried in noise, • measuring frequency distribution of power or energy. CS3291 : Section 8

  42. Spectral analysis of 'power signals' {x[n]} • These exist for all time. • Extract segment {x[n]}0,N-1 to represent {x[n]}. • Energy of {x[n]}0,N-1 is: • Mean square value (MSV) of this segment is (1/N) E • This is “power” of a periodic discrete time signal, of period N samples, for which a single cycle is {x[n]}0,N-1 . • It may be used as an estimate of the power of {x[n]}. CS3291 : Section 8

  43. Parseval's Theorem for the DFT It may be shown that for a real signal segment {x[n]}0,N-1: Proof CS3291 : Section 8

  44. Parseval’s Theorem allows energy & power estimates to be calculated in frequency-domain instead of in time-domain. • Allows us to see how power is distributed in frequency. • Is there more power at high frequencies than at low frequencies or vice-versa? • By Parseval's thm, estimate of power of {x[n]}, obtained by analysing {x[n]}0,N-1 , is: • Usefulness of this estimate illustrated by following example. CS3291 : Section 8

  45. Example: Real periodic signal {x[n]} is rect windowed to give {x[n]}0,39 . 40-point DFT gives magnitude spectrum below. Estimate power of {x[n]} & comment on reliability of estimate. If {x[n]} is passed thro’ ideal digital low-pass filter with cut-off /2 radians/sample how is power likely to be affected? CS3291 : Section 8

  46. Ans: • MSV of {x[n]}0,39 power of {x[n]} • = (1/1600)[2*402 +2*302 +2*202+2*102] = 3.75 watts. • Reduces to 3.125 watts, i.e. by 0.8 dB • Care needed to interpret such results as power estimates. • For periodic or deterministic (non-random) signals: • estimates from segments extracted from different parts of {x[n]} may be similar, & estimates could be fairly reliable. • For random signals : may be considerable variation from estimate to estimate. Averaging may be necessary. CS3291 : Section 8

  47. Power spectral density (PSD) estimate • For N-point DFT, X[k]2 / N2 is estimate of PSD in Watts per “bin”. • A “bin” is a band-width 2/N radians/sample centred on k • ( fS / N Hz centred on k fS / N ) • Instead of |X[k]| often plot 10 log10 (|X[k]|2/N2) dB. against k. • “PSD estimate” graph. • Careful with random signals : each spectral estimate different. • Can average several PSD estimates. CS3291 : Section 8

  48. Spectral analysis of signals: • File: OPERATOR.pcm contains sampled speech. • SNR-12dB.pcm contains sine-wave corrupted with noise. • Sampled at 8 kHz using 12-bit A/D converter. • May be read into "MATLAB" program in Table 5 & spectrally analysed using the FFT. • Meaningless to analyse a large speech file at once. • Divide into blocks of 128 or 256 samples & analyse separately. • Blocks of N (= 512) samples may be read in and displayed. • Programs in notes available on: ~barry/mydocs/CS3291 • In each case, a rectangular window is used. • Exercise, modify programs to have a Hann window. • Notice effectiveness of DFT in locating sine-wave in noise • even when it cannot be seen in time-domain graph. CS3291 : Section 8

More Related