頻率無偏估計的有限沖激響應自適應陷波器方法研究

李 明 涂亞慶 萬 平 肖 瑋 陳 鵬

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

1 引言

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

為此,本文以結構較為簡單的FIR-ANF為基礎,利用誤差函數的性能分析結果,推導出頻率估計遞推表達式,提高其頻率估計方法的收斂速度。在分析遞推表達式頻率估計性能的基礎上,提出偏差補償項,有效提高FIR-ANF頻率估計方法的精度和抗噪性能,獲得無偏的頻率估計結果。

2 方法原理

設正弦輸入信號為:

(1)

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

FIR-ANF的傳遞函數為

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

(2)

其中a=-2cosω,ω為陷波頻率,當x(n)通過式(2)所示的ANF后,其值y(n)如式(3)和式(4)所示。當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)

根據文獻[11],則式(3)和式(4)的誤差函數為

J=y2(n)

(5)

使其最小化

(6)

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

(7)

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

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

(8)

即在穩態條件下,存在式(9)的關系

(9)

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

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頻率估計方法的性能,驗證其偏差性,需要求取式(9)的期望,由于其計算時是按照分子分母分別迭代計算,故其期望值為:

(11)

在式(11)中對其中各項分別計算,可得

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

(12)

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

(13)

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

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

(14)

(15)

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

(16)

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

(17)

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

(18)

將式(18)所示的偏差補償項代入式(15)可得,

(19)

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

(20)

(21)

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

(22)

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

3.2 方差分析

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

(23)

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

(24)

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

(25)

其中,

(26)

CT1為協方差矩陣,其值為

(27)

(28)

(29)

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

(30)

(31)

對式(31)求取方差可得

(32)

其中,

(33)

CT2為協方差矩陣,其值為

(34)

(35)

其中,

(36)

(37)

(38)

(39)

4 計算驗證

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

圖1 A=1,ω=0.1π,θ∈(0, 2π),λ=0.9999,ρ=0.95,μ=10-4,a(0)=0,SNR分別為5 dB和20 dB時,獨立運行50次,不同方法的頻率估計結果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

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

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

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

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

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

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

5 結論

論文針對現有自適應陷波器頻率估計方法存在結構復雜、收斂速度慢、抗噪性弱和頻率估計結果有偏的問題,以結構較為簡單的FIR自適應陷波器為基礎,在分析其誤差函數性能的基礎上,提出了一種頻率無偏估計的FIR-ANF頻率估計方法。研究表明有如下結論:

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

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

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