基于全相位谱分析的多频内插迭代频率估计器
2019-01-23黄翔东黎鸣诗
黄翔东， 黎鸣诗， 罗 蓬， 马 欣
(1. 天津大学 电气自动化与信息工程学院，天津 300072； 2. 国网河北省电力公司电力科学研究院， 石家庄 050021)
频率估计的应用非常广泛(如多普勒效应检测[1]，振动学中的转速测量[2]都可转化为频率估计问题)。相比于基于子空间分解(如多重信号分类法[3]、旋转不变参数估计法[4]等)的频率估计算法，FFT(Fast Fourier Transform)因其具有运算效率高的优势，因而基于FFT的频率估计器一直是学术研究的热点问题。?
Stoica等[5]指出：给定有限长样本，最大似然意义上的频率估计值位于该序列的DTFT(Discrete Time Fourier Transform)的峰值位置，但是DTFT的理想谱峰需要对频率做高分辨率扫描才可搜索得到, 显然频率扫描会增加估计器的计算复杂度。为降低复杂度, 通常的做法是对序列做FFT(而不是DTFT)，进一步对FFT峰值谱附近的谱线采取内插、迭代等措施，对FFT谱峰附近的谱值做进一步修正和细化来提升频率估计精度。从而不同的内插算法就可以产生不同的频率估计器(如Quinn估计器[6]、能量重心估计器[7]、Macleod估计器[8]、Jacobsen估计器[9]、Candan估计器[10]、相位差估计器[11-14]、比值估计器[15]、AM谱估计器[16]、细化谱估计器[17]、Tsui估计器[18]等)。
需指出，对于单频信号的FFT频率估计器，经过不断完善，其精度几乎不存在提升的空间(如文献[16]提出的AM估计器的频率估计方差仅为克拉美-罗限的1.014 7倍，而文献[18]提出的Tsui估计器的估计方差几乎逼近克拉美-罗限[19](Cramer-Rao lower Bound，CRB)。但是当对多频信号做FFT谱分析时，由于FFT固有的谱泄漏，各频率成分之间会严重的谱间干扰，导致Quinn等的估计器精度都会大幅度降低。
因而，提升多频信号的频率估计精度的关键在于：选择比FFT更优良的抑制谱泄漏性能的谱分析方法，基于此设计出具体的内插措施。王兆华等[20-21]指出，全相位FFT(all-phase FFT，apFFT)相比于传统FFT，具有更优良的抑制谱泄漏性能，黄翔东等[22]还把离散的全相位FFT谱分析延伸到连续的全相位DTFT谱分析(all-phase DTFT, apDTFT)领域，从而拓宽了其理论内涵。故本文用全相位FFT替代传统FFT谱分析，通过分析信号频偏值与谱峰附近apDTFT谱的内在联系，设计出一种迭代多频内插估计器。仿真实验表明：给定同样数目的观测样本，本文估计器比Tsui估计器具有更高的估计精度。由于工程应用中的多频信号比单频信号更为广泛，故本文估计器具有较高的应用价值。
1 全相位FFT和全相位DTFT基本原理
1.1 全相位FFT谱分析
王兆华等[23]提出的apFFT谱分析简化流程，如图1所示。
图1 apFFT谱分析简化流程(N=4)
图1的apFFT谱分析分为两个简单步骤：
步骤1全相位数据处理，用长为2N-1的卷积窗wc(n)对输入数据x(n)加权，然后将间隔为N的数据两两叠加而得到y(0),y(1),…,y(N-1)；
步骤2对y(0),y(1),…,y(N-1)做FFT即得全相位离散谱Y(k)。
图1中的卷积窗wc(n)是由长为N的窗f(n)和翻转的窗b(n)通过卷积得到，即
wc(n)=f(n)*b(-n),n∈[-N+1,N-1]
(1)
本文取f(n)=b(n)=RN(n)，其中RN是长为N的矩形窗。
令Δω=2π/N，ω0=(k*+δ)Δω,其中k*∈Z+,δ∈[-0.5,0.5)。黄翔东等[22]已证明，对于单频复指数信号
x(n)=A·ej(ω0n+θ0),n∈[-N+1,N-1]
(2)
其归一化的apFFT谱值为
k=0…,N-1
(3)
式(3)的平方项使得旁谱线Y(k),k≠k*,相对于峰值谱Y(k*)出现大幅度衰减，使得峰值谱更突出，因而apFFT具有比传统FFT更优良的抑制谱泄漏性能。
1.2 全相位DTFT谱分析
黄翔东等提出的apDTFT的简化流程，如图2所示。
图2 apDTFT谱分析的简化流程
图2的apDTFT输出为
(4)
式中:ω为频率扫描变量。因而不同于离散的apFFT谱，apDTFT谱Y(jω)为连续谱。
黄翔东等已证明，对于式(2)中的单频复指数信号，其apDTFT谱的理论值Y(jω)为
(5)
联立式(3)，式(5)，有
Y(k)=Y(jω)|ω=kΔω,k=0,…,N-1
(6)
因而apFFT谱Y(k)可看作是在ω∈[0,2π)内对apDTFT谱Y(jω)的等间隔离散采样结果。两者都是全相位谱分析理论的重要组成部分。
特别地，易证明在距离谱峰的两个对称频点ωL=ω0-0.5Δω、ωR=ω0+0.5Δω处的apDTFT值同为
(7)
2 全相位谱分析的迭代多频内插估计器设计
2.1 apDTFT峰值谱附近的频偏估计机理
ωL=kLΔω=(kc-0.5)Δω
(8)
ωR=kRΔω=(kc+0.5)Δω
(9)
把式(8)和式(9)代入式(5)，可推导出这两个频点处的apDTFT谱幅值分别为
(10)
(11)
联立式(10)和式(11)可得左谱线与右谱线幅值之比为
(12)
因为δ∈[-0.5,0.5)，则式(12)等号右边去掉绝对值符号后可表示为
(13)
当N足够大，根据无穷小关系：sin[(δ-0.5)π/N] ～(δ-0.5)π/N，sin[(δ+0.5)π/N]～(δ+0.5)π/N，则式(13)可近似表示为
由表5可知，安徽省区域旅游发展水平差异明显．黄山市与池州市旅游区位熵得分较高，旅游产业发展水平在省内处于领先地位．黄山市、池州市、安庆市、宣城市、合肥市、芜湖市、六安市旅游区位熵均大于0.74，旅游业发展状况较好，可划分为旅游发展的优势区域;其他城市均小于0.74，旅游业发展状况一般，可划分为一般旅游区域．
(14)
利用上式可解得频偏估计值
(15)
2.2 迭代估计过程
从图3可看出，只在初始化阶段做了1次apFFT谱分析，而且迭代过程中也只涉及两个特殊频点ωL、ωR的apDTFT值的更新，从而绕开了apDTFT连续谱的频率扫描过程，故图3流程复杂度低。
2.3 迭代趋势分析
搜索出apFFT峰值谱位置后，在第1次迭代中，两个频点ωL、ωR上的apDTFT谱不外乎两种情况：|Y(jωL)|≤|Y(jωR)|和|Y(jωR)|<|Y(jωL)|，其分布分别如图4(a)、4(b)所示。
图3 单频迭代估计流程图
(a) |Y(jωL)|<|Y(jωR)|
(b) |Y(jωL)|>|Y(jωR)|
理想情况下，由式可推知|Y(jωL)|=|Y(jωR)|成立。故对于图3迭代流程，可采用两者差异的比例作为迭代终止条件。经验表明，将阈值取为ε∈(0.004，0.3)，只需2次迭代即可获得高精度频率估计。
2.4 多频成分的频率估计
需指出，与单频成分不同的是，多频成分的频率估计误差来源除了外加噪声和估计器本身的系统误差(如式(14)中等价无穷小引入的误差)外，各成分的谱间干扰也是误差源的主要方面。由于apFFT谱分析具有比传统FFT谱分析更优良的抑制谱泄漏性能，故谱间干扰更小，因而其多频估计器精度高于文献[6-18]中的估计器。以上结论将通过第3节中的实验来证明。
2.5 多频信号重构
(16)
(17)
式中,ang(·)表示取相角操作。
从而，包含多个频率成分的实信号的重构式为
n∈[-N+1,N-1]
(18)
3 实 验
例1令N=64,M=3,ω1=7.1Δω,ω2=9.35Δω,ω3=12.2Δω,对信号x(n)=ejω1n+3ejω2n+ejω3n,n=-N+1,…,N-1，分别进行64阶apFFT、apDTFT及其127点的FFT和DTFT，其谱图分别如图5(a)、5(b)所示。并且分别用本文的频率估计器和文献[18]的Tsui估计器(选该估计器的原因是：其误差方差逼近CRB，几乎是精度最高的FFT估计器)分别进行频率估计，表1给出两个估计器的相对误差，其算式为
(19)
从图5明显可看出，apFFT、apDTFT的谱泄漏远小于FFT、DTFT的谱泄漏(即除3个主峰外，前者的带外起伏比后者小得多)，这意味着对于3个成分间的谱干扰，前者也比后者小得多。从表1数据也可看出，本文估计器的相对误差也远小于Tsui估计器情况。
(a) apFFT、apDTFT谱图
(b) FFT、DTFT谱图
η1/%η2/%η3/%本文估计器0.120.010.06Tsui估计器0.330.030.36
例2调研含噪情况下，不同频率间距的各成分的谱间干扰对频率估计器的影响。不妨令N=64,M=2,输入信号为x(n)=ejω1n+ejω2n+ξ(n)，其中ξ(n)为加性高斯复噪声。固定ω1=5.25Δω，ω2分别取值为9.15Δω，8.15Δω，7.15Δω，6.15Δω。图6～图9给出不同信噪比(Signal to Noise Ratio, SNR)条件下，本文估计器、Tsui估计器的均方根误差以及CRB的平方根曲线。
从图6～图9的RMSE曲线，可得出如下结论：
(1) 在低信噪比区域(SNR<10 dB)，本文估计器的SNR阈值普遍比Tsui估计器低，这意味着其抗噪鲁棒性强于Tsui估计器。这是因为对于Tsui估计器，多频信号较大的谱间干扰削弱了抵御强噪声的能力。
(2) 在中、高信噪比区域(SNR>10 dB)，对于Tsui估计器，其RMSE呈现出平坦的形状，与CRB曲线偏离非常远，而本文估计器仅在两频率成分间距较小时(即图8、图9的情况)，才会与CRB曲线产生较大的偏离。这说明，在中、高信噪比情况，影响估计器精度的主导因素，既不是外加噪声，也不是系统误差，而是各成分间的干扰，而全相位谱分析的谱间干扰远小于传统谱分析情况，故其RMSE曲线更贴近于理论下限。
(3) 当频率间距过小时(低于1Δω，即对应图9的情况)，两个估计器的RMSE都变得很大，这是因为无论全相位谱，还是传统FFT谱，其各成分谱的主瓣都相互重叠了，故频率估计失效。
(a) ω1=5.25Δω
(b) ω2=9.15Δω
(a) ω1=5.25Δω
(b) ω2=8.15Δω
(a) ω1=5.25Δω
(b) ω2=7.15Δω
(a) ω1=5.25Δω
(b) ω2=6.15Δω
γ=10lg(Px/Pe)
(20)
式中，Px为原信号x(n)的平均功率，Pe为e(n)的平均功率。表2给出了两种估计器的重构。
表2 机械故障信号的两种估计器重构信噪比
以测量面缺陷信号为例，图10(a)、(b)分别给出了用全相位FFT的振幅谱|Y(k)|、传统FFT的振幅谱|X(k)|，图10(c)给出了原始波形、全相位重构波形和Tsui重构波形。
(a) apFFT振幅谱
(b) FFT振幅谱
(c) 原始波形与两种重构波形
可明显看出，图10(a)的振幅谱|Y(k)|的谱间干扰程度比图10(b)的振幅谱|X(k)|的谱小很多，这保证了前者的各成分的频率、幅值、相位参数估计精度高于后者，从而获得更高的重构信噪比(如表2所示，前者比后者高出2.7 dB)。该结论对于其他两种故障信号的重构也成立。
4 结 论
本文提出基于全相位谱分析的多频内插迭代频率估计器，利用全相位谱分析方法有效地抑制谱泄漏，降低多频成分间的干扰；再结合内插迭代提高多频信号的频率估计精度。由于本文方法在有外界干扰的低信噪比环境仍能高精度地做频率估计，且仿真实验和实际信号重构实验均验证了本文方法具有更高的参数估计精度，故在多频信号估计领域有望得到更广泛的应用。
