李 明， 涂亞慶， 沈廷鰲， 毛育文

(后勤工程學(xué)院信息工程系，重慶 401311)

引 言

振動(dòng)工程中常見的極端頻率信號，其頻率估計(jì)精度具有很重要的理論意義和使用價(jià)值[1,2]。在離散頻譜校正理論中的極端頻率，主要相對于頻率分辨率而言(即采樣頻率和采樣時(shí)間長度的比值)，由于極端頻率頻譜泄露嚴(yán)重、存在柵欄效應(yīng)和負(fù)頻率成分干涉，基于離散頻譜校正的頻率估計(jì)方法存在較大誤差[3,4]，且計(jì)算復(fù)雜，抗噪性較差，精度還有待進(jìn)一步提高[5]，也無法估計(jì)時(shí)變信號的頻率值。

為改變此現(xiàn)狀，提高頻率估計(jì)方法的實(shí)時(shí)性、準(zhǔn)確性和抗噪性，采用自適應(yīng)陷波器(ANF)進(jìn)行頻率估計(jì)。ANF主要是根據(jù)被處理信號的特點(diǎn)，對其進(jìn)行參數(shù)優(yōu)化，自動(dòng)調(diào)節(jié)自身模型參數(shù)，令誤差函數(shù)達(dá)到最小值，使陷波頻率與信號頻率相等，從而估計(jì)出信號頻率值。每估計(jì)一次頻率值只需要3個(gè)采樣點(diǎn)，同離散頻譜校正理論中采樣時(shí)間長度的含義完全不同。該方法從時(shí)域角度進(jìn)行迭代頻率估計(jì)，完全避免了負(fù)頻率成分的干擾，不存在頻譜泄露和柵欄效應(yīng)問題，不僅可以估計(jì)頻率恒定的時(shí)不變信號，還可以估計(jì)頻率發(fā)生變化的時(shí)變信號；同時(shí)可以濾除信號噪聲，提高信噪比，針對噪聲環(huán)境下的弱信號也具有較好的檢測效果，且結(jié)構(gòu)簡單、計(jì)算量小[6]。

由于ANF頻率估計(jì)方法同離散頻譜校正理論的差別較大，故本文所指的極端頻率信號，主要指歸一化頻率接近于0或π兩端的信號。例如，當(dāng)采樣頻率設(shè)為4 096 Hz，型號為6307E型的滾動(dòng)軸承，其內(nèi)圈、外圈的故障特征頻率分別為76.5和98.8 Hz[7]，相對于采樣頻率4 096 Hz則為極低頻信號，其歸一化頻率接近于0。此外，當(dāng)現(xiàn)有的采集裝置的采樣頻率限定為1 024 Hz時(shí)，旋轉(zhuǎn)機(jī)械轉(zhuǎn)子系統(tǒng)松動(dòng)故障產(chǎn)生的10倍頻簡諧分量(基頻50 Hz，50×10 Hz=500 Hz，接近奈奎斯特頻率512 Hz)相對于采樣頻率1 024 Hz則為極高頻信號，其歸一化頻率接近于π。

針對上述信號進(jìn)行頻率估計(jì)時(shí)，由于噪聲干擾增大，其頻率估計(jì)精度將出現(xiàn)明顯下降，且存在算法不穩(wěn)定，收斂速度偏慢的缺點(diǎn)[8,9]。針對此問題，文獻(xiàn)[10,11]對ANF誤差函數(shù)進(jìn)行了改進(jìn)，使該誤差函數(shù)具備一定斜直線的特性，一定程度上加快了ANF收斂速度。文獻(xiàn)[12,13]針對上述文獻(xiàn)結(jié)果有偏的問題進(jìn)行了研究，取得了近似無偏的頻率估計(jì)結(jié)果，但上述方法的缺點(diǎn)是最優(yōu)頻率估計(jì)解為局部極小值點(diǎn)，對頻率初始值的設(shè)定有一定的要求，且為間接頻率估計(jì)，對估計(jì)后參數(shù)的取值有一定的限制；同時(shí)其頻率估計(jì)結(jié)果存在振蕩，針對極端頻率信號的估計(jì)精度和穩(wěn)定性還有待進(jìn)一步提升。

為此，本文為解決現(xiàn)有ANF方法估計(jì)極端頻率信號時(shí)存在收斂速度慢、精度不高、穩(wěn)定性差的問題，通過提出一種新誤差函數(shù)，改善ANF在整個(gè)頻率范圍內(nèi)的迭代收斂曲面，加快ANF收斂速度與算法穩(wěn)定性；同時(shí)，采用偏差補(bǔ)償方式，降低噪聲對ANF的影響，提供近似無偏的極端頻率信號頻率估計(jì)結(jié)果。

