基于迭代更新策略的快速高精度频率估计方法

高志峰, 彭喜元， 彭 宇

1.哈尔滨工业大学 自动化测试与控制研究所，哈尔滨 150080；2. 山东大学 机电与信息工程学院，山东 威海 264209)



基于迭代更新策略的快速高精度频率估计方法

高志峰1,2, 彭喜元1， 彭 宇1

1.哈尔滨工业大学 自动化测试与控制研究所，哈尔滨 150080；2. 山东大学 机电与信息工程学院，山东 威海 264209)

基于滤波器组的谱估计方法用于信号频率估计频率虽分辨率高,但滤波器组中心频率格点划分缺乏先验知识，谱峰搜索过程存在计算复杂、信号不匹配问题。基于MVDR(Minimum Variance Distortionless Response)谱构造指标函数，将谱峰搜索问题等价为标量指标函数局部极小值问题求解，通过构造一组最速下降方向及自适应步长对极值频率迭代更新，实现对信号频率的直接估计。新算法不仅回避中心频率点组划分及传统的谱峰搜索，且有效缓解信号不匹配，估计精度、计算效率更高。对单成分信号频率估计精度与计算量进行新算法与现有算法的性能比较，并用于多成分信号频率估计。

频率估计；迭代算法；自适应步长；MVDR

信号频率估计在机械振动信号监测及故障特征参数提取、电网状态监测、雷达、通讯、语音信号识别等众多领域占据重要地位[1-8]。频率估计方法主要有参数化法及非参数法两类。参数化方法假设信号满足特定的数学模型，频率分辨率更高，但模型失配时性能不佳且数学模型对相关参数初始值有诸多限制；而非参数化法无需建模应用更灵活，且可同时获取频率、幅值及相位等特征信息，但大多通过频域谱峰搜索实现频率估计，普遍存在计算量与估计精度的矛盾。如何快速、准确得到信号成分的频率估计具有重要研究、应用价值。

最大似然估计(Maximum Likehood, ML)方法及非线性最小二乘(Nonlinear Least Square,NLS)方法可获得最优信号频率估计，但计算量巨大，且提高频率分辨率需较长采样序列。基于频域的改进算法大多包含两步，即采取快速离散傅里叶频谱粗搜索提高计算效率，进而采取相位差法(Phase Different,PD)或迭代频域细化搜索或频域插值(frequency-domain interpolation)[9-10]法提高估计精度。在时域中通过信号建模可提高计算效率，包括基于线性预测(Linear Prediction,LP)模型及相位平均(Phase Average,PA)[11-14]。基于采样数据自相关函数(autocorrelation-based)方法与基于迭代自回归滑动(Autoregressive Moving Average，ARMA)模型[15]方法等。信号建模频率分辨率较高，但大多属于次优估计，在低信噪比时频率估计性能不佳且存在模型失配问题。以上算法多适用于单频率信号。

将短采样数据上具有高分辨率的功率谱估计方法用于多频率成分信号频率估计[16-20]。基于协方差矩阵子空间分解的多重信号分类(Multiple Signal Classification, MUSIC)[16]及旋转不变信号参数估计(Estimation of Signal Parameters via Rotational Invariance Technique,ESPRIT)[17]方法具有噪声抑制能力，但其特征空间分解缺乏先验知识，易造成空间维数判别失衡而得到有偏估计。基于匹配滤波组理论[18]的谱估计方法，利用一组均匀分割的频率格点进行功率谱估计，并通过谱峰搜索实现。而预划分频率格点数目有限难以与信号频率成分完全匹配，存在固有信号不匹配(Signal Mismatch Problem,SMP)问题[19]，包括最小方差无失真响应 (Minimum Variance Distortionless Response,MVDR)[20]及幅值与相位估计(Amplitude and Phase Estimation,APES)[21]等方法。理论表明傅里叶谱估计器也属于匹配滤波器组谱估计方法[22]。基于预划分均匀频率格点的谱峰搜索，计算量与估计精度间存在固有矛盾关系。采取两步式搜索实现频率估计[23]：进行谱峰频率粗搜索，采取二分法迭代搜索以提高估计精度缓解计算量，但粗搜索步骤所需计算量不可忽视。

