利用Mathcad进行时间序列的谱分析I.docx
- 文档编号:29666561
- 上传时间:2023-07-26
- 格式:DOCX
- 页数:14
- 大小:404.13KB
利用Mathcad进行时间序列的谱分析I.docx
《利用Mathcad进行时间序列的谱分析I.docx》由会员分享,可在线阅读,更多相关《利用Mathcad进行时间序列的谱分析I.docx(14页珍藏版)》请在冰豆网上搜索。
利用Mathcad进行时间序列的谱分析I
利用Mathcad进行时间序列的谱分析(
)
在这一节中我们采用与“利用Excel进行时间序列的谱分析(
)”一节相同的例子,以便计算结果能够相互参照。
关于Mathcad的基本操作技巧(包括绘图)已在“回归分析”中备述,本节不再详细说明。
【例1】某水文观测站测得一条河流从1979年6月到1980年5月共计12月份的断面平均流量。
试借助功率谱分析判断时间序列是否存在周期性特征,并对周期长度进行估计。
第一步,录入或拷入数据
可以将在Excel中经过中心化的数据直接录入或拷入Mathcad的数据表中,并定义为“data”。
图1录入的中心化数据
需要说明两点:
第一,如果是从Excel中拷入数据,一定记住检查最后一行数据,如果多出一行0数据点,必须将其删除。
第二,也可以在Mathcad中对数据中心化。
中心化的方法非常简单,将原始数据拷入Mathcad的数据表并将其定义为“data”(图2)。
图2原始数据表的局部显示(数据未经中心化)
在被定义过的表下定义变量为
然后利用求均值命令mean,键入公式
回车,立即得到中心化的结果(图3)。
图3中心化的数据
不过,中心化以后需要重新插入一个数据表并重新定义该表(例如定义为biao)。
将中心化的数据拷入“biao”中,然后在左方插入单元格(InsertCells),再在新插入的一列单元格中录入数据序号(图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进行时间序列的谱分析(
)”一节的图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)。
图6Fourier变换结果与功率谱密度
第四步,绘制频谱图,识别周期
利用Mathcad的Graph工具箱容易绘制频谱图(图7)。
在图上点右键,出现一个选择菜单(图8);选中Trace,弹出X-YTrace对话框(图9);用鼠标点击谱线的最大值点,则Trace自动给出谱密度的极大值P(f)max=63157,对应的频率为f1=0.0625;点击CopyX按钮,将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的时间序列变化图(1961年6月-1963年5月)
利用Fourier变换结果(图12)画出频谱图(图13)。
借助Trace功能容易找到最大功率谱密度P(f)max=63157及其对应的频率f1=0.09375(图14);点击CopyX按钮,将0.09375复制到Mathcad的工作表上,取倒数,可以得到周期P=1/f1=10.667。
在最大值附近存在一个与最大值非常接近的功率谱密度值为31107,对应的频率为0.0625。
由于这个次最大点同样代表一个突变点——它与临近值有很大差距,故须考虑这个数值。
根据上例可知,其倒数为16。
由此判断,时间序列的周期大约在10~16之间。
我们知道,河流径流量的周期为12个月。
可见,对于较短的时间序列,功率谱分析会有较大的误差。
图12例2的Fourier变换结果
图13例2的频谱图
图14例2的功率谱密度最大值及其对应的频率
【例3】某海域海面平均高度年际变化的功率谱分析。
本例时间序列长度为100年,即N=100。
将时间序列延长到T=27=128。
图15例3数据的局部显示(已经中心化)
录入数据以后,根据时间序列长度对变量和Fourier变换定义如下:
从原始数据的图像可以看到,时间序列可能具有某种内在节律(图16)。
图16海面高度的年际变化曲线
借助Fourier变换的结果,作出频谱图如下(图17)。
图17海面高度变化的频谱图
借助Mathcad的Trace功能发现功率谱密度最大值对应的频率为f=0.09375,于是得到周期约为P=1/f=10.667(年)。
图18最大功率谱密度值对应的频率
【例4】在例2中,利用时间序列的自相关函数分析其周期特性。
在“利用Excel进行时间序列的谱分析(
)”一节,我们谈到可以定义时间序列的周期图为
,
式中I(fi)为频率fi处的强度。
以fi为横轴,以I(fi)为纵轴,绘制时间序列的周期图,可以在最大值处找到时间序列的周期。
实际上,周期图还可以表作
这里R为时间序列的自协方差函数,即有
这意味着,周期图是对功率谱直接估计的结果。
根据上面几个式子及其内在关系,我们可以利用相关函数估计时间序列的周期。
我们知道,计算相关系数的公式有两种,本例采用下面公式
采用这个公式的原因之一是:
我们的时间序列较短,而这个公式的有效性较好;原因之二是可以直接借助SPSS给出结果。
将结果复制到Excel中间,转换格式以后,再复制到Mathcad的数据表(图19),然后进行Fourier变换。
图19径流量的自相关系数
进行快速Fourier运算的有关定义如下:
根据计算结果容易得到频谱图(图20)。
利用Trace可知,最大功率谱密度对应的频率为f=0.09375,于是得到周期约为P=1/f=10.667(月),与直接采用中心化数据进行Fourier变换得到的结果相同(参见例2)。
图20基于径流量自相关系数的频谱图
前面说过,周期图是对功率谱直接估计的结果,而这种估计通常是不准确的。
但是,在时间序列较短的情况下,利用周期图估计时间序列的周期反而更为可靠。
因此建议,在时间序列较短时运用周期图谱——此时直接采用功率谱反而不准确;在时间序列较长时运用功率谱——此时采用周期图计算工作量较大,而且有误差。
可见虽然理论上二者等价,在实用时却各有适应的时间尺度范围。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 利用 Mathcad 进行 时间 序列 谱分析