基于迭代的多頻頻率估計方法

陳 皓，雷 藝

（1.中國電子科技集團(tuán)第三十八研究所孔徑陣列與空間探測安徽省重點(diǎn)實(shí)驗(yàn)室，合肥，230031；2.合肥工業(yè)大學(xué)計算機(jī)與信息學(xué)院，合肥，230009）

引 言

調(diào)頻連續(xù)波（Frequency modulated continuous wave，F(xiàn)MCW）雷達(dá)是一種常用的雷達(dá)，在軍民產(chǎn)業(yè)中均有廣泛的使用。FMCW 雷達(dá)發(fā)射連續(xù)調(diào)頻信號，利用發(fā)送和接收的頻率之差計算目標(biāo)的距離。因此，雷達(dá)需要對發(fā)送與接收信號的差頻信號進(jìn)行頻率估計，頻率估計的精度直接影響雷達(dá)測距的精度。頻率估計中最簡單有效的方法是利用離散傅里葉變換（Discrete Fourier transform，DFT）把信號從時域變換到頻域，在頻域可以直接得到信號的頻率值。但是利用DFT 進(jìn)行頻率估計容易帶來一些測量誤差。誤差主要來源于兩方面：（1）離散變換的柵欄效應(yīng)，DFT 的結(jié)果是頻譜的采樣值，難以保證頻譜峰值正好和待測頻率值重合；（2）負(fù)頻譜泄漏，實(shí)數(shù)信號的頻譜既包含正頻譜也包含負(fù)頻譜，正負(fù)頻譜疊加后的總頻譜最大值與信號的真實(shí)頻率會出現(xiàn)偏差。

為了解決柵欄效應(yīng)帶來的誤差，該領(lǐng)域研究者們已提出很多有效的算法［1?4］，其中Candon 算法［5］和Orguner 算法［6］是復(fù)指數(shù)信號頻率估計中最經(jīng)典的算法，只需利用DFT 結(jié)果的3 個采樣點(diǎn)即可計算出精確的頻率值。但由于沒有考慮負(fù)頻譜干擾，這兩種方法在實(shí)數(shù)信號的頻率估計中并不適用。為了抑制負(fù)頻譜干擾，學(xué)者們又提出了加窗算法，其中包括一些常用的加窗算法［7］（如Hamming，Hann 和Blackman 窗）以及基于最大旁瓣衰減（Maximum?sidelobe?decay，MSD）窗的高精度算法。使用MSD 窗的算法有經(jīng)典三點(diǎn)法［8］，以及在經(jīng)典三點(diǎn)法基礎(chǔ)上改進(jìn)而來的新三點(diǎn)法［9］。相比經(jīng)典三點(diǎn)法，改進(jìn)三點(diǎn)法增加了糾正因子，能夠比經(jīng)典三點(diǎn)法有更好的負(fù)頻譜抑制效果。同時，Borkowski 等基于MSD 窗提出的算法［10］也能夠達(dá)到非常高的頻率估計精度。實(shí)際應(yīng)用中，F(xiàn)MCW 雷達(dá)測距時往往不止探測一個目標(biāo)，當(dāng)對多個頻率同時進(jìn)行頻率估計時，除了單頻估計遇到的柵欄效應(yīng)和負(fù)頻譜干擾以外，多個頻率頻譜之間的互相干擾也嚴(yán)重影響頻率估計的精度。加窗方法雖然也能對多個頻率之間的互相干擾有一定抑制，但無論在頻率估計精度還是計算復(fù)雜度上都有一定的局限性。本文首先對多頻信號的離散頻域形式進(jìn)行了詳細(xì)推導(dǎo)，得到每個頻率對應(yīng)的最大譜線和其左右兩譜線的數(shù)學(xué)表達(dá)式。進(jìn)而定義了干擾系數(shù)，并采取迭代的手段，利用兩種不同的方法逐步消除負(fù)頻譜和其他頻率帶來的干擾。仿真實(shí)驗(yàn)證明，本文提出方法相比傳統(tǒng)加窗算法在頻率估計精度上有明顯提升，而且由于不需要進(jìn)行加窗處理，也降低了計算復(fù)雜度。

1 多頻頻率估計基礎(chǔ)理論

假設(shè)采集的信號中有M個頻率分量（f1，f2，…，fM），高斯白噪聲下的多頻離散時間信號表示為

式中：fm為信號的頻率；fs為采樣頻率；Am為第m個頻率信號的幅度；φm為第m個頻率信號的初始相位。式（1）也可以表示成復(fù)數(shù)信號相加的形式

進(jìn)而對式（2）進(jìn)行N點(diǎn)DFT 可以得到信號的離散頻域表示為

