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

1海平面高度变化的海平面高度变化的 Fourier 分析（分析（Matlab6.x）以海平面高度变化的时间序列为例进行 Fourier 分析，供有兴趣的同学参考。第一步，存入数据 在 Matlab 的编辑窗口中将时间序列的序号定义为 x 向量，将时间序列定义为 y 向量，然后存为名为 TS1 的 m 文件，以便调用（图 1）。数据可从 Excel 中拷贝、粘贴。注意：数据的保存与执行程序的保存不同，后者可以直接在编辑窗口保存，在命令窗口运行；前者则必须在命令窗口保存。因此，需要将图 1 中的数据文件复制到命令窗口(Command Window)中，回车即可保存（图 2）。图 1 时间序列的向量表示 图 1 时间序列的向量表示 图 2 在命令窗口保存数据文件 图 2 在命令窗口保存数据文件 2第二步，编写 Fourier 变换的 m 文件 根据快速 Fourier 变换的原理，借助 Matlab 的有关命令，不难编写进行 FFT 的计算程序。需要说明的是，由于调用数据，首先需要清屏（clear）并装入（load）数据文件。M 文件中各个语句的意义已在旁边注释，供参考（图 3）。第三步，执行计算命令 将编写成功的 FFT 文件从编辑窗口复制到命令窗口，回车，立即得到计算结果（图 4，图 5）。计算结果中包括最大谱密度 Pm 和时间序列的周期 T（图 4）。图 3 Fourier 变换的 m 文件之一 图 3 Fourier 变换的 m 文件之一 在默认的情况下，Matlab 总是根据最后一个绘图指令给出图像。在我们的计算程序中，最后一个绘图指令是局部放大的频谱图，执行计算命令以后，Matlab 立即弹出该图（图 5），图上加有格网（grid），目的是便于直观估计最大频谱及其对应的频率。在 Matlab 的图形窗口上，可以对图形进行编辑处理，当然也可以放大和缩小。要想通过曲线直观地估计周期长度，可以利用图形工具栏（Figure Toolbar）中的（Zoom in）对曲线进行局部放大，借助放大的结果估计最大功率谱密度 Pm 对应的频率 fp，该频率的倒数就是时间序列的周期,即有 T=1/fp。不过，本程序已设计为直接计算出周期长度。3 图 4 Fourier 文件的执行结果 图 4 Fourier 文件的执行结果 如果需要绘出其他图形，则可在命令窗口中单独执行有关命令。例如，将以下指令“plot(f(1:N/2+1),P(1:N/2+1)%以对称点为界绘制频谱图”拷贝到命令窗口的提示符下：回车，立即得到的以对称点为界的频谱图（图 6）。执行命令 可得对称分布的频谱全图（图 7）从该图可以看出我们进行 FFT 为什么只取一半的结果。执行命令 可得原始时间序列变化的曲线图式（图 8）。400.020.040.060.080.10.120.140.160.180.2020004000600080001000012000140001600018000Freq-Spec Fig 图 5 局部放大并加上格网的频谱图图 5 局部放大并加上格网的频谱图 00.050.10.150.20.250.30.350.40.450.5020004000600080001000012000140001600018000 图 6 以对称点为界的频谱图 图 6 以对称点为界的频谱图 500.10.20.30.40.50.60.70.80.91020004000600080001000012000140001600018000 图 7 对称分布的频谱全图 图 7 对称分布的频谱全图 0102030405060708090100-60-40-20020406080100120 图 8 时间序列的原始信号图 图 8 时间序列的原始信号图 6【附录】我们曾经讲过，计算时间序列的周期至少有三个方面的意义：分析）以进行实践意义的预测发展预测（有了周期可）响中国大陆的洪灾频率果关系（如太阳黑子影揭示地理系统演化的因含规律）律（认识地理现象的隐寻找地理系统的内在节 功率谱分析是寻找时间序列周期最为有效的方法之一，Matlab 提供了 Fourier 分析的全套工具。今后各位可以在不断地应用中逐步掌握。最后，作为参考，给出 FFT 的另一个 m 文件（图 9），下面这种文件的与前面那个文件大同小异。文件中使用了线性等分向量的命令“linspace”，这个命令为绝大多少进行 FFT的编程者使用，但是，使用这个命令进行 Fourier 变换的结果不是十分精确的，原因在于：FFT 虽然是对 2n个数据进行，但结果却有 2n+1 个数据点。但是，应用线性等分向量的命令“linspace”以后，系统自觉的将变换结果视为 2n个数据点，而将最后一个数据舍弃。因此，最后给出的最大功率谱密度虽然正确，但对应的频率却有误差，从而计算的周期也不精确。例如，在命令窗口执行图 9 中的 m 文件，得到的周期为 T=10.5833（图 10）。不过，当时间序列足够长时，误差将会很小。还有一点需要说明的是：我们之所以将数据和 FFT 程序存为不同的文件，主要是想为同一个程序计算多个不同的时间序列提供方便。将时间序列与计算过程的有关指令放在一起也没有什么不妥下面的文件即是如此（图 9）。图 9 Fourier 变换的 m 文件之二 图 9 Fourier 变换的 m 文件之二 7 图 10 第二个 FFT 程序的执行结果图 10 第二个 FFT 程序的执行结果