本文提出基于迭代策略的谱峰频率搜索算法，以回避传统的谱峰搜索方式，在提高估计精度的同时降低计算量。基于MVDR谱定义标量指标函数，由指标函数极小值估计MVDR谱峰频率。新算法基于最速下降法，由标量梯度函数构造一组迭代方向并计算最优更新步长，且迭代更新过程收敛到指标函数局部极小值。因无需预先划分频率格点，可缓解信号不匹配，且在频率估计精度提高的同时降低计算量。新算法适用于单频率信号，若结合初始频率选择策略，可直接推广到多频率信号频率估计。在仿真、试验中从精度、计算量角度对算法性能进行分析比较。

1 基于MVDR谱的频率估计

含噪声的时域采样信号可近似为多成分正弦信号组合，即

(1)

式中：Ai,fi,φi分别为频率成分幅值、频率、初始相位；Ts为采样周期；v(n)为加性高斯白噪声。

则问题等价为通过MVDR谱峰角频率ωi得到对信号时频率fi的估计值fi=ωi/(2πTs)。

长度为N的采样信号通过Hilbert变换转换为复信号进行功率谱估计，其MVDR谱估计记为

(2)

式中：上标H表示共轭转置；R为L维协方差矩阵；f(ω)为定义在角频率ω∈[0,2π)的L维搜索向量，即

f(ω)=[1 e-jω… e-j(L-1)ω]T

(3)

式中：上标T表示矩阵及向量转置运算。

若搜索向量f(ω)与信号x(n)频率相同，则对应MVDR谱中的尖锐谱峰(局部极大值)。取式(2)分母作为指标函数，即

g(ω)=fH(ω)R-1f(ω)

(4)

则所有的MVDR谱峰均对应指标函数g(ω)的局部极小值点。本文采取最速下降方法构造g(ω)局部极小值频率迭代式为

ωk=ωk-1-μkωg(ωk-1)

(5)

式中：μk为第k次迭代中正实数值步长；ωg(ωk-1)为指标函数g(ω)在角频率ωk-1处标量梯度函数。

2 迭代频率估计算法

2.1 梯度函数及自适应步长

搜索向量f(ω)关于角频率ω的梯度函数定义为

[0 … -j(L-1)e-j(L-1)ω]T=-jJf(ω)

(6)

式中：J=diag([0,1,…,L-1])为常系数对角矩阵。

指标函数g(ω)关于角频率ω的梯度函数定义为

jfH(ω)JR-1f(ω)-jfH(ω)R-1Jf(ω)

(7)

由jfH(ω)JR-1f(ω)=[-jfH(ω)R-1Jf(ω)]H简化为

ωg(ω)=-2Im[fH(ω)JR-1f(ω)]

(8)

则频率ωk上的梯度函数ωg(ωk)为实数值标量函数，且在局部极小值频率上满足)=0。

由式(4)、(5)知，最优迭代步长μk对应最小化问题

(9)

随角频率变量ωk迭代更新，函数g(ω)单调下降并收敛到局部极小值点。为获得指标函数g(ω)闭式表达式，将指数函数ejmωk进行一阶泰勒展开并线性化处理

e-jmωk≈e-jmωk-1-jme-jmωk-1(ωk-ωk-1)≈

e-jmωk-1[1+jmμkωg(ωk-1)]

代入式(3)得f(ωk)关于f(ωk-1)的近似表达式为

f(ωk)=(I+jμkωg(ωk-1)J)f(ωk-1)

(10)

式中：I为L维的单位对角阵。

将式(10)代入式(4)，得函数g(ωk)闭式表达式为

g(ωk)=fH(ωk-1)R-1f(ωk-1)-

jμkωg(ωk-1)fH(ωk-1)JR-1f(ωk-1)+

jμkωg(ωk-1)fH(ωk-1)R-1Jf(ωk-1)+

(11)

对步长μk求偏导，得梯度函数表达式为

μg(ωk)=-jωg(ωk-1)fH(ωk-1)JR-1f(ωk-1)+

jωg(ωk-1)fH(ωk-1)R-1Jf(ωk-1)+

(12)

将式(7)代入式(12)，得

(13)

由于Jf(ωk-1)为非零向量及逆矩阵R-1存在，则迭代步长μk均为正实数标量。

