%SUNSPOTS The answer is 11.08, what is the question?
echo off
%	CBM+LS, 8-27-92.
%	Copyright (c) 1984-93 by The MathWorks, Inc.

echo on
clf
clc
% You are probably aware that sunspot activity is cyclical,
% reaching a maximum about every 11 years.  Let's confirm that.
% We will analyze a quantity known as the Wolfer number, which
% measures both number and size of sunspots, and which astronomers
% have tabulated for almost 300 years.
%
% See chapter 11 of "Numerical Methods and Software, by
% D. Kahaner, C. Moler and S. Nash, Prentice-Hall, 1989, 
% for a more detailed discussion.

pause % Press any key ...
clc
% The demo data file, sunspot.dat, contains the data.

load sunspot.dat
year = sunspot(:,1);
wolfer = sunspot(:,2);
whos

pause % Press any key ...
clc
% Plot Wolfer number versus year.

subplot(2,1,1)
plot(year,wolfer)
title('Sunspot Data')

pause % Press any key ...

% And take a closer look at the first 50 years.

subplot(2,1,2)
plot(year(1:50),wolfer(1:50),'-',year(1:50),wolfer(1:50),'go');

pause % Press any key ...
clf
clc
% The fundamental tool of signal processing is the FFT, 
% the fast Finite Fourier Transform.

Y = fft(wolfer);
n = length(Y);

pause % Press any key ...
clc
% Here are the first few Fourier coefficients.  They are complex.

Y(1:5)

% The first component, Y(1), is simply the sum of the data
% and can be removed.

Y(1) = [];

pause % Press any key ...
clc
% A graph of the distribution of the Fourier coefficients
% in the complex plane is pretty, but difficult to interpret.

plot(Y,'co')

pause % Press any key ...
clc
% The complex magnitude squared of Y is called the power,
% and a plot of power versus frequency is a "periodogram".

power = abs(Y(1:n/2)).^2;
nyquist = 1/2;
freq = (1:n/2)/(n/2)*nyquist;
subplot(2,1,1)
plot(freq,power)
xlabel('cycles/year')
title('Periodogram')

pause % Press any key ...

% Zoom in on the first 40 components.

subplot(2,1,2)
plot(freq(1:40),power(1:40))
xlabel('cycles/year')

pause % Press any key ...
clc
% Finally, we use some specialized features of MATLAB's
% "Handle Graphics" to change the labeling on the x-axis
% in the bottom graph from frequency to its reciprocal.

f = get(gca,'xtick');
c = zeros(length(f),5);
for i = 1:length(f)
   c(i,1:5) = sprintf('%5.1f',1/f(i));
end
set(gca,'xticklabels',c)
xlabel('years/cycle')

pause % Press any key ...
clc
% As expected, there is a very prominent cycle with a
% length of about 11 years.

k = find(power == max(power))
cycle = n/k

% End
echo off
set(gcf,'nextplot','replace')
