利用EXCEL“傅理叶分析”工具进行频谱分析

发表于 讨论求助 2023-05-10 14:56:27

Source:声振之家


EXCEL分析工具库中内置了分析工具“傅立叶分析”,其功能是进行离散型快速傅立叶变换(FFT),也可进行傅立叶逆变换。


傅立叶变换是将时间序列数据转换为频率序列数据,以便了解序列的频率构成。


对于时间序列可展开为傅立叶级数:


式中:


N为观测值个数;k为周期分量个数;fj为频率(=j/N),εt为误差项,是由于选取级数前k项所产生的。


时间域序列xk变换到频率域序列ωj的公式如下:


式中:N为序列数据项数,对第j个分量:aj为实部,bj为虚部,其模为:



辐角为:


我们所关心的是序列主要由哪些频率成份构成及其振幅。操作方法如下:


Step1:取得时间序列数据{xk},要求项数N为2的整数次幂,即2、4、8、16、32…,项数最大限制为4096。N选择多大为好,要视频谱分析时要分析的项数。N个数据中,最多能分析N/2+1个频率构成。


Step2:时间序列{xk}各项减其平均数E(xk)得中心化时间序列;


Step3:利用“傅立叶分析”工具进行快速傅立叶变换,得ωj


Step4:利用IMABS()函数求得复数的模。该序列第1项为0,去掉之,从第2项起共奇数项,中间项为常数项,两侧是完全对称的。绘制折线图观察之,通常只观察前半部分;


Step5:更改横坐标,观察频率分布。


需要指出的是:数据系列的周期性,是系统的特性,并不是由采样的时间间隔和样本量的多寡所决定。

应用举例

何先生,在自己所从事的工作中,以每分钟等间隔抽样200次,抽取了168383条记录,下图中只列出前几条:

图1 部分数据资料


以CH#1为例,从中按顺序选择样本量为128的样本,编制频谱分析图如下:


Step1:先按顺序截取128个样本单位的样本(必需是2的整数幂,本例为2⁷)(见图2)

图2 128个样本单位的傅立叶分析


Step2:在C1单元格输入“=AVERAGE(B2:B129)”求得平均数,在C3单元格输入公式,求得观测值与平均值之差,并向下复制到B3:B129。


Step3:工具|数据分析|傅立叶分析,设置对话框(如图3),求得(如图2)中D列的傅立叶变换。

图3 傅立叶分析对话框设置


Step4:在E2单元格输入(如图2)所示函数,求得D2单元格复数的模,向下复制到B3:B129。将B2:B129制成折线图(如图4)。

图4 频率各构成成分的模


由图可见,图形是完成对称的,通常只看前面一半。需要指明的是该频谱图是由系统特性决定的,样本量不同,其频谱是类似的,只是图形密集程度不同和模的大小不同。模是由多个周期样本模的叠加的结果,样本量越大,模越大。但这一点并不影响分析的结果,我们只考虑频率强度从大到小的有限个频率,即考虑主要频率构成。


Step5:确定横轴分类标志:


将图形的横轴先进行编号,编号从0开始,本例选择128个样本单位,编号为0~127,然后再用编号值除以128,得到一个周期,周期的倒数即为频率。


按此方法制作了N=512、N=2048和N=4096的频谱图如下:


图5 N=128的频谱图

图6 N=512的频谱图

图7 N=2048的频谱图

图8 N=4096的频谱图


由图可见,样本量越多频率构成越丰富。但分析频谱时,都集中在峰值附近,不能反映面上的情况。由图可见,模较大的频率成分周期分布在0.13~0.26之间,也就是频率在4~8之间,我们选择N=256项进行分析完全够用。由于图形是对称的,只看前半部分,128项,分析占总数约10%的成份,即分析12个主要频率。

操作

(1) 按顺序选择256项数据,并求平均数,进行中心化平均,使均值为0;


(2) 利用“傅立叶分析”工具求得快速傅立叶变换;


(3) 选择一半的数据(从第2项到第129项)


  • 利用“=IMABS(D3)”求得复数的模。


  • 从第2项开始从1进行编号;以编号值除256得周期序列。


  • 将周期序列求倒数得频率。


  • 以频率为横坐标、模为纵坐标绘制频率分布图:

图9 选择256个样本数据,显示快速傅立叶变换的前半部分


(4) 利用“=IMREAL”函数提取实部,用“=IMAGINARY”函数提取虚部,形成序列值。


(5) 筛选主要成分。


  • 在L列输入第k大的顺序号;


  • 在M列输入“=LARGE($E$3:$E$130,L2)”提取第k大傅立叶变换的模;


  • 在N列输入:“=MATCH(M2,$E$3:$E$130,0)”提取第k大的顺序号;


  • 在O列输入:“=INDEX(I$3:I$130,$N2)”提取第k大的实部;


  • 在P列输入:“=INDEX(J$3:J$130,$N2)”提取第k大的虚部;


  • 在Q列输入:“=INDEX($H$3:$H$130,N2)”提取第k大的频率。


图10 主要频率成分提取


于是得函数主要成分的傅立叶级数中的主成分:

本文来源于新浪博客:关文中工作室,作者:重庆三峡学院关文忠。


发表
26906人 签到看排名