基于FFT和DTFT的指數(shù)衰減復(fù)正弦信號參數(shù)估計(jì)算法

王 鴻 劉春華 齊國清

(大連海事大學(xué)信息科學(xué)技術(shù)學(xué)院 遼寧 大連 116026)

0 引 言

噪聲背景中的指數(shù)衰減復(fù)正弦信號參數(shù)估計(jì)技術(shù)是數(shù)字信號處理領(lǐng)域中的一個(gè)經(jīng)典課題，廣泛應(yīng)用于低頻機(jī)械光譜學(xué)[1]、線性系統(tǒng)識別[2]、核磁共振波譜分析[3]等領(lǐng)域。為提高該類信號頻率和衰減因子的估計(jì)精度，國內(nèi)外眾多學(xué)者做了大量研究，主要分為頻域和時(shí)域兩方面[4]。許多時(shí)域算法都以奇異值分解為基礎(chǔ)，雖然參數(shù)估計(jì)精度較高，但計(jì)算量大，效率不高[5]。為提高指數(shù)衰減復(fù)正弦信號參數(shù)估計(jì)算法的運(yùn)算速度及抗噪聲性能，本文從頻域角度進(jìn)行研究。

由于FFT得到的是離散頻率值，信號頻率可表示為f0=(m+δ)Δf。其中：m是幅度最大譜線對應(yīng)的離散頻率索引值;δ是信號頻率與幅度最大譜線位置的相對頻率偏差，且δ∈[-0.5,0.5];Δf為頻率分辨率[6]?；贔FT插值的指數(shù)衰減復(fù)正弦信號參數(shù)估計(jì)算法可分為粗估計(jì)和精估計(jì)兩個(gè)步驟，而各算法的差異僅體現(xiàn)在精估計(jì)中[7]。文獻(xiàn)[8]提出了利用信號FFT離散頻譜主瓣內(nèi)幅度最大譜線和次大譜線得到估計(jì)頻率和衰減因子的算法，該算法簡單、計(jì)算量小，但當(dāng)信號實(shí)際頻率與幅度最大譜線對應(yīng)的頻率接近(δ接近0)時(shí)，次大譜線幅度很小、易受噪聲影響導(dǎo)致頻率和衰減因子估計(jì)的誤差增大。文獻(xiàn)[3]提出Bertocco-Yoshida一階離散傅里葉插值算法(BY1)，該算法使用了FFT離散頻譜中最大的三根譜線，在低信噪比和FFT點(diǎn)數(shù)較少時(shí)具有一定優(yōu)勢，但整體估計(jì)精度不高。利用FFT得到的信號x(n)的離散頻譜X(k)只在離散頻率索引值k取整數(shù)值時(shí)有值，文獻(xiàn)[9]將k推廣到非整數(shù)值的情況，提出了一種Aboutanios-Mulgrew算法(A&M)。k=m時(shí)，X(m)為幅度最大譜線，該算法利用x(n)的DTFT(連續(xù)時(shí)間傅里葉變換)抽樣得到X(m)位置+0.5Δf和-0.5Δf處的兩根離散譜線來估計(jì)頻率和衰減因子。文獻(xiàn)[10]提出了一種基于指數(shù)衰減窗的改進(jìn)A&M算法，但利用原信號衰減因子的估計(jì)值選擇最優(yōu)窗函數(shù)e-γn/fs的系數(shù)γ的公式是由大量實(shí)驗(yàn)數(shù)據(jù)擬合得到，難以適用于所有的應(yīng)用場合。文獻(xiàn)[11]利用文獻(xiàn)[10]算法的原理估計(jì)多頻率指數(shù)衰減復(fù)正弦信號參數(shù)，但存在相同的缺點(diǎn)。針對指數(shù)衰減復(fù)正弦信號的參數(shù)估計(jì)問題，A&M算法在現(xiàn)有算法中性能最好，但該算法只討論了當(dāng)i=0.5這一特定值時(shí)，原信號的DTFT在X(m)位置+iΔf和-iΔf處的兩根離散譜線對信號參數(shù)估計(jì)性能的影響。對于恒定幅度正弦信號，文獻(xiàn)[12]提出了一種結(jié)合FFT和DTFT利用X(m)及其+iΔf和-iΔf位置處的兩根離散譜線估計(jì)信號頻率的算法。為盡量避免利用幅度很小的離散頻譜，文中約定i∈(0,1]，并發(fā)現(xiàn)當(dāng)i<0.5時(shí)，算法性能接近CRLB，且隨i的增加，性能基本不變；當(dāng)i≥0.5時(shí)，隨著i的增加，性能逐漸變差。與利用信號補(bǔ)零達(dá)到相同效果的FFT類算法[13-14]相比，文獻(xiàn)[12]算法的計(jì)算量與i的取值無關(guān)；而信號補(bǔ)零后長度與i的倒數(shù)成正比，當(dāng)i較小時(shí)，補(bǔ)零類算法的計(jì)算量是難以接受的。對于恒定幅度正弦信號頻率的估計(jì)，文獻(xiàn)[12]提出的方法可達(dá)到最優(yōu)效果。但是，對于指數(shù)衰減復(fù)正弦信號，通過實(shí)驗(yàn)發(fā)現(xiàn)，并不能直接應(yīng)用文獻(xiàn)[12]算法進(jìn)行參數(shù)估計(jì)，尤其當(dāng)衰減因子較大時(shí)估計(jì)性能很差。因此，借鑒文獻(xiàn)[12]算法的思路，本文提出一種結(jié)合FFT和DTFT估計(jì)指數(shù)衰減復(fù)正弦信號頻率和衰減因子的三譜線插值算法，分析了i對算法性能的影響并給出了選擇i的建議。

