白噪聲背景下基于MUSIC算法的信號頻率檢測仿真

秦 陽，劉洪彪，胥勇軍，趙 樺，周陳超，房 敏，鄧朝平

（解放軍95806部隊(duì)，北京 100076）

0 引 言

在實(shí)際應(yīng)用中，常常需要對空間中存在的多個(gè)信號源進(jìn)行分解，以便跟蹤或檢測感興趣的空間信號，抑制那些被認(rèn)為是干擾的空間信號。例如對天線陣列接收的空間信號所進(jìn)行的分析與處理稱為陣列信號處理。而空間譜估計(jì)技術(shù)是在波束形成技術(shù)、零點(diǎn)技術(shù)和時(shí)域譜估計(jì)技術(shù)的基礎(chǔ)上發(fā)展起來的一種技術(shù)。與頻譜表示信號在各個(gè)頻率上的能量分布相對應(yīng)，空間譜則可解釋為信號在空間各個(gè)方向上的能量分布，空間譜估計(jì)技術(shù)的目標(biāo)是研究提高在處理帶寬內(nèi)空間信號角度的估計(jì)精度、角度分辨率和提高運(yùn)算速度的各種算法。經(jīng)過多年的發(fā)展，已經(jīng)產(chǎn)生了大量性能優(yōu)異的測向算法可資利用，典型的有 MUSIC［1］、ESPRIT［2］等。

MUSIC算法是基于特征結(jié)構(gòu)分析的空間譜估計(jì)方法，是空間譜估計(jì)技術(shù)的典型代表。其測向原理是根據(jù)矩陣特征分解的理論［3］對陣列輸出協(xié)方差矩陣進(jìn)行特征分解，將信號空間分解為噪聲子空間D和信號子空間C，利用噪聲子空間D與陣列的方向矩陣A的列矢量正交的性質(zhì)，構(gòu)造空間譜函數(shù)P（w）并進(jìn)行譜峰搜索，從而估計(jì)出信號方向信息。本文通過仿真驗(yàn)證了MUSIC算法可以很好地對信號參數(shù)進(jìn)行估計(jì)。

1 MUSIC算法基本原理

1.1 信號子空間和噪聲子空間［4］

一個(gè)系統(tǒng)的輸入信號可以表示為（有用）信號S及噪聲N之和，即：

式中：X＝ ［x1，x2…xp］T，為輸入信號向量；S＝［s1，s2…sp］T， 為 有 用 信 號 向 量；N＝［n1，n2…np］T，為噪聲信號向量。

于是，可以定義輸入信號的相關(guān)函數(shù)矩陣Rxx為：

式中：“＊”表示復(fù)數(shù)共軛；由定義可知，信號矩陣Rxx為Hermite矩陣，即滿足Rxx＊＝RxxT。

因此，若輸入信號和噪聲信號互不相關(guān)，則相關(guān)矩陣Rxx可寫為：

式中：Rs和Rn分別為信號與噪聲的相關(guān)矩陣。

假如有p個(gè)互不相關(guān)的信號源同時(shí)作用于系統(tǒng)的m個(gè)輸入端，系統(tǒng)的輸入向量X構(gòu)成了p維酉向量空間Up，則第p個(gè)源在第m個(gè)輸入端上的輸入信號可以寫為：

式中：ap（t）為輸入信號的復(fù)包絡(luò)；ω為信號載頻；φpm為信號的相位滯后。

若第p個(gè)信號源的復(fù)包絡(luò)ap（t）為隨機(jī)的，則它的有用信號向量S也為一隨機(jī)向量。可以用它在m個(gè)輸入端的相位滯后來定義一個(gè)信號特征向量Sp，即：

則由p個(gè)信號方向向量所張成的空間Span｛Sp，p＝1，2，…，P｝稱為在全空間UM中的p維信號子空間C（其中p≤m），而將它在全空間UM中的正交補(bǔ)空間C⊥稱為M-P維噪聲子空間D。

信號子空間和噪聲子空間的物理意義是明顯的，若一個(gè)理想的系統(tǒng)輸入無噪聲，則經(jīng)過統(tǒng)計(jì)平均所得的輸入向量必落在信號子空間C內(nèi)，但是若系統(tǒng)輸入有噪聲時(shí)，它就不能完全落在信號子空間內(nèi)，而它中間的噪聲分量則應(yīng)落在噪聲子空間C⊥內(nèi)。

1.2 MUSIC算法［5］

MUSIC算法的基本思想是將陣列輸出數(shù)據(jù)的協(xié)方差矩陣進(jìn)行特征分解，從而得到與信號分量相對應(yīng)的信號子空間和與信號分量相正交的噪聲子空間，利用2個(gè)子空間的正交特性構(gòu)造出“針狀”空間譜峰來估計(jì)信號的方向信息。

