付立叶变换理解及excel实现

傅里叶变换的本质傅里叶变换的公式为F(o)=jf(t)e-网出-s可以把傅里叶变换也成另外一种形式：F(o)=—Yf(t),ejot>2兀可以看出，傅里叶变换的本质是内积，三角函数是完备的正交函数集，不同频率的三角函数的之间的内积为0，只有频率相等的三角函数做内积时，才不为0。YejQ/,e02A=jejq-q2)tdt=2兀8(。—。)下面从公式解释下傅里叶变换的意义因为傅里叶变换的本质是内积，所以f(t)和ejo求内积的时候，只有f(t)中频率为O的分量才会有内积的结果，其余分量的内积为0。可以理解为f(t)在ejot上的投影，积分值是时间从负无穷到正无穷的积分，就是把信号每个时间在o的分量叠加起来，可以理解为f(t)在ejot上的投影的叠加，叠加的结果就是频率为o的分量，也就形成了频谱。傅里叶逆变换的公式为1+sf(t)= jF(o)ejodo2兀—s下面从公式分析下傅里叶逆变换的意义傅里叶逆变换就是傅里叶变换的逆过程，在F(o)和ej求内积的时候，F(o)只有t时刻的分量内积才会有结果，其余时间分量内积结果为0，同样积分值是频率从负无穷到正无穷的积分，就是把信号在每个频率在t时刻上的分量叠加起来，叠加的结果就是f(t)在t时刻的值，这就回到了我们观察信号最初的时域。离散付立叶变换的理解FFT是离散傅立叶变换的快速算法，可以将一个信号变换到频域。有些信号在时域上是很难看出什么特征的，但是如果变换到频域之后，就很容易看出特征了。这就是很多信号分析采^FT变换的原因。另外，FFT可以将一个信号的频谱提取出来，这在频谱分析方面也是经常用的。虽然很多人都知道FFT是什么，可以用来做什么，怎么去做，但是却不知^FT之后的结果是什意思、如何决定要使用多少点来做FFT。现在就根据实际经验来说说FFT结果的具体物理意义。一个模拟信号，经过ADC采样之后，就变成了数字信号。采样定理告诉我们，采样频率要大于信号频率的两倍，这些我就不在此啰嗦了。采样得到的数字信号，就可以做FFT变换了。N个采样点，经过FFT之后，就可以得到N个点的FFT结果。为了方便进行FFT运算，通常N取2的整数次方。假设采样频率为Fs，信号频率F，采样点数为N。那么FFT之后结果就是一个为N点的复数。每一个点就对应着一个频率点。这个点的模值，就是该频率值下的幅度特性。具体跟原始信号的幅度有什么关系呢？假设原始信号的峰值为A，那么FFT的结果的每个点(除了第一个点直流分量之外)的模值就是A的N/2倍。而第一个点就是直流分量，它的模值就是直流分量的N倍。而每个点的相位呢，就是在该频率下的信号的相位。第一个点表示直流分量(即0日2)，而最后一个点N的再下一个点(实际上这个点是不存在的，这里是假设的第N+1个点，也可以看做是将第一个点分做两半分，另一半移到最后)则表示采样频率Fs，这中间被N-1个点平均分成N等份，每个点的频率依次增加。例如某点n所表示的频率为：Fn=(n-1)*Fs/N。由上面的公式可以看出，Fn所能分辨到频率为为Fs/N，如果采样频率Fs为1024Hz，采样点数为1024点，则可以分辨到1Hz。1024Hz的采样率采样1024点，刚好是1秒，也就是说，采样1秒时间的信号并做FFT，则结果可以分析到1Hz，如果采样2秒时间的信号并做FFT，则结果可以分析到0.5Hz。如果要提高频率分辨力，则必须增加采样点数，也即采样时间。频率分辨率和采样时间是倒数关系。这一部分的描述很不清晰这部分的分析很关键，有利于理解excel假设FFT之后某点n用复数a+bi表示，那么这个复数的模就是An=根号a*a+b*b，相位就是Pn=atan2(b,a)。根据以上的结果，就可以计算出n点(n手1且n<=N/2)对应的信号的表达式为：An/(N/2)*cos(2*pi*Fn*t+Pn)，即2*An/N*cos(2*pi*Fn*t+Pn)。对于n=1点的信号，是直流分量，幅度即为A1/N。由于FFT结果的对称性，通常我们只使用前半部分的结果，即小于采样频率一半的结果。下面以一个实际的信号来做说明。假设我们有一个信号，它含有2V的直流分量，频率为50Hz、相位为-30度、幅度为3V的交流信号，以及一个频率为75Hz、相位为90度、幅度为1.5V的交流信号。用数学表达式就是如下：S=2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180)。式中cos参数为弧度，所以-30度和90度要分别换算成弧度。我们以256Hz的采样率对这个信号进行采样，总共采样256点。按照我们上面的分析，Fn=(n-1)*Fs/N，我们可以知道，每两个点之间的间距就是1Hz，第n个点的频率就是n-1。我们的信号有3个频率：0Hz、50Hz、75Hz，应该分别在第1个点、第51个点、第76个点上出现峰值，其它各点应该接近0。实际情况如何呢？我们来看看FFT的结果的模值如图所示。从图中我们可以看到，在第1点、第51点、和第76点附近有比较大的值。我们分别将这三个点附近的数据拿上来细看：点：512+0i点：-2.6195E-14-1.4162E-13i点：-2.8586E-14-1.1898E-13i点：-6.2076E-13-2.1713E-12i点:332.55-192i点：-1.6707E-12-1.5241E-12i点:-2.2199E-13-1.0076E-12i点：3.4315E-12+192i点:-3.0263E-14+7.5609E-13i很明显，1点、51点、76点的值都比较大，它附近的点值都很小，可以认为是0,即在那些频率点上的信号幅度为0。接着，我们来计算各点的幅度值。分别计算这三个点的模值，结果如下：1点：51251点:38476点：192按照公式，可以计算出直流分量为：512/N=512/256=2；50Hz信号的幅度为：384/(N/2)=384/(256/2)=3；75Hz信号的幅度为192/(N/2)=192/(256/2)=1.5。可见，从频谱分析出来的幅度是正确的。然后再来计算相位信息。直流信号没有相位可言，不用管它。先计算50Hz信号的相位，atan2(-192,332.55)=-0.5236,结果是孤度，换算为角度就是180*(-0.5236)/pi=-30.0001。再计算75Hz信号的相位，atan2(192,3.4315E-12)=1.5708孤度，换算成角度就是180*1.5708/pi=90.0002。可见，相位也是对的。根据FFT结果以及上面的分析计算，我们就可以写出信号的表达式了，它就是我们开始提供的信号。利用Excel进行FFT和Fourier分析的基本步骤杭州市2000人口分布密度［根据2000年人口普查的街道数据经环带(rings)平均计算得到的结果，数据由冯健博士处理］。下面的变换实质是一种空间自相关的分析过程。第一步，录入数据在Excel中录入数据不赘述(见表1)。表1原始数据序列 表2补充后的数据序列

