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

1、利用 Excel 进行 FFT 和 Fourier分析的基本步骤实例 ：杭州市 2000 人口分布密度根据 2000 年人口普查的街道数据经环带(rings) 平均计算得到的结果，数据由冯健博士处理。下面的变换实质是一种空间自相关的分析过程。第一步，录入数据在 Excel 中录入数据不赘述（见表 1）。表 1 原始数据序列表 2 补充后的数据序列第二步，补充数据由于 Fourier 变换（FT ）一般是借助快速 Fourier 变换（ Fast Fourier Transformation, FFT ）算法，而这种算法的技术过程涉及到 对称处理 ，故数据序列的长度必须是 2N（N=1,2,3,

2、，）。如果数据序列长度不是 2N，就必须对数据进行补充或者裁减。 现在数据长度是 26，介于 24=16到 25=32 之间，而 26 到 32 更近一些，如果裁减数据，就会损失许多信息。因此，采用补充数据的方式。补充的方法非常简单， 在数据序列后面加 0，直到序列长度为 32=25 为止（表 2）。当然，延续到 64=26 也可以，总之必须是 2 的整数倍。不过，补充的“虚拟数据”越多，变换结果的误差也就越大。1第三步， Fourier变换的选项设置沿着 工具（ Tools） 数据分析 （ Data Analysis ）的路径打开数据分析 复选框（ 图 1）。图 1 数据分析（ Data A

3、nalysis ）的路径在数据分析 选项框中选择 傅立叶分析 （Fourier Analysis ）（图 2）。图 2 数据分析（ Data Analysis ）在 Fourier分析对话框中进行如下设置：在输入区域 中输入数据序列的单元格范围“$B$1:$B$33 ”；选中“标志位于第一行 （L）”；将输出区域 设为“$C$2” 或者“$C$2:$C$33”（图 3a）。a2b图 3 傅立叶分析（ Fourier Analysis）注意 ：如果“输入区域”设为“$B$2:$B$33” ，则不选“标志位于第一行（L） ”（ 图 3b）。表 3 FFT 的结果3第四步，输出FFT 结果选项设置完

4、毕以后，确定（OK ），立即得到FFT 结果（ 表 3）。显然，表3 给出的都是复数（complex numbers ）。假定一个数据序列表为f(t)，则理论上 Fourier 变换的结果为- j t=Ff(t), (-)F () =f (t )edt< < - 表 3 中给出的正是相应于F()的复数，这里 为角频率。第五步，计算功率谱Excel 好像不能自动计算功率谱，这需要我们利用有关函数进行计算。计算公式为P() = 1 F () 2 = 1 (A 2 + B 2 )TT式中 A 为复数的 实部 （ real number）， B 为虚部 （ imaginary number

5、 ），T 为假设的周期长度，实则补充后的数据序列长度。对于本例，T=32 。注意复数的平方乃是一个复数与其共轭（conjugate ）复数的乘积，若 F()=a+bj ，则 |F()|2=(a+bj)*( a-bj)= a2+b2。这样，根据表 3 中的 FFT 结果，我们有( 218701.857 2 + 02 ) / 32 = 1494703196(104459.634 2 +103400 .538 2 ) / 32 = 675108949其余依次类推。显然，这样计算非常繁琐。一个简单的办法是调用Excel 的模数 （modulus ）计算函数ImAbs ，方法是在函数类别中找“其他” ，

6、在其他类中找“工程”类，在工程类中容易找到 ImAbs 函数（ 图 4）。确定以后，弹出一个选项框，选中第一个FFT 结果，确定，得到218701.857 （图 5）。我们知道，复数的模数计算公式为M = (A2 + B2)1/2图 4 模数计算函数4对于第一个FFT 结果，由于虚部为0，模数就是其自身，即(218701.8572 + 0 2 )1 / 2 = 218701.857但对于后面真正的复数，就不一样了。 抓住第一个模数所在的单元格的右下角往下一拉，或者用鼠标双击该单元格的右下角，立即得到全部模数。图 5 计算模数最后，用 模数的 2 次方除以 数据长度32 立即得到全部 功率谱密度