設(shè)空間有p個(gè)互不相關(guān)的信號，以方位角θ1，θ2，…，θp入射到具有m個(gè)接收陣元的接收陣元陣列中，入射信號的數(shù)目p小于陣列的陣元數(shù)m，則此陣列系統(tǒng)的信號模型為：

對由上式描述的陣列信號觀測模型做以下假設(shè)：

假設(shè)1：對于不同的wi值，向量a（wi）相互線性獨(dú)立；

假設(shè)2：加性噪聲向量N（n）的每個(gè)原色都是零均值的復(fù)白噪聲，它們不相關(guān)，并且具有相同的方差σn2；

假設(shè)3：矩陣P＝E｛S（n）SH（n）｝非奇異，即rank（P）＝p。

上述3個(gè)假設(shè)都只是一般的假設(shè)，在實(shí)際中容易得到滿足。可見：

式中：AH為A的共軛轉(zhuǎn)置；I為單位矩陣。

于是，Rxx是個(gè)對稱矩陣，令其特征值分解為：

由于A列滿秩，故rank（APAH）＝rank（P）＝p，且這里假定p＜m，于是有：

式中：a12，…，ap2為無加性噪聲時(shí)的觀測信號AX（n）的自相關(guān)矩陣APAH的特征值。

同時(shí)用UH左乘和U右乘上式，得：

這表明，自相關(guān)矩陣Rxx的特征值為：

根據(jù)信號特征值和噪聲特征值，將特征矩陣U的列向量分為兩部分，即：

式中：信號子空間C＝ ［S1，…，SP］＝ ［U1，…，UP］，由信號特征向量組成；噪聲子空間D＝ ［g1，…，

又由Rxx＝APAH＋σn2I，有RxxD＝APAHD＋σn2D，利用上式的結(jié)果，得到：

上式成立的充分必要條件是AHD＝O。

將A＝ ［a（w1），…，a（wp）］代入式（15）得到：

顯然，當(dāng)w≠w1，…，wp時(shí)，aH（w）D≠OT。

將上式改成標(biāo)量形式，可以定義一種類似于功率譜的函數(shù)：gm-p］＝ ［Up＋1，…，Um］，由噪聲特征向量組成。

考察：

式（17）取峰值的p個(gè)w值給出p個(gè)信號的波達(dá)方向θ1，θ2，…，θp。

這樣定義的函數(shù)描述了空間參數(shù)（即波達(dá)方向）的分布，因此成為空間譜。它能夠?qū)Χ鄠€(gè)信號進(jìn)行識別。可以看出，此算法的關(guān)鍵是要對自相關(guān)矩陣進(jìn)行特征分解。

2 計(jì)算機(jī)仿真研究

2.1 仿真數(shù)據(jù)模型

通過對前面MUSIC算法基本原理的分析，N個(gè)正弦信號加白噪聲的自相關(guān)函數(shù)可以表示為：

式中：rxx（n）為輸入信號的自相關(guān)函數(shù)；rss（n）為有用信號的自相關(guān)函數(shù)；ruu（n）為白噪聲信號的自相關(guān)函數(shù)；wi為正弦波的頻率；Ai為正弦波的幅度；σn2為白噪聲的方差；δ（k）為沖擊函數(shù)。

本文對這種正弦波加白噪聲的情形進(jìn)行仿真，假定仿真的觀測數(shù)據(jù)分別由下面兩種情況產(chǎn)生：

（1）單個(gè)正弦信號檢測情況

式中：u（n）為一高斯白噪聲，其均值為0，方差為1，采用混合同余法產(chǎn)生。

用MUSIC方法估計(jì)觀測數(shù)據(jù)中正弦波的頻率，并給出白噪聲方差σn2與復(fù)正弦波的振幅A的估計(jì)值。

2.2 仿真圖形結(jié)果

westimate1＝0.250 0π（估計(jì)頻率），S＝1.405 8（噪聲估計(jì)功率），A＝3.983 5（信號振幅估計(jì)），仿真結(jié)果如圖1所示。

（2）多個(gè)正弦信號情況

圖1 單個(gè)正弦信號短數(shù)據(jù)（N＝64）時(shí)的估計(jì)結(jié)果

westimate1＝0.250 0π（估計(jì)頻率），S＝0.964 6（噪聲估計(jì)功率），A＝4.008 3（信號振幅估計(jì)），仿真結(jié)果如圖2所示。

圖2 單個(gè)正弦信號長數(shù)據(jù)（N＝1 024）時(shí)的估計(jì)結(jié)果

