sunspot.m

来自「利用1700年以来持续300年的太阳黑子活动数据进行周期估算」· M 代码 · 共 38 行

load sunspot.dat
year=sunspot(:,1);
wolfer=sunspot(:,2);
plot(year,wolfer)
title('太阳黑子的数据')
figure;
plot(year(1:50),wolfer(1:50),'b.-');
Y = fft(wolfer);
Y(1)=[];
figure;
plot(Y,'ro')
title('复平面上的Fourier系数');
xlabel('实轴');
ylabel('虚轴');
n=length(Y);
power = abs(Y(1:floor(n/2))).^2;
nyquist = 1/2;
freq = (1:n/2)/(n/2)*nyquist;
figure;
plot(freq,power)
xlabel('周期/年')
title('周期图')
figure;
plot(freq(1:40),power(1:40))
xlabel('周期/年')
period=1./freq;
figure;
plot(period,power);
axis([0 40 0 2e+7]);
ylabel('功率');
xlabel('周期 (年/周期)');
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),['周期 = ',mainPeriodStr]);
hold off;