1 算法原理

加性高斯白噪聲背景下，觀測信號序列表示為：

x(n)=s(n)+u(n)

(1)

s(n)=b0e-αn/fsexp[j(2πf0n/fs+θ0)]

(2)

式中：n=0,1,…,N-1;s(n)為指數(shù)衰減復(fù)正弦信號;f0、α、b0、θ0分別為頻率、衰減因子、幅度和初始相位；u(n)為復(fù)高斯白噪聲，均值為0,方差為σ2；fs為采樣頻率。

對s(n)做N點(diǎn)FFT，搜索幅度最大譜線的位置m，F(xiàn)FT幅度最大譜線可表示為：

(3)

s(n)經(jīng)DTFT后，計(jì)算s(n)的DTFT在(m+i)Δf和(m-i)Δf處離散抽樣，并用S±i表示，約定i∈(0,1]。

S±i=S(ejω)|ω0=DTFT[s(n)]|ω0=

(4)

式中：ω=2πf/fs為圓周頻率;ω0為抽樣點(diǎn)處的離散頻率：

ω0=2π(m±i)Δf/fs=2π(m±i)/N

根據(jù)式(4)可得Si和S-i的表達(dá)式為：

(5)

(6)

用復(fù)數(shù)變量Z表示相對頻率偏差δ和衰減因子α，即Z=e-α/fsej2πδ/N，并計(jì)算以下比值：

(7)

(8)

聯(lián)立式(7)、式(8)，得到Z與S0、Si、S-i關(guān)系為：

(9)

根據(jù)Z的定義，可得δ和α的估計(jì)公式為：

(10)

(11)

式中：argZ為取Z的相角；ln|Z|為取|Z|的幅值。因此，信號頻率f0的估計(jì)公式可表示為：

(12)

2 算法精度分析

在噪聲環(huán)境中，x(n)的N點(diǎn)FFT為：

X(k)=S(k)+U(k)k=0,1,…,N-1

(13)

當(dāng)k=m時(shí)，幅度最大譜線為X(m)，記為X0=S0+U0；x(n)的DTFT抽樣后得到X(m+i)和X(m-i)，分別記為Xi=Si+Ui和X-i=S-i+U-i。

2.1 噪聲系數(shù)的統(tǒng)計(jì)性質(zhì)

復(fù)高斯白噪聲頻域系數(shù)U0、Ui、U-i的方差皆為Nσ2，但當(dāng)i為非整數(shù)時(shí)，U0、Ui、U-i兩兩間的協(xié)方差不再為0[10]。按協(xié)方差定義整理得到公式如下：

(14)

(15)

同理，可得到其他兩兩系數(shù)間的協(xié)方差公式。

2.2 頻率和衰減因子估計(jì)精度分析

令a=e-j2πi-1，b=j2sin2πi，c=ej2πi-1，d=ej2πi/N，e1=e-j2πi/N。式(9)可重寫為：

(16)

(17)

替換相關(guān)變量后，可得：

(18)