7、 结果（ 表 4）。表 4 功率谱密度下表是利用Mathcad2000 计算的功率谱密度（表 5）。利用 Mathcad 进行 FFT，过程要简单得多，只要调用FFT 命令，可以直接给出各种结果（包括图表）。但 Mathcad 的计算不求精度，有一定误差。将Mathcad 的变换结果copy 到 Excel 中进行比较，可以看到，如果5不计误差，二者是一致的（表 4）。表 5 借助 Mathcad2000 进行 FFT的结果010 1.495?0 91 6.751?0 82 2.948?0 83 1.026?0 84 3.188?0 75 2.476?0 76 2.985?0 7Power =

8、72.446?0 78 1.172?0 79 6.238?0 610 8.908?0 611 1.073?0 712 1.042?0 713 9.42?0 614 6.494?0 615 4.79?0 616 4.697?0 6第六步，功率谱分析功率谱分析目前主要用于两个方面， 一是侦测系统变化的某种 周期 或者 节律，据此寻找因果关系 （解释）或者进行某种 发展预测 （应用）；二是寻找周期以外的 某些规律 ，据此对系统的时空结构特征进行解释。表 6 以对称点（ f=0.5 ）为界，从完整的数据序列中截取一半6上面基于杭州人口密度数据的FFT，实际上是一种空间自相关分析过程，属于FT 的第二类

9、应用。这种过程不以寻找周期为目标，实际上也不存在任何周期。不论目标是什么，都必须借助频谱图 （频率功率谱密度图）进行分析和解释。下面第一步就是绘制频谱图。首先要计算频率，线频 或角频 都可以，因为二者相差常数倍（2）。一个简单的办法是，用0 到 T=32 的自然数列除以T=32（表 6）。如果采用的频率变化范围01 ，则绘制的频谱图是对称的（图 6）。实际上，另一半是多余的， Mathcad2000 自动生成的频谱图就没有考虑另外一半儿（图 7）。因此，我们可以以对称点 f=0.5 为界，截取前面一半的数据，在Excel 上绘制频谱图（图 8）。160000000014000000001200

10、000000度 1000000000密谱 800000000率功 600000000400000000200000000000.20.40.60.81频率图 6 对称的频谱图（基于完整的数据序列）991.5 .101.495×10.91 10度密谱Power j率功8.5 104.697×106000.10.20.30.40.50freq j0.5频率杭州人口密度衰减的频谱图（2000）图 7 Mathcad2000 生成的频谱图下图是常用的频谱图形式， 如果存在周期， 则在尖峰突出的最大点可以找到。 这个图中是没有显示任何周期的， 但并不意味着没有重要信息。 在理论上，

11、如果人口密度分布服从负指数模型，则其频率与功率谱之间应该满足如下关系P ( f ) f - 27为了检验这种推断，不妨用下式进行拟合P( f ) f这正是 噪声（ -noise）表达式。160000000014000000001200000000度 1000000000密谱 800000000率功 6000000004000000002000000000- 00.10.20.30.40.5频率图 8 利用 Excel 绘制的频谱图（常用形式）为了拟合幂指数模型，去掉0 频率点，结果得到P ( f ) = 1280514 f - 1 .7983 , R2=0.9494多种模型比较的结果，发现幂指

12、数模型的拟合效果最好（图 9）。将 图 9 转换成对数刻度，拟合效果就尤其明确（图 10）。显然， =1.7983 2。800000000700000000-1.7983P (f ) = 1280514.1795 f600000000R2 = 0.9494度 500000000密谱 400000000率功 300000000200000000100000000000.10.20.30.40.5频率图 9 频谱图的模型拟合结果（去掉0 频点）81.E+091.E+08度密谱率功1.E+07-1.7983P (f ) = 1280514.1795 f2R = 0.94941.E+060.010.11频率图 10 双对数频谱图利用模型及其参数，我们可以对杭州市人口分布特征及其变化进行系统分析。但是， 深入的分析仅仅借助一个参数是不够的。具体的分析过程将用专门的文章进行论述。最后说明一点：前面的公式P()