利用Mathcad进行时间序列的谱分析Read

1、利用Mathcad进行时间序列的谱分析（I）在这一节中我们采用与“利用Excel进行时间序列的谱分析（I）” 一节相同的例子，以便计算结果能够相互参照。关于Mathcad的基本操作技巧（包括绘图）已在“回归分析”中备述，本节不再详细说明。【例1】某水文观测站测得一条河流从1979年6月到1980年5月共计12月份的断面平均流量。试借助功率谱分析判断时间序列是否存在周期性特征，并对周期长度进行估计。第一步，录入或拷入数据可以将在Excel中经过中心化的数据直接录入或拷入Mathcad的数据表中，并定义为“ data”。data010121.141279.1423384.1434117.14453

2、5.1456-62.8667-103.4678-121.9689-135.06910-117.061011-90.461112-5.86图1录入的中心化数据需要说明两点：第一，如果是从Excel中拷入数据，一定记住检查最后一行数据，如果多出一行0数据点，必须将其删除。第二，也可以在 Mathcad中对数据中心化。中心化的方法非常简单，将原始数据拷入 Mathcad的数据表并将其定义为“ data”（图2 ”。data010119012248图2原始数据表的局部显示（数据未经中心化）在被定义过的表下定义变量为t data0 x data 1 然后利用求均值命令 mea n,键入公式x mea n

3、（x）回车，立即得到中心化的结果（图3）。01021.142179.1422384.1423117.142435.1425-62.8586-103.4587-121.9588-135.0589-117.05810-90.45811-5.858x mean(x)图3中心化的数据不过，中心化以后需要重新插入一个数据表并重新定义该表（例如定义为biao）。将中心化的数据拷入"biao”中，然后在左方插入单元格（Insert Cells），再在新插入的一列单元格中 录入数据序号（图4）。从这个意义上讲，在Mathcad中进行数据处理不如直接将经过Excel预处理即中心化的数据拷贝过来。bia

4、o010121.141279.14图4新建的数据表第二步，绘制数据序列的变化图假定录入的是 图1所示的中心化的数据，将表定义为data以后，在数据表下对变量进行如下定义：xk0t data0 x data 1k 12 15然后利用Graph中的曲线图工具容易绘制曲线图（ 图5）。时序图5中心化时间序列的变化曲线（1979年6月一1980年5月）化儿中量济径绘图的主要目的有两个：其一，可以通过曲线图直观地判断数据变化是否存在周期，对于本例，当然达不到这种效果；其二，如果拷贝的数据多出0数据行，可以及时将其删除，以免计算出错。第三步，进行快速 Fourier变换首先是定义变量。在 Mathcad中

5、，系统默认的数据排列方式是从0到K。我们有12个月份的数据，也就相当于拥有12个样本，即N=12，从而K=12-仁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将给出对x进行快速Fou

6、rier变换的结果；定义 M为最后一个变换结果的序号，则可以通过last自动从对称点（参见“利用Excel进行时间序列的谱分析（I）” 一节的图15、图16）截取前半部结果，M = N/2=8;定义功率谱为 Poweri，则有2门YPower T注意：如果使用命令 FFT，则上式的分母应为 T2。定义频率变化的步长为 w=1/2 ;则频率为 freqj=j/（2M）=j/T，并取j=0,1,2,M。全部定义的表达式如下：Y fft(x)M last(Y)M 8i 0 M211Power.Y.w -j 0 Mfreq. w j jMii2定义完成以后，键入“ 丫=”，回车，得到x的FFT结果；键

7、入“ Poweri= ”或“丨Yi | "= 回车，得到功率谱密度（图6 ）。Y. iPower.iYi 2 i5.425?0-8i2.943?0-15I2.943?0-15232.503+95.388i6.316?0 46.316?0 4-74.087+144.926i2.649?0 42.649?0 4-75.626+68.218i1.037?0 41.037?0 4-67.25-22.525i5.03?0 35.03?0 3-50.648-45.43i4.629?0 34.629?0 3-0.442-53.645i2.878?0 32.878?0 349.971-53.401i

8、5.349?0 35.349?0 355.7253.105?0 33.105?0 3图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/fi=i6。当然，由于时间序列太短，这个周期是不准确的。6.316 106

9、Powerj42152.943 10h7图7例1的频谱图Powerj4 104囹 Prjoperlies.Fermata.PP43X1000 1ZsQiriii <Diab l& Evaluatio nTroe.,电 CopyPaste6 1042 if? if0.5图8右键菜单与Trace的位置rower-10"褂图9借助Trace显示功率谱密度最大值对应的频率【例2】对同一条河流的连续两年的平均月径流量（1961年6月一1963年5月）进行Fourier分析，并识别其周期。步骤与例1相同，不详述。在录入中心化数据以后（图10），需要重新定义数据的范围。data010

