CN105865616B - 基于fft的调制谱快速细化方法 - Google Patents

Abstract

本发明公开一种基于FFT的调制谱快速细化方法，该方法包括如下步骤：（1）由分析频率范围初步确定采样频率；（2）设定观察窗口宽度以得到预期的采样数据个数；（3）数据分组以得到实际采样数据个数；（4）调整采样频率使观察窗口与分析频率范围对应；（5）采样数据并分组做FFT；（6）计算起始谱线，求细化谱。具有如下有益效果及优点：①本发明采用基于直接FFT与任意基FFT相结合的算法，对调制谱细微结构的分辨率与直接FFT一致、比ZFFT更高；②本发明方法的运算速度比直接FFT、ZFFT快2倍以上，针对不同应用要求通过合理选取观察窗口宽度2^M、半带宽扩展系数α或采样系数β可提高到直接FFT的4倍以上。

Description

基于FFT的调制谱快速细化方法

技术领域

本发明属于细化谱分析技术，具体地说是一种基于FFT的调制谱快速细化方法。

背景技术

调制谱是模拟信号调制产生的连续已调制波的频谱。信号调制就是用一个信号（调制信号）去控制另一作为载体的信号（载波信号），让后者的某一参数（幅值、频率、相位、脉冲宽度等）按前者的值变化；通常以一个高频正弦信号作为载波信号，常用的调制方法有调幅、调频、调相、脉冲调宽。在实际工程信号监测中，由于存在干扰、固有频率等调制因素，信号形式往往如附图2所示。对于附图2中信号的采集，满足Nyquist采样定理的采样角频率f _s≥2(f _c+αf _z)，而采样周期T≥1/f _z；因此，FFT谱分析所需的采样数据点数N=2ⁿ≥2(f _c+αf _z)/f _z。以旋转机械振动监测为例，载波信号通常是振动监测传感器的固有频率信号，其频率f _c为几千或几十千赫兹；调制信号通常是旋转机械零部件的故障振动冲击信号，其频率f _z为在零点几到几十赫兹之间；这会导致采样数据过多，使FFT的时间加长，无法满足实时性要求；且造成了大量时间和资源的浪费，因所关心的频谱仅在f _c-f _z到f _c+f _z之间。

从FFT谱分析方法中已知，频率分辨率Δf为采样频率f _s与采样点数N 的比值。在对具有密集型频谱的信号进行分析时，我们希望频率分辨率Δf越小越好，同时，采样频率f _s必须保持较大值以包含所有感兴趣的频率成分。要两者兼顾，必须增大采样点数。这就意味着运算速度和分析效率将大大地降低。为了解决这一矛盾，上世纪70年代后相继出现了基于不同原理的各种频率细化分析方法，其主要目的是识别谱图上的细微结构。例如：扫频窄带分析法、基于复调制的Zoom-FFT法、直接选抽法、相位补偿细化和最大频谱的局部表示法等。最常用的有复调制Zoom-FFT、相位补偿细化等方法。然而在计算效率、精度和灵活性等方面都比较理想的方法还是基于复调制的Zoom-FFT，因此得到了较多的应用。以基于复调制的ZFFT方法为例，很多人对ZFFT有个误区，事实上，ZFFT所谓的提高分辨率是针对作同样点数的离散傅立叶变换而言的，即作细化D倍的Ñ点细化谱实际上原始数据必须采N=D*Ñ点数据，这时候它比用Ñ点原始数据直接做FFT而言频率分辨率提高了D倍，但如果把N=D*Ñ点原始数据全部做FFT，那它和细化谱的分辨率是一样的。也就是说，ZFFT本质上可以说是一种快速算法，并不能提高分辨率。它通过滤波重采样来降低采样频率，这样就可以用较少点数的FFT来实现较高的频率分辨率。当然，提高速度的代价就是只能对局部频带进行细化；而如果将ZFFT利用的所有原始数据全部直接做FFT的话，它做出的是整个频域的，而且频率分辨率和细化后的一样，甚至如果考虑细化时滤波所需去掉的点，直接FFT的频率分辨率可以更高。

