极端频率直接估计的新自适应陷波器方法

李 明， 涂亚庆， 沈廷鳌， 毛育文

(后勤工程学院信息工程系，重庆 401311)

引 言

振动工程中常见的极端频率信号，其频率估计精度具有很重要的理论意义和使用价值[1,2]。在离散频谱校正理论中的极端频率，主要相对于频率分辨率而言(即采样频率和采样时间长度的比值)，由于极端频率频谱泄露严重、存在栅栏效应和负频率成分干涉，基于离散频谱校正的频率估计方法存在较大误差[3,4]，且计算复杂，抗噪性较差，精度还有待进一步提高[5]，也无法估计时变信号的频率值。

为改变此现状，提高频率估计方法的实时性、准确性和抗噪性，采用自适应陷波器(ANF)进行频率估计。ANF主要是根据被处理信号的特点，对其进行参数优化，自动调节自身模型参数，令误差函数达到最小值，使陷波频率与信号频率相等，从而估计出信号频率值。每估计一次频率值只需要3个采样点，同离散频谱校正理论中采样时间长度的含义完全不同。该方法从时域角度进行迭代频率估计，完全避免了负频率成分的干扰，不存在频谱泄露和栅栏效应问题，不仅可以估计频率恒定的时不变信号，还可以估计频率发生变化的时变信号；同时可以滤除信号噪声，提高信噪比，针对噪声环境下的弱信号也具有较好的检测效果，且结构简单、计算量小[6]。

由于ANF频率估计方法同离散频谱校正理论的差别较大，故本文所指的极端频率信号，主要指归一化频率接近于0或π两端的信号。例如，当采样频率设为4 096 Hz，型号为6307E型的滚动轴承，其内圈、外圈的故障特征频率分别为76.5和98.8 Hz[7]，相对于采样频率4 096 Hz则为极低频信号，其归一化频率接近于0。此外，当现有的采集装置的采样频率限定为1 024 Hz时，旋转机械转子系统松动故障产生的10倍频简谐分量(基频50 Hz，50×10 Hz=500 Hz，接近奈奎斯特频率512 Hz)相对于采样频率1 024 Hz则为极高频信号，其归一化频率接近于π。

针对上述信号进行频率估计时，由于噪声干扰增大，其频率估计精度将出现明显下降，且存在算法不稳定，收敛速度偏慢的缺点[8,9]。针对此问题，文献[10,11]对ANF误差函数进行了改进，使该误差函数具备一定斜直线的特性，一定程度上加快了ANF收敛速度。文献[12,13]针对上述文献结果有偏的问题进行了研究，取得了近似无偏的频率估计结果，但上述方法的缺点是最优频率估计解为局部极小值点，对频率初始值的设定有一定的要求，且为间接频率估计，对估计后参数的取值有一定的限制；同时其频率估计结果存在振荡，针对极端频率信号的估计精度和稳定性还有待进一步提升。

为此，本文为解决现有ANF方法估计极端频率信号时存在收敛速度慢、精度不高、稳定性差的问题，通过提出一种新误差函数，改善ANF在整个频率范围内的迭代收敛曲面，加快ANF收敛速度与算法稳定性；同时，采用偏差补偿方式，降低噪声对ANF的影响，提供近似无偏的极端频率信号频率估计结果。

综上所述，本文提出了一种ANF极端频率信号直接频率估计新方法，并与原有算法进行了有关性能的比较分析。

1 极端频率信号

设正弦信号为

(1)

根据采样定理，fs≥2f0，故可知信号频率的取值范围为0<ω0≤π；同时现有的ANF研究表明[14]，当信号频率0<ω0≤0.1π或0.9π<ω0≤π时，ANF的性能出现明显下降，所以本文讨论的极频信号主要指频率0<ω0≤0.1π的极低频信号或0.9π<ω0≤π的极高频信号。

2 新自适应陷波器方法

2.1 新误差函数