綜上所述，本文提出了一種ANF極端頻率信號直接頻率估計(jì)新方法，并與原有算法進(jìn)行了有關(guān)性能的比較分析。

1 極端頻率信號

設(shè)正弦信號為

(1)

根據(jù)采樣定理，fs≥2f0，故可知信號頻率的取值范圍為0<ω0≤π；同時(shí)現(xiàn)有的ANF研究表明[14]，當(dāng)信號頻率0<ω0≤0.1π或0.9π<ω0≤π時(shí)，ANF的性能出現(xiàn)明顯下降，所以本文討論的極頻信號主要指頻率0<ω0≤0.1π的極低頻信號或0.9π<ω0≤π的極高頻信號。

2 新自適應(yīng)陷波器方法

2.1 新誤差函數(shù)

ANF傳遞函數(shù)為

(2)

式中ρ為極半徑，控制陷波寬度，且0?ρ<1。ANF參數(shù)a=-2cosω，ω為ANF陷波頻率，在頻率估計(jì)過程中，ω將趨近于ω0，故ω可看作為ω0的估計(jì)值，此時(shí)，a→a0=-2cosω0。

在頻率估計(jì)過程中，若估計(jì)a值，則為間接頻率估計(jì)，此時(shí)需保證-2≤a≤2，以使ω=arccos(-a/2)有解。但若估計(jì)ω，則為直接頻率估計(jì)，不需保證ω的取值范圍，步驟較為簡化。

為此，將a=-2cosω代入式(2)可得

(3)

式中N(z,ω)和H(z,ω)分別為ANF的FIR和IIR結(jié)構(gòu)。則信號x(k)通過N(z,ω)和H(z,ω)的信號e1(k)和e2(k)分別為

(4)

為了獲得最佳頻率估計(jì)值，ANF通過調(diào)整陷波頻率ω，使下式的新誤差函數(shù)最小化

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

(5)

式中E[·]表示求取期望，其保留了傳統(tǒng)IIR自適應(yīng)陷波器頻率估計(jì)精度高，在最優(yōu)頻率值附近下降速度快的優(yōu)點(diǎn)；同時(shí)使誤差函數(shù)在遠(yuǎn)離最優(yōu)頻率值時(shí)仍能夠保持一定的梯度值，使其可保持較快的收斂速度。在實(shí)際計(jì)算中可按下式進(jìn)行計(jì)算

(6)

由于輸入信號為正弦信號，故式(3)中N(z,ω)和H(z,ω)在輸入信號頻率ω0處的幅值為

(7)

相角為

(8)

結(jié)合式(4)，(7)和(8)，可得

(9)

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

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

(10)

為說明所提新誤差函數(shù)的優(yōu)越性，故比較文獻(xiàn)[10]所采用的如下式所示的原有誤差函數(shù)為

(11)

實(shí)際計(jì)算時(shí)按下式計(jì)算

(12)

為說明本文所提如式(5)所示新誤差函數(shù)的優(yōu)越性，故所取參數(shù)同文獻(xiàn)[10]保持一致，A=1，θ=π/6，ρ=0.95，σ2=0，輸入信號頻率分別取3個(gè)不同的值，π/3，π/2和2π/3，則文獻(xiàn)[10]和本文方法所提誤差函數(shù)如圖1所示。又由于極高頻和極低頻具有對稱性，故只分析頻率ω0為0.01π時(shí)的極低頻信號，其他參數(shù)設(shè)置保持不變，則極端頻率下相應(yīng)的誤差函數(shù)如圖2所示。

圖1 不同方法所提誤差函數(shù)的理論值與計(jì)算值

圖2 ω0=0.01π時(shí)不同方法所提誤差函數(shù)

由圖1，2可知，兩種不同誤差函數(shù)的理論計(jì)算值同實(shí)際計(jì)算值基本保持一致。同時(shí)，本文所提誤差函數(shù)J(ω)相較J0(ω)的梯度值較大，在迭代計(jì)算過程中的收斂速度將更快，且最優(yōu)頻率值的選取較為明顯。特別當(dāng)輸入信號頻率為極低頻，ω0=0.01π時(shí)，文獻(xiàn)[10]所提誤差函數(shù)在這一區(qū)域附近基本為平直線，喪失了對最優(yōu)解的收斂性能，且容易在最優(yōu)解附近產(chǎn)生收斂振蕩，導(dǎo)致其算法不穩(wěn)定。而本文所提誤差函數(shù)仍然具有較大的梯度值，仍可收斂至最優(yōu)頻率值，提升了算法的穩(wěn)定性。

圖3 ω0=π/3時(shí)本文方法所提誤差函數(shù)曲面圖

圖4 ω0=π/3時(shí)原有誤差函數(shù)曲面圖

圖5 ω0=0.01π時(shí)本文方法所提誤差函數(shù)曲面圖