本文通過DFT 結(jié)果中每個頻率附近的3 個采樣點(diǎn)對多個頻率值進(jìn)行估計。待測頻率fm可以表示為fm=fs?(km+δm)/N，其中km為整數(shù)部分，δm(-0.5 ≤δm<0.5)為分?jǐn)?shù)部分。km可以從 DFT 的結(jié)果直接得到，精確頻率估計的目的即為估計δm的值。首先在DFT 結(jié)果中的前N/2 個點(diǎn)里尋找第m個頻率對應(yīng)的最大譜線和其左右兩譜線。其對應(yīng)的數(shù)學(xué)表達(dá)式為

為了后續(xù)的推導(dǎo)方便，本文把R[km]、R[km-1]、R[km+1]表示成理想的無干擾項(xiàng)和干擾項(xiàng)的相加（干擾項(xiàng)包括負(fù)頻譜干擾和其他頻率帶來的干擾），于是本文定義R[km]、R[km-1]、R[km+1]所對應(yīng)的理想譜線值為am、bm、cm（即為無負(fù)頻譜干擾和其他頻率頻譜干擾時的理想譜線值），并定義：

頻率i對頻率m最大譜線的干擾系數(shù)為

頻率i對頻率m最大譜線左側(cè)譜線的干擾系數(shù)為

頻率i對頻率m最大譜線右側(cè)譜線的干擾系數(shù)為

頻率i的負(fù)頻譜對頻率m最大譜線的干擾系數(shù)為

頻率i的負(fù)頻譜對頻率m最大譜線左側(cè)譜線的干擾系數(shù)為

頻率i的負(fù)頻譜對頻率m最大譜線右側(cè)譜線的干擾系數(shù)為

定義了上述參量以后，式（4—6）即可表示為

2 基于迭代的多頻頻率估計算法

本文目的是根據(jù)已知的R[km]、R[km-1]、R[km+1]求得am、bm、cm，進(jìn)而利用 Orguner 算法進(jìn)行精確的頻率估計。Orguner 算法在無噪聲無負(fù)頻率干擾時對單頻信號進(jìn)行無偏估計，是對復(fù)指數(shù)信號進(jìn)行頻率估計的最佳算法之一，其計算公式為［6］

由式（7—9）可以看出，干擾系數(shù)是與頻率估計結(jié)果相關(guān)的值，方程沒有解析解。因此，為了得到am、bm、cm的值，本文提出了迭代的方法。通過數(shù)次迭代得到am、bm、cm的近似值后，再利用 Orguner 算法進(jìn)行精確的頻率估計。

迭代法在一些單頻頻率估計方法中也有應(yīng)用［11?13］。在這些方法中，迭代目的都是為了盡可能多地去除負(fù)頻譜帶來的干擾，因?yàn)樨?fù)頻譜的干擾越小，頻率估計就越精確。但是這些方法都只關(guān)心單頻情形下負(fù)頻譜的干擾，對于多頻情形并沒有考慮，而實(shí)際FMCW 雷達(dá)應(yīng)用中，多目標(biāo)即多頻的情況更為普遍，而頻率之間的干擾也比負(fù)頻譜的干擾更大。因此本文迭代的思路為：首先假設(shè)沒有負(fù)頻譜和其他頻率的干擾，得到初步的頻率估計值，進(jìn)而根據(jù)初步的頻率估計值計算出負(fù)頻率和其他頻率的干擾值并進(jìn)行剔除，從而得到更精確的頻率估計值。這樣循環(huán)迭代下去，頻率估計的結(jié)果也會逐漸收斂到一個相對較為精確的值。在每個循環(huán)中進(jìn)行干擾剔除時，為了求得第q次迭代的值，本文又提出了兩種不同的計算方法：（1）利用上一次迭代中的計算出的干擾系數(shù)；（2）只利用由計算出的干擾系數(shù)，不用上一次迭代中的值。具體的詳細(xì)步驟如表1 和表2 所示。

表1 多頻頻率估計算法1 步驟Table 1 Procedure of Algorithm 1

表2 多頻頻率估計算法2 步驟Table 2 Procedure of Algorithm 2

3 仿真結(jié)果與分析

本文研究了兩種多頻迭代算法的收斂速度。仿真模擬了無噪聲情況下同時對2 個頻率進(jìn)行頻率估計的場景，其中1 個頻率固定為f1=20 kHz，另一個頻率f2的變化范圍為40～80 kHz。仿真中用兩個不同頻率的正弦波模擬兩個頻率的信號，采樣速率設(shè)置為200 kHz，采樣點(diǎn)數(shù)N為64。仿真共進(jìn)行5 次迭代，實(shí)驗(yàn)中記錄每一次迭代后的結(jié)果。圖1，2 為算法1 和算法2 對2 個頻率（f1和f2）進(jìn)行頻率估計的結(jié)果。從圖中可以看出，算法1 和算法2 均在3 次迭代后達(dá)到穩(wěn)定狀態(tài)，因此基于計算復(fù)雜度考慮，后續(xù)仿真均使用3 次迭代。另外，比較圖1 和2 可以看出，算法1 相比算法2 的最終收斂誤差更小，因此后文與其他算法比較時均使用效果更好的算法1。

