sunspots_demo.m

数字信号处理（机械工业出版社）的源码

% This demonstration uses the FFT function to analyze the variations in sunspot
% activity over the last 300 years.
% Sunspot activity is cyclical, reaching a maximum about every 11 years.  Let's
% confirm that.  Here is a plot of a quantity called the Wolfer number, which
% measures both number and size of sunspots.  Astronomers have tabulated this
% number for almost 300 years.

load sunspot.dat
year=sunspot(:,1); 
wolfer=sunspot(:,2);
figure
plot(year,wolfer)
xlabel('Years'); ylabel(' Sunspot Data '); title('Sunspot Data')
pause

figure
plot(year(1:50),wolfer(1:50),'b.-');
xlabel('Years');ylabel(' Sunspot Data '); title('At the first 50 years')
Y = fft(wolfer);
Y(1)=[];

pause
figure
plot(Y,'ro')
title('Fourier Coefficients in the Complex Plane');
xlabel('Real Axis');
ylabel('Imaginary Axis');
pause

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

figure
plot(freq(1:40),power(1:40))
xlabel('cycles/year')
pause

figure
period=1./freq;
plot(period,power);
axis([0 40 0 2e+7]);
ylabel('Power');
xlabel('Period (Years/Cycle)');
pause

figure
hold on;
index=find(power==max(power));
mainPeriodStr=num2str(period(index));
plot(period(index),power(index),'r.', 'MarkerSize',25);
text(period(index)+2,power(index),['Period = ',mainPeriodStr]);
hold off;