实际工程信号监测中，为了快速识别调制谱中的细微结构，就必须要求信号分析既要具有高的频率分辨率、又要有较宽的频率范围，且满足系统实时性；然而，常规的FFT谱、频率细化方法都难以同时兼顾。能否发明一种快速识别调制谱中细微结构的方法，目前国内外科技界一直没有解决，尚无相关的研究成果报道。

发明内容

针对目前调制谱分析中存在的不足，本发明为解决该问题提供了一种基于FFT的调制谱快速细化方法。

本发明采取如下技术方案：

基于FFT的调制谱快速细化方法包括如下步骤：

步骤一，由分析频率范围初步确定采样频率；

步骤二，设定观察窗口宽度以得到预期的采样数据个数；

步骤三，数据分组以得到实际采样数据个数；

步骤四，调整采样频率使观察窗口与分析频率范围对应；

步骤五，采样数据并分组做FFT；

步骤六，计算起始谱线，求细化谱。

在本发明中，对于调制频率f _z、载波频率f _c的信号，取半带宽扩展系数为α，则所需的分析频率范围是：f _c-αf _z到f _c+αf _z；初步确定满足采样定理的采样频率fi＝β（f _c+αf _z），其中采样系数β≥2。

在本发明中，根据实际应用需要设定观察窗口宽度为2^M条谱线，可知满足要求的采样数据个数Ñ=fi/Δfi，其中观察窗口预定频谱分辨率Δfi＝αf _z /2^M-1。

在本发明中，为能观察到2^M条有效谱线，据FFT对称性取Q＝2^M＋1，按任意基FFT原则将数据分为P=INT（Ñ/Q+0.5）组且每组Q点，则实际采样数据个数N=P×Q。

在本发明中，为使2^M条有效谱线能表示所需的分析频率范围：f _c-αf _z到f _c+αf _z，需调整采样频率f _s＝f _c*N/λ使载波频率f _c的谱线位于观察窗口中心；若f _s<2(f _c+αf _z)，则返回第一步重新确定采样系数β值。观察窗口中心谱线号λ与预定频谱分辨率Δfi、数据分组大小Q的关系为λ＝(INT( f _c /(Δfi*Q)+0.5))*Q+2^M-1。

在本发明中，按采样频率f _s采样N个数据，并在采样同时分为P组Q点的形式存放，采样结束后对已采集P组Q点数据分别做直接FFT。

在本发明中，由采样频率f _s和数据个数N计算观察窗起始谱线L₀＝INT(N*(f _c-αf _z)/ f _s+0.5)，从L₀处开始按下式计算2^M条谱线即为所需的细化谱

式中k＝L₀, L₀+1,……,L₀+2^M -1；l＝0, 1,……,P-1；G _l (k)= G _l (Z*Q+y)= G _l (y)，G _l (y)为P组FFT结果中的一个值，Z为正整数，0≤y≤Q-1。

本发明具有如下有益效果及优点：（1）本发明采用基于直接FFT与任意基FFT相结合的算法，对调制谱细微结构的分辨率与直接FFT一致、比ZFFT更高；（2）本发明方法的运算速度比直接FFT、ZFFT快2倍以上，针对不同应用要求通过合理选取观察窗口宽度2^M、半带宽扩展系数α或采样系数β可提高到直接FFT的4倍以上。

附图说明

图1是本发明的运算模块流程图；

图2是调制谱形成示意图；

图3是本发明的实施例细化谱图。

具体实施方式

下面结合附图和实施例对本发明的方案作进一步详细说明。

结合附图1、附图2、附图3，基于FFT的调制谱快速细化方法包括如下步骤：（1）由分析频率范围初步确定采样频率；（2）设定观察窗口宽度以得到预期的采样数据个数；（3）数据分组以得到实际采样数据个数；（4）调整采样频率使观察窗口与分析频率范围对应；（5）采样数据并分组做FFT；（6）计算起始谱线，求细化谱。

步骤一，由分析频率范围初步确定采样频率

