home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 3.ddi / DEMOS.DI$ / FFTDEMO.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  2.1 KB  |  86 lines
  1. %FFTDEMO Demonstrate use of FFT for spectral analysis.
  2.  
  3. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  4.  
  5. echo on
  6. clc
  7. %    This example shows the use of the FFT function for spectral 
  8. %    analysis.  A common use of FFT's is to find the frequency
  9. %    components of a signal buried in a noisy time domain signal.
  10. %
  11.  
  12. pause % Strike any key to continue.
  13. clc
  14. %    First we need to create some data.  Consider data sampled at
  15. %    1000 Hz.  We start by forming a time axis for our data, running
  16. %    from t=0 until t=.25 in steps of 1 millisecond:
  17.  
  18. t = 0:.001:.25;
  19.  
  20. pause % Strike any key to continue.
  21.  
  22. %    Next we can form a signal containing 50 Hz and 120 Hz:
  23.  
  24. x = sin(2*pi*50*t) + sin(2*pi*120*t);
  25.  
  26. pause % Strike any key to continue.
  27.  
  28. %    and add some random noise with a standard deviation of 2 to 
  29. %    produce a noisy signal y:
  30.  
  31. y = x + 2*randn(size(t));
  32.  
  33. pause % Strike any key to continue.
  34. x=[];t=[];
  35. clc
  36. %    Let's take a look at our noisy signal y by plotting it.
  37.  
  38. pause % Strike any key for plot.
  39.  
  40. plot(y(1:50)), title('Noisy time domain signal'), pause
  41. clc
  42. %    Clearly, it is difficult to identify the frequency components
  43. %    from looking at the original signal; that's why spectral analysis
  44. %    is so popular.
  45. %
  46. %    Finding the discrete Fourier transform of the noisy signal y
  47. %    is easy; we just take the fast-Fourier transform (FFT) :
  48.  
  49. Y = fft(y,256);
  50.  
  51. pause % Strike any key to continue.
  52.  
  53. %    The power spectral density, a measurement of the energy at
  54. %    various frequencies, is found with:
  55.  
  56. Pyy = Y.*conj(Y)/256;
  57.  
  58. pause % Strike any key to continue.
  59. Y=[];y=[];
  60. clc
  61. %    To plot the power spectral density, we must first form a 
  62. %    frequency axis:
  63.  
  64. f = 1000/256*(0:127);
  65.  
  66. %    which we do for the first 127 points. (The remainder of the 256
  67. %    points are symmetric.)  We can now plot the power spectral
  68. %    density:
  69.  
  70. pause % Strike any key for plot.
  71.  
  72. plot(f,Pyy(1:128)), title('Power spectral density'), ...
  73. xlabel('Frequency (Hz)'), pause
  74. clc
  75. %    Let's zoom in and plot only up to 200 Hz:
  76.  
  77. pause % Strike any key for plot.
  78.  
  79. plot(f(1:50),Pyy(1:50)), title('Power spectral density'), ...
  80. xlabel('Frequency (Hz)'), pause
  81.  
  82. clc
  83. echo off
  84.  
  85. disp('End')
  86.