HAHAB□1距离人」密度20.3281S430.92682141.52462152.12317662.71891013.31960183.91694594.51032^105.17282115.T6200126.3S&44136.94297147.53306158.13153168.T2683179.32354189.920281910.51S282011.11S512111.715812212.314902312.914您2413.512782514.110332614.79582715.3882183.91694F94.51082^105.17282115.T6200126.35&44136.94297147.53806158.13153168.72683179.323541S9.920281910.513282011.11&512111.715812212.314902312.914652413.512782514.110332614.795S2715.38822S029030031032033034第二步，补充数据由于Fourier变换（FT）一般是借助快速Fourier变换（FastFourierTransformation,FFT）算法，而这种算法的技术过程涉及到对称处理，故数据序列的长度必须是2n（N=1,2,3,・・・，）。如果数据序列长度不是2n,就必须对数据进行补充或者裁减。现在数据长度是26,介于24=16到25=32之间，而26到32更近一些，如果裁减数据，就会损失许多信息。因此，采用补充数据的方式。补充的方法非常简单，在数据序列后面加0，直到序列长度为32=25为止（表2）。当然，延续到64=26也可以，总之必须是2的整数倍。不过，补充的“虚拟数据”越多，变换结果的误差也就越大。第三步，Fourier变换的选项设置沿着工具（Tools）一数据分析（DataAnalysis）的路径打开数据分析复选框（图1）。

