《利用Excel进行FFT和Fourier分析的基本步骤》由会员分享,可在线阅读,更多相关《利用Excel进行FFT和Fourier分析的基本步骤(9页珍藏版)》请在金锄头文库上搜索。
1、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 的整数倍。不过,补充的“虚拟数据”越多,变换结果的误差也就越大。2第三步,Fourier 变换的选项设置沿着工具(Tools)数据分析 (Data Analysis)的路径打开数据分析复选框(图 1) 。图 1 数据分析(Data Analys
3、is)的路径在数据分析选项框中选择傅立叶分析(Fourier Analysis) (图 2) 。图 2 数据分析(Data Analysis)在 Fourier 分析对话框中进行如下设置:在输入区域中输入数据序列的单元格范围“$B$1:$B$33”;选中“标志位于第一行(L) ”;将输出区域设为“$C$2”或者“$C$2:$C$33”(图 3a) 。a 3b图 3 傅立叶分析(Fourier Analysis)注意:如果“输入区域”设为“$B$2:$B$33” ,则不选“标志位于第一行(L) ”(图 3b) 。表 3 FFT的结果4第四步,输出 FFT 结果选项设置完毕以后,确定(OK) ,立
4、即得到 FFT 结果(表 3) 。显然,表 3 给出的都是复数(complex numbers) 。假定一个数据序列表为 f(t),则理论上 Fourier 变换的结果为=Ff(t), ( )detfFtj)()( 表 3 中给出的正是相应于 F()的复数,这里 为角频率。第五步,计算功率谱Excel 好像不能自动计算功率谱,这需要我们利用有关函数进行计算。计算公式为 )(1)()( 22BATFP式中 A 为复数的实部(real number) ,B 为虚部(imaginary number) ,T 为假设的周期长度,实则补充后的数据序列长度。对于本例,T=32。注意复数的平方乃是一个复数与
5、其共轭(conjugate)复数的乘积,若 F()=a+bj,则|F()| 2=(a+bj)*(a-bj)=a2+b2。这样,根据表 3中的 FFT 结果,我们有 1497036/)0857.21(2582/.364.09其余依次类推。显然,这样计算非常繁琐。一个简单的办法是调用 Excel 的模数(modulus)计算函数ImAbs,方法是在函数类别中找 “其他” ,在其他类中找“ 工程”类,在工程类中容易找到ImAbs 函数(图 4) 。确定以后,弹出一个选项框,选中第一个 FFT 结果,确定,得到 218701.857(图 5) 。我们知道,复数的模数计算公式为 2/12)(BAM图 4
6、 模数计算函数5对于第一个 FFT 结果,由于虚部为 0,模数就是其自身,即 857.210)857.21(/2但对于后面真正的复数,就不一样了。抓住第一个模数所在的单元格的右下角往下一拉,或者用鼠标双击该单元格的右下角,立即得到全部模数。图 5 计算模数最后,用模数的 2 次方除以数据长度 32 立即得到全部功率谱密度结果(表 4) 。表 4 功率谱密度下表是利用 Mathcad2000 计算的功率谱密度(表 5) 。利用 Mathcad 进行 FFT,过程要简单得多,只要调用 FFT 命令,可以直接给出各种结果(包括图表) 。但 Mathcad 的计算不求精度,有一定误差。将 Mathca
7、d 的变换结果 copy 到 Excel 中进行比较,可以看到,如6果不计误差,二者是一致的(表 4) 。表 5 借助 Mathcad2000进行 FFT的结果Power0 10123456789101112131415161.49510 96.75110 82.94810 81.02610 83.18810 72.47610 72.98510 72.44610 71.17210 76.23810 68.90810 61.07310 71.04210 79.4210 66.49410 64.7910 64.69710 6第六步,功率谱分析功率谱分析目前主要用于两个方面,一是侦测系统变化的某种周
8、期或者节律,据此寻找因果关系(解释)或者进行某种发展预测(应用) ;二是寻找周期以外的某些规律,据此对系统的时空结构特征进行解释。表 6 以对称点(f=0.5)为界,从完整的数据序列中截取一半7上面基于杭州人口密度数据的 FFT,实际上是一种空间自相关分析过程,属于 FT 的第二类应用。这种过程不以寻找周期为目标,实际上也不存在任何周期。不论目标是什么,都必须借助频谱图(频率功率谱密度图)进行分析和解释。下面第一步就是绘制频谱图。首先要计算频率,线频或角频都可以,因为二者相差常数倍(2) 。一个简单的办法是,用 0 到 T=32 的自然数列除以 T=32( 表 6) 。如果采用的频率变化范围
9、01,则绘制的频谱图是对称的(图 6) 。实际上,另一半是多余的,Mathcad2000 自动生成的频谱图就没有考虑另外一半儿(图 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.50510811091.510920001.495 1094.69
10、7 106Powerj0.50 freqj图 7 Mathcad2000生成的频谱图下图是常用的频谱图形式,如果存在周期,则在尖峰突出的最大点可以找到。这个图中是没有显示任何周期的,但并不意味着没有重要信息。在理论上,如果人口密度分布服从负指数模型,则其频率与功率谱之间应该满足如下关系 2)(fP8为了检验这种推断,不妨用下式进行拟合 fP)(这正是 噪声( -noise)表达式。020000000040000000060000000080000000010000000001200000000140000000016000000000 0.1 0.2 0.3 0.4 0.5频 率功率谱密度图
11、8 利用 Excel绘制的频谱图(常用形式)为了拟合幂指数模型,去掉 0 频率点,结果得到, R2=0.94947983.1254)(ffP多种模型比较的结果,发现幂指数模型的拟合效果最好(图 9) 。将图 9 转换成对数刻度,拟合效果就尤其明确(图 10) 。显然,=1.79832。P(f) = 1280514.1795 f-1.7983R2 = 0.949401000000002000000003000000004000000005000000006000000007000000008000000000 0.1 0.2 0.3 0.4 0.5频 率功率谱密度图 9 频谱图的模型拟合结果(去
12、掉 0频点)9P(f) = 1280514.1795 f-1.7983R2 = 0.94941.E+061.E+071.E+081.E+090.01 0.1 1频 率功率谱密度图 10 双对数频谱图利用模型及其参数,我们可以对杭州市人口分布特征及其变化进行系统分析。但是,深入的分析仅仅借助一个参数是不够的。具体的分析过程将用专门的文章进行论述。最后说明一点:前面的公式 2)(1)(FTP给出的是功率谱。有时在理论上进行讨论时,采用下式,2)(S这里给出的是能量谱。能量谱的计算假定数据序列无穷长,积分范围一般从负无穷到正无穷;功率谱主要用于对实际遇到的有限长度的数据。二者在数值上相差常数倍。因此,在理论讨论时,采用能量谱公式比较方便;在具体应用时采用功率谱公式便于比较。二者的数理本质是一致的,故一般行文过程中无需澄清二者的关系。