付立叶变换理解及excel实现.docx

付立叶变换理解及excel实现.docx

《付立叶变换理解及excel实现.docx》由会员分享，可在线阅读，更多相关《付立叶变换理解及excel实现.docx（16页珍藏版）》请在冰点文库上搜索。

付立叶变换理解及excel实现

傅里叶变换的本质

傅里叶变换的公式为

可以把傅里叶变换也成另外一种形式：

可以看出，傅里叶变换的本质是内积，三角函数是完备的正交函数集，不同频率的三角函数的之间的内积为0，只有频率相等的三角函数做内积时，才不为0。

下面从公式解释下傅里叶变换的意义

因为傅里叶变换的本质是内积，所以f（t）和

求内积的时候，只有f（t）中频率为

的分量才会有内积的结果，其余分量的内积为0。

可以理解为f（t）在

上的投影，积分值是时间从负无穷到正无穷的积分，就是把信号每个时间在

的分量叠加起来，可以理解为f（t）在

上的投影的叠加，叠加的结果就是频率为

的分量，也就形成了频谱。

傅里叶逆变换的公式为

下面从公式分析下傅里叶逆变换的意义

傅里叶逆变换就是傅里叶变换的逆过程，在

和

求内积的时候，

只有t时刻的分量内积才会有结果，其余时间分量内积结果为0，同样积分值是频率从负无穷到正无穷的积分，就是把信号在每个频率在t时刻上的分量叠加起来，叠加的结果就是f（t）在t时刻的值，这就回到了我们观察信号最初的时域。

离散付立叶变换的理解

FFT是离散傅立叶变换的快速算法，可以将一个信号变换到频域。

有些信号在时域上是很难看出什么特征的，但是如果变换到频域之后，就很容易看出特征了。

这就是很多信号分析采用FFT变换的原因。

另外，FFT可以将一个信号的频谱提取出来，这在频谱分析方面也是经常用的。

虽然很多人都知道FFT是什么，可以用来做什么，怎么去做，但是却不知道FFT之后的结果是什意思、如何决定要使用多少点来做FFT。

现在就根据实际经验来说说FFT结果的具体物理意义。

一个模拟信号，经过ADC采样之后，就变成了数字信号。

采样定理告诉我们，采样频率要大于信号频率的两倍，这些我就不在此啰嗦了。

采样得到的数字信号，就可以做FFT变换了。

N个采样点，经过FFT之后，就可以得到N个点的FFT结果。

为了方便进行FFT运算，通常N取2的整数次方。

假设采样频率为Fs，信号频率F，采样点数为N。

那么FFT之后结果就是一个为N点的复数。

每一个点就对应着一个频率点。

这个点的模值，就是该频率值下的幅度特性。

具体跟原始信号的幅度有什么关系呢？

假设原始信号的峰值为A，那么FFT的结果的每个点（除了第一个点直流分量之外）的模值就是A的N/2倍。

而第一个点就是直流分量，它的模值就是直流分量的N倍。

而每个点的相位呢，就是在该频率下的信号的相位。

第一个点表示直流分量（即0Hz），而最后一个点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点附近有比较大的值。

我们分别将这三个点附近的数据拿上来细看：

1点：

512+0i 

2点：

-2.6195E-14-1.4162E-13i 

3点：

-2.8586E-14-1.1898E-13i 

50点：

-6.2076E-13-2.1713E-12i 

51点：

332.55-192i 

52点：

-1.6707E-12-1.5241E-12i 

75点：

-2.2199E-13-1.0076E-12i 

76点：

3.4315E-12+192i 

77点：

-3.0263E-14+7.5609E-13i 

很明显，1点、51点、76点的值都比较大，它附近的点值都很小，可以认为是0，即在那些频率点上的信号幅度为0。

接着，我们来计算各点的幅度值。

分别计算这三个点的模值，结果如下：

1点：

512 

51点：

384 

76点：

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补充后的数据序列

第二步，补充数据

由于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）。

图1数据分析（DataAnalysis）的路径

在数据分析选项框中选择傅立叶分析（FourierAnalysis）（图2）。

图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的结果

第四步，输出FFT结果

选项设置完毕以后，确定（OK），立即得到FFT结果（表3）。

显然，表3给出的都是复数（complexnumbers）。

假定一个数据序列表为f（t），则理论上Fourier变换的结果为

=F[f（t）],（

）

表3中给出的正是相应于F（ω）的复数，这里ω为角频率。

第五步，计算功率谱