ANF传递函数为

(2)

式中ρ为极半径，控制陷波宽度，且0≪ρ<1。ANF参数a=-2cosω，ω为ANF陷波频率，在频率估计过程中，ω将趋近于ω0，故ω可看作为ω0的估计值，此时，a→a0=-2cosω0。

在频率估计过程中，若估计a值，则为间接频率估计，此时需保证-2≤a≤2，以使ω=arccos(-a/2)有解。但若估计ω，则为直接频率估计，不需保证ω的取值范围，步骤较为简化。

为此，将a=-2cosω代入式(2)可得

(3)

式中N(z,ω)和H(z,ω)分别为ANF的FIR和IIR结构。则信号x(k)通过N(z,ω)和H(z,ω)的信号e1(k)和e2(k)分别为

(4)

为了获得最佳频率估计值，ANF通过调整陷波频率ω，使下式的新误差函数最小化

J(ω)=E[e2(k)(e1(k)+e2(k))]

(5)

式中E[·]表示求取期望，其保留了传统IIR自适应陷波器频率估计精度高，在最优频率值附近下降速度快的优点；同时使误差函数在远离最优频率值时仍能够保持一定的梯度值，使其可保持较快的收敛速度。在实际计算中可按下式进行计算

(6)

由于输入信号为正弦信号，故式(3)中N(z,ω)和H(z,ω)在输入信号频率ω0处的幅值为

(7)

相角为

(8)

结合式(4)，(7)和(8)，可得

(9)

将式(7)～(9)代入式(5)，可得

J(ω)=E[e2(k)(e1(k)+e2(k))]=

(10)

为说明所提新误差函数的优越性，故比较文献[10]所采用的如下式所示的原有误差函数为

(11)

实际计算时按下式计算

(12)

为说明本文所提如式(5)所示新误差函数的优越性，故所取参数同文献[10]保持一致，A=1，θ=π/6，ρ=0.95，σ2=0，输入信号频率分别取3个不同的值，π/3，π/2和2π/3，则文献[10]和本文方法所提误差函数如图1所示。又由于极高频和极低频具有对称性，故只分析频率ω0为0.01π时的极低频信号，其他参数设置保持不变，则极端频率下相应的误差函数如图2所示。

图1 不同方法所提误差函数的理论值与计算值

图2 ω0=0.01π时不同方法所提误差函数

由图1，2可知，两种不同误差函数的理论计算值同实际计算值基本保持一致。同时，本文所提误差函数J(ω)相较J0(ω)的梯度值较大，在迭代计算过程中的收敛速度将更快，且最优频率值的选取较为明显。特别当输入信号频率为极低频，ω0=0.01π时，文献[10]所提误差函数在这一区域附近基本为平直线，丧失了对最优解的收敛性能，且容易在最优解附近产生收敛振荡，导致其算法不稳定。而本文所提误差函数仍然具有较大的梯度值，仍可收敛至最优频率值，提升了算法的稳定性。

图3 ω0=π/3时本文方法所提误差函数曲面图

图4 ω0=π/3时原有误差函数曲面图

图5 ω0=0.01π时本文方法所提误差函数曲面图

图6 ω0=0.01π时原有误差函数曲面图

2.2 频率估计偏差问题分析

结合式(5)所示新误差函数，可知ANF直接频率估计方法为

ω(k)-μ[(g1(k)+2g2(k))e2(k)+e1(k)g2(k)]

(13)

式(13)中的e1(k)g2(k)相对于其他项较小，为简化计算，可略去e1(k)g2(k)[10]，于是可得

ω(k+1)=ω(k)-μ[g1(k)+2g2(k)]e2(k)

(14)

令

(15)

可知g1(k)和g2(k)分别为x(k)通过G1(z,ω)和G2(z,ω)后的信号。

一般而言，式(14)所示频率估计方法结果有偏，为消除该偏差，对式(14)进行稳态条件下的偏差分析，此时ω→ω0，故可得[15]

