李 明 涂亞慶 萬 平 肖 瑋 陳 鵬

(1. 陸軍勤務(wù)學(xué)院軍事物流系, 重慶 401311; 2. 中國空氣動力研究與發(fā)展中心設(shè)備設(shè)計(jì)及測試技術(shù)研究所, 四川綿陽 621000)

1 引言

自適應(yīng)陷波器(ANF)頻率估計(jì)方法在衛(wèi)星導(dǎo)航定位、工業(yè)控制與儀表以及醫(yī)學(xué)工程中有著較為廣泛的應(yīng)用[1-5]。該頻率估計(jì)方法同基于FFT離散頻譜校正的頻率估計(jì)方法[6-7]相比,具有無頻譜泄漏、計(jì)算簡單的優(yōu)點(diǎn),但也存在易受噪聲影響,頻率估計(jì)結(jié)果有偏的問題。ANF頻率估計(jì)方法按照陷波器傳遞函數(shù)的不同,大致可分為兩類[8],一類是基于有限沖激響應(yīng)自適應(yīng)陷波器(FIR-ANF)的頻率估計(jì)方法,該類方法的結(jié)構(gòu)較為簡單,但易受噪聲干擾,估計(jì)結(jié)果往往存在較大的偏差;另一類是基于無限沖激響應(yīng)自適應(yīng)陷波器(IIR-ANF)的頻率估計(jì)方法,該類方法雖然抗噪性較好,頻率估計(jì)精度有一定的提升,但其結(jié)構(gòu)較為復(fù)雜,且存在收斂速度慢的缺點(diǎn)。針對上述問題,相關(guān)的文獻(xiàn)開展了一些研究。如文獻(xiàn)[9]針對FIR-ANF的頻率估計(jì)問題進(jìn)行了研究,一定程度上提高了頻率估計(jì)精度,具備了一定的抗噪性,但其頻率估計(jì)有偏和收斂速度不快的問題還有待解決。文獻(xiàn)[10]則主要針對IIR-ANF的頻率估計(jì)問題開展研究,雖然提高了頻率估計(jì)精度,但其頻率估計(jì)偏差依然存在,且收斂速度慢的問題仍較為嚴(yán)重,亟待進(jìn)一步改善。

為此,本文以結(jié)構(gòu)較為簡單的FIR-ANF為基礎(chǔ),利用誤差函數(shù)的性能分析結(jié)果,推導(dǎo)出頻率估計(jì)遞推表達(dá)式,提高其頻率估計(jì)方法的收斂速度。在分析遞推表達(dá)式頻率估計(jì)性能的基礎(chǔ)上,提出偏差補(bǔ)償項(xiàng),有效提高FIR-ANF頻率估計(jì)方法的精度和抗噪性能,獲得無偏的頻率估計(jì)結(jié)果。

2 方法原理

設(shè)正弦輸入信號為:

(1)

式中,A為信號幅值,θ為信號相位;ω0為信號頻率,單位rad,f0為信號實(shí)際頻率,單位Hz,fs為采樣頻率;υ0(n)為高斯(正態(tài))白噪聲,其均值為0,方差為σ2。

FIR-ANF的傳遞函數(shù)為

H(z)=1+az-1+z-2

(2)

其中a=-2cosω,ω為陷波頻率,當(dāng)x(n)通過式(2)所示的ANF后,其值y(n)如式(3)和式(4)所示。當(dāng)ANF收斂后,此時a→a0=-2cosω0,即ω→ω0,陷波頻率等于輸入信號頻率。

x(n)→(1+az-1+z-2)→y(n)

(3)

y(n)=x(n)+ax(n-1)+x(n-2)

(4)

根據(jù)文獻(xiàn)[11],則式(3)和式(4)的誤差函數(shù)為

J=y2(n)

(5)

使其最小化

(6)

常見的基于FIR-ANF的頻率估計(jì)是在式(6)的基礎(chǔ)上形成如下的基于梯度下降的頻率迭代遞推計(jì)算式,

(7)

式(7)所示的梯度下降頻率估計(jì)方法,在y(n)中包含有待估計(jì)參數(shù)a值,即梯度下降的速度與估計(jì)參數(shù)a值的關(guān)聯(lián)性較大,會導(dǎo)致整個頻率估計(jì)方法出現(xiàn)收斂速度偏慢的問題。

為此,將式(6)展開后取0,可得

(8)

即在穩(wěn)態(tài)條件下,存在式(9)的關(guān)系

(9)

式(9)中,r1(n)=x(n-1)(x(n)+x(n-2)),r0(n)=x2(n-1),與估計(jì)參數(shù)a完全沒有關(guān)聯(lián),只與輸入信號有關(guān),而輸入信號一般來說相對固定,可以有效的提高整個頻率估計(jì)方法的收斂性,在實(shí)際計(jì)算中,為保證計(jì)算結(jié)果平滑,可以采用式(10)所示的迭代計(jì)算方式。

