一、 引言太阳黑子是人们最早发现也是人们最熟悉的一种太阳表面活动。
因为太阳内部磁场发生变化，太阳黑子的数量并不是固定的，它会随着时间的变化而上下波动，每隔一定时间会达到一个最高点，这段时间就被称之为一个太阳黑子周期。
太阳黑子的活动呈现周期性变化是由施瓦贝首次发现的。
沃尔夫 (R.Wolfer)继而推算出11年的周期规律。
实际上，太阳黑子的活动不仅呈11年的周期变化，还有海耳在研究太阳黑子磁场分布时发现的22年周期；格莱斯堡等人发现的80年周期以及蒙德极小期等。
由于太阳黑子的活动规律极其复杂，时至今日科学家们仍在努力研究其内在的规律和特性。
事实上，对太阳黑子活动规律的研究不仅具有理论意义，而且具有直接的应用需求。
太阳黑子的活动呈现周期性变化的，沃尔夫(R.Wolfer)根据在过去的288 年(1700年~1987 年)间每年太阳黑子出现的数量和大小的观测数据推算出11 年的周期规律。
我们利用Matlab 强大的数据处理与仿真功能，对Wolfer 数进行功率谱密度分析从而可以得到对太阳黑子活动周期的结论。
二、 实验原理在该实验中，对Wolfer 数序列做FFT 变换后得到Y （长度为n ），只取其前n2个数据的功率谱密度的估计值2|Y |。
原因是时域为离散的实序列的傅立叶变换对应于具有周期性且偶对称的频域特性，因此Y 的前n2个数据已经包含了Wolfer 数的全部信息。
根据DFT 的频域单位k 与DTFT 的频域单位ϖ的表达关系式2k Nπϖ=⋅以及ϖ与f 对应关系=2s f f ϖπ⋅，可以看出k 与f 呈线性关系s f f k N=⋅；同样地，因为Y 的前n2个数据已经包含了Wolfer 数的全部信息，只取前n2个数据分析功率-频率图时，对应的横坐标时应取2s f f k N =⋅。
（注：s f 为采样频率）三、方案组成与流程图1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111 1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111四、实现步骤及内容(1) 读取数据并绘制Wolfer 图，观察太阳黑子的活动现象。
程序如下： load sunspot_sun.dat %读取数据 year=sunspot(:,1); %读取年份信息wolfer=sunspot(:,2); %读取黑子活动数据 figure%新建图像plot(year,wolfer) %画出时域图xlabel('Years'); ylabel(' Sunspot Data '); title('Sunspot Data') %标注横纵坐标和标题pause1700175018001850190019502000020406080100120140160180200YearsS u n s p o t D a t aSunspot Data为清楚起见，画出部分时间段(如实验中采取开始50年)的Wolfer 图进行观察。
程序如下：figureplot(year(1:50),wolfer(1:50),'b.-'); %用前50年数据做图xlabel('Years');ylabel(' Sunspot Data '); title('At the first 50 years') %标注横纵坐标和标题YearsS u n s p o t D a t aAt the first 50 years(2) 对已经得到的Wolfer 数应用FFT 技术分析它的变化规律。
程序如下：Y = fft(wolfer); %对全部数据做FFTY(1)=[]; %由于算法问题，Y(1)为所有fft 数值之和，因此需舍弃第一个点 pause观察Wolfer 数在FFT 分析后得到的在复平面中的分布图。
程序如下： figureplot(Y,'ro') %在复平面做图，空心点title('Fourier Coefficients in the Complex Plane'); xlabel('Real Axis'); %标注横坐标为’实部’ylabel('Imaginary Axis'); %标注纵坐标为’虚部’ pause-4000-3000-2000-10000100020003000-2500-2000-1500-1000-50005001000150020002500Fourier Coefficients in the Complex PlaneReal AxisI m a g i n a r y A x i s(3) 功率与频率的关系曲线则被定义为周期图。
用周期图法对 Wolfer 数进行功率谱密度分析。
程序如下：figuren=length(Y);%取FFT 结果长度为npower = abs(Y(1:n/2)).^2; %取前n/2个数据求其功率 nyquist = 1/2; %取最大频率为0.5freq = (1:n/2)/(n/2)*nyquist; %将FFT 转换为DTFT plot(freq,power) %画周期图xlabel('cycles/year');title('Periodogram') pause00.050.10.150.20.250.30.350.40.450.50.20.40.60.811.21.41.61.827cycles/yearPeriodogram为清楚起见，对功率和频率的前50个分量作它的周期图。
程序如下： figureplot(freq(1:50),power(1:50)) %为观察方便，取前50年数据做图 xlabel('cycles/year') pause00.020.040.060.080.10.120.140.160.180.20.40.60.811.21.41.61.82cycles/year(4) 确定出太阳黑子的活动周期。
为清楚起见，画出功率与周期（频率的倒数）的关系曲线图。
程序如下：figureperiod=1./freq; %将横坐标单位转换为年/周期 plot(period,power);%画周期图 axis([0 50 0 2e+7]); ylabel('Power');xlabel('Period (Years/Cycle)');%单位为年/周期 pause510152025303540455000.20.40.60.811.21.41.61.82P o w e rPeriod (Years/Cycle)(5) 在功率与周期关系曲线图中标出功率的最高点，该位置对应的周期则为太阳黑子活动周期。
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;510152025303540455000.20.40.60.811.21.41.61.82P o w e rPeriod (Years/Cycle)五、 实验结果分析与总结实验得出的太阳黑子的活动周期为11.0385年,与沃尔夫得出的11年的周期规律一致，说明实验结果正确。