10、139.5512168.55图10例2的数据表（局部，已中心化）考虑到这次N=24,可以将其延长到T=25=32。于是Fourier变换定义如下:t data0data 12431xk但不能确画出时间序列的变化曲线图,Y fft(x)M last(Y)M 16i 0 M211Power.Y.w -j 0 Mfreq. w j 一JMii2从图上可以直观地判断径流量变化具有某种周期性,定（图11）。2001000213.546 3584.354ioo 01is m101520t时序2524图11例2的时间序列变化图（1961年6月一1963年5月）利用Fourier变换结果（图12）画出频谱图（

11、图13）。借助Trace功能容易找到最大功 率谱密度P（f）max=63157及其对应的频率 20.09375 （图14）；点击Copy X按钮，将0.09375 复制到Mathcad的工作表上，取倒数，可以得到周期P=1/f1=10.667。在最大值附近存在一个与最大值非常接近的功率谱密度值为3110 7,对应的频率为0.0625。由于这个次最大点同样代表一个突变点一一它与临近值有很大差距，故须考虑这个 数值。根据上例可知，其倒数为16。由此判断，时间序列的周期大约在1016之间。我们知道，河流径流量的周期为12个月。可见，对于较短的时间序列，功率谱分析会有较大的误差。Yi-3.005?0

12、刑31.059+32.314i154.717-85.201i-23.963+177.338i22.68-20.025i17.391+59.409i-32.591+12.514i21.468+16.304i-2.775+32.173i-3.39-6.952i11.502+45.611i-17.362+2.902i1.68+37.675i-32.942-6.489i9.594+28.757i-51.658-3.068i图12例Power.iYi 2 i0102.009?0 32.009?0 33.12?0 43.12?0 43.202?0 43.202?0 4915.377915.3773.832

13、?0 33.832?0 31.219?0 31.219?0 3726.679726.6791.043?0 31.043?0 359.81659.8162.213?0 32.213?0 3309.87309.871.422?0 31.422?0 31.127?0 31.127?0 3918.989918.9892.678?0 32.678?0 32的Fourier 变换结果.3.202 10 _3Power j210_图13例2的频谱图K-Y Trace2S1y_ValueSy 'L |丽 Trtok ToiktX-|o. 0625T-Valne |311STTrack Foint图14

14、例2的功率谱密度最大值及其对应的频率【例3】某海域海面平均高度年际变化的功率谱分析。本例时间序列长度为 100年，即N=100。将时间序列延长到 T=27=128。data0101-41.0312-35.03图15例3数据的局部显示（已经中心化）t data 0ydata 1k 100 127yk 0丫 fft (y)Mlast(Y)M 64i 0 M211Power.Y.w _j 0 Mfreq j wii2jM从原始数据的图像可以看到，时间序列可能具有某种内在节律（图 16）。录入数据以后，根据时间序列长度对变量和Fourier变换定义如下:度高均平面海108.27446.026图16海面

15、高度的年际变化曲线借助Fourier变换的结果，作出频谱图如下（ 图17）。4 21.683 101041.51044Powerj 1 1050000 000.10.20.30.40.50freqj0.5图17海面高度变化的频谱图f=0.09375，于是得到借助Mathcad的Trace功能发现功率谱密度最大值对应的频率为 周期约为P=1/f=10.667 （年）。1.5 10Poiverj 1 10"Track Data Point0.1|0.093T5Y-7alue |16&3302丄百33迅心0.40.5图18最大功率谱密度值对应的频率【例4】 在例2中，利用时间序列的

16、自相关函数分析其周期特性。在“利用Excel进行时间序列的谱分析(I)” 一节，我们谈到可以定义时间序列的周期图为N 22I (fi)2 bi )， i 1,2, k式中I(fi)为频率fi处的强度。以fi为横轴，以l(fi)为纵轴，绘制时间序列的周期图，可以在 最大值处找到时间序列的周期。实际上，周期图还可以表作N 1l(fi) R?(k)ej2fkk 1 N这里R为时间序列的自协方差函数，即有R(k)1 N ikNt1XtXfk'这意味着，周期图是对功率谱直接估计的结果。 根据上面几个式子及其内在关系， 我们可以 利用相关函数估计时间序列的周期。 我们知道，计算相关系数的公式有两种

17、，本例采用下面 公式n krk(Xt X)(Xt k X)t 1n(Xt X)2t 1采用这个公式的原因之一是： 我们的时间序列较短， 而这个公式的有效性较好； 原因之二是 可以直接借助SPSS给出结果。将结果复制到 Excel中间，转换格式以后，再复制到Mathcad 的数据表(图19)，然后进行Fourier变换。01010.66120.42230.0334-0.2945-0.5156-0.5767-0.5378-0.3989-0.159100.0910110.311120.412130.4513140.2214150.131516-0.02图19径流量的自相关系数data进行快速Fourier运算的有关定义如下:t data0xdata1Y fft(x)Mlast(Y)Power.iYi21w -2根据计算结果容易得到频谱图（图20）。利用k 16 31xk0M 16i 0Mj 0 M鬥.1 w jMTrace可知，最大功率谱密度对应的