基于期望最大化算法的雙音頻率估計方法研究*

袁 恩，劉 鵬，陳隆亮，張文宇，常 青

(陸軍工程大學，江蘇 南京 210007)

0 引言

雙音信號在多頻連續波測距雷達[1]、無線傳感器網絡定位[2-3]、無線電信號導航[4]等方面具有廣泛應用。文獻[4]針對當前衛星導航信號易受干擾的問題，提出一種應用于航空導航的抗干擾測距導航方法。該方法在高速跳頻信號的每一跳發送不同的雙音信號，構建時頻矩陣用于解算地面站與飛機之間的偽距。在上述應用中，都需要估計雙音信號的頻率，但是，文獻[4]對雙音信號頻率估計方法提出了新的挑戰：(1)高速跳頻信號駐留時間短，使得信號采樣時間較短；(2)為了提高偽距解算的性能，需要提高解算帶寬，那么就需要時頻矩陣中第一行的雙音頻率的間隔盡可能小。此外，在基于雙音信號的無線傳感器網絡定位中，為了定位需要多個節點同時發送雙音信號，不同節點發送的雙音信號要保持一定間隔以防止相互影響。當頻譜資源緊張時，雙音信號頻率間隔可能會比較小。

雙音信號的頻率估計問題是一個經典的問題，具有廣泛的研究[5-9]。文獻[5]給出了多音信號參數的近似的最大似然估計量，并利用基于離散傅里葉變換(Discrete Fourier Transform,DFT)的方法解算。但是，當音數大于1時，基于DFT的方法會使估計值存在偏差。雖然加窗可以降低這種偏差，但是卻會降低頻率的分辨率。基于旋轉因子不變法(Estimating Signal Parameters Viarotational Invariance Techniques,ESPRIT)、多重信號分類(Multiple Signal Classification,MUSIC)的頻率估計方法[7-8]具有高分辨率的特性，但是當采樣樣本較少時，協方差的估計誤差會比較大；另外，在多音情況下ESPRIT等方法對頻率估計的均方誤差(Mean-Square Error,MSE)也無法達到克拉美羅下界(Cramer-Rao Lower Bound,CRLB)[7]。文獻[7]提出NS-ESPRIT雖然能夠使頻率估計的MSE達到CRLB，但是，頻率子帶范圍的選取對結果影響很大，而文獻[7]并沒有給出頻率子帶范圍的有效選取方法。

文獻[9]提出了一種高效的迭代算法，它利用期望最大化(Expectation Maximization,EM)算法結合IEEE 1057標準實現對多音實信號的參數估計。該方法首先基于DFT的方法得到多音信號參數的初估計；然后將多音信號的采樣數據分解為多個單音采樣數據；最后利用IEEE 1057標準采用最小二乘擬合的方法迭代計算各個單音的參數。上述過程中，多音信號采樣數據的分解和各個單音參數的估計迭代進行，得到最終的估計值。當雙音信號頻率間隔足夠大、采樣次數足夠多時，該方法頻率估計的MSE可以達到CRLB。但是，文獻[9]并沒有考慮雙音間隔較小、采樣數據較短的情況，同時文獻[9]中的方法是針對實信號的。

針對雙音信號頻率間隔小、采樣數據少情況下估計雙音信號頻率面臨的問題，本文在文獻[9]的基礎上提出基于EM算法的雙音頻率估計算法，主要改進在于以下幾點：(1)EM算法中初值對結果有重要影響，針對雙音間隔小的情況，選擇具有高分辨率特性的MUSIC算法作為雙音頻率初值的獲取方法；(2)針對樣本數據少的問題，在MUSIC算法中采用修正的前向-后向線性預測(Modified Forward-Backward Linear Prediction,MFBLP)方法估計協方差[10]；(3)對于單音頻率的估計，采用高效的二分搜索方法[11]。

本文提出的算法在EM算法的初估計時采用MUSIC算法，因此，在文中簡稱MEM算法。

1 信號模型

復雙音信號可表示為：

(1)

其中，ai為雙音信號的幅度，fi為雙音信號頻率，θi為雙音信號相位。那么，雙音信號的未知參數為：

(2)

雙音信號經過采樣頻率為fs的模數轉換器采樣后，得到信號樣本為：

