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

海平面高度变化的 Fourier 分析（Matlab6.x） 以海平面高度变化的时间序列为例进行Fourier 分析，供有兴趣的同学参考。 第一步，存入数据 在 Matlab 的编辑窗口中将时间序列的序号定义为 x 向量，将时间序列定义为y 向量， 然后存为名为 TS1 的m 文件，以便调用（图 1）。数据可从Excel 中拷贝、粘贴。 注意：数据的保存与执行程序的保存不同，后者可以直接在编辑窗口保存，在命令窗口 运行；前者则必须在命令窗口保存。因此，需要将图 1 中的数据文件复制到命令窗口 (Command Window)中，回车即可保存（图2 ）。 图 1 时间序列的向量表示 图 2 在命令窗口保存数据文件 1 第二步，编写 Fourier 变换的 m 文件 根据快速 Fourier 变换的原理，借助 Matlab 的有关命令，不难编写进行 FFT 的计算程序。 需要说明的是，由于调用数据，首先需要清屏（clear ）并装入（load ）数据文件。M 文件中 各个语句的意义已在旁边注释，供参考（图3 ）。 第三步，执行计算命令 将编写成功的 FFT 文件从编辑窗口复制到命令窗口，回车，立即得到计算结果（图4 ， 图5 ）。计算结果中包括最大谱密度Pm 和时间序列的周期 T （图4 ）。 图 3 Fourier 变换的 m 文件之一 在默认的情况下，Matlab 总是根据最后一个绘图指令给出图像。在我们的计算程序中， 最后一个绘图指令是局部放大的频谱图，执行计算命令以后，Matlab 立即弹出该图（图5 ）， 图上加有格网（grid ），目的是便于直观估计最大频谱及其对应的频率。 在 Matlab 的图形窗口上，可以对图形进行编辑处理，当然也可以放大和缩小。要想通 过曲线直观地估计周期长度，可以利用图形工具栏（Figure Toolbar） 中的 （Zoom in ）对曲线进行局部放大，借助放大的结果估计最大功率谱密度 Pm 对应的 频率 fp，该频率的倒数就是时间序列的周期,即有 T=1/fp。不过，本程序已设计为直接计算 出周期长度。 2 图 4 Fourier 文件的执行结果 如果需要绘出其他图形，则可在命令窗口中单独执行有关命令。例如，将以下指令 “plot(f(1:N/2+1),P(1:N/2+1)) % 以对称点为界绘制频谱图” 拷贝到命令窗口的提示符下： 回车，立即得到的以对称点为界的频谱图（图6 ）。 执行命令 可得对称分布的频谱全图（图7 ）——从该图可以看出我们进行 FFT 为什么只取一半的结果。 执行命令 可得原始时间序列变化的曲线图式（图 8 ）。 3 18000 16000 14000 12000 10000 8000 6000 4000 2000 0 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 Freq-Spec Fig 图 5 局部放大并加上格网的频谱图 18000 16000 14000 12000 10000 8000