太阳黑子MATLAB

附参考程序：
装载太阳黑子的数据程序如下：
load sunspot.dat; %这组数据记录了在过去的288年（1700-1987年）间每年太
阳黑子出现的数量和大小的观测数据。

绘制Wolfer图的程序如下：
year=sunspot(:,1); wolfer=sunspot(:,2);
plot(year,wolfer);
xlabel('Years');ylabel(' Sunspot Data ')
title('Sunspot Data')
求Wolfer数的FFT程序如下：
Y = fft(wolfer);
在复平面上绘制由Y给出的傅立叶系数的分布图程序如下：
Y(1)=[ ];
plot(Y,'ro')
title('Fourier Coefficients in the Complex Plane');
xlabel('Real Axis');
ylabel('Imaginary Axis');
set(gca,'YTick',[],'XTick',[]);
Wolfer数的周期图程序如下：
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')
% 加载.据
load sunspot.dat
year=sunspot(:,1);
relNums=sunspot(:,2);
plot(year,relNums)
title(‘太阳黑子wolfer图’)
xlabel(‘年份’);
ylabel(‘太阳黑子数量’ );
year=sunspot(:,1);
%从数据中读取年份的相关信息
wolfer=sunspot(:,2);
%从数据中读取关于黑子出现数量的相关信息plot(year,wolfer,'-c');
%绘制太阳黑子Wolfer图
title('太阳黑子Wolfer图');
grid on;
hold on;
% 前五⼗十年的数据
figure(1);
plot(year(1:50),relNums(1:50),'c.-');
xlabel('Year')
hold on;
% 将头⼗一个Y的值省略
Y = fft(relNums);
Y(1)=[];
%在复平⼗面中进.行复平⼗面
figure(2);
plot(Y,'bo')
title('在复平⼗面傅里叶系数’);
xlabel('Real Axis');
ylabel('Imaginary Axis');
hold on;
Y = fft(wolfer);
Y(1)=[ ];
%去掉Y的第一个数据
plot(Y,'yo')
%在复平面上绘制由Y给出的傅立叶系数的分布图 title('傅里叶系数分布图');
xlabel('实数轴');
ylabel('虚数轴');
set(gca,'YTick',[],'XTick',[]);
%将功率⼗比频率定义为”preiodogram”,也就是周期图，并显示。

n=length(Y);
power = abs(Y(1:floor(n/2))).^2;
nyquist = 1/2;
freq = (1:n/2)/(n/2)*nyquist;
figure(3);
plot(freq,power,'k')
xlabel('cycles/year')
title('Periodogram-周期图')
hold on;
Y = fft(wolfer);
Y(1)=[ ];
%去掉Y的第一个数据
plot(Y,'yo')
%在复平面上绘制由Y给出的傅立叶系数的分布图 title('傅里叶系数分布图');
xlabel('实数轴');
ylabel('虚数轴');
set(gca,'YTick',[],'XTick',[]);
%数据大小略大，调成40年的数据
figure(4);
plot(freq(1:30),power(1:30),'m')
xlabel('周期/年')
title('前30年的周期图')
hold on;
%再画出功率⼗周期（频率倒数）的图来指明周期大概是11年
period=1./freq;
figure(5);
plot(period,power);
axis([0 30 0 2e+7]);
ylabel('太阳黑子爆发能量');
xlabel('爆发周期/年');
title('黑子爆发周期图’)
period=1./freq;
plot(period,power,'-k');
axis([0 50 0 2.0e+7]);
ylabel('太阳黑子爆发能量')
xlabel('爆发周期/年份')
title('黑子爆发周期图')
grid on;
hold on;
%通过星形来指明11年的具体值，得出结论
index=find(power==max(power));
mainPeriodStr=num2str(period(index));
plot(period(index),power(index),'r*', 'MarkerSize',20);
text(period(index)+2,power(index),['周期= ',mainPeriodStr]); hold off;。