海平面高度变化的Fourier分析(Matlab).pdf

《海平面高度变化的Fourier分析(Matlab).pdf》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。

1、海平面高度变化的Fourier分析（Matlab6.x）以海平面高度变化的时间序列为例进行Fourier分析，供有兴趣的同学参考。第一步，存入数据在Matlab的编辑窗口中将时间序列的序号定义为x向量，将时间序列定义为y向量，然后存为名为TS1的m文件，以便调用（图1）。数据可从Excel中拷贝、粘贴。注意：数据的保存与执行程序的保存不同，后者可以直接在编辑窗口保存，在命令窗口运行；前者则必须在命令窗口保存。因此，需要将图1中的数据文件复制到命令窗口(CommandWindow)中，回车即可保存（图2）。图1时间序列的向量表示图2在命令窗口保存数据文件1第二步，编写Four

2、ier变换的m文件根据快速Fourier变换的原理，借助Matlab的有关命令，不难编写进行FFT的计算程序。需要说明的是，由于调用数据，首先需要清屏（clear）并装入（load）数据文件。M文件中各个语句的意义已在旁边注释，供参考（图3）。第三步，执行计算命令将编写成功的FFT文件从编辑窗口复制到命令窗口，回车，立即得到计算结果（图4，图5）。计算结果中包括最大谱密度Pm和时间序列的周期T（图4）。图3Fourier变换的m文件之一在默认的情况下，Matlab总是根据最后一个绘图指令给出图像。在我们的计算程序中，最后一个绘图指令是局部放大的频谱图，执行计算命令以后，Ma

3、tlab立即弹出该图（图5），图上加有格网（grid），目的是便于直观估计最大频谱及其对应的频率。在Matlab的图形窗口上，可以对图形进行编辑处理，当然也可以放大和缩小。要想通过曲线直观地估计周期长度，可以利用图形工具栏（FigureToolbar）中的（Zoomin）对曲线进行局部放大，借助放大的结果估计最大功率谱密度Pm对应的频率fp，该频率的倒数就是时间序列的周期,即有T=1/fp。不过，本程序已设计为直接计算出周期长度。2图4Fourier文件的执行结果如果需要绘出其他图形，则可在命令窗口中单独执行有关命令。例如，将以下指令“plot(f(1:N/2+1),P(1

4、:N/2+1))%以对称点为界绘制频谱图”拷贝到命令窗口的提示符下：回车，立即得到的以对称点为界的频谱图（图6）。执行命令可得对称分布的频谱全图（图7）——从该图可以看出我们进行FFT为什么只取一半的结果。执行命令可得原始时间序列变化的曲线图式（图8）。318000160001400012000100008000600040002000000.020.040.060.080.10.120.140.160.180.2Freq-SpecFig图5局部放大并加上格网的频谱图18000160001400012000100008000600040002000000.050.10.15

5、0.20.250.30.350.40.450.5图6以对称点为界的频谱图418000160001400012000100008000600040002000000.10.20.30.40.50.60.70.80.91图7对称分布的频谱全图120100806040200-20-40-600102030405060708090100图8时间序列的原始信号图5【附录】我们曾经讲过，计算时间序列的周期至少有三个方面的意义：⎧寻找地理系统的内在节律（认识地理现象的隐含规律）⎪⎨揭示地理系统演化的因果关系（如太阳黑子影响中国大陆的洪灾频率）⎪⎩发展预测（有了周期可以进行实践意义的预测分

6、析）功率谱分析是寻找时间序列周期最为有效的方法之一，Matlab提供了Fourier分析的全套工具。今后各位可以在不断地应用中逐步掌握。最后，作为参考，给出FFT的另一个m文件（图9），下面这种文件的与前面那个文件大同小异。文件中使用了线性等分向量的命令“linspace”，这个命令为绝大多少进行FFT的编程者使用，但是，使用这个命令进行Fourier变换的结果不是十分精确的，原因在于：nnFFT虽然是对2个数据进行，但结果却有2+1个数据点。但是，应用线性等分向量的命令n“linspace”以后，系统自觉的将变换结果视为2个数据点，而将最后一个数据舍弃。因此，最后给出的最

7、大功率谱密度虽然正确，但对应的频率却有误差，从而计算的周期也不精确。例如，在命令窗口执行图9中的m文件，得到的周期为T=10.5833（图10）。不过，当时间序列足够长时，误差将会很小。还有一点需要说明的是：我们之所以将数据和FFT程序存为不同的文件，主要是想为同一个程序计算多个不同的时间序列提供方便。将时间序列与计算过程的有关指令放在一起也没有什么不妥——下面的文件即是如此（图9）。图9Fourier变换的m文件之二6图10第二个FFT程序的执行结果7