0≤n≤N-1

(3)

s=[s(0)s(1)…s(N-1)]T

(4)

其中，w(n)為方差為σ2的高斯白噪聲，Ts=1/fs為信號的采樣間隔，N為樣本數，NTs為總采樣時間。文獻[4]中，高速跳頻信號單跳駐留時間短，令駐留時間為Tr，那么總采樣時間NTs≤Tr，這使得在采樣頻率一定的條件下，樣本數N較小。s為得到的樣本向量。

雙音信號可以表示成兩個單音信號的和，則公式(3)可以表示為：

(5)

xi(n)=aiexp(j(2πfiTsn+θi))+wi(n)

(6)

xi=[xi(0)xi(1)…xi(N-1)]T

(7)

其中，xi表示單音信號的樣本序列，wi(n)是噪聲w(n)的分量，可表示為：

wi(n)=βiw(n)

(8)

文獻[9]指出βi可以是任意值，但是必須滿足β1+β2=1。此外，各個單音信號的未知參數可表示為Ai=[aifiθi]T。

2 MEM雙音頻率估計算法

2.1 算法描述

MEM雙音頻率估計算法大致可以分為3個步驟：

(1)估計未知參數A的初值；

第(2)步和第(3)步是一個迭代過程，即第(3)步的估計結果用于下一次迭代時第(2)步的估計。這種迭代算法稱為EM算法[9]。

初值的估計對最終的結果具有重要的影響，文獻[9]提出可用DFT計算初值，但是要求雙音信號滿足一定的條件，即雙音頻率的間隔Δf≥3/T，其中T表示信號采樣時間。文獻[9]并沒有考慮雙音頻率間隔較小的情況，而文獻[4]面臨雙音頻率間隔較小的情況，因此需要選擇合適的初值估計方法。本文在2.2節中對初值估計方法進行詳細描述。

由于雙音信號采樣時引入了噪聲，在式(3)中，高斯白噪聲w(n)是未知的。在雙音未知參數A的初估計后，令A的初估計為A(0)= [a1(0)f1(0)θ1(0)a2(0)f2(0)θ2(0)]T，那么，在第(2)步對于第j次迭代，可以得到噪聲的估計為：

(9)

那么，此時單音信號的樣本的估計值可表示為：

(10)

(11)

此時，得到兩個單音信號的樣本的估計值為：

(12)

通過上述方法將雙音信號的樣本值分解成兩個單音信號的樣本值，那么，在第(3)步可以采用單音信號的估計方法估計單音信號的參數。采用二分搜索的方法[11]估計單音信號的頻率，單音信號相位和幅度可利用文獻[12]給出的估計方法進行估計。

將第(3)步得到的兩個單音信號的參數的估計值作為下一次迭代時的輸入，在迭代時設置門限Δf，當滿足下列條件時，迭代結束：

(13)

2.2 初值估計方法

算法第(1)步中的初估計值對算法的最終結果有重要的影響，特別是當雙音信號的頻率間隔比較小時體現得更為明顯。文獻[9]提出采用DFT的方法進行初值估計。文獻[11]提出的二分搜索方法是在DFT的基礎上進行迭代計算，其針對單音信號能夠獲得較高的估計精度，可以作為一種EM算法中估計初值的方法。但是，基于DFT的方法會使估計值存在偏移，特別是當雙音頻率間隔較小時，這種偏移會增大。

MUSIC算法具有高分辨率的特性。令雙音信號的采樣數據s的協方差為R，對R進行特征分解，特征值較小的N-2個特征向量構成的矩陣為UN，那么雙音信號的頻率估計值為下列函數的兩個最大的峰值：

(14)

(15)

其中，[]H表示共軛轉置。通常協方差R可通過計算樣本的自相關矩陣進行估計，但是，當采樣數據較少時協方差矩陣的估計會存在較大的誤差。借助MFBLP方法對協方差的估計進行修正，可提高頻率估計的精度。此時，協方差R的估計值為：

(16)

(17)

其中，m為MFBLP方法的階數，[]*表示共軛。

2.3 雙音頻率估計的CRLB

將CRLB作為衡量估計方法性能的標準。單音信號頻率估計的CRLB由文獻[12]給出。如果不考慮雙音間的影響，可以將雙音信號看成兩個獨立的單音信號的疊加，那么，雙音信號頻率估計的CRLB為：

