home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 3.ddi / DEMOS.DI$ / SUNSPOTS.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  2.6 KB  |  116 lines
  1. %SUNSPOTS The answer is 11.08, what is the question?
  2. echo off
  3. %    CBM+LS, 8-27-92.
  4. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  5.  
  6. echo on
  7. clf
  8. clc
  9. % You are probably aware that sunspot activity is cyclical,
  10. % reaching a maximum about every 11 years.  Let's confirm that.
  11. % We will analyze a quantity known as the Wolfer number, which
  12. % measures both number and size of sunspots, and which astronomers
  13. % have tabulated for almost 300 years.
  14. %
  15. % See chapter 11 of "Numerical Methods and Software, by
  16. % D. Kahaner, C. Moler and S. Nash, Prentice-Hall, 1989, 
  17. % for a more detailed discussion.
  18.  
  19. pause % Press any key ...
  20. clc
  21. % The demo data file, sunspot.dat, contains the data.
  22.  
  23. load sunspot.dat
  24. year = sunspot(:,1);
  25. wolfer = sunspot(:,2);
  26. whos
  27.  
  28. pause % Press any key ...
  29. clc
  30. % Plot Wolfer number versus year.
  31.  
  32. subplot(2,1,1)
  33. plot(year,wolfer)
  34. title('Sunspot Data')
  35.  
  36. pause % Press any key ...
  37.  
  38. % And take a closer look at the first 50 years.
  39.  
  40. subplot(2,1,2)
  41. plot(year(1:50),wolfer(1:50),'-',year(1:50),wolfer(1:50),'go');
  42.  
  43. pause % Press any key ...
  44. clf
  45. clc
  46. % The fundamental tool of signal processing is the FFT, 
  47. % the fast Finite Fourier Transform.
  48.  
  49. Y = fft(wolfer);
  50. n = length(Y);
  51.  
  52. pause % Press any key ...
  53. clc
  54. % Here are the first few Fourier coefficients.  They are complex.
  55.  
  56. Y(1:5)
  57.  
  58. % The first component, Y(1), is simply the sum of the data
  59. % and can be removed.
  60.  
  61. Y(1) = [];
  62.  
  63. pause % Press any key ...
  64. clc
  65. % A graph of the distribution of the Fourier coefficients
  66. % in the complex plane is pretty, but difficult to interpret.
  67.  
  68. plot(Y,'co')
  69.  
  70. pause % Press any key ...
  71. clc
  72. % The complex magnitude squared of Y is called the power,
  73. % and a plot of power versus frequency is a "periodogram".
  74.  
  75. power = abs(Y(1:n/2)).^2;
  76. nyquist = 1/2;
  77. freq = (1:n/2)/(n/2)*nyquist;
  78. subplot(2,1,1)
  79. plot(freq,power)
  80. xlabel('cycles/year')
  81. title('Periodogram')
  82.  
  83. pause % Press any key ...
  84.  
  85. % Zoom in on the first 40 components.
  86.  
  87. subplot(2,1,2)
  88. plot(freq(1:40),power(1:40))
  89. xlabel('cycles/year')
  90.  
  91. pause % Press any key ...
  92. clc
  93. % Finally, we use some specialized features of MATLAB's
  94. % "Handle Graphics" to change the labeling on the x-axis
  95. % in the bottom graph from frequency to its reciprocal.
  96.  
  97. f = get(gca,'xtick');
  98. c = zeros(length(f),5);
  99. for i = 1:length(f)
  100.    c(i,1:5) = sprintf('%5.1f',1/f(i));
  101. end
  102. set(gca,'xticklabels',c)
  103. xlabel('years/cycle')
  104.  
  105. pause % Press any key ...
  106. clc
  107. % As expected, there is a very prominent cycle with a
  108. % length of about 11 years.
  109.  
  110. k = find(power == max(power))
  111. cycle = n/k
  112.  
  113. % End
  114. echo off
  115. set(gcf,'nextplot','replace')
  116.