1、 1 Excel 进行 FFT 和 Fourier 分析 实例 ： 杭州市 2000 人口分布密度 根据 2000 年人口普查的街道数据经环带 (rings)平均计算得到的结果，数据由冯健博士处理 。下面的变换实质是一种空间自相关的分析过程。 第一步，录入数据 在 Excel 中录入数据不赘述（见 表 1）。 表 1 原始数据序列 表 2 补充后的数据序列 第二步，补充数据 由于 Fourier 变换（ FT）一般是借助快速 Fourier 变换（ Fast Fourier Transformation, FFT）算法，而这种算法的技术过程涉及到 对称处理 ，故数据序列的长度必须是 2N（ N
2、=1,2,3,，）。如果数据序列长度不是 2N，就必须对数据进行补充或者裁减。现在数据长度是 26，介于 24=16到 25=32 之间，而 26 到 32 更近一些，如果裁减数据，就会损失许多信息。因此，采用补充数据的方式。 补充的方法非常简单，在数据序列后面加 0，直到序列长度为 32=25为止（ 表 2）。当然，延续到 64=26也可以，总之必须是 2 的整数倍。不过，补充的“虚拟数据”越多，变换结果的误差也就越大。 2 第三步， Fourier 变换的选项设置 沿着 工具 （ Tools） 数据分析 （ Data Analysis）的路径打开 数据分析 复选框（ 图 1）。 图 1 数
3、据分析（ Data Analysis）的路径 在 数据分析 选项框中选择 傅立叶分析 （ Fourier Analysis）（ 图 2）。 图 2 数据分析（ Data Analysis） 在 Fourier 分析对话框中进行如下设置：在 输入区域 中输入数据序列的单元格范围“ $B$1:$B$33”；选中“标志位于第一行（ L）”；将 输出区域 设为“ $C$2”或 者“ $C$2:$C$33”（ 图 3a）。 a 3 b 图 3 傅立叶分析（ Fourier Analysis） 注意 ： 如果“输入区域”设为 “ $B$2:$B$33”，则不选“标志位于第一行（ L）” （ 图 3b）。
4、表 3 FFT的结果 4 第四步，输出 FFT 结果 选项设置完毕以后，确定（ OK），立即得到 FFT 结果（ 表 3）。 显然，表 3 给出的都是复数（ complex numbers）。假定一个数据序列表为 f(t)，则理论上 Fourier 变换的结果为 dtetfF tj )()( =Ff(t), ( ) 表 3 中给出的正是相应于 F()的复数，这里 为角频率。 第五步，计算功率谱 Excel 好像不能自动计算功率谱，这需要我们利用有关函数进行计算。计算公式为 )(1)(1)( 222 BATFTP 式中 A为复数的 实部 （ real number）， B为 虚部 （ imagi
5、nary number）， T为假设的周期长度，实则补充后的 数据序列长度 。对于本例， T=32。注意复数的平方乃是一个复数与其 共轭（ conjugate）复数的乘积，若 F()=a+bj，则 |F()|2=(a+bj)*(a-bj)=a2+b2。这样，根据表 3 中的 FFT 结果，我们有 149470319632/)0857.218701( 22 67510894932/)538.103400634.104459( 22 其余依次类推。 显然，这样计算非常繁琐。一个简单的办法是调用 Excel 的 模数 （ modulus） 计算函数ImAbs，方法是在函数类别中找“其他”，在其他类中
6、找“工程”类，在工程类中容易找到ImAbs 函数（ 图 4）。 确定以后，弹出一个选项框，选中第一个 FFT 结果，确定，得到 218701.857（ 图 5）。我们知道，复数的模数计算公式 为 2/122 )( BAM 图 4 模数计算函数 5 对于第一个 FFT 结果，由于虚部为 0，模数就是其自身，即 857.218701)0857.218701( 2/122 但对于后面真正的复数，就不一样了。抓住第一个模数所在的单元格的右下角往下一拉，或者用鼠标双击该单元格的右下角，立即得到全部模数。 图 5 计算模数 最后，用 模数的 2 次方 除以 数据长度 32 立即得到全部 功率谱密度 结果（
7、 表 4）。 表 4 功率谱密度 下表是利用 Mathcad2000 计算的功率谱密度（ 表 5）。利用 Mathcad 进行 FFT，过程要简单得多，只要调用 FFT 命令，可以直接给出各种结果（包括图表）。但 Mathcad 的计算不求精度，有一定误差。将 Mathcad 的变换结果 copy 到 Excel 中进行比较，可以看到，如果6 不计误差，二者是一致的（ 表 4）。 表 5 借助 Mathcad2000进行 FFT的结果 P o w e r0 10123456789101112131415161 . 4 9 5 ? 0 96 . 7 5 1 ? 0 82 . 9 4 8 ? 0
8、81 . 0 2 6 ? 0 83 . 1 8 8 ? 0 72 . 4 7 6 ? 0 72 . 9 8 5 ? 0 72 . 4 4 6 ? 0 71 . 1 7 2 ? 0 76 . 2 3 8 ? 0 68 . 9 0 8 ? 0 61 . 0 7 3 ? 0 71 . 0 4 2 ? 0 79 . 4 2 ? 0 66 . 4 9 4 ? 0 64 . 7 9 ? 0 64 . 6 9 7 ? 0 6第六步，功率谱分析 功率谱分析目前主要用于两个方面，一是侦测系统变化的某种 周期 或者 节律 ，据此寻找因果关系 （解释）或者进行某种 发展预测 （应用）；二是寻找周期以外的 某些规律
9、，据此对系统的时空结构特征进行解释。 表 6 以对称点（ f=0.5）为界，从完整的数据序列中截取一半 7 上面基于杭州人口密度数据的 FFT，实际上是一种空间自相关分析过程，属于 FT 的第二类应用。这种过程不以寻找周期为目标，实际上也不存在任何周期。 不论目标是什么，都必须借助 频谱图 （频率功率谱密度图）进行分析和解释。下面第一步就是绘制频谱图。首先要计算频率， 线频 或 角频 都可以，因为二者相差常数倍（ 2）。一个简单的办法是，用 0 到 T=32 的自然数列除以 T=32（ 表 6）。 如果采用的频率变化范围 01，则绘制的频谱图是对称的（ 图 6）。实际上，另一半是多余的， Ma
10、thcad2000 自动生成的频谱图就没有考虑另外一半儿（ 图 7）。因此，我们可以以 对称点 f=0.5 为界，截取前面一半的数据，在 Excel 上绘制频谱图（ 图 8）。 020000000040000000060000000080000000010000000001200000000140000000016000000000 0.2 0.4 0.6 0.8 1频率功率谱密度图 6 对称的频谱图（基于完整的数据序列） 0 0 . 1 0 . 2 0 . 3 0 . 4 0 . 505 1081 1091 . 5 109杭 州 人 口 密 度 衰 减 的 频 谱 图 （ 2 0 0 0 ）
11、频率功率谱密度1 . 4 9 5 1094 . 6 9 7 106P o w e r j0 . 50 f r e q j图 7 Mathcad2000生成的频谱图 下图是常用的频谱图形式，如果存在周期，则在尖峰突出的最大点可以找到。这个图中是没有显示任何周期的，但并不意味着没有重要信息。在理论上，如果人口密度分布服从负指数模型，则其频率与 功率谱之间应该满足如下关系 2)( ffP 8 为了检验这种推断，不妨用下式进行拟合 ffP )( 这正是 噪声（ -noise）表达式。 02000000004000000006000000008000000001000000000120000000014
12、0000000016000000000 0.1 0.2 0.3 0.4 0.5频率功率谱密度图 8 利用 Excel绘制的频谱图（常用形式） 为了拟合幂指数模型，去掉 0 频率点，结果得到 7983.11280514)( ffP , R2=0.9494 多种模型比较的结果，发现幂指数模型的拟合效果最好（ 图 9）。将 图 9 转换成对数刻度，拟合效果就 尤其明确（ 图 10）。显然， =1.7983 2。 P ( f ) = 1 2 8 0 5 1 4 . 1 7 9 5 f- 1. 7983R2= 0 . 9 4 9 4010000000020000000030000000040000000
13、05000000006000000007000000008000000000 0.1 0.2 0.3 0.4 0.5频率功率谱密度图 9 频谱图的模型拟合结果（去掉 0频点） 9 P ( f ) = 1 2 8 0 5 1 4 . 1 7 9 5 f- 1. 79 83R2= 0 . 9 4 9 41.E + 061.E + 071.E + 081.E + 090.01 0.1 1频率功率谱密度图 10 双对数频谱图 利用模型及其参数，我们可以对杭州市人口分布特征及其变化进行系统分析。但是，深入的分析仅仅借助一个参数是不够的。具体的分析过程将用专门的文章进行论述。 最后说明一点：前面的公式 2)(1)( FTP 给出的是功率谱。有时在理论上进行讨论时，采用下式 2)()( FS , 这里给出的是能量谱。能量谱的计算假定数据序列无穷长，积分范围一般从 负无穷 到 正无穷 ；功率谱主要用于对实际遇到的有限长度的数据。二者在数值上相差常数倍。因此，在理论讨论时，采用能量谱公式比较方便；在具体应用时采用功率谱公式便于比较。二者的数理本质是一致的，故一般行文过程中无需澄清二者的关系。