Excel好像不能自动计算功率谱，这需要我们利用有关函数进行计算。

计算公式为

式中A为复数的实部（realnumber），B为虚部（imaginarynumber），T为假设的周期长度，实则补充后的数据序列长度。

对于本例，T=32。

注意复数的平方乃是一个复数与其共轭（conjugate）复数的乘积，若F（ω）=a+bj，则|F（ω）|2=（a+bj）*（a-bj）=a2+b2。

这样，根据表3中的FFT结果，我们有

其余依次类推。

显然，这样计算非常繁琐。

一个简单的办法是调用Excel的模数（modulus）计算函数ImAbs，方法是在函数类别中找“其他”，在其他类中找“工程”类，在工程类中容易找到ImAbs函数（图4）。

确定以后，弹出一个选项框，选中第一个FFT结果，确定，得到218701.857（图5）。

我们知道，复数的模数计算公式为

图4模数计算函数

对于第一个FFT结果，由于虚部为0，模数就是其自身，即

但对于后面真正的复数，就不一样了。

抓住第一个模数所在的单元格的右下角往下一拉，或者用鼠标双击该单元格的右下角，立即得到全部模数。

图5计算模数

最后，用模数的2次方除以数据长度32立即得到全部功率谱密度结果（表4）。

表4功率谱密度

下表是利用Mathcad2000计算的功率谱密度（表5）。

利用Mathcad进行FFT，过程要简单得多，只要调用FFT命令，可以直接给出各种结果（包括图表）。

但Mathcad的计算不求精度，有一定误差。

将Mathcad的变换结果copy到Excel中进行比较，可以看到，如果不计误差，二者是一致的（表4）。

表5借助Mathcad2000进行FFT的结果

第六步，功率谱分析

功率谱分析目前主要用于两个方面，一是侦测系统变化的某种周期或者节律，据此寻找因果关系（解释）或者进行某种发展预测（应用）；二是寻找周期以外的某些规律，据此对系统的时空结构特征进行解释。

表6以对称点（f=0.5）为界，从完整的数据序列中截取一半

上面基于杭州人口密度数据的FFT，实际上是一种空间自相关分析过程，属于FT的第二类应用。

这种过程不以寻找周期为目标，实际上也不存在任何周期。

不论目标是什么，都必须借助频谱图（频率－功率谱密度图）进行分析和解释。

下面第一步就是绘制频谱图。

首先要计算频率，线频或角频都可以，因为二者相差常数倍（2π）。

一个简单的办法是，用0到T=32的自然数列除以T=32（表6）。

如果采用的频率变化范围0~1，则绘制的频谱图是对称的（图6）。

实际上，另一半是多余的，Mathcad2000自动生成的频谱图就没有考虑另外一半儿（图7）。

因此，我们可以以对称点f=0.5为界，截取前面一半的数据，在Excel上绘制频谱图（图8）。

图6对称的频谱图（基于完整的数据序列）

图7Mathcad2000生成的频谱图

下图是常用的频谱图形式，如果存在周期，则在尖峰突出的最大点可以找到。

这个图中是没有显示任何周期的，但并不意味着没有重要信息。

在理论上，如果人口密度分布服从负指数模型，则其频率与功率谱之间应该满足如下关系

为了检验这种推断，不妨用下式进行拟合

这正是β噪声（β-noise）表达式。

图8利用Excel绘制的频谱图（常用形式）

为了拟合幂指数模型，去掉0频率点，结果得到

R2=0.9494

多种模型比较的结果，发现幂指数模型的拟合效果最好（图9）。

将图9转换成对数刻度，拟合效果就尤其明确（图10）。

显然，β=1.7983≠2。

图9频谱图的模型拟合结果（去掉0频点）

图10双对数频谱图

利用模型及其参数，我们可以对杭州市人口分布特征及其变化进行系统分析。

但是，深入的分析仅仅借助一个参数是不够的。

具体的分析过程将用专门的文章进行论述。

最后说明一点：

前面的公式

给出的是功率谱。

有时在理论上进行讨论时，采用下式

这里给出的是能量谱。

能量谱的计算假定数据序列无穷长，积分范围一般从负无穷到正无穷；功率谱主要用于对实际遇到的有限长度的数据。

二者在数值上相差常数倍。

因此，在理论讨论时，采用能量谱公式比较方便；在具体应用时采用功率谱公式便于比较。

二者的数理本质是一致的，故一般行文过程中无需澄清二者的关系。