基于全相位滤波的密集成分的频率估计

杨俊东， 余 江， 黄 铭， 张 果， 葛孚华

(1. 云南大学 信息学院，昆明 650504； 2. 昆明理工大学 信息工程与自动化学院，昆明 650500)

频率测量是振动工程中最常见的问题之一。许多物理量的测量(如大气动力学中的风速测量[1]、阵列信号处理中的波达方向估计[2]、振动分析中的转速测量[3]等)都可转化为频率估计问题；频率估计也是频谱分析理论中最基本的问题，而其中基于FFT(Fast Fourier Transform)的经典谱分析方法，相比于各种现代谱估计方法(如AR模型法[4]、多重信号分类法[5]、旋转不变参数估计法[6]等)，具有其所不能比拟的运算效率高、耗费内存小的优势，故在仪表设计、电力谐波分析、旋转机械故障诊断等工程领域中得到非常广泛的应用，基于FFT的频率估计器一直是学术研究的热点问题。

影响FFT频率估计精度主要有两方面因素：频谱泄漏和栅栏效应。频谱泄漏可通过加窗来改善，对于单复指数成分，理想的傅里叶谱是冲激函数(即对应单根谱线)，但因存在谱泄漏，其直接FFT的可能结果是产生五根以上的大能量值谱线，加窗FFT可起到抑制泄漏的效果，表现在其主要能量通常集中在三根谱线上；而栅栏效应则可以通过频谱校正的途径来解决，具体说来，就是借助内插措施，对峰值谱附近的谱线做进一步修正和细化来提升估计精度，设计不同的内插算法就产生多种频率校正器(如Quinn校正器、能量重心校正器、Macleod校正器、Candan校正器、比值校正器、AM谱估计器、Tsui估计器等)。

文献[7]指出，频谱校正误差来源于三个方面：系统固有误差、谱间干扰误差和噪声干扰误差。需指出，经过文献[8-15]在内插算法设计方面的不断完善，其系统误差做得越来越小，精度越来越高，可以说几乎不存在提升的空间(如AM估计器的频率估计均方误差仅为克拉美-罗限的1.014 7倍，再降低几乎很困难)。然而，需指出的是，文献[8-15]的频谱校正算法都是针对单频复指数信号而推导出的，即没有考虑存在谱间干扰误差的情况，故只适用于单频信号或者间隔较远的多频信号情况，对于频率间隔较近的密集谱信号，由于存在严重的谱间干扰，其频率估计精度会急剧降低。因而如何提升密集谱校正精度是亟待解决的问题。

密集谱校正具有很高的工程意义，因为密集谱实际上是和短样本情况紧密联系在一起的：在采样速率fs不变的情况下，其FFT频率分辨单位Δf=fs/N(N为样本长度)，因而当工程上无法采集到足够样本时(如地震波测量和爆破勘探中，N只能取小数值)， 则Δf变大， 这时会导致两个频率成分的间距|f1-f2|/Δf减小，从而短样本情况趋于变为频率估计难度更高的密集谱情况。

本文提出基于全相位滤波的密集谱频率估计算法，该算法以FFT为中心的三根谱线的相位谱值的线性偏离程度作为密集谱的判据，借助全相位点通滤波器依次去除所关注的成分的带外干扰，进而结合比值频谱校正器，可以大大减小系统误差，提升频率估计精度。具有较高的实用价值。

1 密集谱频率估计问题

1.1 单频信号的直接FFT谱和加窗FFT谱特性

若以采样速率fs对频率为f1的单频复指数信号采集N个样点， 其离散序列可表示为

(1)

令FFT的频率分辨率Δω=2π/N， 相应地，x1(n)对应的数字角频率ω1可表示为

(2)

令矩形窗傅里叶幅度谱函数为RN(ω)， 不难证明，x1(n)的直接FFT谱X1(k)表达式为

(3)

令某对称窗序列w(n)的傅里叶幅度谱W(ω)， 则对x1(n)做加窗FFT得到的离散谱表达式为

(4)

令信号频率f1=5.2 Hz，幅值a1=1， 初相θ1=60°，采样速率fs=32 Hz，样本长度N=32，则Δf=1 Hz，其直接FFT和加汉宁窗FFT的振幅谱，如图1所示。

图1 直接FFT和加窗FFT的振幅谱图Fig.1 Amplitude spectra for direct FFT and windowed FFT

从图1可知，加窗FFT相比于直接FFT， 其谱泄漏范围大为减少(由于W(ω)的旁瓣衰减高于RN(ω)的缘故)。进而对于密集谱估计场合，为尽量减小谱间干扰误差，只能用加窗FFT。

