1 / 84

UNIVERSITY of MANCHESTER School of Computer Science Comp30291 Digital Media Processing 2009-10

UNIVERSITY of MANCHESTER School of Computer Science Comp30291 Digital Media Processing 2009-10 Section 7 Discrete Fourier Transform. Introduction. Given an analogue signal x a (t) with Fourier Transform:. Assume x a (t) is band-limited between 0 & F s/2.

gemma-ware
Download Presentation

UNIVERSITY of MANCHESTER School of Computer Science Comp30291 Digital Media Processing 2009-10

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 School of Computer Science Comp30291 Digital Media Processing 2009-10 Section 7 Discrete Fourier Transform Comp30291 : Section 7

  2. Introduction • Given an analogue signal xa(t) with Fourier Transform: • Assume xa(t) is band-limited between 0 & Fs/2. • If xa(t) is sampled at Fs = 1/T to give {x[n]} with DTFT: it may be shown that: DTFT is convenient way of calculating Xa(j) by computer. Comp30291 : Section 7

  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 samples of X(ej) Comp30291 : Section 7

  4. Sampling X(ej) • Storing M equally spaced samples of X(ej) over the range = - to  would be possible. • For real signals, only the range 0 to  is of interest. • In practice, the range is made 0 to 2 rather than - to . • So we generally disregard the range from  to 2. Comp30291 : Section 7

  5. With M equally spaced samples over 0  < 2 { X[k] } 0, M-1  { X[0], X[1],…, X[M-1] } |X(ej)| |X[3]| |X[1]| |X[M - 1]| |X[0]| W 2p W 0 W1 W2 W3 p M - 1 Comp30291 : Section 7

  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. Comp30291 : Section 7

  7. Discrete Fourier Transform (DFT) • Transforms: {x[n]} 0, N-1to {X[k]} 0, N-1 •  • (complex) (complex) • DFT transforms a finite sequence to another finite sequence. • DTFT transforms infinite sequence to continuous functn of  Comp30291 : Section 7

  8. Inverse DFT Inverse DFT: {X[k]}0, N-1 {x[n]}0, N-1 Note Similarity with DFT: Comp30291 : Section 7

  9. 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. Comp30291 : Section 7

  10. 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] Comp30291 : Section 7

  11. 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;} } } Comp30291 : Section 7

  12. Fast Fourier Transform (FFT) • An FFT algorithm in ‘C ‘ is given on next slide. • Gives ‘exactly’ same results as DFT only much faster. • Its detail & how speed is achieved is outside our syllabus. • Works best when N is a power of 2, e.g. 32, 64, 512, 1024, etc. • 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! • ‘C’ version of FFT for reference on next slide. • Direct DFT programs of academic interest only. Comp30291 : Section 7

  13. void fft(void) { int N1,N2,NM,N2M,i,j,k,l; float e,c,s,a,xt,yt; N2=N ; NM=N-1; if(Invers==1) e= -PI/(float)N; else e= PI/(float)N; for ( k=2; k<N+1; k=k+k) { N1=N2; N2=N2/2; N2M=N2-1 ; e=e*2.0; a=0.0; for (j=0; j < N2M+1; j++) { c=cos(a); s=sin(a); a=a+e; for ( i=j; i<NM+1; i=i+N1) { l=i+N2; xt=x[i]-x[l]; x[i]=x[i]+x[l]; yt=y[i]-y[l]; y[i]=y[i]+y[l]; x[l]=xt*c+yt*s; y[l]=yt*c-xt*s;} } //end of j loop } //end of k loop if(Invers==1) for (k=0; k<(NM+1); k++) {x[k]=x[k]/(float)N; y[k]=y[k]/(float)N ;} j=0; for(i=0; i<N-1; i++) { if (i<j) {xt=x[j]; x[j]=x[i]; x[i]=xt; yt=y[j]; y[j]=y[i]; y[i]=yt;} k=N/2; while((k-1)<j) { j=j-k; k=k/2;} j=j+k; } // end of i loop } // end of procedure FFT Comp30291 : Section 7

  14. DFT (FFT) for spectral analysis • Given {x[n]}0,N-1 we get {X[k]}0,N-1 • N complex spectral samples from 0 to fs. • When {x[n]} is real, plot magnitudes of X[k] for k=0 to N/2. • When N=512, k=256 corresponds to fs/2. • MATLAB arrays cannot start from zero, so x[n] is stored as x(1+n) • Similarly {X[k]}0,N-1 are stored as X(1) ... X(N). • We analyse some sinusoids on next slide. Comp30291 : Section 7

  15. FFT analysis of sinusoids for n=0:63 x(1+n)=100*cos(0.25*pi*n) + 100*sin(0.5*pi*n); end; X=fft(x); plot(abs(X(1:32))); grid on; There are 64 pts but we only plot 32. The horiz axis goes from 1 to 32 for frequencies 0 to fs/2 (). The vertical axis can tell us the amplitudes of the sinusoids. Divide by N/2 i.e. 32 here. It can’t tell whether they are sines or cos as we plot only magnitude. Comp30291 : Section 7

  16. Scaled FFT analysis of sinusoids N=64; for n=0:N-1 x(1+n) = 100*cos(0.25*pi*n)+ 100*sin(0.5*pi*n); end; X=fft(x)/(N/2); plot(abs(X(1:N/2)),'*-'); grid on; Now it seems we can read the amplitudes directly from the graph. Comp30291 : Section 7

  17. Increased order FFT analysis of sinusoids clear all; N=128; for n=0:N-1 x(1+n) = 50*cos(0.25*pi*n)+ 100*sin(0.5*pi*n); end; X=fft(x)/(N/2); plot(abs(X(1:N/2)),'*-'); grid on; Now have more pts & we can read amplitudes & frequencies more accurately. Comp30291 : Section 7

  18. Trouble!! clear all; N=64; for n=0:N-1 x(1+n) = 100*cos((pi/4)*n) + 100*sin(1.522*n); end; X=fft(x)/(N/2); plot(abs(X(1:N/2)),'*-'); grid on; Changed frequency of 2nd sine-wave from /4 to 15.5/32=1.522. Now we don’t get the correct amplitude as 1.522 does not correspond to one of our frequency sampling pts. 35% error in amplitude reading. Any suggestions?. Comp30291 : Section 7

  19. Increase N to 128 clear all; N=128; for n=0:N-1 x(1+n) = 100*cos((pi/4)*n) + 100*sin(1.522*n); end; X=fft(x)/(N/2); plot(abs(X(1:N/2)),'*-'); grid on; This works as 15.5/32 = 31/64 But it may not be an option. We may only have 64 samples available as signal may be changing rapidly. Any other suggestions?. Comp30291 : Section 7

  20. Help from ‘zero-padding’ clear all; N=64; for n=0:N-1 x(1+n) = 100*cos((pi/4)*n) + 100*sin(1.522*n); end; for n = N:2*N-1, x(1+n)=0; end; % Length of x is now 128 X=fft(x)/(N/2); figure(1); plot(abs(X(1:N)),'*-'); grid on; Double length of ‘x’ by ‘zero-padding’. This doubles the no. of freq sampling pts & 1.522 now hits one of these. We are seeing the same graph with extra pts inserted by interpolation. The ‘ripples’ were there before but we did not sample them. Comp30291 : Section 7

  21. Zero-padding • Take a simple example with N=4 and zero-pad to N=8:: • {7 1 2 4} becomes {7 1 2 4 0 0 0 0} or, if you wish, {0 0 0 0 7 1 2 4} • Simple as that. • Now have twice as many samples • So we get twice as many points in freq-domain . • Can interpret FFT spectral graph more clearly. • (Note that scaling in ‘X=fft(x)/(N/2)’ has not changed) Comp30291 : Section 7

  22. Extend by symmetry • Another way of increasing no. of freq-domain points without increasing no. of time-domain samples is to extend by even symmetry: {6 1 2 4} becomes {6 1 2 4 4 2 1 6 } or {4 2 1 6 6 1 2 4 } • Both are OK but we normally take the second. Comp30291 : Section 7

  23. Extend by ‘even symmetry’ (2nd way) clear all; N=64; for n=0:N-1 x(1+n) = 100*cos((pi/4)*n) + 100*sin(1.522*n); end; for n = 0 : N-1 y(1+N+n)=x(1+n); y(N-n)=x(1+n); end; Y=fft(y) / N; plot(abs((Y(1:N))),'*-'); grid on; Doubles length of ‘x’ by ‘even symmetry’. Doubles no. of freq sampling pts About 30% error in amplitude of higher freq sin-wave. Comp30291 : Section 7

  24. Extending by even symmetry (cont) {6 1 2 4} becomes { 4 2 1 6 6 1 2 4} DFT of the extended sequence is: = e-3.5jk{4e3.5jk+2e2.5jk + e1.5jk + 6e0.5jk + 6e-0.5jk+e-1.5jk+ 2e-2.5jk+ 4e-3.5jk } = 2e-3.5jk{4cos(3.5k) + 2cos(2.5k) +cos(1.5k)+6cos(0.5k)} = 2e-3.5jkDCT[k] where DCT[k] = 6cos(0.5k)+cos(1.5k)+ 2cos(2.5k)+4cos(3.5k) Comp30291 : Section 7

  25. Discrete Cosine Transform (DCT) • The DCT has many forms: we consider one • Given {x[n]}0,N-1 its DCT is • Note that k = 2k/(2N) • Means that freq range 0 to 2 is sampled at 2N points. • If {x[n]}0,N-1 is real, we only need DCT[k] for k = 0 to N-1. • Takes us from 0 to  ( fS/2) in N steps. • (FFT of same order gives only N/2 pts in range 0 to fS/2). Comp30291 : Section 7

  26. DCT (Cont) • Given {x[n]}0,N-1 the DFT (FFT) of this signal extended symmetrically (front-wise) is: For real signal, we might only go from k=0 to k=N-1, and we may evaluate only the modulus or modulus squared. Then, DCT gives us modulus spectrum or ESD without complex numbers. Comp30291 : Section 7

  27. Evaluating the DCT Direct implementation summing cosines will be slow. Pre-extending & performing an FFT may be faster There are faster algorithms. MATLAB has X=dct(x); Also inverse: x= dct(X); Comp30291 : Section 7

  28. 2-Dimensional DCT • DCT is widely used in sound & image processing. • For the latter, a 2D version is needed. • Performed by ‘X=dct2(x)’ in MATLAB • Inverse is available ‘x = idct2(X);’ • More about the 2-d DCT later. Comp30291 : Section 7

  29. 1-dimensional DCT by computed direct method clear all; N=64; for n = 0:N-1 x(1+n) = 100*cos((pi/4)*n) + 100*sin(1.522*n); end; for k=0:N-1 X(k+1)=0; Wk = -k*pi/(2*N); for n=0:N-1 X(k+1) = X(k+1) + x(n+1)*cos((2*n+1)*Wk); end; end; plot(abs(X/(N/2)),'*-' ); grid on; Gives exactly same graph as for ‘even symmetry’ FFT method. This is slower. . But no complex numbers! Comp30291 : Section 7

  30. 1-d DCT by MATLAB clear all; N=64; for n = 0:N-1 x(1+n) = 100*cos((pi/4)*n) + 100*sin(1.522*n); end; X = dct(x); plot(abs(X/(sqrt(N/2)),'*-' ); grid on; Gives same graph, but scaling is different. See MATLAB ‘help’ & documentation. X[k] scaled by [2/N] when k>1 &[2/N] when k=1. A small complication! Comp30291 : Section 7

  31. Applications of DFT (FFT) & DCT • 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. Comp30291 : Section 7

  32. Other applications • Many other uses in signal processing, e.g. filtering. • To filter out a range of frequencies, • perform DFT or DCT, • set the unwanted frequency samples to zero • perform inverse DFT or DCT. • DCT used in MP3 compression to remove components we can’t hear anyway, to save bits. • Illustrate by using the DCT & IDCT to remove the higher frequency sine-wave. Comp30291 : Section 7

  33. DCT to filter out sine-wave clear all; N=64; for n=0:N-1 x(1+n)=100*cos((pi/4)*n)+100*sin(1.522*n); end; X = dct(x); for k=25:40 X(k) = 0; end; y=idct(X); plot(y); grid on; There’s no sign of the higher frequency sine-wave, but the ‘edge’ effects at the beginning & end need to be improved. As we are using ‘dct’ & corresponding ‘idct’ we don’t have to worry abt the funny scaling. Comp30291 : Section 7

  34. Use of FFT to spectrally analyse music Next slide (Table 4 in notes) 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. Comp30291 : Section 7

  35. % Table 4: MATLAB program (dmp7t4.m) to analyse a music file clear all; N = 512 ; % DFT order; [music, fs, nbits] = wavread('cap4th.wav'); L=length(music); for frame = 10 : 20 for n =0 : N-1 x(1+n) = music( 1+ (frame-1)*N + n) ; end; figure (1); plot( x); grid on; X=fft(x); figure(2); plot( abs(X(1:N/2)) ); grid on; SOUND([x x x],fs,nbits); % Play it 3 times. pause; % Press return for next frame end; Comp30291 : Section 7

  36. Musical (violin) note & its mag-spectrum Magnitude spectrum shows fundamental & harmonics Time-domain: (512 samples) Comp30291 : Section 7

  37. Spectral analysis of speech • 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 (next slide) & spectrally analysed using the FFT. • Meaningless to analyse a large speech file at once. • Divide into blocks of samples & analyse separately. • Blocks of N (= 512) samples may be read in and displayed. • Programs in notes available on: ~barry/mydocs/Comp30291 Comp30291 : Section 7

  38. %MATLAB program (dmp7t5.m)for spectrally analysing speech. clear all; N = 512 ; % DFT order IFid=fopen('operator.pcm','rb'); speech = fread(IFid, 'int16'); H=hann(N); % samples of a Hann window of order N for frame = 25 : 200 for n =0 : N-1 x(1+n) = speech( 1+ (frame-1)*N + n); winx(1+n) = x(1+n)*H(1+n) ; end; figure(1); plot(x); X=fft(winx); figure(2); plot( 20*log10(abs(X(1:N/2) )) ); ylim([40 100]); SOUND(x/4000,8000,16); % listen to the frame pause; % Press return for next frame end; fclose('all'); Comp30291 : Section 7

  39. Spectrum of a segment of voiced male speech Freq-domain (dB against freq point) Time-domain(Volts against sample no. Comp30291 : Section 7

  40. Spectrum of another segment of voiced male speech formants Freq-domain (dB against freq point) Time-domain(Volts against sample no. Comp30291 : Section 7

  41. Spectrum of a segment of voiced female speech Freq-domain (dB against freq point) Time-domain(Volts against sample no. Comp30291 : Section 7

  42. Comments on speech graphs • Periodicity seen in time-domain for voiced speech (vowels) • Mag-spectrum has fundamental & many harmonics • Can measure fundamental to determine pitch of the voice. • Male has lower fundamental than female speech • Unvoiced speech (consonants) has no or less periodicity. • ‘Formants’ seen as peaks in spectral envelope. • Caused by vocal tract resonance. • Can determine vowel sound (a e i o u etc) from formants. • In principle, can do speech recognition this way. • Bit-rate compression based on same understanding. Comp30291 : Section 7

  43. Apply 2-d DCT to a picture clear all; rgbpic = imread('autumn.tif'); imview(rgbpic); % input image & display bwpic = rgb2gray(rgbpic); % Convert to gray scale imview(bwpic); imwrite(bwpic,'bwaut.tif','tif'); % Display & store image BWspectrum = dct2(bwpic); % Apply dct figure(1); imshow(log(abs(BWspectrum)),[]), colormap(jet), colorbar; BWspectrum(abs(BWspectrum)<10) = 0.001; % Make zero if <10 figure(2); imshow(log(abs(BWspectrum)),[]), colormap(jet), colorbar; reconpic = idct2(BWspectrum); % Apply inverse dct imview(reconpic,[0 255]); imwrite(bwpic,'bwreconaut.tif','tif'); % Display & store Comp30291 : Section 7

  44. Original picture ‘autumn.tif’ Comp30291 : Section 7

  45. Original & reconstructed images with DCT spectra Reconstructed Original Comp30291 : Section 7

  46. Comments on pictures • Previous program takes a coloured picture & converts it to gray scale. • Then takes DCT & plots 3-D mag-spectrum with colour scale. • Notice a concentration of energy in top corner. • We set to zero any energy values <10. • Useful for image compression as we don’t have to code these. • Then go back to an image via an inverse dct. • We can see the reconstructed image & its modified spectrum (with lots of blue). • Any perceivable loss of quality? Comp30291 : Section 7

  47. n 0 N • Windowing • Given {x[n]}0,N-1 {x[0], x[1], ..., x[N-1] } • Assumed obtained by multiplying infinite sequence {x[n]} by non-symmetric ‘rectangular window’ {r[n]} r[n] Comp30291 : Section 7

  48. r[n] n 0 N • Windowing: symmetric & non-symmetric • Defn of rect window changed from that used in Section 4. • Now, non-symmetric window: In Section 4, symmetric window: rM+1[n] n -M/2 M/2 Has N non-zero samples Has M+1 non-zero samples Comp30291 : Section 7

  49. DTFT of non-symmetric rectangular window (Details of calculation omitted) Comp30291 : Section 7

  50. DTFT of non-symmetric rect window DTFT is now not purely real because of e-j(N-1)/2 term. Mostly interested in magnitude (modulus) which is: Comp30291 : Section 7

More Related