對式(18)進(jìn)行一階泰勒級數(shù)展開并省去高次項(xiàng)：

(19)

(20)

ln(Z)=-α/fs+j2πδ/N

(21)

(22)

(23)

聯(lián)立式(22)、式(23)整理得到：

(24)

(25)

利用式(14)和式(15)推出E[HH*]的表達(dá)式為：

E[HH*]=(|a(1-Zd)|2+|b(1-Z)|2+

|c(1-Ze1)|2)Nσ2+

(26)

根據(jù)ξ和θ的表達(dá)式可推出：

(27)

(28)

3 迭代策略的應(yīng)用

3.1 非迭代情況下的算法性能分析

本節(jié)分析i分別為0.1、0.3、0.5、0.7、0.9時(shí)，f0和α估計(jì)的歸一化均方根誤差隨相對頻率偏差δ變化的情況。參數(shù)設(shè)為：α=5，fs=51.2 kHz，N=512，f0=(N/4+δ)fs/N，則f0∈[12.75,12.85] kHz，θ0隨機(jī)變化，SNR=37 dB，結(jié)果如圖1和圖2所示。

圖1 f0估計(jì)的歸一化均方根誤差隨δ變化的仿真結(jié)果

圖2 α估計(jì)的歸一化均方根誤差隨δ變化的仿真結(jié)果

仿真結(jié)果表明，當(dāng)i一定時(shí)，f0和α估計(jì)的歸一化均方根誤差隨|δ|增加而增加，故可借鑒文獻(xiàn)[12]的迭代策略，提高參數(shù)的估計(jì)精度。

3.2 迭代步驟

從圖1和圖2中也可以看出當(dāng)0

2) 對x(n)進(jìn)行FFT變換，找到m，并得到X0。

3) 信號x(n)經(jīng)DTFT后，在連續(xù)頻譜位置m+p處離散抽樣，并計(jì)算X0.1和X-0.1。

(29)

(30)

4 仿真實(shí)驗(yàn)

為驗(yàn)證本文算法性能，本節(jié)通過計(jì)算機(jī)Monte Carlo模擬實(shí)驗(yàn)進(jìn)行仿真。

文獻(xiàn)[12]算法、BY1算法和A&M算法估計(jì)δ和α的表達(dá)式如下所示(信號模型以式(1)為準(zhǔn)，BY1算法和A&M算法的公式有細(xì)微調(diào)整)。

文獻(xiàn)[12]算法:

BY1算法[3]：

A&M的迭代算法[9]:

并依次計(jì)算：

(1) 仿真實(shí)驗(yàn)1。圖3給出了本文算法和文獻(xiàn)[12]算法的f0估計(jì)性能隨α變化的關(guān)系曲線。參數(shù)設(shè)為：fs=51.2 kHz，N=512，i=0.2，δ=0.2，f0=(N/32+δ)fs/N，則f0=1.62 kHz，SNR=37 dB，θ0隨機(jī)變化，α在[0,20]內(nèi)變化，步進(jìn)1，迭代次數(shù)均為1。

圖3 f0估計(jì)的歸一化均方根誤差隨α變化的仿真結(jié)果對比

仿真結(jié)果表明，α較大時(shí)，直接應(yīng)用文獻(xiàn)[12]算法估計(jì)指數(shù)衰減正弦信號參數(shù)的性能很差，而本文算法是基于指數(shù)衰減的信號模型進(jìn)行推導(dǎo)的，即使在α較高的情況下，算法性能幾乎不受影響。

(2) 仿真實(shí)驗(yàn)2。圖4和圖5給出了δ分別為0.05、0.15、0.25、0.35、0.45以及-0.05、-0.15、-0.25、-0.35、-0.45時(shí)，本文算法的f0和α估計(jì)性能隨i變化的關(guān)系曲線。參數(shù)設(shè)為：α=5，fs=51.2 kHz，N=512，SNR=37 dB，f0=(N/4+δ)fs/N，則f0∈[12.75,12.85] kHz，θ0隨機(jī)變化，i在[0.01,1]內(nèi)變化，步進(jìn)0.01，迭代次數(shù)為1。

圖4 f0估計(jì)的歸一化均方根誤差隨i變化的仿真結(jié)果

圖5 α估計(jì)的歸一化均方根誤差隨i變化的仿真結(jié)果

