基于FFT和DTFT的指數衰減復正弦信號參數估計算法

王 鴻 劉春華 齊國清

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

0 引 言

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

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

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

(3)

s(n)經DTFT后，計算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為抽樣點處的離散頻率：

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

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

(5)

(6)

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

(7)

(8)

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

(9)

根據Z的定義，可得δ和α的估計公式為：

(10)

(11)

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

(12)

2 算法精度分析

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

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

(13)

當k=m時，幅度最大譜線為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 噪聲系數的統計性質

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

(14)

(15)

同理，可得到其他兩兩系數間的協方差公式。

2.2 頻率和衰減因子估計精度分析

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

(16)

(17)

替換相關變量后，可得：

(18)

對式(18)進行一階泰勒級數展開并省去高次項：

(19)

(20)

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

(21)

(22)

(23)

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

(24)

(25)

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

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

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

(26)

根據ξ和θ的表達式可推出：

(27)

(28)

3 迭代策略的應用

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

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

圖1 f0估計的歸一化均方根誤差隨δ變化的仿真結果

圖2 α估計的歸一化均方根誤差隨δ變化的仿真結果

仿真結果表明，當i一定時，f0和α估計的歸一化均方根誤差隨|δ|增加而增加，故可借鑒文獻[12]的迭代策略，提高參數的估計精度。

3.2 迭代步驟

從圖1和圖2中也可以看出當0

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

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

(29)

(30)

4 仿真實驗

為驗證本文算法性能，本節通過計算機Monte Carlo模擬實驗進行仿真。

文獻[12]算法、BY1算法和A&M算法估計δ和α的表達式如下所示(信號模型以式(1)為準，BY1算法和A&M算法的公式有細微調整)。

文獻[12]算法:

BY1算法[3]：

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

并依次計算：

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

圖3 f0估計的歸一化均方根誤差隨α變化的仿真結果對比

仿真結果表明，α較大時，直接應用文獻[12]算法估計指數衰減正弦信號參數的性能很差，而本文算法是基于指數衰減的信號模型進行推導的，即使在α較高的情況下，算法性能幾乎不受影響。

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

圖4 f0估計的歸一化均方根誤差隨i變化的仿真結果

圖5 α估計的歸一化均方根誤差隨i變化的仿真結果

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

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

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

圖6 三種算法的f0估計性能隨δ變化的仿真結果對比

圖7 三種算法的α估計性能隨δ變化的仿真結果對比

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

圖8 兩種算法的f0估計性能隨δ變化的仿真結果對比

圖9 兩種算法的α估計性能隨δ變化的仿真結果對比

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

實驗結果如表1和表2所示。在表1中，本文算法和A&M算法迭代1次；在表2中，本文算法和A&M算法迭代2次。實驗結果表明，在不同信噪比條件下，本文算法頻率和衰減因子估計的歸一化均方根誤差比A&M算法和BY1算法的更小；經過兩次迭代后，本文算法的f0和α估計性能更接近CRLB。

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

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

5 結 語

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