利用Excel进行时间序列的谱分析(I)

60000
50000
功率P(f)
40000
30000
20000
10000
0
0
0.1
0.2
0.3
0.4
0.5
0.6
频率f
图 17 实用的频谱图 第五步，功率谱分析
频域分析的主要目标之一是判断时间序列是否存在周期性。从图 17 可以看出，功率最 大点对应的频率是 f1=0.0625，该频率对应的周期长度为 16。可见，在时间序列较短的情况 下之间用功率谱图寻找时间序列的周期不如周期图准确。
120000
100000
频率强度
80000
60000
40000
20000
0
0
0.1
0.2
0.3
0.4
0.5
0.6
频率
图 8 某河流径流量的周期图（1961 年 6 月－1963 年 5 月）
2 时间序列的频谱图
【例 3】首先考虑对例 1 的数据进行功率谱分析。例 1 的时间序列较短，分析的效果不 佳，但计算过程简短。给出这个例子，主要是帮助大家理解 Fourier 变换过程和方法。为了 进行 Fourier 分析，需要对数据进行预处理。第一，将数据中心化，即用原始数据减去其平 均值。中心化的数据均值为 0，我们对中心化的数据进行变换，其周期更为明显。第二，由 于 Fourier 分析通常采用所谓快速 Fourier 变换（Fast Fourier Transformation，FFT），而 FFT 要求数据必须为 2n 个，这里 n 为正整数（1,2,3,…），而我们的样本为 N=12，它不是 2 的某
将上面公式中的余弦值换成正弦值，即可得到 bi 值（见下表）。上面的计算过程相当于
2
采用式（2）进行逐步计算。
第四步，计算频率强度
利用式（3），非常容易算出 I(fi)值。例如
I( f1)
=
N 2
(a12
+ b12 )
=
6 * (6.5972
+ 170.2132 )
= 174096.914
其余依此类推（见图 4）。
图 10 在数据分析选项框中选择 Fourier 分析 第二步，定义变量和输出区域
确定之后，弹出傅立叶分析对话框，根据数据在工作表中的分布情况进行如下设置： 将光标置于“输入区域”对应的空白栏，然后用鼠标选中单元格 C1-C17，这时空白栏 中自动以绝对单元格的形式定义中心化数据的区域范围（即$C$1-$C$17）。 选中“标志位于第一行”。 选中输出区域，定义范围为 D2-D17（图 11）。 注意：如果输入区域的数据范围定义为 C2-C17，则不要选中“标志位于第一行”，这与 回归分析中的原始变量定义是一样的（图 12）。如果不定义输出区域范围，则变换结果将会 自动给在新的工作表组上。这一点也与回归分析一样。
最大值处找到时间序列的周期。对于本例，N=12，t=1,2,…,N，fi=i/N，下面借助 Excel，利
用上述公式，计算有关参数并分析时间序列的周期特性。
第一步，录入数据，并将数据标准化或中心化（图 1）。
图 1 录入的数据及其中心化结果
1
中心化与标准化的区别在于，只需将原始数据减去均值，而不必再除以标准差。不难想 到，中心化的数据均值为 0，但方差与原始数据相同（未必为 1）。 第二步，计算三角函数值
S(ω) = F (ω)F (ω) = F (ω) 2
(4)
相应地，有了 XT(f)也就容易计算功率谱（密度）
P( f ) = XT ( f ) 2
(5)
T
式中 f 表示线频率，与角频率 ω 的转换关系是 ω=2π/T，这里 T 为数据区间长度。 如果将 XT(f)表作 XT(f)=A+jB（这里 A 为实部，B 为虚部），则有
图 14 功率谱密度的计算结果
8
第四步，绘制频谱图 线频率 fi 可以表作
fi = i / T = i / N ， i = 0,1,2,Λ , N -1
显然 f0=0/16=0，f1=1/16=0.0625，f2=2/16=0.125，…，f15=15/16=0.9375。在 Excel 中，容易 计算频率的数值。将频率与功率谱对应起来（图 15），就可以画出频谱图。如果补上最后一 个频率数值 f16=1 及其对应的功率 0，则可画出完全对称的谱图（图 16）。
图 3 计算正弦值的表格 第三步，计算参数 ai、bi
利用中心化的数据（仍然表作 xt）计算参数 ai、bi。首先算出 xtcos2πfit 和 xtsins2πfit。在 D9 单元格中输入公式“=D1*D3”，回车得到 18.309；按住单元格的右下角右拉至 O9 单元 格，得到 f=1/12=0.083，t=1,2,…,12 时的全部 xtcos2πfit 值；加和得 39.584，再除以 6，即得 a1=6.597。在 D10 单元格中输入公式“=D1*D4”，回车得到 10.571；按住单元格的右下角右 拉至 O10 单元格，得到 f=2/12=0.083，t=1,2,…,12 时的全部 xtcos2πfit 值；加和得-365.25，再 除以 6，得到 a2=-60.875。其余依此类推。
功率谱密度P(f)
70000
图 15 功率谱密度与频率的对应关系
60000
50000
40000
30000
20000
10000
0
0
0.2
0.4
0.6
0.8
1
频率f
图 16 对称分布的频谱图
由于功率谱图的对称性，画出整个谱图实在没有必要，因此，在实际工作中，通常只画 出左半边（图 17）。
9
70000
为了借助式（1）计算参数 ai、bi，首先需要计算正弦值和余弦值。
取 i = 1,2,Λ ,6 ，则频率为 fi = i / N = 1/12,2 /12,Λ ,6 /12 （图 1）。
将频率写在单元格 C3-C14 中（根据对称性，我们只用前 6 个），将中心化的数据转置 粘贴于第一行的单元格 D1-O1 中，月份的序号写在单元格 D2-O2 中（与中心化数据对齐）。
图 4 计算频率强度 第五步，绘制时间序列周期图
利用图 4 中的数据，不难画出周期图（图 5）。
频率强度
200000 180000 160000 140000 120000 100000 80000 60000 40000 20000
0 0
0.1
0.2
0.3
0.4
0.5
0.6
频率
图 5 某河流径流量的周期图（1979 年 6 月－1980 年 5 月） 第六步，周期识别
P( fi ) =
Ai2 + Bi2 T
（7）
考虑到本例中 T=N=16，在 J2 单元格中输入公式“=(H2^2+I2^2)/16”，回车得到第一个功率 谱密度 0；下拉至 J17，得到全部谱密度数值（图 14）。基于 FFT 的谱密度分布是对称的， 可以看出，以 J10 所在的 3105.28 为对称点，上下的数值完全对称。
图 2 计算余弦值的表格 在 D2 单元格中输入公式“=COS($B$1*$D$2*C3)”，回车得到 0.866；按住单元格的右 下角右拉至 O3 单元格，得到 f=1/12=0.083，t=1,2,…,12 时的全部余弦值。在 D2 单元格中输 入公式“=COS($B$1*$D$2*C4)”，回车得到 0.5；按住单元格的右下角右拉至 O4 单元格， 得到 f=2/12=0.167，t=1,2,…,12 时的全部余弦值。依次类推，可以算出全部所要的余弦值（在 D3-O8 区域中）。根据对称性，我们的计算到 k=6 为止（图 2）。注意，这里 B1 单元格是 2π=6.28319（图中未能显示）。 在上面的计算中，只要将公式中的“COS”换成“SIN”，即可得到正弦值，不过为了 计算过程清楚明白，最好在另外一个区域给出结算结果（在 D17-O22 区域中，参见图 3）。
5
个 n 次方。因此，在中心化的数据后面加上 4 个 0，这样新的样本数为 N′=12＋4＝16＝24 个，这才符合 FFT 的需要（图 9）。下面，我们对延长后的中心化数据进行 Fourier 变换。
图 9 数据的中心化与“延长” 第一步，打开 Foureir 分析对话框
沿着主菜单的“工具（Tools）”→“数据分析（Data Analysis）”路径打开数据分析选项 框（图 10），从中选择“傅立叶分析（Fourier Analysis）”。
图 7 参数和频率强度的计算结果
从图 8 中可以看出，频率强度的最大值（极值点）对应于频率 f1=1/12=0.0833，故时间 序列的周期判断为 T=1/f1=12。这与用 12 月的数据进行估计的结果是一致的，但由于例 2 的 时间序列比例 1 的时间序列长 1 倍，故判断结果更为可靠。
140000
分析：假定将时间序列 xt 展开为 Fourier 级数，则可表示为
Hale Waihona Puke Baidu
k
∑ xt = (ai cos 2πfit + bi sin 2πfit) + ε t
(1)
i =1
式中 fi 为频率，t 为时间序号，k 为周期分量的个数即主周期（基波）及其谐波的个数，εt 为标准误差（白噪声序列）。当频率 fi 给定时，式（1）可以视为多元线性回归模型，可以证 明，待定系数 ai、bi 的最小二乘估计为
X T ( fi ) 2 = X T ( fi ) X T ( fi ) = ( Ai + jBi )( Ai − jBi ) = Ai2 + Bi2
（6）
因此这一步是要分离变换结果的实部与虚部。逐个手动提取是非常麻烦而且容易出错 的，可以利用 Excel 有关复数计算的命令。提取实部的 Excel 命令是 imreal。在 H2 单元格 中输入命令“=IMREAL(D2)”（这里 D2 为变换结果的第一个复数所在的单元格），回车得到 第一个复数的实部 0；点中 H2 单元格的右下角，揿住鼠标左键下拉至 H17，得到全部的实 部数值。提取虚部的命令是 imaginary。在 I2 单元格中输入公式“=IMAGINARY(D2)”，回 车得到第一个复数的虚部 0；下拉至 I17，得到全部的虚部数值。根据式（5）、（6），功率谱 密度的计算公式为
图 6 原始数据及部分处理结果
将原始数据回车时间序列变化图，可以初步估计具有 12 月变化周期，但不能肯定（图 6）。
350
300 250
月平均径流量
200
150
100
50
0
0
5
10
15
20
25
30
时间（月份序号）
图 6 径流量的月变化图（1961 年 6 月－1963 年 5 月）
4
按照例 1 给出的计算步骤，计算参数数 ai、bi，进而计算频率强度（结果将图 7）。然后 绘制时间序列的周期图（图 8）。注意这里，N=24，我们取 k=12。
∑ aˆi
=
2 N
N t =1
xt
cos 2πfit
∑ bˆi
=
2 N
N
xt sin 2πfit
t =1
（2）
这里 N 为观测值的个数。定义时间序列的周期图为
I( fi)
=
N 2
(ai2
+ bi2 ) ， i
= 1,2,Λ
,k
(3)
式中 I(fi)为频率 fi 处的强度。以 fi 为横轴，以 I(fi)为纵轴，绘制时间序列的周期图，可以在
6
图 11 选中“标志位于第一行”与数据输入范围的定义
图 12 不选中“标志位于第一行”与数据输入范围的定义 第三步，结果转换
定义数据输入－输出区域完成之后，确定，立即得到 Fourier 变换的结果（图 13）。
图 13 傅立叶变换的结果
7
变换的结果为一组复数，相当于将 f(t)变成了 F(ω)，实际上是将 xt 变成了 XT(f)。我们知 道，有了 f(t)的象函数 F(ω)就可以计算能量谱密度函数 S(ω)，即有
关键是寻找频率的极值点或突变点。在本例中，没有极值点，但在 f1=1/12=0.0833 处， 频率强度突然增加（陡增），而此时 T=1/f1=12，故可判断时间序列可能存在一个 12 月的周 期，即 1 年周期。
3
【例 2】为了映证上述判断，我们借助同一条河流的连续两年的平均月径流量（1961 年 6 月－1963 年 5 月）。原始数据见下图（图 6）。
利用 Excel 进行时间序列的谱分析（I）
在频域分析中，功率谱是揭示时间序列周期特性的最为有力的工具之一。下面列举几个 例子，分别从不同的角度识别时间序列的周期。
1 时间序列的周期图
【例 1】某水文观测站测得一条河流从 1979 年 6 月到 1980 年 5 月共计 12 月份的断面平均 流量。试判断该河流的径流量变化是否具有周期性，周期长度大约为多少？