圖6 ω0=0.01π時(shí)原有誤差函數(shù)曲面圖

2.2 頻率估計(jì)偏差問題分析

結(jié)合式(5)所示新誤差函數(shù)，可知ANF直接頻率估計(jì)方法為

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

(13)

式(13)中的e1(k)g2(k)相對于其他項(xiàng)較小，為簡化計(jì)算，可略去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)所示頻率估計(jì)方法結(jié)果有偏，為消除該偏差，對式(14)進(jìn)行穩(wěn)態(tài)條件下的偏差分析，此時(shí)ω→ω0，故可得[15]

(16)

穩(wěn)態(tài)下，由于ω→ω0，故可用ω0替代式(15)中的ω(k)，由此可得e1(k)，e2(k)，g1(k)和g2(k)的表達(dá)式

(17)

其中，δω(k)=ω(k)-ω0；v1(k)，v2(k)，v3(k)和v3(k)為噪聲v0(k)分別通過N(z,ω)，G1(z,ω)，H(z,ω)和G2(z,ω)后所產(chǎn)生。

由此，將式(14)的兩邊同時(shí)減去ω0，可得

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

(18)

對式(18)兩邊同取期望，并將式(17)代入。同時(shí)，令Ri,j=Rj,i=E[vi(k)vj(k)]表示噪聲vi(k)和vj(k)的相關(guān)值，計(jì)算時(shí)由于ρ→1，為簡化計(jì)算，設(shè)(1-ρ)2≈0，則可得

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

(19)

由式(19)可知，其右邊最后一項(xiàng)，即輸出信號間噪聲的相關(guān)性，是導(dǎo)致如式(14)所示頻率估計(jì)算法結(jié)果有偏的根本原因。同時(shí)，該項(xiàng)是輸入信號噪聲σ2、極半徑ρ和頻率估計(jì)結(jié)果ω的函數(shù)，且無論如何調(diào)整上述參數(shù)，都無法使該項(xiàng)為0。為此，需消除此項(xiàng)的影響，提高其頻率估計(jì)精度。

2.3 無偏頻率估計(jì)方法

為得到無偏頻率估計(jì)結(jié)果，需對輸入信號噪聲σ2進(jìn)行估計(jì)。為此，令c(k)=4(ρ-1)sin2ω(k)，分析可知

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

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

(20)

式(20)中右邊最后一項(xiàng)包含了2R3,4+R2,3的值，于是，將式(20)中的c(k)x(k)e1(k)代入式(18)，兩邊同取期望。雖然引入了多余項(xiàng)c(k)A2sinω0cosφ1·E[δω(k)]，但是該項(xiàng)可與式(19)中的M1合并，對頻率估計(jì)精度的影響較小，由此可知基于新誤差函數(shù)的ANF極端頻率無偏直接估計(jì)新方法為

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

(21)

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

該算法流程為：

Step1. 設(shè)置頻率初始值ω(0)，算法步長μ和ANF參數(shù)ρ；

Step2. 在k時(shí)刻，將信號x(k)分別通過N(z,ω)和H(z,ω)，得到e1(k)和e2(k)，并計(jì)算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. 計(jì)算k時(shí)刻的信號g1(k)和g2(k)，同時(shí)計(jì)算G(k)的值；

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

Step4. 按式(21)更新頻率估計(jì)值，得到

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

Step5. 重復(fù)步驟Step2～Step4，直至算法收斂。

3 性能分析

3.1 穩(wěn)態(tài)偏差分析

對式(21)進(jìn)行穩(wěn)態(tài)條件下的偏差分析，式(21)兩邊同減去ω0并取期望，得

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

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

(22)

穩(wěn)態(tài)條件下，可知

(23)

式(22)可變化為

(24)

3.2 穩(wěn)態(tài)MSE分析

對式(21)兩邊同減去ω0并平方，求取期望可得

(25)

式(25)左邊第二項(xiàng)的計(jì)算可利用文獻(xiàn)[13]所用方法，得

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

(26)

則結(jié)合式(26)，將式(25)展開，經(jīng)過一系列復(fù)雜繁瑣的推導(dǎo)后，可得