1.2 密集成分的直接FFT谱和加窗FFT谱特性

在单频成分基础上，假定信号包含两个余弦成分

x1(n)=a1cos(ω1n+θ1), 0≤n≤N-1

(5)

x2(n)=a2cos(ω2n+θ2), 0≤n≤N-1

(6)

从而复合信号表示为

x(n)=x1(n)+x2(n)

(7)

令x1(n)的参数与图1情况相同， 对x2(n)做如下参数设置：f2=6.9 Hz，幅值a2=1，初相θ2=45°， 则单频信号x1(n)，x2(n)和复合信号x(n)的直接FFT和加汉宁窗FFT的振幅谱和相位谱，如图2所示。

图2 密集谱成分及复合信号的振幅谱和相位谱Fig.2 Amplitude spectra & phase spectra for individual dense components and their composite signal

从图2(a)的振幅谱可看出， |X1(k)|主峰位于k=5处， 在左右相邻的两根谱线(即k=4，k=6)处产生大幅值旁谱泄漏；而|X2(k)|主峰位于k=7处，在相邻的两根谱线(即k=6，k=8)处产生大幅值旁谱泄漏。很明显，由于两频率成分相隔太近，两者在共同位置k=6处产生了谱重叠，这必然引起谱间干扰而降低谱校正精度。因而，降低谱间干扰是提升密集谱校正精度的关键。

2 基于全相位滤波的密集谱频率估计算法

2.1 基于相位谱线性度偏离检测的密集谱识别判据

显然，密集谱校正的前提是能识别出密集成分。这单单依据振幅谱是不够的：以图2(a)为例， 从其复合振幅谱|X(k)|中仍能判断出k=5和k=6的两处谱峰，但却不能推断出两者的相互影响程度。因而，除了振幅谱检测外，引入相位谱检测是必要的，故本文提出基于相位谱线性度偏离检测的密集谱识别方法，原理如下。

从式(4)可知，加窗后的单频复指数信号的FFT相位谱的表达式为

(8)

具体说来，图2(b)给出了两个单频信号的相位谱φ1(k)，φ2(k)和复合信号的相位谱φ(k)(为尽量保持数值连续性，对这三个相位谱值都做了解缠绕处理)，可看出：单频相位谱φ1(k)在k=4,k=5,k=6三根谱线上(即在虚框内)可连成直线；类似地，单频相位谱φ2(k)则在k=6,k=7,k=8三根谱线上(即在虚框内)也可连成直线。

具体说来，先计算如下两个差分相位

(9)

再求取其均值，即

(10)

进而定义如下相位线性偏离因子

(11)

2.2 基于全相位滤波的谱成分分离

为提升密集谱校正精度，不妨采用滤波方法分离各个密集成分，这对滤波器性能提出如下要求：①为分离出密集成分，滤波器带宽必须足够窄；②为精确定位所关注频率的位置，滤波器的中心频点必须可以精确控制；③为不改变各成分之间的相对时延，滤波器需具有线性相位特性；④滤波器设计应足够简单。

为满足以上需求，本文选取Huang等提出的全相位窄带线性相位FIR滤波器来实现密集谱分离。该滤波器的优势在于：一个N阶的全相位滤波器蕴含了N个子滤波器的相互补偿的过程，因而传输特性好；而且全相位滤波器系数可推出解析公式[16-19]，经过Huang等的改进，对于其窄带滤波情况，其滤波器设计非常简单，仅依据下式即可确定所有滤波器系数

(12)

(13)

其中，式(13)的归一化因子C为

(14)

仍以上述例子说明对复合密集谱信号的滤波过程：图3(a)为前述的复合信号加窗FFT谱，为保留第1个谱成分，滤除k=6处的谱间干扰，将ωp=5Δω代入式(12)可得到全相位滤波器系数g(n)， 其对应的传输曲线|G(jw)|如图3(b)所示，从图3(b)可知， |G(jw)|不仅将中心频点精确定位在ωp=5Δω处，而且还可以抑制k=6处的谱间干扰。图3(c)给出了滤波前复合信号的32个样点波形, 图3(d)给出了滤后的的94个样点波形。

图3 全相位滤波器传输曲线及滤波波形Fig.3 Transfer curve of the all-phase filter and the waveforms