2.2 迭代算法和收敛性分析

由式(5)、(7)、(13)获得指标函数局部极值的迭代求解算法(Single local minimum Iterative algorithm，SITER )，主要步骤为

步骤1: 起始角频率ω0，Q=JR-1,U=JR-1J，k=1

步骤2: 由式(7)计算频率点ωk-1对应的梯度函数

ωg(ωk-1)=-2Im[fH(ωk-1)Qf(ωk-1)]

步骤3: 由式(13)计算迭代步长μk

μk=1/[2fH(ωk-1)Uf(ωk-1)]

步骤4: 更新频率变量ωk=ωk-1-μkωg(ωk-1)，k=k+1，跳转到Step2继续执行。

证明：由式(5)及迭代步长μk>0，性质1自然成立。

证明：由式(4)，对可逆矩阵R, 指标函数g(ω)≥0为正值标量函数，存在下界。由式(9)获得自适应步长满足g(ωk)≤g(ωk-1)，即指标函数值单调下降。由式(11)得，当且仅当ωg(ωk-1)=0时等式成立g(ωk)=g(ωk-1)，即得到局部极小值点。

2.3 频率估计解析解法

若信号中除噪声外仅含单频率成分，取矩阵维数L=2即可区分。由于此时指标函数g(ω)仅存在一个极小值点，可推导获得频率估计得解析表达式。记L=2维满秩酉对称协方差矩阵R及逆矩阵R-1分别为

(14)

式中：a22>0为正实数；a12=-r12/det(R),det(R)>0为矩阵R的行列式。

式(7)进一步简化为

ωg(ω)=-2Im(fH(ω)JR-1f(ω))=

-2Im(a12e-jω+a22)=

-2{cosωIm(a12)-Re(a12)sin(ω)}

(15)

(16)

由a12=-r12/det(R), 可回避矩阵求逆运算，得

(17)

则单成分信号频率可由式(17)直接求得，计算时间几乎可忽略。式(17)为SITER算法在L=2的特例，不采取迭代过程，直接由(17)式获得信号成分频率估计 (Direct solution of the Iterative algorithm，DITER)。该算法仅适用于单频率成分信号频率估计，简单易行，其频率估计值常作为单成分信号时SITER算法迭代频率初始值。

3 仿真实验

3.1 单频率成分信号的频率估计

式(1)中频率f=203.1 Hz，A=1，采样长度32，采样频率fs=1 500 Hz，采样周期Ts=1/fs。

DITER算法与Rife-ML方法及未加窗的线性预测(Kay-ULP)、相位平均(Kay-UPA)方法的均方误差(Mean Square Error, MSE)随输入信噪比(Signal Noise Ratio, SNR)变化，并与Cramer-Rao Lower Bound (CRLB)比较见图1。由图1看出，Rife-ML方法基于周期图谱峰搜索实现，验证ML估计的最优性，但计算量繁重。而DITER、Kay-ULP、Kay-UPA算法性能接近，与CRLB存在约6 dB阈值，因数据截取时均未加窗导致频谱泄漏，可加窗修正。DITER算法计算量少为优点所在。

图1 不同信噪比对应的均方误差Fig.1 Mean square error versus different SNR

取L=2, SNR=40 dB，在K=200，2 000，3 000频率格点的MVDR谱估计结果见图2。由图2看出，K越大频率间隔越小，频率估计精度越高，但计算量呈级数增加。图2(d)为由式(17)直接所得DITER算法的频率估计值，并由式(2)计算对应MVDR谱估计，其频率估计最接近真实频率203.1 Hz。

图2 高信噪比时单频率估计结果Fig.2 The frequency estimation with high SNR

不同信噪比及最大迭代次数对应的均方误差见图3。将矩阵维数增大L=8以获取更高频率分辨率及更窄匹配滤波器主瓣，提高估计精度。以DITER算法(L=2)的频率估计值作为初始频率, 用SITER算法进行迭代频率更新。由图3看出，随SITER算法迭代次数增大(分别为3,5,10,100,200)，其估计值的均方误差MSE越接近CLRB，频率估计精度提高。

图3 不同信噪比、最大迭代次数对应的均方误差Fig.3 Mean square error versus different SNR and different maximum iterations