对于调制频率f _z、载波频率f _c的信号，取半带宽扩展系数为α，则所需的分析频率范围是：f _c-αf _z到f _c+αf _z；初步确定满足采样定理的采样频率fi＝β（f _c+αf _z），其中采样系数β≥2。以某旋转机械的振动监测为例，振动监测传感器的固有频率为12KHz，而监测点的故障特征频率不大于11Hz，则设载波频率f _c＝12KHz、调制频率f _z＝11Hz。振动监测所关心的频谱仅在f _c-f _z到f _c+f _z之间，考虑到转速波动等因素的影响，取半带宽扩展系数为α＝1.1，则所需的分析频率范围是：11987.9Hz~12012.1Hz。对该监测点信号的采集，根据奈魁斯特（Nyquist）采样定理，采样系数β≥2，取β＝2.1；这样，可初步确定该监测点满足采样定理的采样频率fi＝β（f _c+αf _z）＝2.1×12012.1Hz＝25225.41 Hz。

步骤二，设定观察窗口宽度以得到预期的采样数据个数

根据实际应用需要设定观察窗口宽度为2^M条谱线，可知满足要求的采样数据个数Ñ=fi/Δfi，其中观察窗口预定频谱分辨率Δfi＝αf _z /2^M-1。对于步骤一中所述的旋转机械振动监测例，取M=8，即采用宽度为2⁸＝256条谱线的观察窗口；由此得观察窗口预定频谱分辨率Δfi＝2αf _z /2^M＝αf _z /2^M-1＝12.1 /128＝0.094531Hz，可知满足此频谱分辨率要求的采样数据个数Ñ=fi/Δfi=2^M-1 fi/αf _z＝25225.41/0.094531＝266847.31个。

步骤三，数据分组以得到实际采样数据个数

为能观察到2^M条有效谱线，据FFT对称性取Q＝2^M＋1，按任意基FFT原则将数据分为P=INT（Ñ/Q+0.5）组且每组Q点，则实际采样数据个数N=P×Q。根据FFT的对称性，FFT谱中的前半部与后半部对称一致。对于步骤一中所述的旋转机械振动监测例，为能观察到2^M＝256条有效谱线，所需的数据个数应为Q＝2^M＋1＝512；这样，以每组数据为Q点，按任意基FFT原则采用四舍五入取整方式将数据Ñ分为P组，P=INT（Ñ/Q+0.5）=INT（266847.31/512+0.5）＝521，则得到实际所需的采样数据个数N=P×Q＝521×512＝266752个。显然，本实施例的采样数据个数N不是2的整数幂；若要进行频谱变换，有以下两种处理方法：① 通过对采样数据序列x(n)补0，使x(n)的点数成为2的整数幂以便直接FFT；② 采用任意基FFT算法，其基本思路仍然是将长点的DFT尽可能分解为短点的DFT运算。为介绍方便，以15点的数据序列x(n)说明任意基FFT的数据分组原则、运算特点，设N=P×Q＝3×5＝15点数据可分为3组、每组5点

一般项为x(Pr+l)；r＝0, 1,……,Q-1；l＝0, 1,……,P-1；其DFT为

也就是说，可将分解P＝3组、Q＝5点的DFT，再利用的周期性将3个5点的DFT合成一个15点的DFT。

步骤四，调整采样频率使观察窗口与分析频率范围对应

为使2^M条有效谱线能表示所需的分析频率范围：f _c-αf _z到f _c+αf _z，需调整采样频率f _s＝f _c*N/λ使载波频率f _c的谱线位于观察窗口中心；若f _s<2(f _c+αf _z)，则返回第一步重新确定采样系数β值；观察窗口中心谱线号λ与预定频谱分辨率Δfi、数据分组大小Q的关系为λ＝(INT( f _c /(Δfi*Q)+0.5))*Q+2^M-1。对于步骤一中的旋转机械振动监测例，采样频率f _s＝f _c*N/λ＝f _c*N/[(INT( f _c/(Δfi*Q)+0.5))*Q+2^M-1] ＝12000*266752/[(INT(12000/512*0.094531)+0.5)*512+128]＝25184.29Hz，因f _s＝25184.29Hz >2(f _c+αf _z)＝24024.2Hz，则本例不需重新确定采样系数β值，频谱分辨率Δf＝f _s/N＝25184.29/266752=0.094411 Hz。