(16)

稳态下，由于ω→ω0，故可用ω0替代式(15)中的ω(k)，由此可得e1(k)，e2(k)，g1(k)和g2(k)的表达式

(17)

其中，δω(k)=ω(k)-ω0；v1(k)，v2(k)，v3(k)和v3(k)为噪声v0(k)分别通过N(z,ω)，G1(z,ω)，H(z,ω)和G2(z,ω)后所产生。

由此，将式(14)的两边同时减去ω0，可得

δω(k+1)=δω(k)-μ[g1(k)+2g2(k)]e2(k)

(18)

对式(18)两边同取期望，并将式(17)代入。同时，令Ri,j=Rj,i=E[vi(k)vj(k)]表示噪声vi(k)和vj(k)的相关值，计算时由于ρ→1，为简化计算，设(1-ρ)2≈0，则可得

E[δω(k+1)]=E[δω(k)]-μE[g1(k)+2g2(k)]e2(k)=(1-μM1)E[δω(k)]+

(19)

由式(19)可知，其右边最后一项，即输出信号间噪声的相关性，是导致如式(14)所示频率估计算法结果有偏的根本原因。同时，该项是输入信号噪声σ2、极半径ρ和频率估计结果ω的函数，且无论如何调整上述参数，都无法使该项为0。为此，需消除此项的影响，提高其频率估计精度。

2.3 无偏频率估计方法

为得到无偏频率估计结果，需对输入信号噪声σ2进行估计。为此，令c(k)=4(ρ-1)sin2ω(k)，分析可知

E[c(k)x(k)e1(k)]=

c(k)A2sinω0cosφ1E[δω(k)]+c(k)σ2

(20)

式(20)中右边最后一项包含了2R3,4+R2,3的值，于是，将式(20)中的c(k)x(k)e1(k)代入式(18)，两边同取期望。虽然引入了多余项c(k)A2sinω0cosφ1·E[δω(k)]，但是该项可与式(19)中的M1合并，对频率估计精度的影响较小，由此可知基于新误差函数的ANF极端频率无偏直接估计新方法为

ω(k+1)=ω(k)-μG(k)

(21)

其中，G(k)=[g1(k)+2g2(k)]e2(k)-c(k)x(k)e1(k)。

该算法流程为：

Step1. 设置频率初始值ω(0)，算法步长μ和ANF参数ρ；

Step2. 在k时刻，将信号x(k)分别通过N(z,ω)和H(z,ω)，得到e1(k)和e2(k)，并计算c(k)；e1(k)=x(k)-2cosω(k)x(k-1)+x(k-2),e2(k)=e1(k)+2ρcosω(k)e2(k-1)-ρ2e2(k-2),c(k)=4(ρ-1)sin2ω(k)。

Step3. 计算k时刻的信号g1(k)和g2(k)，同时计算G(k)的值；

G(k)=[g1(k)+2g2(k)]e2(k)-c(k)x(k)e1(k)

Step4. 按式(21)更新频率估计值，得到

ω(k+1)=ω(k)-μG(k)

Step5. 重复步骤Step2～Step4，直至算法收敛。

3 性能分析

3.1 稳态偏差分析

对式(21)进行稳态条件下的偏差分析，式(21)两边同减去ω0并取期望，得

E[δω(k+1)]=E[δω(k)]-μE[G(k)]=

[1-μ(M1-c(k)A2sinω0cosφ1)]E[δω(k)]+

(22)

稳态条件下，可知

(23)

式(22)可变化为

(24)

3.2 稳态MSE分析

对式(21)两边同减去ω0并平方，求取期望可得

(25)

式(25)左边第二项的计算可利用文献[13]所用方法，得

2μE[G(k)δω(k)]=2μE[δω(k-2)G(k)]-

(26)

则结合式(26)，将式(25)展开，经过一系列复杂繁琐的推导后，可得

(27)