S＝2.1686（噪聲估計(jì)功率），A1＝3.887 2，A2＝2.187 4，A3＝2.910 4（信號振幅估計(jì)），仿真結(jié)果如圖3所示。

圖3 多個(gè)正弦信號短數(shù)據(jù)（N＝64）時(shí)的估計(jì)結(jié)果

S＝1.038 5（噪聲估計(jì)功率），A1＝4.009 8，A2＝2.021 8，A3＝2.991 7（信號振幅估計(jì)），仿真結(jié)果如圖4所示。

圖4 多個(gè)正弦信號長數(shù)據(jù)（N＝1 024）時(shí)的估計(jì)結(jié)果

2.3 仿真結(jié)果分析

（2）在信號數(shù)據(jù)較少，即信號長度較短的情況下（N＝64），噪聲功率和信號幅度的估計(jì)值有一定的偏差，這主要是因?yàn)樾盘柖蹋哟靶?yīng)明顯，自相關(guān)矩陣的估計(jì)誤差比較大，因此特征值分解存在誤差。增加信號長度可以減少信號自相關(guān)矩陣的估計(jì)誤差，能夠得到比較理想的結(jié)果，這可以從長數(shù)據(jù)（N＝1 024）中看出。

（3）在MUSIC方法中，前面幾個(gè)噪聲特征值對應(yīng)的特征矢量都是相同的，于是嘗試用單個(gè)特征矢量對單個(gè)正弦信號進(jìn)行頻率估計(jì)，結(jié)果出現(xiàn)了虛假峰，即在不等于信號頻率的地方多出來一些峰值。研究特征矢量的正交性可知，特征矢量是兩兩正交的。對矩陣特征分解得到的特征矢量組成的矩陣U進(jìn)行處理，程序如下：

計(jì)算得到的結(jié)果為：

圖5 出現(xiàn)虛假峰的結(jié)果

3 結(jié)束語

通過對MUSIC算法的仿真研究討論，可以看出該算法具有高分辨特性及較小的估計(jì)方差，有著優(yōu)異的性能。在信號處理中幾個(gè)重要問題，如本文討論的噪聲中幾個(gè)疊加信號的參數(shù)估計(jì)，還有窄帶傳感器陣列的方向搜索，以及重疊回波的分辨率等，都可以采取此算法進(jìn)行處理。但采樣協(xié)方差矩陣的特征分解涉及大量不規(guī)則矩陣運(yùn)算，運(yùn)算量較大，限制了算法的實(shí)時(shí)運(yùn)用。而且這里假定的信號源都是不相關(guān)的，對相關(guān)信號源運(yùn)用MUSIC算法還需進(jìn)行去相關(guān)處理［6-7］。因此，對算法的改進(jìn)研究成為近年來許多學(xué)者研究的課題［8-10］。

［1］ Schmid T R O.Multiple emitter location and signal parameter estimation［J］.IEEE Transactions on AP，1986，34（3）：276-280.

［2］ ROY R，Kailath T.ESPRIT-subspace rotation approach to estimation of parameters of cissoids in noise［J］.IEEE Transactions on ASSP，1986，34（10）：1340-1343.

［3］ 程云鵬，張凱院，徐仲.矩陣論［M］.第3版.西安：西北工業(yè)大學(xué)出版社，2006.

［4］ 沈福民.自適應(yīng)信號處理［M］.西安：西安電子科技大學(xué)出版社，2001.

［5］ 王永良，陳輝，彭應(yīng)寧，等.空間譜估計(jì)理論與算法［M］.北京：清華大學(xué)出版社，2004.

［6］ Lin J D，F(xiàn)ang W H，Wang Y Y.FSF MUSIC for joint DOA and frequency estimation and its performance analysis［J］.IEEE Transactions on Signal Processing，2006，52（12）：4529-4542.

［7］ Mestr X，Langunas M A.Modified subspace algorithms for DOA estimation with large arrays［J］.IEEE Transactions on Signal Processing，2008，56（2）：598-614.

［8］ 黃磊，吳順君，張林讓，等.快速子空間分解方法及其維數(shù)的快速估計(jì)［J］.電子學(xué)報(bào)，2005，33（6）：977-981.

［9］ 張立杰，黃建國，史文濤，等.正交投影陣列信號子空間估計(jì)方法［J］.系統(tǒng)工程與電子技術(shù)學(xué)報(bào)，2009（9）：2063-2066.

［10］余繼周，陳定昌.一種DOA估計(jì)的快速子空間算法［J］.現(xiàn)代電子技術(shù)，2005，28（12）：90-92.