步骤五，采样数据并分组做FFT

按采样频率f _s采样N个数据，并在采样同时分为P组Q点的形式存放，采样结束后对已采集P组Q点数据分别做直接FFT。本发明采用Q为2的整数幂的任意基FFT算法，而对Q点数据则直接FFT运算，充分利用了任意基FFT和直接FFT运算的优点；对于采样数据序列x(n)的分组可以在采样数据保存时或做直接FFT运算取数时进行，这运用CPU的间接寻址、变址寻址方式很容易快速实现。

步骤六，计算起始谱线，求细化谱

由采样频率f _s和数据个数N计算观察窗起始谱线L₀＝INT(N*(f _c-αf _z)/ f _s +0.5)，从L₀处开始按下式计算2^M条谱线即为所需的细化谱

式中k＝L₀, L₀+1,……,L₀+2^M -1；l＝0, 1,……,P-1；由DFT的周期性可知G _l (k)=G _l (Z*Q+y)= G _l (y)，G _l (y)为P组FFT结果中的一个复数值，Z为正整数，0≤y≤Q-1。

对于采样数据序列x(n) 的DFT计算量为N²次复数乘和N(N-1)次复数加，而FFT计算量为0.5N*(LOG(N,2))次复数乘和N*(LOG(N,2))次复数加。以N＝1024为例，DFT运算需要1024²＝1048576次复乘，而按时间抽取的FFT算法仅需512*10=5120次复乘，两者相差200多倍。

本发明是P个Q点直接FFT与P组Q点任意基FFT的有机结合，对调制谱细微结构的分辨率与直接FFT完全一致；关于计算量的考查，以其中的复数乘计算为例，对于步骤一中的旋转机械振动监测例，其复数乘计算量为P*0.5Q*(LOG(Q,2))+ 0.5Q *(P-1)＝521*256*(LOG(512,2))+ 256 *(521-1)＝1402880次，由上式可知：通过减少Q可降低前半部直接FFT的运算量、但后半部任意基FFT的运算量增加，通过增加Q可降低后半部任意基FFT的运算量、但前半部直接FFT的运算量增加，根据DFT的周期性，Q的大小与分析窗口宽度2^M并无直接关系，本发明取Q＝2^M+1是运算量最优的情况；若对同样数据个数N=P×Q＝521×512＝266752直接FFT，因2¹⁸<N<2¹⁹，则需补0到2¹⁹个数，复数乘为2¹⁸*(LOG(2¹⁹,2))＝4980736次；在该实施例中，直接FFT的复数乘是本发明方法的3.7倍。附图3是步骤一中的旋转机械振动监测例，图中上半部分为测点监测的时域信号，从中难以判断被监测机械的运行状态；下半部分为细化谱图，分析频率范围是：11987.9154Hz~12012.0846Hz，图中可明显观察到3.3Hz的机械运行状态特征谱线，若仅为重复该特征谱线，可观察窗口缩小为原来的1/3、其它不变，则直接FFT的复数乘将是本发明方法的4.0倍。因此，针对不同应用要求，通过合理选取观察窗口宽度2^M、半带宽扩展系数α或采样系数β可使本发明方法的运算速度比直接FFT快4倍以上。

