nyquist= 1/2;
freq = (1:n/2)/(n/2)*nyquist;%%求频率
plot(freq,power)%%作出功率与频率的关系曲线
xlabel('cycles/year')%%标注横坐标
title('Periodogram')
（2）对功率和频率的前100个分量作它的周期图，方便观察。
Y=fft(wolfer); %%对全部数据做FFT
Y(1)=[ ]; %%舍弃第一个点，因为Y(1)为所有FFT数值之和
plot(Y,'ro') %%在复平面作图，且为空心点
>> figure%%作图
plot(year(1:100),wolfer(1:100),'b.-');%%取年份1700-1800
xlabel('Years');ylabel(' Sunspot Data ');%%标注横纵坐标
title('At the first 100 years ')
分析：由上图可以更加清楚地看到周期大约为10到12年之间。
[2]孙燕. 用数字信号处理方法分析太阳黑子活动的周期性[J]. 青海师范大学学报（自然科学版）, 2006, (4): 1-3
[3]唐洁. 功率谱分析方法在周期分析中的应用[J]. 陕西理工学院学报, 2013, 29(5): 1-5
七、MATLAB程序源代码
>> load sunspotyear.csv %%装载太阳黑子的数据
用pwd命令来查找当前工作路径。
>>pwd
ans=
C:\Program Files\MATLAB\R2012b\bin
下面的分析选择1700年到2013年期间，共313年的数据。首先用如下程序语句装载太阳黑子的年度数据：
>> load sunspotyear.csv%%装载太阳黑子的数据
2.以横坐标表示年份，纵坐标表示黑子出现的数量，绘制Wolfer图。
（1）由1700年到2013年期间Wolfer图。
>> year=sunspotyear(:,1);%%读取年份信息
wolfer=sunspotyear(:,2);%%读取太阳黑子活动数据
plot(year,wolfer);%%作以year为横坐标，wolfer为纵坐标的图象
xlabel('Years');ylabel(' Sunspot Data ')%%标注横坐标为年份，纵坐标为黑子数量
title('Sunspot Data')
分析：由上图可见，太阳黑子数据具有周期性，但是仅根据上图周期的准确值难以确定。不过可以确定上图中任意两个峰值之间的年份之差，从而获得一个周期的估计值。比如第3个峰值出现在1728年，第4个峰值出现在1739年，因此周期的近似值是11年。
（2）由1700到1800年期间Wolfer图。
>>figure
plot(freq(1:100),power(1:100))
xlabel('cycles/year')
pause
分析：图象中出现的其它尖峰可能是因为太阳黑子数据中的噪声引起的；或者是因为泄露产生的。
（3）确定太阳黑子活动周期
确定出太阳黑子的活动周期。为清楚起见，画出功率与周期（频率的倒数）的关系曲线图。在功率与周期关系曲线图中标出功率的最高点.
程序如下：
figure
period=1./freq;%%求周期
plot(period,power);%%作功率与周期的关系
axis([0 50 0 2e+7]);%%确定横纵坐标的范围
ylabel('Power');%%标注纵坐标为功率
xlabel('Period (Years/Cycle)');%%标注横坐标为周期
text(period(index)+2,power(index),['Period = ',mainPeriodStr]); %%文字标注该点
holdoff;
由上图可以近似得出太阳黑子活动周期为11年。
四、结论及优缺点改进
实验得出的太阳黑子的活动周期为11.0385年,与沃尔夫得出的11年的周期规律一致，
pause
index=find(power==max(power)); %%找到频率最大点，该点横坐标即为太阳黑子周期
mainPeriodStr=num2str(period(index));
plot(period(index),power(index),'r.', 'MarkerSize',25); %%用实心点指出该点
22年模型：美国天文学家乔治•海尔【George Ellery Hale(1868 - 1938)】发现在每一个施瓦贝周期，太阳的磁场都会扭转，因此磁极要两次扭转之后才会回到相同磁极的状态，他因此发现了22年周期。
其他曾经发现的模型：87年模型、210年模型、哈尔斯塔周期等。
在太空运行的人造卫星、飞行器、宇宙飞船和高空飞行的飞机，受到来自太阳的高能带电粒子的袭击，会造成有些零部件损坏导致整个设备系统不能正常工作。乘坐在飞机和飞行器上的飞行员和宇航员的身体也会受到来自太阳的高能带电粒子伤害。受太阳活动影响的还有长距离的高压输电系统和输油管道
则可将DFT化为
利用系数 的可约性，即 ，上式可表示成
（3）功率谱密度： ，
则
3.流程图
三、实验过程及分析
1.下载太阳黑子的数据monthssv.dat，下载数据的时间段从1700年1月到2013年12月，共314年。这个数据文件的第1列是年和月，第2列是该月太阳黑子的平均数，sunspotyear.csv同时保存到matlab工作路径。
根据DFT的频域单位K与DTFT的频域单位w的关系表达式 以及W与f对应
关系 （ 为采样频率），可以看出K与f成线性关系 。又因为前 个数据已经包含了Wolfer数的全部信息，只取前 个数据分析功率-频率图时，对应的横坐标应取 。
2.函数名称
（1）离散时间傅里叶正变换：
离散时间傅里叶逆变换：
（2）按时间抽选的FFT算法：将N= 的序列 先按n的奇偶分成两组：
xlabel('Real Axis');%%标注横坐标为实轴
ylabel('Imaginary Axis&#模的平方被定义为功率，功率与频率的关系曲线则被定义为周期图）。
（1）绘制周期图
>> n=length(Y);
power= abs(Y(1:n/2)).^2;%%求功率
我们利用Matlab强大的数据处理与仿真功能，对Wolfer数进行功率谱密度分析从而可以得到对太阳黑子活动周期的结论。
二、分析方法
1.实验原理
在对太阳黑子活动周期的分析中，对Wolfer数序列作FFT变换后得到Y（长度为n），只
取前 个数据的功率谱密度的估计值 。原因是时序为离散的实序列的傅里叶变换对应于具有周期性且偶对称的频域特性，因此Y的前 个数据已经包含了Wolfer数的全部信息。
year=sunspotyear(:,1); %%读取年份信息
wolfer=sunspotyear(:,2); %%读取太阳黑子活动数据
plot(year,wolfer); %%作以year为横坐标，wolfer为纵坐标的图象
xlabel('Years');ylabel(' Sunspot Data ') %%标注横坐标为年份，纵坐标为黑子数量
作业1太阳黑子活动周期的分析
一、问题描述
太阳黑子的活动是周期的，大约每11年达到一个爆发高峰。
太阳周期是太阳行为上的循环变化，人们已经构建了许多可能的太阳黑子活动模型，
但在观测上只有11年和22年的周期容易被清除地观察到。
11年模型：1843年，塞瑟尔•海因里希•施瓦布【Samuel HeinrichSchwabe(1789–1875)】发现了太阳活动周期。从1826年到1843年，施瓦贝每天仔细观察看太阳表面，记录太阳上的黑子数，经过17年间的长期艰辛观测，他整理了观测资料，于1843年发表了一篇题为《1843年间的太阳观测》的论文，文章指出：“太阳的年平均黑子数具有周期性的变化，变化的周期约十年”。伯尔尼天文台台长的鲁道夫沃尔夫在搜集整理太阳黑子数观测资料的过程中，为使不同观测台站以及不同人的太阳黑子观测资料具有可比性，于1848年提出了太阳黑子相对数的概念。沃尔夫经过几年的仔细观测和精心的资料整理，返现太阳黑子数变化周期平均为11.1年。观测到的最短黑子周期为9年，最长黑子周期为14年。不同周期之间和字数的变化也非常明显。到1852年他还发现地磁活动和极光与太阳活动有关。沃尔夫提出将太阳黑子数从一个极小到另一个极小之间的事件定为一个周期，并将1755年之1766年的周期定为第一个太阳活动周。根据连续的观测记录推算下来，2008是第24个太阳活动周。
title('Sunspot Data')
figure %%作图
plot(year(1:100),wolfer(1:100),'b.-'); %%取年份1700-1800
xlabel('Years');ylabel(' Sunspot Data '); %%标注横纵坐标
title('At the first 100 years ')
总体而言，这次实验虽只是对信号的简单处理，但是在程序设计中遇到很多问题，这些都是在平时学习中所不曾注意到的，面对这些问题我通过查阅资料、网络、与同组同学的探讨，得到很大程度的解决。