www.pudn.com > NumericalComputingwithMatlabCode.zip > sunspotstx.m, change:2003-12-30,size:3112b

% SUNSPOTSTX  Updated version of toolbox/matlab/demos/sunspots.m 
% For centuries people have noted that the face of the sun is not 
% constant or uniform in appearance, but that darker regions appear 
% at random locations on a cyclical basis.  This activity is correlated 
% with weather and other economically significant terrestrial phenomena. 
% In 1848, Rudolf Wolfer proposed a rule that combined the number and 
% size of these sunspots into a single index.  Using archival records, 
% astronomers have applied Wolfer's rule to determine sunspot activity 
% back to the year 1700.  Today the sunspot index is measured by many 
% astronomers and the worldwide distribution of the data is coordinated 
% by the Solar Influences Data Center at the Royal Observatory of Belgium. 
% The text file sunspot.dat in the matlab/demos directory has two columns 
% of numbers.  The first column is the years from 1700 to 1987 and the 
% second column is the average Wolfer sunspot number for each year. 
    load sunspot.dat 
    t = sunspot(:,1)'; 
    wolfer = sunspot(:,2)'; 
    n = length(wolfer); 
% There is a slight upward trend to the data.  A least squares fit 
% gives the trend line. 
    c = polyfit(t,wolfer,1); 
    trend = polyval(c,t); 
    plot(t,[wolfer; trend],'-',t,wolfer,'k.') 
    ylabel('Wolfer index') 
    title('Sunspot index with linear trend') 
% You can definitely see the cyclic nature of the phenomenon. 
% The peaks and valleys are a little more than 10 years apart. 
% Now, subtract off the linear trend and take the finite Fourier transform. 
    y = wolfer - trend; 
    Y = fft(y); 
% The vector |Y|.^2^2 is the power in the signal. 
% A plot of power versus frequency is a periodogram. 
% We prefer to plot |Y|, rather than |Y|.^2, because the 
% scaling is not so exaggerated. 
% The sample rate for this data is one observation per year, 
% so the frequency #f# has units of cycles per year. 
    Fs = 1;  % Sample rate 
    f = (0:n/2)*Fs/n; 
    pow = abs(Y(1:n/2+1)); 
    pmax = 5000; 
    plot([f; f],[0*pow; pow],'c-', f,pow,'b.', ... 
    axis([0 .5 0 pmax]) 
% The maximum power occurs near frequency = 0.09 cycles/year. 
% We would like to know the corresponding period in years/cycle. 
% Let's zoom in on the plot and use the reciprocal of frequency 
% to label the x-axis. 
    k = 0:44; 
    f = k/n; 
    pow = pow(k+1); 
    plot([f; f],[0*pow; pow],'c-',f,pow,'b.', ... 
    axis([0 max(f) 0 pmax]) 
    k = 2:3:41; 
    f = k/n; 
    period = 1./f; 
    periods = sprintf('%5.1f|',period); 
    title('Periodogram detail') 
% As expected, there is a very prominent cycle with a length of 
% about 11.1 years.  This shows that over the last 300 years, 
% the period of the sunspot cycle has been slightly over 11 years.