需指出：由于输入信号长度为N, 式(12)的全相位滤波器的系数长度为2N-1， 而FIR滤波输入输出之间为线性卷积关系，故图3(d)的滤波输出为N+(2N-1)-1=3N-2=94个样点，文献[20]指出，全相位滤波器存在(N-1)个样点的群延时，故应取从n=N=32～n=2N-1=63的N个样本段做为有效的滤波输出y(n)(如虚框所界定)。对比图3(d)和图3(c)可知，滤波输出的有效样本段y(n)相比于输入x(n)波形更为平稳，有利于提升后续的谱校正精度。进而对y(n)做汉宁窗FFT，即可得到频谱|Y(k)|。

类似地，令ωp=7Δω，重复以上过程，也可分离出第2个频率成分。

2.3 比值频谱校正法

针对密集谱情况，显然应选取使用谱线数少的校正算法。而比值校正法因仅需两根谱线(峰值谱和幅度次大的旁谱线)故成为首选。其谱校正过程如下

(1) 算出加窗FFT峰值谱和次高谱的幅度比值

(15)

u=(2-v)/(1+v)

(16)

(17)

2.4 算法总结

基于以上分析，将提出的算法流程总结如图4所示。

图4 密集谱频率估计算法流程Fig.4 Dataflow of frequency estimator for dense components

3 仿真试验

3.1 无噪试验

分别用原有比值校正法、Tsui校正法(因其单频情况做频率估计噪声估计方差接近克拉美罗限而选作对比ξ=0.5)和本文提出的基于全相位滤波的密集谱校正法(设置阈值)对图2的两个密集成分的信号进行频率校正，其校正结果如表1所示。

表1 三种校正法的频率估计结果

从表1可知，由于引入了相位线性偏离度判决和全相位滤波，本文提出的校正器比原比值校正器和Tsui估计器的精度有所提升，第1个成分的校正误差分别从0.255 9 Hz，0.152 1 Hz降低到-0.105 7 Hz，第2个成分的校正误差分别从-0.336 9 Hz，-0.176 2 Hz降低到-0.105 7 Hz。

这验证了以上提及的结论：对于现有的校正器，当未非密集谱情况时，其单频成分的频率校正均具有很高的精度；当所关注的成分周围存在干扰谱时，因没有考虑去除干扰的滤波措施，原有各校正器的内插公式准确度会降低。

3.2 加噪测量试验

(18)

表2给出了三种信噪比情况下的RMSE(Root Mean Square Error)统计结果。

表2 三种校正法的均方根误差

从表2的RMSE统计结果，可得出如下结论：

(1) 总体上，无论在哪种信噪比情况下，本文方法在估计密集成分时(即f1=5.1 Hz与f2=6.8 Hz为密集成分，f3=15.1 Hz与f4=17.2 Hz也为密集成分)，其均方根误差都低于原比值法和Tsui方法，这验证了本文提出的相位线性偏离度因子在辨识密集成分和全相位滤波在分离密集成分方面是有效的。

(2) 进一步观察，可发现SNR越高，本文方法均方误差相比于原比值法和Tsui方法的均方根误差的相对比例就越小，即精度改善越明显。这是因为，如前所述，校正器的误差包括系统误差、谱间干扰误差和噪声误差三部分，由于都是基于比值法校正，两者的系统误差可近似认为相等，而由于本文引入了全相位滤波，故谱间干扰误差相比于传统比值法要减小。故随着信噪比升高，在同样的低噪声干扰情况下，这就越加凸显出全相位滤波在降低谱间干扰误差方面对精度改善的贡献。

(3) 另外可看出，对于离密集谱较远的成分f5=22.3 Hz，本文方法的RMSE值反而高于原比值法和Tsui方法。这也很容易得到解释：一方面，对于这种非密集谱，谱间干扰其实可以忽略，如前所述，在没有谱间干扰情况下，对于现有的频谱校正器来说，几乎不存在精度提升的空间；另一方面，由于噪声干扰是随机的，总体说来，相位谱相比于振幅谱而言，对噪声的敏感度会更高些，因而对于非密集谱情况，理论上根本不需要引入相位线性偏离度判据和全相位窄带滤波的场合，反而出现小概率的判据出错情况而增大总体的RMSE误差。

4 结 论

本文提出基于全相位窄带滤波的密集谱频率估计法，该方法通过所提出的谱峰附近的相位线性偏离度因子，可以辨识密集谱是否存在，进而借助全相位窄带滤波，可实现各密集频率成分的有效分离，从而降低谱间干扰误差，再结合现有的频谱校正法(如比值校正法)，即可提升密集谱的频率估计精度。

鉴于频率估计在雷达、光学、声纳、电力、地震勘探等领域的应用广泛性，本文提出的基于全相位滤波的频率估计器具有较广阔的应用前景。