r1(n)=λr1(n-1)+
(1-λ)[x(n-1)(x(n)+x(n-2))]
r0(n)=λr0(n-1)+(1-λ)x2(n-1)

(10)

式(10)中λ是遺忘因子。

3 無偏性與方差分析

3.1 無偏性分析

為分析所提FIR-ANF頻率估計(jì)方法的性能,驗(yàn)證其偏差性,需要求取式(9)的期望,由于其計(jì)算時是按照分子分母分別迭代計(jì)算,故其期望值為:

(11)

在式(11)中對其中各項(xiàng)分別計(jì)算,可得

E[r1(n)]=E[x(n-1)(x(n)+x(n-2))]

(12)

E[r0(n)]=E[x2(n-1)]

(13)

為計(jì)算式(11),需要分別計(jì)算出E[r1(n)]、E[r0(n)],為此將式(1)代入式(12)～式(13),在計(jì)算過程中注意噪聲及其在不同時刻的自相關(guān)性和互相關(guān)性,可得

E[r1(n)]=E[x(n-1)(x(n)+x(n-2))]=A2cosω0

(14)

(15)

將式(14)～式(15)代入式(11)可得,

(16)

式(16)中,影響頻率估計(jì)均值的關(guān)鍵因素是其中的噪聲σ2,如果在無噪聲的情況下,即σ2=0時,則式(16)可簡化為:

(17)

為獲得無偏的頻率估計(jì)結(jié)果,需要將r0(n)中的σ2消除。由文獻(xiàn)[12]可知,當(dāng)FIR-ANF處于收斂的穩(wěn)定狀態(tài)下時,即a0→-2cosω0,可得

(18)

將式(18)所示的偏差補(bǔ)償項(xiàng)代入式(15)可得,

(19)

由此可將于r0(n)中的σ2消掉,則對式(9)和式(10)所示方法進(jìn)行偏差補(bǔ)償后的頻率估計(jì)方法為

(20)

(21)

對式(20)所示的頻率估計(jì)方法求取期望,可得

(22)

由此可見所提的頻率估計(jì)方法為無偏的,為進(jìn)一步分析所提方法的方差性能,對式(9)和式(20)進(jìn)行方差分析。

3.2 方差分析

針對式(9),為求取其方差,需要按照文獻(xiàn)[13]對其進(jìn)行處理,

(23)

將f1在點(diǎn)T10=E[T1]=μ1=[E[r1(n)]E[r0(n)]]T附近用一階泰勒級數(shù)展開[13],可得

(24)

對式(24)求取方差可得,

(25)

其中,

(26)

CT1為協(xié)方差矩陣,其值為

(27)

(28)

(29)

同理,對式(22)進(jìn)行處理,可得

(30)

(31)

對式(31)求取方差可得

(32)

其中,

(33)

CT2為協(xié)方差矩陣,其值為

(34)

(35)

其中,

(36)

(37)

(38)

(39)

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

為了驗(yàn)證式(9)和式(20)所提頻率估計(jì)方法之間的差異和有效性,在MatLab R2018環(huán)境下,給出不同信號參數(shù)下的頻率估計(jì)結(jié)果。令A(yù)=1,ω=0.1π,θ在區(qū)間(0, 2π)上服從均勻分布,λ=0.9999,頻率初值a(0)=0,信噪比SNR=5 dB和20 dB,獨(dú)立運(yùn)行50次,則式(9)、式(20)和文獻(xiàn)[9]、[10]和[12](參數(shù)設(shè)置ρ=0.95,μ=10-4)所提的頻率估計(jì)結(jié)果如圖1所示。由圖1可知,當(dāng)信噪比較高時,式(9)、式(20)和文獻(xiàn)[9]、[10]和[12]所提方法都可以較準(zhǔn)確的給出頻率估計(jì)結(jié)果,但式(20)方法的收斂性能最佳;當(dāng)信噪比較低時,則只有式(20)和文獻(xiàn)[9]、[10]可以給出較滿意的頻率估計(jì)結(jié)果,而式(9)和文獻(xiàn)[12]所給頻率估計(jì)結(jié)果都已經(jīng)偏離了頻率真值。

圖1 A=1,ω=0.1π,θ∈(0, 2π),λ=0.9999,ρ=0.95,μ=10-4,a(0)=0,SNR分別為5 dB和20 dB時,獨(dú)立運(yùn)行50次,不同方法的頻率估計(jì)結(jié)果Fig.1 A=1, ω=0.1π, θ∈(0, 2π), λ=0.9999, ρ=0.95, μ=10-4, a(0)=0, SNR is 5 dB or 20 dB,frequency estimation results by different methods, 50 runs

