利用Excel进行FFT和Fourier分析的基本步骤

表 6 以对称点（f=0.5）为界，从完整的数据序列中截取一半
6
上面基于杭州人口密度数据的 FFT，实际上是一种空间自相关分析过程，属于 FT 的第 二类应用。这种过程不以寻找周期为目标，实际上也不存在任何周期。
不论目标是什么，都必须借助频谱图（频率－功率谱密度图）进行分析和解释。下面第 一步就是绘制频谱图。首先要计算频率，线频或角频都可以，因为二者相差常数倍（2π）。 一个简单的办法是，用 0 到 T=32 的自然数列除以 T=32（表 6）。
如果采用的频率变化范围 0~1，则绘制的频谱图是对称的（图 6）。实际上，另一半是多 余的，Mathcad2000 自动生成的频谱图就没有考虑另外一半儿（图 7）。因此，我们可以以对 称点 f=0.5 为界，截取前面一半的数据，在 Excel 上绘制频谱图（图 8）。
功率谱密度
1600000000
1400000000
9 6.238?06
10 8.908?06
11 1.073?07
12 1.042?07
13 9.42?06
14 6.494?06
15 4.79?06
16 4.697?06
第六步，功率谱分析 功率谱分析目前主要用于两个方面，一是侦测系统变化的某种周期或者节律，据此寻找
因果关系（解释）或者进行某种发展预测（应用）；二是寻找周期以外的某些规律，据此对 系统的时空结构特征进行解释。
a
2
b 图 3 傅立叶分析（Fourier Analysis） 注意：如果“输入区域”设为“$B$2:$B$33”，则不选“标志位于第一行（L）”（图 3b）。
表 3 FFT 的结果
3
第四步，输出 FFT 结果 选项设置完毕以后，确定（OK），立即得到 FFT 结果（表 3）。 显然，表 3 给出的都是复数（complex numbers）。假定一个数据序列表为 f(t)，则理论
0.2
0.3
0.4
0.5
0
freq j
0.5
频率
杭州人口密度衰减的频谱图（2000）
图 7 Mathcad2000 生成的频谱图
下图是常用的频谱图形式，如果存在周期，则在尖峰突出的最大点可以找到。这个图中 是没有显示任何周期的，但并不意味着没有重要信息。在理论上，如果人口密度分布服从负 指数模型，则其频率与功率谱之间应该满足如下关系
最后说明一点：前面的公式
P(ω) = 1 F (ω) 2 T
给出的是功率谱。有时在理论上进行讨论时，采用下式
S(ω ) = F (ω) 2 ,
这里给出的是能量谱。能量谱的计算假定数据序列无穷长，积分范围一般从负无穷到正无穷； 功率谱主要用于对实际遇到的有限长度的数据。二者在数值上相差常数倍。因此，在理论讨 论时，采用能量谱公式比较方便；在具体应用时采用功率谱公式便于比较。二者的数理本质 是一致的，故一般行文过程中无需澄清二者的关系。
图 1 数据分析（Data Analysis）的路径 在数据分析选项框中选择傅立叶分析（Fourier Analysis）（图 2）。
图 2 数据分析（Data Analysis） 在 Fourier 分析对话框中进行如下设置：在输入区域中输入数据序列的单元格范围 “$B$1:$B$33”；选中“标志位于第一行（L）”；将输出区域设为“$C$2”或者“$C$2:$C$33” （图 3a）。
确定以后，弹出一个选项框，选中第一个 FFT 结果，确定，得到 218701.857（图 5）。 我们知道，复数的模数计算公式为
M = (A2 + B 2 )1/ 2
图 4 模数计算函数
4
对于第一个 FFT 结果，由于虚部为 0，模数就是其自身，即 (218701.857 2 + 0 2 )1/ 2 = 218701.857
0 0
0.1
0.2
0.3
0.4
0.5
频率
图 8 利用 Excel 绘制的频谱图（常用形式）
为了拟合幂指数模型，去掉 0 频率点，结果得到 P( f ) = 1280514 f −1.7983 , R2=0.9494
多种模型比较的结果，发现幂指数模型的拟合效果最好（图 9）。将图 9 转换成对数刻度， 拟合效果就尤其明确（图 10）。显然，β=1.7983≠2。
补充的方法非常简单，在数据序列后面加 0，直到序列长度为 32=25 为止（表 2）。当然， 延续到 64=26 也可以，总之必须是 2 的整数倍。不过，补充的“虚拟数据”越多，变换结果 的误差也就越大。
1
第三步，Fourier 变换的选项设置 沿着工具（Tools）→数据分析（Data Analysis）的路径打开数据分析复选框（图 1）。
的 FFT 结果，我们有
(218701.857 2 + 02 ) / 32 = 1494703196
(104459.6342 +103400.5382 ) / 32 = 675108949
其余依次类推。 显然，这样计算非常繁琐。一个简单的办法是调用 Excel 的模数（modulus）计算函数
ImAbs，方法是在函数类别中找“其他”，在其他类中找“工程”类，在工程类中容易找到 ImAbs 函数（图 4）。
P(ω) = 1 F (ω) 2 = 1 ( A2 + B 2 )
T
T
式中 A 为复数的实部（real number），B 为虚部（imaginary number），T 为假设的周期长度，
实则补充后的数据序列长度。对于本例，T=32。注意复数的平方乃是一个复数与其共轭
（conjugate）复数的乘积，若 F(ω)=a+bj，则|F(ω)|2=(a+bj)*(a-bj)=a2+b2。这样，根据表 3 中
表 2 补充后的数据序列
第二步，补充数据 由于 Fourier 变换（FT）一般是借助快速 Fourier 变换（Fast Fourier Transformation, FFT）
算法，而这种算法的技术过程涉及到对称处理，故数据序列的长度必须是 2N（N=1,2,3,…，）。 如果数据序列长度不是 2N，就必须对数据进行补充或者裁减。现在数据长度是 26，介于 24=16 到 25=32 之间，而 26 到 32 更近一些，如果裁减数据，就会损失许多信息。因此，采用补充 数据的方式。
但对于后面真正的复数，就不一样了。抓住第一个模数所在的单元格的右下角往下一拉，或 者用鼠标双击该单元格的右下角，立即得到全部模数。
图 5 计算模数 最后，用模数的 2 次方除以数据长度 32 立即得到全部功率谱密度结果（表 4）。
表 4 功率谱密度
下表是利用 Mathcad2000 计算的功率谱密度（表 5）。利用 Mathcad 进行 FFT，过程要 简单得多，只要调用 FFT 命令，可以直接给出各种结果（包括图表）。但 Mathcad 的计算不 求精度，有一定误差。将 Mathcad 的变换结果 copy 到 Excel 中进行比较，可以看到，如果
9
功率谱密度
800000000 700000000 600000000 500000000 400000000 300000000 200000000 100000000
0 0
P (f ) = 1280514.1795 f -1.7983 R2 = 0.9494
0.1
0.2
0.3
0.4
0.5
频率
图 9 频谱图的模型拟合结果（去掉 0 频点）
8
1.E+09
1.E+08
功率谱密度
P (f ) = 1280514.1795 f -1.7983 R2 = 0.9494
0.01
0.1
频率
图 10 双对数频谱图
1.E+07
1.E+06 1
利用模型及其参数，我们可以对杭州市人口分布特征及其变化进行系统分析。但是，深 入的分析仅仅借助一个参数是不够的。具体的分析过程将用专门的文章进行论述。
1200000000
1000000000
800000000
600000000
400000000
200000000
0
0
0.2
0.4
0.6
0.8
1
频率
图 6 对称的频谱图（基于完整的数据序列）
1.5 .109 1.495×109
功率谱密度
1 .109 Powerj
5 .108
4.697×106
0
0
0.1
利用 Excel 进行 FFT 和 Fourier 分析的基本步骤
实例：杭州市 2000 人口分布密度 [根据 2000 年人口普查的街道数据经环带(rings)平均 计算得到的结果，数据由冯健博士处理]。下面的变换实质是一种空间自相关的分析过程。
第一步，录入数据 在 Excel 中录入数据不赘述（见表 1）。 表 1 原始数据序列
5
不计误差，二者是一致的（表 4）。 表 பைடு நூலகம் 借助 Mathcad2000 进行 FFT 的结果
0
1
0 1.495?09
1 6.751?08
2 2.948?08
3 1.026?08
4 3.188?07
5 2.476?07
6 2.985?07
Power = 7 2.446?07 8 1.172?07
P( f ) ∝ f −2
7
为了检验这种推断，不妨用下式进行拟合 P( f ) ∝ f −β
这正是 β 噪声（β-noise）表达式。
功率谱密度
1600000000 1400000000 1200000000 1000000000
800000000 600000000 400000000 200000000