可以看出f0和α估計(jì)的理論曲線和仿真曲線吻合，歸一化均方根誤差隨|δ|的增大而增大。若δ為定值，當(dāng)i<0.5時(shí)，f0和α估計(jì)的歸一化均方根誤差基本不變；當(dāng)i≥0.5時(shí)，歸一化均方根誤差隨i的增加而增加。原因是當(dāng)i<0.5時(shí)，S0的幅值大于Si、S-i的幅值，且Si、S-i的幅值相對較大，不易受噪聲影響。因此，在實(shí)際應(yīng)用中，輔助譜線應(yīng)當(dāng)選擇位于幅度最大譜線位置+(i<0.5)Δf和-(i<0.5)Δf處的兩根離散譜線。

(3) 仿真實(shí)驗(yàn)3。參數(shù)設(shè)置為：α=5，fs=51.2 kHz，N=512，θ0隨機(jī)，f0=(N/4+δ)fs/N，則f0∈[12.75,12.85] kHz，SNR=37 dB。

圖6和圖7給出了本文算法、BY1算法和A&M算法的f0和α估計(jì)性能隨δ變化的關(guān)系曲線。當(dāng)i=0.5時(shí)，根據(jù)式(9)可知，本文算法只用了S0.5、S-0.5兩根離散譜線，從圖中可以看出，本文算法(i=0.5)性能和A&M算法基本一致。當(dāng)i=1時(shí)，本文算法采用了S0、S1、S-1三根離散譜線，本文算法(i=1)性能與BY1算法基本一致。這里i的取值只是為了方便與另外兩種算法進(jìn)行比較。

圖6 三種算法的f0估計(jì)性能隨δ變化的仿真結(jié)果對比

圖7 三種算法的α估計(jì)性能隨δ變化的仿真結(jié)果對比

圖8和圖9給出了本文算法(i=0.1)和A&M算法的f0和α估計(jì)性能隨δ變化的關(guān)系曲線。實(shí)際上，當(dāng)0

圖8 兩種算法的f0估計(jì)性能隨δ變化的仿真結(jié)果對比

圖9 兩種算法的α估計(jì)性能隨δ變化的仿真結(jié)果對比

(4) 仿真實(shí)驗(yàn)4。實(shí)驗(yàn)4研究了不同信噪比下，本文算法、BY1算法和A&M算法的f0和α估計(jì)性能。參數(shù)設(shè)為：α=5，N=512，θ0隨機(jī)，fs=51.2 kHz，δ=0.2，i=0.1，f0=(N/4+δ)fs/N，則f0=12.82 kHz，信噪比步進(jìn)1 dB。

實(shí)驗(yàn)結(jié)果如表1和表2所示。在表1中，本文算法和A&M算法迭代1次；在表2中，本文算法和A&M算法迭代2次。實(shí)驗(yàn)結(jié)果表明，在不同信噪比條件下，本文算法頻率和衰減因子估計(jì)的歸一化均方根誤差比A&M算法和BY1算法的更?。唤?jīng)過兩次迭代后，本文算法的f0和α估計(jì)性能更接近CRLB。

表1 δ=0.2時(shí)三種算法的f0和α估計(jì)性能對比(迭代1次)

表2 δ=0.2時(shí)兩種算法的f0和α估計(jì)性能對比(迭代2次)

5 結(jié) 語

本文探討了加性高斯白噪聲環(huán)境中，指數(shù)衰減復(fù)正弦信號頻率和衰減因子的估計(jì)問題，提出一種結(jié)合FFT和DTFT的三譜線插值估計(jì)算法。算法將輔助譜線的選擇推廣到了幅度最大譜線位置+iΔf和-iΔf處兩根離散譜線；推導(dǎo)了頻率和衰減因子的估計(jì)公式和理論均方根誤差公式；分析了i取值對算法性能的影響。仿真結(jié)果表明，當(dāng)i=0.5時(shí)，本文算法與A&M算法性能基本一致；當(dāng)i=1時(shí)，本文算法(迭代1次)與BY1算法性能基本一致；當(dāng)i的取值小于0.5時(shí)，本文算法性能優(yōu)于現(xiàn)有性能最好的方法，接近克拉美羅下界。在實(shí)際應(yīng)用中，為取得更好的參數(shù)估計(jì)性能，應(yīng)選擇最大幅度譜線位置+(i<0.5)Δf和-(i<0.5)Δf處的兩根輔助譜線參與算法估計(jì)。