為進(jìn)一步分析式(9)和式(20)方法在ω∈(0, π)范圍內(nèi)的頻率估計(jì)性能,其頻率估計(jì)期望與均方差如圖2所示。其中A=1,θ∈(0, 2π),λ=0.9999,a(0)=0,SNR=5 dB,獨(dú)立運(yùn)行50次。由圖2可知,在整個頻率范圍內(nèi),式(20)比式(9)的頻率估計(jì)精度要高,特別是在頻譜的兩端,即ω→0或ω→π時,此時式(20)的頻率估計(jì)精度要優(yōu)于式(9)的頻率估計(jì)精度,但式(20)的方差要略大于式(9)的方差,這是式(20)的不足之處,有待深入研究加以改善。值得注意的是,當(dāng)ω處于頻率的中段,即ω在0.5π附近時,此時式(9)和式(20)頻率估計(jì)精度相當(dāng)。

圖2 A=1,θ∈(0, 2π),λ=0.9999,a(0)=0,SNR=5 dB,ω在(0, π)的范圍內(nèi)時,獨(dú)立運(yùn)行50次的頻率估計(jì)結(jié)果Fig.2 A=1, θ∈(0, 2π), λ=0.9999, a(0)=0, SNR=5 dB,ω in (0, π), frequency estimation results, 50 runs

為進(jìn)一步分析式(9)和式(20)方法在λ∈[0.96, 1)范圍內(nèi)的頻率估計(jì)性能,其頻率估計(jì)期望與均方差如圖3所示。其中A=1,ω=0.1π,θ∈(0, 2π),a(0)=0,SNR=5 dB,獨(dú)立運(yùn)行50次。由圖3可知,式(9)對λ值的選取不敏感,受λ值的影響較小;但式(20)則對λ值的選取非常敏感,存在不穩(wěn)定的現(xiàn)象,其選擇的區(qū)間相對較小,且越接近于1時其頻率估計(jì)精度越高,但過于接近1會導(dǎo)致收斂速度變慢,在實(shí)際應(yīng)用中需綜合考慮。

圖3 A=1,ω=0.1π,θ∈(0, 2π),a(0)=0,SNR=5 dB,λ在[0.96, 1)的范圍內(nèi)時,獨(dú)立運(yùn)行50次的頻率估計(jì)結(jié)果Fig.3 A=1, ω=0.1π, θ∈(0, 2π), a(0)=0, SNR=5 dB,λ in [0.96, 1), frequency estimation results, 50 runs

為進(jìn)一步分析式(9)和式(20)方法在信噪比SNR∈[-10, 20]時的頻率估計(jì)性能,其頻率估計(jì)期望與均方差如圖4所示。其中A=1,ω=0.1π,θ∈(0, 2π),λ=0.9999,a(0)=0,獨(dú)立運(yùn)行50次。由圖4可知,式(20)和式(9)的期望與實(shí)際計(jì)算值吻合較好,但方差的分析實(shí)際值與理論值有一定的差距,需下一步深入研究。式(20)相比式(9)方差要大,這與前面的分析基本一致,但式(20)的頻率估計(jì)方法的估計(jì)準(zhǔn)確度較好,具有無偏的特性。綜上所述,本文所提的FIR-ANF頻率估計(jì)方法具有收斂速度快、抗噪性能好的特點(diǎn),且其頻率估計(jì)結(jié)果無偏。

圖4 A=1,ω=0.1π,θ∈(0, 2π),λ=0.9999,a(0)=0,SNR在[-10, 20]dB的范圍內(nèi)時,獨(dú)立運(yùn)行50次的頻率估計(jì)結(jié)果Fig.4 A=1, ω=0.1π, θ∈(0, 2π), λ=0.9999, a(0)=0, SNR in [-10, 20], frequency estimation results, 50 runs

5 結(jié)論

論文針對現(xiàn)有自適應(yīng)陷波器頻率估計(jì)方法存在結(jié)構(gòu)復(fù)雜、收斂速度慢、抗噪性弱和頻率估計(jì)結(jié)果有偏的問題,以結(jié)構(gòu)較為簡單的FIR自適應(yīng)陷波器為基礎(chǔ),在分析其誤差函數(shù)性能的基礎(chǔ)上,提出了一種頻率無偏估計(jì)的FIR-ANF頻率估計(jì)方法。研究表明有如下結(jié)論:

(1)收斂性能好。在保證頻率估計(jì)精度前提下,所提方法的收斂速度較快,具備較好的實(shí)時性。

(2)抗噪性能好。在不同信噪比條件下,都可以取得滿意的頻率估計(jì)結(jié)果,對噪聲不敏感。

(3)頻率估計(jì)結(jié)果無偏。所提方法可以給出無偏的頻率估計(jì)結(jié)果,且結(jié)構(gòu)簡單,具備較為廣泛的應(yīng)用前景。下一步將針對調(diào)制后的正弦信號頻率估計(jì)開展研究,進(jìn)一步提升所提頻率估計(jì)方法的應(yīng)用范圍。