对于步骤一中的旋转机械振动监测例，若采用ZFFT来实现频谱细化。同样取采样频率f _s＝25184.29Hz，设频谱窗口宽度2^M＝256，需要Q＝512点数据来自采样原始数据，则其FFT的频谱分辨率Δf _a＝f _s/Q＝25184.29/512=49.1880664 Hz，要使谱分辨率提高到Δf＝0.094411 Hz，细化倍数D=Δf _a/Δf≈521，所需的原始采样数为N=D×Q＝521×512＝266752点数据，可见采样数据与本发明方法一致。计算量以基于复解析滤波的ZFFT方法为例，D=521为细化倍数，令M= 4D=2084为滤波器的半阶数，Q＝512为FFT分析数据点数，则细化需要的采样序列为x(n) ，其中n ＝0,1,2,……,DQ+2M-1。从采样序列x(n)中选抽g(r)的方法是：①频移，为了将感兴趣的频段的下限频率移至原来的零频率位置，以便有可能将感兴趣频段放大到整个频率显示范围上，需首先对信号进行频率调制。一般采用的是复数调制法，如果欲将某一频率f _x＝f _c-f _z移至原来的零频处，则以原信号x _i 与e- ^{j2π×fx×k/N×Δf} 相调制得：实部为x _i ×cos(2π×f _x×k/N×Δf)，虚部为-x _i ×sin(2π×f _x×k/N×Δf)；②选抽，确定选抽位置为[r＝2M,2M+D,2M+2D,……,2M+(Q-1)D]，取每个选抽位置以前的2M个采样数据；③波器，用设计好的复解析滤波器对2M个采样数据滤波后得到g(r)中一个数据的实部、虚部；④ 对②、③步重复Q次，即可选抽得到g(r)；⑤ 对g(r)作FFT。可见，频移需作N=D×Q次复乘，滤波需作2M×Q＝8DQ次复乘，做Q点的直接FFT需0.5Q*(LOG(Q,2))次复乘，共需复乘次数9DQ＋0.5Q*(LOG(Q,2))；对于步骤一中的旋转机械振动监测例，D=521，Q＝512，其复乘次数为2403072次，是本发明方法复乘次数的1.713倍；考虑复数加、移频与滤波系数计算等，本发明方法整体可比ZFFT快2倍以上；且ZFFT要通过复调制移频、经复解析滤波减少参与FFT的数据点，不可避免的会损失一些频率成分，所得到频谱的分辨率没有直接FFT的高。

Claims (2)

1.一种基于FFT的调制谱快速细化方法，包括由分析频率范围初步确定采样频率的步骤一，设定观察窗口宽度以得到预期采样数据个数的步骤二，数据分组以得到实际采样数据个数的步骤三，调整采样频率使观察窗口与分析频率范围对应的步骤四，采样数据并分组做FFT的步骤五，计算起始谱线并求细化谱的步骤六；其特征在于：

步骤一，对于调制频率f _z、载波频率f _c的信号，取半带宽扩展系数为α，则所需的分析频率范围是：f _c-αf _z到f _c+αf _z，初步确定满足采样定理的采样频率fi=β(f _c+αf _z)，其中采样系数β≥2；

步骤二，根据实际应用需要设定观察窗口宽度为2^M条谱线，可知预期采样数据Ñ与采样频率fi、观察窗口预定频谱分辨率Δfi的关系为Ñ=fi/Δfi，其中Δfi＝αf _z/2^M-1，M为使2^M满足观察谱线数要求的正整数；

步骤三，为能观察到2^M条有效谱线，据FFT对称性取Q＝2^M＋1，按任意基FFT原则将数据分为P=INT（Ñ/Q+0.5）组且每组Q点，则实际采样数据个数N=P×Q；

步骤四，为使2^M条有效谱线能表示所需的分析频率范围：f _c-αf _z到f _c+αf _z，需调整采样频率f _s＝f _c*N/λ使载波频率f _c的谱线位于观察窗口中心，若f _s<2(f _c+αf _z)，则返回第一步重新确定采样系数β值，其中λ为观察窗口中心谱线号；

步骤五，按采样频率f _s采样N个数据，并在采样同时分为P组Q点的形式存放，采样结束后对已采集P组Q点数据分别做直接FFT；

步骤六，由采样频率f _s和数据个数N计算观察窗起始谱线L₀＝INT(N*(f _c-αf _z)/f _s+0.5)，从L₀处开始按下式计算2^M条谱线即为所需的细化谱

式中k＝L₀, L₀+1,……,L₀+2^M -1；l＝0, 1,……,P-1；G _l (k)= G _l (Z*Q+y)= G _l (y)，G _l (y)为P组FFT结果中的一个值，Z为正整数，0≤y≤Q-1。

2.根据权利要求1所述的基于FFT的调制谱快速细化方法，其特征在于：观察窗口中心谱线号λ与预定频谱分辨率Δfi、数据分组大小Q的关系为λ＝(INT(f _c/(Δfi*Q)+0.5))*Q+2^M-1。