(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。

結(jié)合式(23)可知，穩(wěn)態(tài)下式(27)可化為

(28)

4 計(jì)算驗(yàn)證

4.1 頻率估計(jì)精度分析

圖7 ω0=0.01π時(shí)ANF頻率估計(jì)效果圖

圖8 ANF頻率估計(jì)偏差E[δω(k)]和均方差

圖9 ANF在全頻段頻率估計(jì)偏差E[δω(∞)]和均方差

圖10 SNR=5 dB時(shí)ANF在0<ω0≤0.1π的頻率估計(jì)偏差E[δω(∞)]和均方差

圖11 無噪聲時(shí)ANF在0<ω0≤0.1π的頻率估計(jì)E[δω(∞)]和

4.2 時(shí)變信號頻率估計(jì)精度分析

圖12 ANF時(shí)變信號頻率估計(jì)效果圖

保持參數(shù)設(shè)置不變，待估頻率值ω0開始設(shè)為0.05π，在第5×104個(gè)采樣點(diǎn)后變?yōu)?.3π。則ANF估計(jì)時(shí)變信號頻率的效果如圖12所示。由圖12可知，本文方法能夠較好地跟蹤時(shí)變信號，但文獻(xiàn)[13]所提方法在信號頻率由0.05π變?yōu)?.3π時(shí)，喪失了跟蹤信號頻率的能力，這主要是由于當(dāng)頻率發(fā)生變化時(shí)，誤差函數(shù)曲面也發(fā)生了變化，此時(shí)ANF的頻率估計(jì)值剛好位于變化后誤差函數(shù)的全局極小值點(diǎn)，而最優(yōu)頻率解恰恰位于遠(yuǎn)離此處的局部極小值點(diǎn)，從而導(dǎo)致其無法收斂至最優(yōu)頻率解，反而繼續(xù)收斂至全局極小值點(diǎn)，導(dǎo)致頻率估計(jì)產(chǎn)生偏差，無法實(shí)時(shí)跟蹤時(shí)變信號的頻率。

4.3 步長μ對頻率估計(jì)精度的影響分析

圖13 不同μ值下的ANF頻率估計(jì)偏差E[δω(∞)]和均方差

4.4 極半徑ρ對頻率估計(jì)精度的影響分析

圖14 不同ρ值下的ANF頻率估計(jì)偏差E[δω(∞)]和均方差

4.5 不同信噪比(SNR)對頻率估計(jì)精度的影響分析

保持參數(shù)設(shè)置不變，圖15是不同信噪比下ANF的頻率估計(jì)精度。由圖15可知，當(dāng)信噪比較高時(shí)，本文方法與文獻(xiàn)[13]方法的頻率估計(jì)精度相當(dāng)，而當(dāng)信噪比較低時(shí)，本文方法明顯優(yōu)于文獻(xiàn)[13]方法。

圖15 不同信噪比下的ANF頻率估計(jì)偏差E[δω(∞)]和均方差

5 結(jié) 論

針對ANF極端頻率估計(jì)存在的問題，提出了一種極端頻率估計(jì)ANF新方法。通過新誤差函數(shù)，改善ANF在整個(gè)頻率范圍內(nèi)的迭代收斂曲面，加快了ANF收斂速度，增強(qiáng)了ANF頻率估計(jì)的穩(wěn)定性；通過噪聲σ2估計(jì)，利用偏差補(bǔ)償技術(shù)，有效降低了ANF頻率估計(jì)的偏差與均方差，提高了頻率估計(jì)精度，獲得了近似無偏頻率估計(jì)結(jié)果；分析了所提方法的穩(wěn)態(tài)性能。本文方法彌補(bǔ)了ANF對極端頻率信號進(jìn)行頻率估計(jì)時(shí)的缺陷，拓展了ANF進(jìn)行頻率估計(jì)的適用范圍。研究表明有如下結(jié)論：

(1) 算法精度較高。特別針對極端頻率信號時(shí)，較原方法有明顯提高。

(2) 算法抗噪性較好。在信噪比為-10～20 dB時(shí)，本文方法的頻率估計(jì)精度變化較原方法更小，且MSE值基本保持不變。

(3) 算法實(shí)現(xiàn)簡單。本文算法相對于原算法增加的計(jì)算量較少，且為時(shí)域遞推算法，實(shí)時(shí)性可以得到滿足。

參考文獻(xiàn)：

[1] 毛育文，涂亞慶，張海濤，等. 極低頻信號的一種離散頻譜校正新方法[J]. 振動(dòng)工程學(xué)報(bào)，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ì)及負(fù)頻率的極高頻信號離散頻譜校正新方法[J].振動(dòng)、測試與診斷，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] 陳奎孚, 王建立, 張森文. 頻譜校正的復(fù)比值法[J]. 振動(dòng)工程學(xué)報(bào),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] 楊志堅(jiān), 丁康. 調(diào)制FFT 及其在離散頻譜校正技術(shù)中的應(yīng)用[J] .振動(dòng)工程學(xué)報(bào), 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] 毛育文，涂亞慶，肖瑋，等. 離散密集頻譜細(xì)化分析與校正方法研究進(jìn)展[J]. 振動(dòng)與沖擊, 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] 陳向民，于德介，羅潔思. 基于信號共振稀疏分解的包絡(luò)解調(diào)方法及其在軸承故障診斷中的應(yīng)用[J]. 振動(dòng)工程學(xué)報(bào)， 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.