图4为L=8及SNR=20 dB时，SITER算法不同迭代次数对应的计算量与MSE值，并与DITER算法(L=2)、标准MVDR算法(频率格点数K=200, 2 000,3 000)进行比较(所有结果均为运行1 000次的平均值)。图4(a)中SITER算法耗时与DITER算法相当，略少于K=200点的MVDR算法进行频率估计及谱峰搜索所需时间；但图4(b)中SITER算法的MSE远小于其它结果，对应高估计精度。SITER算法打破了估计精度及复杂计算量间固有矛盾关系。

图4 不同迭代次数对应的时间和均方误差Fig.4 Time and mean square error versus iteration number

3.2 多频率成分信号频率估计

对含随机噪声的多正弦成分信号，有

式中：fi分别为各成分频率{100,254,312,410,503}Hz；信噪比取20 dB, 连续采样256点(采样频率1 500 Hz)，初始相位随机φi∈[0,2π)，取协方差矩阵维数L=30。

MVDR谱分别在100，500，4 000个频率格点的计算见图5(a)、(b)、(c)。图5(d)中共需独立5次SITER算法以获得5个信号频率成分，对应的初始角频率可采取傅里叶频谱谱峰粗搜索获得，或用MVDR算法在K=100频率点上进行谱峰粗搜索获得。

图5 MVDR与SIter算法频率估计Fig.5 The frequency estimations of MVDR and SITER algorithm

不同算法对应的频率估计结果见表1。由表1看出，由SITER算法所得频率估计精度优于MVDR算法在4 000个频率点上所得频率估计结果。

表1 不同算法对应的频率估计结果

4 结 论

(1) 相对于低分辨率的傅里叶频谱，匹配滤波器组谱估计方法在较短采样数据上具有高分辨率。由于信号频率成分数目、位置未知，传统谱峰搜索方法普遍存在谱线不匹配(SMP)问题，估计精度及复杂计算量之间存在固有矛盾关系，影响估计精度的提升。

(2) 为提高频率估计精度、降低计算量，提高计算效率，本文将谱峰搜索问题转化为标量指标函数的局部极值计算问题，获得局部极小值迭代算法SITER。最速下降搜索方向及最优步长的选取满足指标函数单调下降，收敛到局部极小值点。针对单成分信号频率估计，SITER算法简化为解析算法DITER，无需迭代而直接计算信号频率估计值。SITER算法推广到多成分信号频率估计应用中，可对多成分信号频率估计算法有效补充。

(3) 在SITER算法中采取最速下降法，主要考虑其迭代收敛过程中指标函数值及频率更新的单调性，但该方法收敛速度较慢，尤其在极值点附近。

[1] Stoica P, Moses R L.Spectral analysis of signals[M]. Pearson Prentice Hall, 2005.

[2] Li J, Stoica P. Roubust adaptive beamforming[M]. New York:Wiley, 2005.

[3] 胡爱军,朱瑜.基于改进峰值搜索法的旋转机械瞬时频率估计[J].振动与冲击,2013,32(7):113-117. HU Ai-jun, ZHU Yu. Instantaneous frequency estimation of a rotating machinery based on an improved peak search method [J]. Journal of Vibration and Shock,2013, 32(7):113-117.

[4] 沈廷鳌,涂亚庆,张海涛,等.一种改进的自适应格型陷波频率估计算法及其收敛性分析[J].振动与冲击,2013, 32(24): 28-32. SHEN Ting-ao,TU Ya-qing, ZHANG Hai-tao, et al.A modified frequency estimation method of adaptive lattice notch filter and its convergence analysis[J].Journal of Vibration and Shock,2013, 32(24):28-32.

[5] Rife C D, Boorstyn R R. Single-tone parameter estimation from discrete-time observations[J]. IEEE Trans. On Information Theory, 1974,20(5): 591-598.

[6] Stoica P, Nehoral A.Statistical analysis of two non-linear least squares estimators of sine waves parameters in the colored noise [J]. Proceddings of the ICASSP, 1998, 4:2408-2411.