其中，ψ1=A2sinω0[6Bcos(ω0-φ2)-2ccosφ1]，ψ2=A4sin2ω0[9B2+9B2cos(2φ2-2ω0)/2-

3cBcos(φ1-φ2-ω0)-6cBcos(φ2-ω0)cosφ1]，

2cos2(iω0)]+6Bc[cos(φ1-φ2-ω0)cos(2iω0)+

2cos(ω0-φ2)cosφ1]}，ψ4=4sin2ω0(1-ρ)(13-32cos2ω0)σ4。

结合式(23)可知，稳态下式(27)可化为

(28)

4 计算验证

4.1 频率估计精度分析

图7 ω0=0.01π时ANF频率估计效果图

图8 ANF频率估计偏差E[δω(k)]和均方差

图9 ANF在全频段频率估计偏差E[δω(∞)]和均方差

图10 SNR=5 dB时ANF在0<ω0≤0.1π的频率估计偏差E[δω(∞)]和均方差

图11 无噪声时ANF在0<ω0≤0.1π的频率估计E[δω(∞)]和

4.2 时变信号频率估计精度分析

图12 ANF时变信号频率估计效果图

保持参数设置不变，待估频率值ω0开始设为0.05π，在第5×104个采样点后变为0.3π。则ANF估计时变信号频率的效果如图12所示。由图12可知，本文方法能够较好地跟踪时变信号，但文献[13]所提方法在信号频率由0.05π变为0.3π时，丧失了跟踪信号频率的能力，这主要是由于当频率发生变化时，误差函数曲面也发生了变化，此时ANF的频率估计值刚好位于变化后误差函数的全局极小值点，而最优频率解恰恰位于远离此处的局部极小值点，从而导致其无法收敛至最优频率解，反而继续收敛至全局极小值点，导致频率估计产生偏差，无法实时跟踪时变信号的频率。

4.3 步长μ对频率估计精度的影响分析

图13 不同μ值下的ANF频率估计偏差E[δω(∞)]和均方差

4.4 极半径ρ对频率估计精度的影响分析

图14 不同ρ值下的ANF频率估计偏差E[δω(∞)]和均方差

4.5 不同信噪比(SNR)对频率估计精度的影响分析

保持参数设置不变，图15是不同信噪比下ANF的频率估计精度。由图15可知，当信噪比较高时，本文方法与文献[13]方法的频率估计精度相当，而当信噪比较低时，本文方法明显优于文献[13]方法。

图15 不同信噪比下的ANF频率估计偏差E[δω(∞)]和均方差

5 结 论

针对ANF极端频率估计存在的问题，提出了一种极端频率估计ANF新方法。通过新误差函数，改善ANF在整个频率范围内的迭代收敛曲面，加快了ANF收敛速度，增强了ANF频率估计的稳定性；通过噪声σ2估计，利用偏差补偿技术，有效降低了ANF频率估计的偏差与均方差，提高了频率估计精度，获得了近似无偏频率估计结果；分析了所提方法的稳态性能。本文方法弥补了ANF对极端频率信号进行频率估计时的缺陷，拓展了ANF进行频率估计的适用范围。研究表明有如下结论：

(1) 算法精度较高。特别针对极端频率信号时，较原方法有明显提高。

(2) 算法抗噪性较好。在信噪比为-10～20 dB时，本文方法的频率估计精度变化较原方法更小，且MSE值基本保持不变。

(3) 算法实现简单。本文算法相对于原算法增加的计算量较少，且为时域递推算法，实时性可以得到满足。

参考文献：

[1] 毛育文，涂亚庆，张海涛，等. 极低频信号的一种离散频谱校正新方法[J]. 振动工程学报，2012，25(4)：474—480.Mao Yuwen, Tu Yaqing, Zhang Haitao, et al. A discrete spectrum correction method for ultra low frequency signals [J]. Journal of Vibration Engineering, 2012, 25(4): 474—480.