(18)

3 仿真分析

根據文獻[4]中仿真的高速跳頻信號的單跳駐留時間，在仿真中信號采樣時間T為50 μs，基于此，文獻[9]要求雙音間的間隔大于3/T=60 kHz，而本文提出的MEM算法則不受這個限制。在仿真中，采樣頻率fs為250 kHz，那么樣本數N為12，雙音信號的幅度a1=a2=1，雙音信號的初相θ1和θ2隨機生成，雙音信號頻率f1=10 kHz，f2則選擇不同的頻率。迭代時的門限值Δf設置為1 Hz。

在仿真中，為了分析MEM算法的性能，對比三種不同方法的雙音頻率估計性能。第一種方法是本文提出的MEM算法。文獻[10]指出當MFBLP方法中階數m為3N/4時性能最佳，因此取m=9。第二種方法初值采用基于二分搜索的DFT方法進行估計[11]，其他步驟與第一種方法相同，稱這種方法為初值采用二分搜索的EM算法。第三種方法采用MUSIC算法直接估計雙音頻率。

門特卡羅仿真次數為104次，結果如圖1所示。在圖1中，在不同信噪比(Signal-to-Noise Ratio,SNR)下，仿真三種頻率估計算法的均方根誤差(Root-Mean-Square Error,RMSE)，而將CRLB的平方根CRLB1/2作為衡量估計性能的標準。圖1(a)中f2為30 kHz，雙音的頻率間隔為20 kHz。頻率估計的RMSE存在門限效應，初值采用二分搜索的EM算法在這種情況下已經失效，其頻率估計的RMSE都處于無信息區。二分搜索方法是基于DFT的，當雙音頻率間隔較小時，DFT方法可能會無法分辨雙音信號。MUSIC算法具有高分辨特性，在上述情況下，它雖然有效并且門限效應的門限與MEM算法相當，但是，此時頻率估計的RMSE與CRLB1/2相差較大。MEM算法頻率估計的RMSE在高信噪比下明顯優于MUSIC算法，并且與CRLB1/2非常接近。MEM算法頻率估計的RMSE沒有達到CRLB1/2的一個重要原因是當雙音非常接近時，公式(18)不能準確表達CRLB，由于雙音間的相互影響，實際的CRLB值比公式(18)的計算結果大。

圖1 雙音頻率估計的均方根誤差

圖1(b) 中f2為40 kHz，雙音的頻率間隔為30 kHz。在這種情況下，初值采用二分搜索的EM算法頻率估計的RMSE的門限明顯高于其他兩種方法。MUSIC算法與MEM算法相比，頻率估計的RMSE略高于后者。MEM算法頻率估計的RMSE與CRLB1/2更加接近。圖1(c)進一步提高了f2的值將其調整為50 kHz，此時雙音的頻率間隔為40 kHz。這時，初值采用二分搜索的EM算法和MEM算法頻率估計的RMSE幾乎相當，而MUSIC算法頻率估計的RMSE比前兩者略高。MEM算法頻率估計的RMSE已經達到CRLB1/2。可見，隨著雙音頻率間隔的提高，三種算法的估計性能逐漸接近，但是在雙音間隔較小的情況下，MEM算法在三者中的估計性能最好。

4 結論

本文針對雙音信號間隔小、采樣時間短時估計雙音頻率時面臨的問題，提出基于EM算法的雙音頻率估計方法，稱該方法為MEM算法。MEM算法在頻率初估計階段采用MUSIC算法，并且針對采樣時間短的情況采用MFBLP方法估計協方差以提高MUSIC算法估計性能。通過仿真分析了MEM算法的性能，并與初值采用二分搜索的EM算法、MUSIC算法進行了比較，主要結果有以下幾點：(1)當頻率間隔為1/T時，即頻率間隔較小時，MEM算法具有較好的頻率估計性能，MUSIC算法頻率估計的RMSE與MEM算法有較大差距，而初值采用二分搜索的EM算法完全失效；(2)隨著雙音頻率間隔的增大，三者頻率估計的RMSE逐漸接近。可見，在雙音頻率間隔較小的情況下，采用MEM可有效估計雙音信號頻率。