快速傅立叶变换及应用_本科毕业论文.doc

毕业(设计)论文题 目 快速傅里叶变换及其应用 学生姓名 辛鹏宇 专业班级 R计算081班 所在院系 理 学 院 指导教师 职 称副教授 所在单位 理 学 院 教研室主任 完成日期 2013 年6月18日摘 要快速傅氏变换（FFT），是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。它对傅氏变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。傅里叶变换的理论与方法在“数理方程”、“线性系统分析”、“信号处理、仿真”等很多学科领域都有着广泛应用,由于计算机只能处理有限长度的离散的序列,所以真正在计算机上运算的是一种离散傅里叶变换。虽然傅里叶运算在各方面计算中有着重要的作用，但是它的计算过于复杂，大量的计算对于系统的运算负担过于庞大，使得一些对于耗电量少，运算速度慢的系统对其敬而远之，然而，快速傅里叶变换的产生，使得傅里叶变换大为简化，在不牺牲耗电量的条件下提高了系统的运算速度，增强了系统的综合能力，提高了运算速度，因此快速傅里叶变换在生产和生活中都有着非常重要的作用，对于学习掌握都有着非常大的意义。关键字： 快速傅氏变换；快速算法；简化；广泛应用ABSTRACTFast Fourier Transform (FFT), is a discrete fast Fourier transform algorithm, which is based on the Discrete Fourier Transform of odd and even, false, false, and other characteristics of the Discrete Fourier Transform algorithms improvements obtained. Its Fourier transform theory has not found a new, but in the computer system or the application of digital systems Discrete Fourier Transform can be said to be a big step into. Fourier transform theory and methods in the mathematical equation and linear systems analysis and signal processing, simulation, and many other areas have a wide range of applications, as the computer can only handle a limited length of the sequence of discrete, so true On the computers operation is a discrete Fourier transform. Fourier Although all aspects of computing in the calculation has an important role, but its calculation was too complicated, a lot of computing system for calculating the burden is too large for some Less power consumption, the slow speed of operation of its system at arms length, however, have the fast Fourier transform, Fourier transform greatly simplifying the making, not in power at the expense of the conditions to increase the speed of computing systems, and enhance the system The comprehensive ability to improve the speed of operation, the Fast Fourier Transform in the production and life have a very important role in learning to master all have great significance.Key Words： Fast Fourier Transform; fast algorithm; simplified; widely used大连交通大学2013届本科生毕业论文目 录一、快速傅里叶变换原理及性质1（一）快速傅里叶变换原理1（二）快速傅里叶变换的优越性1（三）快速傅里叶变换的意义2二、快速傅里叶变换的算法4（一）快速傅里叶变换算法4三、快速傅里叶变换的应用6（一）利用FFT计算连续时间信号的傅里叶变换6（二）利用FFT计算离散信号的线性卷积9（三）利用FFT进行离散信号压缩11（四）利用FFT对离散信号进行滤波14（五）利用FFT提取离散信号中的最强正弦分量17谢 辞22参考文献23一、快速傅里叶变换原理及性质数字信号的傅里叶变换,通常采用离散傅里叶变换(DFT)方法。DFT 存在的不足是计算量太大,很难进行实时处理。计算一个N 点的DFT ,一般需要次复数乘法和N(N-1)次复数加法运算.因此,当N较大或要求对信号进行实时处理时,往往难以实现所需的运算速度。1965年,J.W.Cooly和J.W.Tukey发现了DFT的一种快速算法,经其他学者进一步改进, 很快形成了一套高效运算方法,这就是现在通用的快速傅里叶变换, 简称FFT( The Fast Fourier Transform)。快速傅里叶变换的实质是利用式(1)中的权函数的对称性和周期性,把N点DFT进行一系列分解和组合,使整个DFT的计算过程变成一系列叠代运算过程,使DFT的运算量大大简化,为DFT及数字信号的实时处理和应用创造了良好的条件。（一）快速傅里叶变换原理快速傅里叶变换原理1. 将长序列DFT分解为短序列的DFT2. 利用旋转因子的周期性、对称性、可约性。将时域序列逐次分解为一组子序列，利用旋转因子的特性，由子序列的DFT来实现整个序列的DFT。其中：快速傅里叶变换分为两种，分为基2时间抽取算法和基2频率抽取算法基2时间抽取(Decimation in time)FFT算法其中:r=0,1,2基2频率抽取(Decimation in frequency)FFT算法 （二）快速傅里叶变换的优越性设为项的复数序列，由DFT变换，任一的计算都需要次复数乘法和次复数加法，而一次复数乘法等于四次实数乘法和两次实数加法，一次复数加法等于两次实数加法，即使把一次复数乘法和一次复数加法定义成一次“运算”（四次实数乘法和四次实数加法），那么求出项复数序列的, 即N点DFT变换大约就需要次运算。当点甚至更多的时候，需要次运算，在FFT中，利用WN的周期性和对称性，把一个N项序列（设,为正整数），分为两个项的子序列，每个点DFT变换需要次运算，再用N次运算把两个点的DFT 变换组合成一个N点的DFT变换。这样变换以后，总的运算次数就变成。继续上面的例子，时，总的运算次数就变成了525312次，节省了大约50%的运算量。而如果我们将这种“一分为二”的思想不断进行下去，直到分成两两一组的DFT运算单元，那么N点的 DFT变换就只需要次的运算，在1024点时，运算量仅有10240次，是先前的直接算法的1%，点数越多，运算量的节约就越大，这就是 FFT的优越性.（三）快速傅里叶变换的意义1.傅立叶变换的物理意义傅立叶变换是数字信号处理领域一种很重要的算法。要知道傅立叶变换算法的意义，首先要了解傅立叶原理的意义。傅立叶原理表明：任何连续测量的时序或信号，都可以表示为不同频率的正弦波信号的无限叠加。而根据该原理创立的傅立叶变换算法利用直接测量到的原始信号，以累加方式来计算该信号中不同正弦波信号的频率、振幅和相位。和傅立叶变换算法对应的是反傅立叶变换算法。该反变换从本质上说也是一种累加处理，这样就可以将单独改变的正弦波信号转换成一个信号。因此，可以说，傅立叶变换将原来难以处理的时域信号转换成了易于分析的频域信号（信号的频谱），可以利用一些工具对这些频域信号进行处理、加工。最后还可以利用傅立叶反变换将这些频域信号转换成时域信号。从现代数学的眼光来看，傅里叶变换是一种特殊的积分变换。它能将满足一定条件的某个函数表示成正弦基函数的线性组合或者积分。在不同的研究领域，傅里叶变换具有多种不同的变体形式，如连续傅里叶变换和离散傅里叶变换。在数学领域，尽管最初傅立叶分析是作为热过程的解析分析的工具，但是其思想方法仍然具有典型的还原论和分析主义的特征。任意的函数通过一定的分解，都能够表示为正弦函数的线性组合的形式，而正弦函数在物理上是被充分研究而相对简单的函数类：1. 傅立叶变换是线性算子,若赋予适当的范数,它还是酉算子;2. 傅立叶变换的逆变换容易求出,而且形式与正变换非常类似;3. 正弦基函数是微分运算的本征函数,从而使得线性微分方程的求解可以转化为常系数的代数方程的求解.在线性时不变杂的卷积运算为简单的乘积运算,从而提供了计算卷积的一种简单手段;4. 离散形式的傅立叶的物理系统内,频率是个不变的性质,从而系统对于复杂激励的响应可以通过组合其对不同频率正弦信号的响应来获取;5. 著名的卷积定理指出:傅立叶变换可以化复变换可以利用数字计算机快速的算出(其算法称为快速傅立叶变换算法(FFT)。正是由于上述的良好性质,傅里叶变换在物理学、数论、组合数学、信号处理、概率、统计、密码学、声学、光学等领域都有着广泛的应用。2.图像傅立叶变换的物理意义图像的频率是表征图像中灰度变化剧烈程度的指标，是灰度在平面空间上的梯度。如：大面积的沙漠在图像中是一片灰度变化缓慢的区域，对应的频率值很低；而对于地表属性变换剧烈的边缘区域在图像中是一片灰度变化剧烈的区域，对应的频率值较高。傅立叶变换在实际中有非常明显的物理意义，设f是一个能量有限的模拟信号，则其傅立叶变换就表示f的谱。从纯粹的数学意义上看，傅立叶变换是将一个函数转换为一系列周期函数来处理的。从物理效果看，傅立叶变换是将图像从空间域转换到频率域，其逆变换是将图像从频率域转换到空间域。换句话说，傅立叶变换的物理意义是将图像的灰度分布函数变换为图像的频率分布函数，傅立叶逆变换是将图像的频率分布函数变换为灰度分布函数。傅立叶变换以前，图像（未压缩的位图）是由对在连续空间（现实空间）上的采样得到一系列点的集合，我们习惯用一个二维矩阵表示空间上各点，则图像可由z=f(x,y)来表示。由于空间是三维的，图像是二维的，因此空间中物体在另一个维度上的关系就由梯度来表示，这样我们可以通过观察图像得知物体在三维空间中的对应关系。为什么要提梯度？因为实际上对图像进行二维傅立叶变换得到频谱图，就是图像梯度的分布图，当然频谱图上的各点与图像上各点并不存在一一对应的关系，即使在不移频的情况下也是没有。傅立叶频谱图上我们看到的明暗不一的亮点，实际上图像上某一点与邻域点差异的强弱，即梯度的大小，也即该点的频率的大小（可以这么理解，图像中的低频部分指低梯度的点，高频部分相反）。一般来讲，梯度大则该点的亮度强，否则该点亮度弱。这样通过观察傅立叶变换后的频谱图，也叫功率图，我们首先就可以看出，图像的能量分布，如果频谱图中暗的点数更多，那么实际图像是比较柔和的（因为各点与邻域差异都不大，梯度相对较小），反之，如果频谱图中亮的点数多，那么实际图像一定是尖锐的，边界分明且边界两边像素差异较大的。对频谱移频到原点以后，可以看出图像的频率分布是以原点为圆心，对称分布的。将频谱移频到圆心除了可以清晰地看出图像频率分布以外，还有一个好处，它可以分离出有周期性规律的干扰信号，比如正弦干扰，一副带有正弦干扰，移频到原点的频谱图上可以看出除了中心以外还存在以某一点为中心，对称分布的亮点集合，这个集合就是干扰噪音产生的，这时可以很直观的通过在该位置放置带阻滤波器消除干扰。二、快速傅里叶变换的算法（一）快速傅里叶变换算法快速傅里叶变换算法如下： (1)由(1)式可知，对每一个n，计算X()须作N次复数乘法及N-1次复数加法，要完成这组变换共需次乘法及N(N-1)次复数加法。但以下介绍的快速傅里叶变换的算法，可大大减少运算次数，提高工作效率。当时，n和k可用二进制数表示：又记 ,则（1）式可改写为 （2）式中： （3）因为所以（2）可改成 （4） （5）则式（5）即为式（4）的分解形式。将初始数据代入式(5的第一个等式，可得每一组计算数据，一般将痗L-1组计算数据代入式(5的第L个等式，计算后可得第L组计算数据（L1，2，），计算公式也可表示为= （6）式中 （7） 根据式(5)，第L个数组中每个 的计算只依赖于上一个数组的两个数据这两个数据的标号相差，即，而且这两个数据只用于计算第L个数组中标号的数据（等号右端为二进制数）。当分别取0和1时，分别有。因此，用上一组的两个数据计算所得的两个新数据仍可储存在原来位置，计算过程中只需要N个存储器。将与称为第L个数组中的对偶结点对。计算每个对偶结点对只需一次乘法，事实上由式（6）可得（8）式中： ；别为式（7）中取0，1时对应的P值。因，于是对偶结点的有如下关系：,因此式（6）可表示为（9）P的求法：在中，i写成二进制数右移位，就成为颠倒位序得式(5)吕，前面的个等式，每个等式均对应一组数据进行计算，每组数据都有N/对结点，根据式(9)，每对结点只需作1次乘法和2次加法，因此，每组数据只需N/2次乘法和N次加法，因而完成组数据的计算共需N/2次乘法和N次加法。三、快速傅里叶变换的应用（一） 利用FFT计算连续时间信号的傅里叶变换设是连续时间信号，并假设时，则其傅里叶变换由下式给出令是一个固定的正实数，是一个固定的正整数。当时，利用FFT算法可计算。已知一个固定的时间间隔，选择足够小，使得每一个秒的间隔内，的变化很小，则式中积分可近似为 （10）假设足够大，对于所有的整数，幅值很小，则式（10）变为 （11）当时，式（10）两边的值为 （12）其中代表抽样信号的点。最后令，则上式变为 （13）首先用FFT算法求出，然后可用上式求出时的。应该强调的是，式（11）只是一个近似表示，计算得到的只是一个近似值。通过取更小的抽样间隔，或者增加点数，可以得到更精确的值。如果时，幅度谱很小，对应于奈奎斯特抽样频率，抽样间隔选择比较合适。如果已知信号只在时间区间内存在，可以通过对时的抽样信号补零，使足够大。例1 利用FFT计算傅里叶变换如图1所示的信号其傅里叶变换为：利用下面的命令，可得到的近似值和准确值。 图1 连续时间信号x(t) N=input(Input N:);T=input(Input T:);%计算X(w)近似值t=0:T:2;x=t-1 zeros(1,N-length(t);X=fft(x);gamma=2*pi/(N*T);k=0:10/gamma;Xapp=(1-exp(-i*k*gamma*T)/(i*k*gamma)*X;%计算真实值X(w)w=0.05:0.05:10;Xact=exp(-i*w)*2*i.*(w.*cos(w)-sin(w)./(w.*w);plot(k*gamma,abs(Xapp(1:length(k),o,w,abs(Xact);legend(近似值,真实值);xlabel(频率(rad/s);ylabel(|X|)运行程序后输入N=128，T=0.1，此时，得到实际的和近似的傅里叶变换的幅度谱如图2所示，此时近似值已经相当准确。通过增加NT可以增加更多的细节，减少T使得到的值更精确。再次运行程序后输入N=512，T=0.05，此时，得到实际的和近似的傅里叶变换的幅度谱如图3所示。图2 N=128，T=0.1时的幅度谱图3 N=512，T=0.05时的幅度谱（二）利用FFT计算离散信号的线性卷积已知两个离散时间信号与，取，对和右端补零，使得 （14）利用FFT算法可以求得和的L点DFT，分别是和，利用DTFT卷积性质，卷积等于乘积的L点DFT反变换，这也可以通过FFT算法得到。例2 利用FFT计算线性卷积已知，其中为单位阶跃序列，信号如图4所示。由于当时，很小，故可以取为17；N取10，。利用下面的Matlab命令，可得到、的卷积图形如图4所示。subplot(3,1,1);n=0:16;x=0.8.n;stem(n,x);xlabel(n);ylabel(xn);subplot(3,1,2);n=0:15;y=ones(1,10) zeros(1,6);stem(n,y);xlabel(n);ylabel(yn)subplot(3,1,3);L=26;n=0:L-1;X=fft(x,L);Y=fft(y,L);Z=X.*Y;z=ifft(Z,L);stem(n,z);xlabel(n);ylabel(zn)图4 信号xn、yn及其卷积zn=xn*yn利用下面的Matlab命令，可得到信号xn、yn的幅度谱与相位谱如图5所示。subplot(2,2,1);L=26;k=0:L-1;n=0:16;x=0.8.n;X=fft(x,L);stem(k,abs(X);axis(0 25 0 5);xlabel(k);ylabel(|Xk|)subplot(2,2,2);stem(k,angle(X);axis(0 25 -1 1);xlabel(k);ylabel(Angle(Xk)(弧度)subplot(2,2,3);y=ones(1,10);Y=fft(y,L);stem(k,abs(Y);axis(0 25 0 10);xlabel(k);ylabel(|Yk|)subplot(2,2,4);stem(k,angle(Y);axis(0 25 -3 3);xlabel(k);ylabel(Angle(Yk)(弧度)图5 信号xn、yn的幅度谱与相位谱（3） 利用FFT进行离散信号压缩利用FFT算法对离散信号进行压缩的步骤如下：1）通过采样将信号离散化；2）对离散化信号进行傅里叶变换；3）对变换后的系数进行处理，将绝对值小于某一阈值的系数置为0，保留剩余的系数；4）利用IFFT算法对处理后的信号进行逆傅里叶变换。例3 对单位区间上的下列连续信号以采样频率进行采样，将其离散化为个采样值用FFT分解信号，对信号进行小波压缩，然后重构信号。令绝对值最小的80%系数为0，得到重构信号图形如图6a)所示，均方差为0.0429，相对误差为0.0449；令绝对值最小的90%系数为0，得到重构信号图形如图6 b)所示，均方差为0.0610，相对误差为0.0638。 a) 绝对值最小的80%系数为0的重构信号（FFT） b) 绝对值最小的90%系数为0的重构信号（FFT）图6 用FFT压缩后的重构信号相关Matlab程序如下function wc=compress(w,r)%压缩函数compress.m%输入信号数据w,压缩率r%输出压缩后的信号数据if(r1) error(r 应该介于0和1之间!);end;N=length(w);Nr=floor(N*r);ww=sort(abs(w);tol=abs(ww(Nr+1);wc=(abs(w)=tol).*w;function unbiased_variance,error=fftcomp(t,y,r)%利用FFT做离散信号压缩%输入时间t,原信号y,以及压缩率r%输出原信号和压缩后重构信号的图像,以及重构均方差和相对l2误差if(r1) error(r 应该介于0和1之间!);end;fy=fft(y);fyc=compress(fy,r); %调用压缩函数compress.myc=ifft(fyc);plot(t,y,r,t,yc,b);legend(原信号,重构信号);unbiased_variance=norm(y-yc)/sqrt(length(t);error=norm(y-yc)/norm(y);输入以下Matlab命令：t=(0:255)/256;f=t+cos(4*pi*t)+1/2*sin(8*pi*t);unbiased_variance,error=fftcomp(t,f,0.8)unbiased_variance = 0.0429error =0.0449如果用Harr尺度函数和Harr小波分解信号，对信号进行小波压缩，然后重构信号。令绝对值最小的80%系数为0，得到重构信号图形如图7 a)所示，均方差为0.0584，相对误差为0.0611；令绝对值最小的90%系数为0，得到重构信号图形如图7 b)所示，均方差为0.1136，相对误差为0.1190。 a) 绝对值最小的80%系数为0的重构信号（Harr） b) 绝对值最小的90%系数为0的重构信号（Harr）图7 用Harr小波压缩后的重构信号相关Matlab程序如下function unbiased_variance,error=daubcomp(t,y,n,r)%利用Daubechies系列小波做离散信号压缩%输入时间t,原信号y,分解层数n,以及压缩率r%输出原信号和压缩后重构信号的图像,以及重构均方差和相对l2误差if(r1) error(r应该介于0和1之间!);end;c,l=wavedec(y,n,db1);cc=compress(c,r); %调用压缩函数compress.myc=waverec(cc,l,db1);plot(t,y,r,t,yc,b);legend(原信号,重构信号);unbiased_variance=norm(y-yc)/sqrt(length(t);error=norm(y-yc)/norm(y);输入以下Matlab命令：t=(0:255)/256;f=t+cos(4*pi*t)+1/2*sin(8*pi*t);unbiased_variance,error=daubcomp(t,f,8,0.8)unbiased_variance = 0.0584error = 0.0611结论：在信号没有突变、快变化或者大致上具有周期性的信号，用FFT可以处理得很好（甚至比小波还要好）。（四）利用FFT对离散信号进行滤波利用FFT算法对信号进行滤波的步骤如下：1）通过采样将信号离散化；2）对离散化信号进行傅里叶变换；3）对变换后的系数进行处理，将绝对值大于某一阈值的系数置为0，保留剩余的系数；4）利用IFFT算法对处理后的信号进行逆傅里叶变换。例4 股票价格分析首先进入网址/q/hp?s=YHOO，点击网页底部位置的Download To Spreadsheet按钮，即可把以Excel表格格式存储的价格数据下载到本地计算机。表格从1列至第6列分别给出了从1996年4月12日至2007年5月30日的交易期里每天的开盘价、最高价、最低价、收盘价、成交量以及趋势。数据下载完成后，需要颠倒顺序，使得最早时间的数据首先显示。然后另存到Matlab所在的目录中，并重新命名为“yhoodata.csv”。 本次分析选择开盘价，时间是从2007年1月1日至2007年5月30日的=102个交易日期。令代表一支股票的开盘价。为了便于分析，可以先从中减去跃变，得到，即 （15）输入以下命令，得到的频谱如图8所示。o=csvread(yhoodata.csv,2700,1,2700 1 2801 1);N=102;for n=1:N x(n)=o(n)-o(1)-(o(N)-o(1)/(N-1)*(n-1);endX=fft(x);k=0:N-1;stem(k,abs(X);axis(0 101 0 300);xlabel(k);ylabel(|Xk|)图8 xn的幅度谱可以看出上图中有5个较强谱分量，频率（）分别对应和。保留这5个频率分量的系数，将其他频率分量的系数置为0，然后再进行逆傅里叶变换，得到滤波后的近似值。输入如下Matlab程序，得到真实值与滤波后的近似值，如图9所示。plot(x);hold on;fliterX=X(1:2) 0 0 X(5) zeros(1,102-9) X(99) 0 0 X(102);fliterx=ifft(fliterX);plot(real(fliterx),r);axis(1 102 0 7);xlabel(n);ylabel(xn的值和它的近似值);legend(xn真实值,xn近似值)图9 xn的真实值与滤波后的近似值从上图可以看出，滤波后的近似值既大致上保留了真实值的变化趋势，而且与其十分接近。与滤波前比较，滤波后的图形要比滤波前平滑得多。再由式（15）即可求得 （16）输入如下Matlab程序，画出真实开盘价与近似开盘价的图形。如图21所示，可以看出是近似基础上的平滑值。plot(o);hold on;for n=1:N oapp(n)=fliterx(n)+o(1)+(o(N)-o(1)/(N-1)*(n-1);endplot(oapp,r);axis(1 102 25 34);xlabel(n);ylabel(on的值和它的近似值);legend(on真实值,on近似值)图10 on的真实值与滤波后的近似值（五）利用FFT提取离散信号中的最强正弦分量这里最强是指在信号中某个正弦分量的振幅远大于其它正弦分量的振幅。可以对求点来确定信号中是否有最强的正弦分量。如果信号是连续时间形式的，首先还要对其进行抽样，得到离散时间形式的信号，根据Nyquist定理，抽样间隔应满，其中是中的最大频率分量。要判断信号中是否包含最强正弦分量，采样数据至少要包含该分量一个完整周期的数据，例5 太阳耀斑数据的分析太阳耀斑活动的周期是11年，这个事实可以通过提取耀斑数据的最强正弦分量加以证实。耀斑数据可以从比利时皇家天文台(Royal Observatory of Belgium)太阳耀斑数据索引中心（Sunspot Index Data Center, SIDC）网站下载。网址是：http:/sidc.oma.be/index.php3下载后的数据存放在文件“monthssn.dat”中，里面有四列数据，第一年是日期，第三列是太阳耀斑的平均数，第四列平滑后太阳耀斑的平均数，可以得到从1749年到当前年月（2006年4月）的耀斑数据。本次分析选择第三列1981年1月作为开始日期，2005年12月作为结束日期，共25年300个月份的数据。为此先把相关数据复制到Excel表格的第一列中，然后保存到Matlab所在目录下，并命名为“sunspotdata.csv”。然后输入以下命令，得到耀斑曲线如图11所示。spd=csvread(sunspotdata.csv,0,0,0 0 299 0);plot(spd);grid;xlabel(月数);ylabel(耀斑平均数);axis(0 300 0 200)图11 1981年1月至2005年12月太阳耀斑的平均数由上图可见，太阳耀斑的活动确实具有周期性，但周期的准确值不明显。可以通过数第一个峰值和第二个峰值之间的月份来估计周期的值。查验表中的数据得，第一个峰值为200.3，出现在第116个月（1990年8月），第二个峰值为170.1，出现在第235个月（2000年7月），所以周期是235-116=119月，和实际值132月比较接近。下面利用来分析。首先，从图11中可以看出从到整个区间的平均值不为0，为了更方便地分析，需要减去该均值，得到结果如下： （17）然后对进行傅里叶变换，得到。在Matlab中输入以下命令，得到的幅度谱的图形如图12所示。x=spd-sum(spd)/300;X=fft(x);k=0:299;plot(k,abs(X);grid;xlabel(k);ylabel(|Xk|)图12 的幅度谱由图12可见，的幅度谱中有一个清晰的尖峰，这就证实中确实包含一个最强正弦分量。为了确定尖峰所对应的频率，用火柴棒图重画当时的图形，在Matlab中输入以下命令，得到图形如图13所示。k=0:10;X=X(1:11);stem(k,abs(X);xlabel(k);ylabel(|Xk|)图13 当时的火柴棒图可以看出上图中有两个强谱分量，频率分别为弧度/月和，周期分别为150月和100月。由于的数据长度不是其中强正弦分量的整数倍，谱分量在出现了泄露。要消除泄露，需要使数据的长度正好是其中强正弦分量的周期（太阳耀斑活动的周期，也即）的整数倍。为此，重新取中从1至264之间的数据进行分析，令 （18）然后对进行傅里叶变换，得到。输入以下