圖3 仿真比較了本文提出多頻頻率估計算法和 Orguner 算法［6］、Borkowski 算法［10］、新三點(diǎn)法［9］在無噪聲情況下的測量誤差。Orguner 算法是用于復(fù)指數(shù)頻譜估計的最精確的方法之一；新三點(diǎn)法是Bele?ga 等人在MSD 窗經(jīng)典三點(diǎn)法基礎(chǔ)上提出來的，專門用于負(fù)頻譜的干擾抑制；Borkowski 算法是電力系統(tǒng)頻率估計中提出的一種性能較好的方法。仿真模擬了無噪聲情況下同時對2 個頻率進(jìn)行頻率估計的場景，其中一個頻率固定為f1=20 kHz，另一個頻率f2的變化范圍為40～80 kHz。仿真中用2 個不同頻率的正弦波模擬2 個頻率的信號，采樣速率設(shè)置為200 kHz，采樣點(diǎn)數(shù)N為64。從圖中可以看出，相比上述幾種算法，本文提出的算法有更高的測量精度。

圖1 算法1 中測量誤差隨頻率2 變化情況（N=64）Fig.1 Frequency estimation error versus the frequency of f2 for different numbers of iteration by Algorithm 1 (N=64)

圖2 算法2 中測量誤差隨頻率2 變化情況（N=64）Fig.2 Frequency estimation error versus the frequency of f2 for different numbers of iteration for Algorithm 2 (N=64)

圖3 無噪聲時使用不同算法測量誤差隨頻率2 變化情況（N=64）Fig.3 Frequency estimation error versus the frequency of f2 for different algorithms without noise (N=64)

本文接下來研究了在有噪聲情況下各方法頻率估計的結(jié)果（假設(shè)噪聲為高斯白噪聲）。仿真模擬了40 dB 信噪比下同時對2 個頻率進(jìn)行頻率估計的場景，其中一個頻率固定為f1=20 kHz，另一個頻率f2的變化范圍為40～80 kHz。仿真中用2 個不同頻率的正弦波模擬2 個頻率的信號，噪聲為高斯白噪聲，采樣速率設(shè)置為200 kHz，采樣點(diǎn)數(shù)N為64。圖4 比較了各算法的均方根誤差（Root mean square error，RMSE）。從圖中可以看出，本文提出算法在有噪聲情況下相比其他算法依然有優(yōu)勢。這是因?yàn)楸疚姆椒]有使用窗函數(shù)，避免了窗函數(shù)引起的信噪比損失［14］，而且本文方法不需要進(jìn)行加窗處理，因此也減輕了計算復(fù)雜度的負(fù)擔(dān)。

圖4 信噪比40 dB 時測量誤差RMSE 隨頻率2 變換情況（N=64）Fig.4 RMSE versus the frequency of f2 for different algorithms with SNR=40 dB (N=64)

本文算法在精度上有提升主要有兩個原因。（1）該算法充分考慮了頻率估計誤差的來源（負(fù)頻率干擾和其他頻率干擾），并通過詳細(xì)的數(shù)學(xué)推導(dǎo)得到干擾項(xiàng)的數(shù)學(xué)表達(dá)式，進(jìn)而通過迭代的方法進(jìn)行干擾的剔除。相比之下，傳統(tǒng)的加窗算法只是通過窗函數(shù)盡量抑制頻譜的副瓣，干擾的抑制效果不如本文算法。（2）窗函數(shù)會帶來信噪比的損失，本文算法沒有用到窗函數(shù)，也因此避免了信噪比損失帶來的精度下降。本文算法在計算效率方面有提升主要體現(xiàn)在少了與窗函數(shù)相乘這一環(huán)節(jié)，減少了計算量。代價是迭代的次數(shù)如果太多，也會增加計算量。但是本文證明，算法經(jīng)過3 次迭代即可收斂，迭代增加的計算量很小，因此總體來說本文算法在計算效率方面有所提升。

4 結(jié)束語

針對FMCW 雷達(dá)測距中遇到的多頻頻率估計問題，本文在Orguner 算法的基礎(chǔ)上利用迭代的手段，逐步消除負(fù)頻譜和其他頻率帶來的干擾，最終實(shí)現(xiàn)了較為精確的頻率估計。本文設(shè)計了2 種迭代方法并通過比較選擇了效果更佳的方法，該方法與傳統(tǒng)頻率估計方法相比，無論在估計精度上還是計算復(fù)雜度上都有明顯優(yōu)勢。