在数据分析选项框中选择傅立叶分析(在数据分析选项框中选择傅立叶分析(FourierAnalysis)(图2)。陞］文件(已编撤日视图(心插入①格式。工具⑴数据E窗口(如帮助皿□芽19药•B争龟窿9*拳■排写⑤,.， F7二自动保存心.，B31v S0ABC共享工作懑(E，,,保护(E) ，联机协作(11) 卜F1距离」人口密度20.32818430.926321规划求解毋•,41.52462152.123176加载宏①,,.自定义©心选项(皿,，62.71391073.319601E3.916945数据尝折回,,g4.510829105.17282图1数据分析(DataAnalysis)的路径图2数据分析(DataAnalysis)在Fourier分析对话框中进行如下设置：在 输入区域中输入数据序列的单元格范围“$B$1:$B$33”；选中“标志位于第一行(L)”；将输出区域设为“$C$2”或者“$C$2:$C$33”图3a)。a

b图3傅立叶分析(FourierAnalysis)注意：如果“输入区域”设为“$B$2:$B$33”，则不选“标志位于第一行(L)(图3b)。表3FFT的结果

A■_Bc1距离大口密度FFT结果20.328184218701.857013486 30.926821104459.634221S27-103400.5376153S6141.52462137684.6丑4056542-39512.0012034132152.12317614197.250796317-55505.05387560S8162.71E91011894.2167827967-29641.7572S0&S23173.31960118938.3527863987-20821.94263532651E3.91694520635.720934404-23004.8596457225194.51082914487.9855812052-23932.2738&S&426I105.172829835.411667338-16677.82583265551115.7S20011919.46805S0233-7585.65666S7729S1126.3564416112.0959220229-E044.183634845861136.9429717113.44821S299-7100.569507992161147.5380616598.4374913769-7622.17758&S9195I158.1315315868.1982894284-7046.263462401091IS8.726S313454.062810094-5175.724568&12521179.3235412259.0748348707-1721.407440893981189.9202812261.30462533861910.5182812259.074834870S+1721.4074403940312011.1165113454.062810094+5175.7245686125212111.715S115868.1982894284+7046.26346240109122IS.3149016598.4374913769+7622.1775868919412312.9146517113.448218299+7100.5695079921412413.5127816112.0959220223+5044.1836348458412514.1103311919.4680580232+7585.6566687729912614.79589S35.41166733S02+16677.825832&555i2715.388214487.9855812052+23932.2738636426123020635.7209344041+23004.859645T225129018938.3527863987+20S21.9426353265130011894.216782796S+29641.7572S0SS23i31014197.2507963171+55505.0538T560S91MM*|R原馆数—FF；/业韭"/第四步，输出FFT结果选项设置完毕以后，确定(OK),立即得到FFT结果(表3)。显然，表3给出的都是复数(complexnumbers)o假定一个数据序列表为f(t)，则理论上Fourier变换的结果为F(w)=Mfq)e-j^tdt=Ff(t)],(一8<W<8)-3表3中给出的正是相应于F(/)的复数，这里①为角频率。第五步，计算功率谱Excel好像不能自动计算功率谱，这需要我们利用有关函数进行计算。计算公式为P(w)=T\f(w)|2=*(A2+B2)式中A为复数的实部(realnumber),B为虚部(imaginarynumber),T为假设的周期长度，实则补充后的数据序列长度。对于本例，T=32。注意复数的平方乃是一个复数与其共轭(conjugate)复数的乘积，若F(g+bj，则IF(刃)l2=(a+bj)*(a-bj)=a2+b2。这样，根据表3中的FFT结果，我们有(218701.8572+02)/32=1494703196(104459.6342+103400.5382)/32=675108949其余依次类推。显然，这样计算非常繁琐。一个简单的办法是调用Excel的模数(modulus)计算函数ImAbs，方法是在函数类别中找“其他”，在其他类中找“工程”类，在工程类中容易找到ImAbs函数(图4)。确定以后，弹出一个选项框，选中第一个FFT结果，确定，得到218701.857(图5)。我们知道，复数的模数计算公式为M=(A2+B2)1/2图4模数计算函数对于第一个FFT结果，由于虚部为0，模数就是其自身，即(218701.8572+02)1/2=218701.857但对于后面真正的复数，就不一样了。抓住第一个模数所在的单元格的右下角往下一拉，或者用鼠标双击该单元格的右下角，立即得到全部模数。ABCD距离0.3人口密度2S1S4FFT结果218701.857013486 :功率谱密度=ImAbs（C2）IniAl'EiRumberKJ="218701.85701348(=218701.B57FieiwnstheEtljEolu+ev：alne(jnodulus.lofaconiplextluhiIier氐sber回 计算结果=218T01.857 |确定|取稍|图5计算模数最后，用模数的2次方除以数据长度32立即得到全部功率谱密度结果（表4）。表4功率谱密度D E FFFT结果 模数 功率请密度阪血M的结果218701.857013486218701.85714947031961495000000104459.634221S27-103400.53761533611469S1.244967S10894967510000037634.6714056542-39512.0012034132197121.227432铜76665129430000014197.250796317-55505.053S75608Si57291.9971410257415410260000011894.2167827967-29641.7572806823131939.1009231878317.73188000013933.3527363937-20321.9426353265i2S146.3052924756703.22476000020635.720934404-23004.8596457225i30903.9891629S45517.12985000014487.9855812052-23932.2738686426i27975.9800324457983.1244600009835.411667338-16677.8258326555i19361.9523111715162.41172000011919.4680580233-7585.65666877298i14128.549326237997.06623800016112.0959220229-5044.lS3634845S6i16SS3.229068907606.99890800017113.448218299-7100.56950799216118528.0381410727756.2107300001&598.4374913769-7622.17758&8919511S264.S766310425178.71042000015868.1982894284-7046.26346240109117362.302439420298.3942000013454.062310094-5175.72456361252i14415.26729649374T.84649400012259.0748348707-1721.40744089398i12379.344064789004.98479000012261.304625338612261.3046346^3112.224697000下表是利用Mathcad2000计算的功率谱密度（表5）。利用Mathcad进行FFT，过程要简单得多，只要调用FFT命令，可以直接给出各种结果（包括图表）。但Mathcad的计算不求精度，有一定误差。将Mathcad的变换结果copy到Excel中进行比较，可以看到，如果不计误差，二者是一致的（表4）。表5借助Mathcad2000进行FFT的结果

0101.495?0916.751?0822.948?0831.026?0843.188?0752.476?0762.985?0772.446?0781.172?0796.238?06108.908?06111.073?07121.042?07139.42?06146.494?06154.79?06164.697?06Power第六步，功率谱分析功率谱分析目前主要用于两个方面，一是侦测系统变化的某种周期或者节律，据此寻找因果关系（解释）或者进行某种发展预测（应用）；二是寻找周期以外的某些规律，据此对系统的时空结构特征进行解释。表6以对称点（f=0.5）为界，从完整的数据序列中截取一半1自然数列频率助率谱密度百然数列频率助率借辑技2001494703196001494703196310.0312567510894910.03125675108949420.062529476665120.0625294766651530.0937510257415430.09375102574154640.12531878317.740.12531878317.7750.1562S24756703.250.1562524756703.2860.187529845517.160.187529845517.1970.2187524457983.170.21875244579S3.110S0.2511715162.4S0.2511715162.41190.2S1256237997.0690.281256237997.0612100.31258907606.99100.31253907606.9913110.3437510727756.2110.3437510727756.214120.37510425178.7120.37510425178.715130.406259420298.3130.406259420293.316140.43756493747.84140.43756493747.8417150.4&S754789004.98150.468754739004.9318160.54698112.22160.54698112.2219170.531254789004.98上面基于杭州人口密度数据的FFT，实际上是一种空间自相关分析过程，属于FT的第二类

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