[7] 胡文彪,夏立,向东阳,等. 一种改进的基于相位差法的频谱校正方法[J]. 振动与冲击,2012, 31(1): 162-166. HU Wen-biao, XIA Li, XIANG Dong-yang, et al. An improved frequency spectrum correction method based on phase difference correction method[J]. Journal of Vibration and Shock, 2012, 31(1): 162-166.

[8] Lin H, Ding K. Energy based signal parameter estimation method and a comparative study of different frequency estimators [J]. Mechanical Systems and Signal Processing, 2011, 25:452-464.

[9] Fu H, Kam P. Sample autocorrelation function based frequency estimation of a single sinusoid in AWGN[C]// Vehicular Technology Conference, IEEE 75th, 2012:1-5.

[10] Aboutanios E, Mulgrew B. Iterative frequency estimation by interpolation on Fourier coefficients[J].IEEE Trans. Signal Processing, 2005, 53:1237-1242.

[11] Jackson L, Tufts D, Soong F, et al. Frequency estimation by linear prediction[J].IEEE International Coference on Acoustics,Speech and Signal Processing,1978, 3:352-356.

[12] Kay S. A fast and accurate single frequency estimator[J]. IEEE Transaction on Acoustics,Speech and Signal Processing, 1989,37(12):1987-1990.

[13] Lui K, So K. Two-stage autocorrelation approach for accurate-single sinusoidal frequency estimation[J]. Signal Processing, 2008, 88(7):1852-1857.

[14] Cao Y,Wei G, Chen F. A closed-form expanded autocorrelation method for frequency estimation of a sinusoid[J]. Signal Processing, 2012, 92:885-892.

[15] Quinn B G, Fernandes J M. A fast efficient technique for the estimation of frequency[J]. Biometrika, 1991,78(3): 489-497.

[16] Schmidt R. Multiple emitter location and signal parameter estimation[C]// Proc. RADC spectrum estimation Workshop, 1979:243-258.

[17] Roy R,Kailath T. Esprit-estimation of signal parameters via rotational invariance techniques[J]. IEEE Trans. Acoust. Speech Signal Process., 1989, 37: 988-995.

[18] Stoica P, Jakobsson A, Li J. Matched-filter bank interpretation of some spectral estimators[J]. IEEE Trans. Signal Processing, 1998,66: 45-59.

[19] Cox H. Resolving power and sensitivity to mismatch of optimum array processors[J]. Journal of the Acoustic Society of America,1973, 54: 771-785.

[20] Benesty J, Chen J, Huang Y. A generalized MVDR spectrum [J]. IEEE Signal Process, 2005, 12(12): 827-830.

[21] Stocia P, Li H, Li J.A new derivation of the APES filter[J]. IEEE Signal Process, 1999, 6(8):205-206.

[22] Zheng C, Zhou M, Li X. On the relationship of non-parametric methods for coherence function estimation[J]. Signal Processing, 2008, 88:2863-2867.

[23] Peng Y, Gao Z, Peng X. MVDR spectral estimation by spectral peak dichotomous search[C]//I2MTC, 2012:1692-1696.

Fast and accurate frequency estimation with a gradient-based iterative algorithm

GAO Zhi-feng1, 2, PENG Xi-yuan1, PENG Yu1

1.Harbin Institute of Technology, Automatic Test and Control Institute, Harbin 150001, China;2.Shandong University, School of Mechanical, Electrical and Information, Weihai 264209, China)

Nonparametric spectral estimation algorithms were applied to frequency estimation for their significant performance. To avoid the signal mismatch problem and to improve the frequency estimation accuracy, a new iterative algorithm was presented based on the minimum variance distortionless response (MVDR) spectrum. With given initial frequency, searching directions and adaptive steps were derived to update the frequency sequence, which converges to a local spectral peak as the scalar gradient function goes to zero. Without the spectral peak searching on predefined analysis frequency grids, the computation is saved. The proposed algorithm was also applied to multiple component frequency estimation with carefully selected initial frequencies.

frequency estimation; iterative algorithm; adaptive step; MVDR

国家自然科学基金(61304142，61305130)

2014-04-17 修改稿收到日期：2014-06-03

高志峰 男，讲师，1979年5月生

彭喜元 男，教授，1961年12月生

TP202.4

A

10.13465/j.cnki.jvs.2015.14.004