利用Mathcad进行时间序列的谱分析（I）.doc

利用Mathcad进行时间序列的谱分析（I） 在这一节中我们采用与“利用Excel进行时间序列的谱分析（I）”一节相同的例子，以便计算结果能够相互参照。关于Mathcad的基本操作技巧（包括绘图）已在“回归分析”中备述，本节不再详细说明。 【例1】某水文观测站测得一条河流从1979年6月到1980年5月共计12月份的断面平均流量。试借助功率谱分析判断时间序列是否存在周期性特征，并对周期长度进行估计。 第一步，录入或拷入数据 可以将在Excel中经过中心化的数据直接录入或拷入Mathcad的数据表中，并定义为“data”。 图1 录入的中心化数据 需要说明两点：第一，如果是从Excel中拷入数据，一定记住检查最后一行数据，如果多出一行0数据点，必须将其删除。 第二，也可以在Mathcad中对数据中心化。中心化的方法非常简单，将原始数据拷入Mathcad的数据表并将其定义为“data”（图2）。 图2 原始数据表的局部显示（数据未经中心化） 在被定义过的表下定义变量为 然后利用求均值命令mean，键入公式 回车，立即得到中心化的结果（图3）。 图3 中心化的数据 不过，中心化以后需要重新插入一个数据表并重新定义该表（例如定义为biao）。将中心化的数据拷入“biao”中，然后在左方插入单元格（Insert Cells），再在新插入的一列单元格中录入数据序号（图4）。从这个意义上讲，在Mathcad中进行数据处理不如直接将经过Excel预处理即中心化的数据拷贝过来。 图4 新建的数据表 第二步，绘制数据序列的变化图 假定录入的是图1所示的中心化的数据，将表定义为data以后，在数据表下对变量进行如下定义： 然后利用Graph中的曲线图工具容易绘制曲线图（图5）。 图5 中心化时间序列的变化曲线（1979年6月－1980年5月） 绘图的主要目的有两个：其一，可以通过曲线图直观地判断数据变化是否存在周期，对于本例，当然达不到这种效果；其二，如果拷贝的数据多出0数据行，可以及时将其删除，以免计算出错。 第三步，进行快速Fourier变换 首先是定义变量。在Mathcad中，系统默认的数据排列方式是从0到K。我们有12个月份的数据，也就相当于拥有12个样本，即N=12，从而K=12-1=11。但是，Fourier变换要求时间序列的长度必须是T=2n个（n为正整数），因此，可以考虑在序列末尾加0，将其延长到T=24个。方法是定义序号k:=12..15，并设xk=0，这相当于在Excel加上4个0元素。这种补充定义已在上文给出。 定义快速Fourier变换，命令为fft（注意：在Mathcad中，大写的FFT与小写的fft变换结果相差T倍，为了与Excel的结果对应，本例采用fft）。我们是要对x进行fft，令变换的结果为Y，Y将给出对x进行快速Fourier变换的结果；定义M为最后一个变换结果的序号，则可以通过last自动从对称点(参见“利用Excel进行时间序列的谱分析（I）”一节的图15、图16)截取前半部结果，M=N/2=8；定义功率谱为Poweri，则有 注意：如果使用命令FFT，则上式的分母应为T2。定义频率变化的步长为w=1/2；则频率为freqj=j/(2M)=j/T，并取j=0,1,2,…,M。全部定义的表达式如下： 定义完成以后，键入“Yi=”，回车，得到x的FFT结果；键入“Poweri=”或“︱Yi︱2＝”，回车，得到功率谱密度（图6）。 图6 Fourier变换结果与功率谱密度 第四步，绘制频谱图，识别周期 利用Mathcad的Graph工具箱容易绘制频谱图（图7）。在图上点右键，出现一个选择菜单（图8）；选中Trace，弹出X-Y Trace对话框（图9）；用鼠标点击谱线的最大值点，则Trace自动给出谱密度的极大值P(f)max=63157，对应的频率为f1=0.0625；点击Copy X按钮，将0.0625复制到Mathcad的工作表上，取倒数，立即到周期P=1/f1=16。当然，由于时间序列太短，这个周期是不准确的。 图7 例1的频谱图 图8 右键菜单与Trace的位置 图9 借助Trace显示功率谱密度最大值对应的频率 【例2】对同一条河流的连续两年的平均月径流量（1961年6月－1963年5月）进行Fourier分析，并识别其周期。 步骤与例1相同，不详述。在录入中心化数据以后（图10），需要重新定义数据的范围。 图10 例2的数据表（局部，已中心化） 考虑到这次N=24，可以将其延长到T=25=32。于是Fourier变换定义如下： 画出时间序列的变化曲线图，从图上可以直观地判断径流量变化具有某种周期性，但不能确定（图11）。 图11 例2的时间序列变化图（1