[2] 毛育文，涂亚庆，张海涛，等. 计及负频率的极高频信号离散频谱校正新方法[J].振动、测试与诊断，2012，32(3)：477— 482.Mao Yuwen, Tu Yaqing, Zhang Haitao, et al. New discrete spectrum correction method with negative frequency contribution for ultra high frequency signals[J]. Journal of Vibration , Measurement & Diagnosis, 2012, 32(3): 477—482.

[3] 陈奎孚, 王建立, 张森文. 频谱校正的复比值法[J]. 振动工程学报,2008,21(3) :314—318.Chen Kuifu, Wang Jianli, Zhang Senwen. Spectrum correction based on the complex ratio of discrete spectrum around the main-lobe [J]. Journal of Vibration Engineering, 2008, 21(3): 314—318.

[4] 杨志坚, 丁康. 调制FFT 及其在离散频谱校正技术中的应用[J] .振动工程学报, 2009,22(6): 623—637.Yang Zhijian, Ding Kang. Modulated FFT and its application in the discrete spectrum correction[J]. Journal of Vibration Engineering, 2009, 22(6): 623—637.

[5] 毛育文，涂亚庆，肖玮，等. 离散密集频谱细化分析与校正方法研究进展[J]. 振动与冲击, 2012, 31(21): 112—119.Mao Yuwen, Tu Yaqing, Xiao Wei, et al. New discrete spectrum correction method with negative frequency contribution for ultra high frequency signals [J]. Journal of Vibration and Shock, 2012, 31(21)：112— 119.

[6] Prioakis J G, Manolakis D G. Digital Signal Processing: Principle, Algorithms, and Application [M]. New Jersey: Prentice Hall, 2006: 112—114.

[7] 陈向民，于德介，罗洁思. 基于信号共振稀疏分解的包络解调方法及其在轴承故障诊断中的应用[J]. 振动工程学报， 2012, 25(6): 628—636.Chen Xiangmin, Yu Dejie, Luo Jiesi. Envelope demodulation method based on resonance-based sparse signal decomposition and its application in roller bearing fault diagnosis [J]. Journal of Vibration Engineering, 2012, 25(6): 628—636.

[8] Johansson A T, White P R. Instantaneous frequency estimation at low signal-to noise ratios using time-varying notch filters [J]. Signal Processing, 2008, 88(5): 1 271—1 288.

[9] Minh T, Victor D B. Robust notch filtering by combining adaptation in both time and frequency [A]. Conference Record of the Forty-First Asilomar Conference on Signals, Systems and Computers[C]． 2007: 1 633—1 637.

[10] Punchalard R, Lorsawatsiri A, Koseeyaporn J, et al. Adaptive IIR notch filters based on new error criteria[J].Signal Processing, 2008, 88 (3):685—703.

[11] Punchalard R, Koseeyaporn J, Wardkein P. Adaptive IIR notch filter using a modified sign algorithm [J]. Signal Processing, 2009, 89(2):239—243.

[12] Loetwassana W, Punchalard R, Koseeyaporn J, et al. Unbiased plain gradient algorithm for a second-order adaptive IIR notch filter with constrained poles and zeros [J]. Signal Processing, 2010, 90(8): 2 513—2 520.

[13] Punchalard R. Mean square error analysis of unbiased modified plain gradient algorithm for second-order adaptive IIR notch filter [J]. Signal Processing, 2012, 92(11): 2 815—2 820.

[14] Minh T, Hieu T, Victor D B. Stochastic search methods to improve the convergence of adaptive notch filters[A]. IEEE 13th 2009 Digital Signal Processing Workshop and 5th IEEE Signal Processing Education Workshop[C]． DSP/SPE 2009,2009: 78—83．

[15] Zhou J, Li G. Plain gradient-based direct frequency estimation using second-order constrained adaptive IIR notch filter [J]. Electronics Letters, 2004, 40(5):351—352.