付立叶变换理解及excel实现.docx
- 文档编号:26487953
- 上传时间:2023-06-19
- 格式:DOCX
- 页数:16
- 大小:258.72KB
付立叶变换理解及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双对数频谱图
利用模型及其参数,我们可以对杭州市人口分布特征及其变化进行系统分析。
但是,深入的分析仅仅借助一个参数是不够的。
具体的分析过程将用专门的文章进行论述。
最后说明一点:
前面的公式
给出的是功率谱。
有时在理论上进行讨论时,采用下式
这里给出的是能量谱。
能量谱的计算假定数据序列无穷长,积分范围一般从负无穷到正无穷;功率谱主要用于对实际遇到的有限长度的数据。
二者在数值上相差常数倍。
因此,在理论讨论时,采用能量谱公式比较方便;在具体应用时采用功率谱公式便于比较。
二者的数理本质是一致的,故一般行文过程中无需澄清二者的关系。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 付立叶 变换 理解 excel 实现
![提示](https://static.bdocx.com/images/bang_tan.gif)