孫永梅，趙維，闞園園

(大連交通大學(xué) 電氣信息學(xué)院，遼寧 大連 116028)＊

0 引言

譜估計(jì)是數(shù)字信號處理的十分重要研究領(lǐng)域，人們先后提出了各種譜估計(jì)的理論和方法，MUSIC譜估計(jì)因其具有分辨率高、穩(wěn)定性好等優(yōu)點(diǎn)，在生物醫(yī)學(xué)工程、雷達(dá)、水聲信號等領(lǐng)域得到廣泛應(yīng)用.MUSIC譜估計(jì)方法基本思想是直接對估計(jì)的隨機(jī)過程相關(guān)矩陣進(jìn)行特征分解，分別生成信號子空間和噪聲子空間，利用信號子空間和噪聲子空間的正交性，構(gòu)造空間譜函數(shù)，從而進(jìn)行譜峰搜索來估計(jì)信號頻率.傳統(tǒng)的MUSIC譜估計(jì)方法假定背景噪聲滿足高斯分布，并采用二階統(tǒng)計(jì)量方法處理，但在地球物理、雷達(dá)、水聲信號等領(lǐng)域中所涉及的噪聲，往往不滿足高斯分布且具有顯著的尖峰脈沖和較厚的拖尾.因此在此種情況下高斯分布模型不再適用，α穩(wěn)定分布模型卻能很好的描述這類伴有顯著尖峰脈沖和較厚拖尾的噪聲［1－5］.α 穩(wěn)定分布由于不存在二階及二階以上統(tǒng)計(jì)量，所以基于自相關(guān)矩陣特征分解的傳統(tǒng)MUSIC譜估計(jì)方法已不再適用，需要研究適用于上述脈沖噪聲環(huán)境的新的譜估計(jì)方法.本文首先介紹了α穩(wěn)定分布信號的數(shù)學(xué)模型和傳統(tǒng)的MUSIC譜估計(jì)方法，然后依據(jù)α穩(wěn)定分布信號的特性，提出了兩種改進(jìn)的MUSIC譜估計(jì)方法.

1 α穩(wěn)定分布模型

如果隨機(jī)變量X具有式(1)所示的特征函數(shù)，則隨機(jī)變量X服從穩(wěn)定分布.

式中，sgn(·)為符號函數(shù);α為特征指數(shù)0＜α≤2，α值越小表明分布的拖尾越厚，尖峰脈沖越顯著，α=2時(shí)滿足高斯分布;β為對稱參數(shù)－1≤β≤1，當(dāng)β=0稱為對稱α穩(wěn)定分布(SαS);γ為分散系數(shù)γ≥0，類似高斯分布中的方差;a為位置參數(shù)，對于SαS分布當(dāng)0＜α≤1時(shí)，a為中值;當(dāng)1＜α≤2 時(shí)，a 為均值［2］.

2 傳統(tǒng)的MUSIC譜估計(jì)方法

假設(shè)信號x(n)是復(fù)正弦信號加白噪聲，為

式中，αk=|αk|ejφk和ωk分別表示信號的復(fù)幅度和角頻率.初始相位φk是在［0，2π］均勻分布的隨機(jī)變量，并且當(dāng)i≠k時(shí)，φi和φk是相互獨(dú)立的;v(n)是均值為0、方差為的白噪聲，且與信號相互獨(dú)立.

定義信號向量

則由式(2)有

式中，

向量a(ω)、S(n)和V(n)分別定義為

向量X(n)的自相關(guān)矩陣為

式中，I是一個(gè)M×M維的單位矩陣;P是信號S(n)的自相關(guān)矩陣.

對R進(jìn)行特征值分解，有

式中，US是由大特征值對應(yīng)的特征向量張成的子空間，即是信號子空間，UN是由小特征值對應(yīng)的特征向量張成的子空間，即是噪聲子空間.在理想條件下，信號頻率向量a(ωk)與噪聲子空間的特征向量正交，即

由噪聲子空間的向量構(gòu)成矩陣

信號角頻率的估計(jì)可以由函數(shù)PMUSIC(ω)的K個(gè)峰值位置確定.

譜函數(shù)PMUSIC(ω)的峰值的位置反映了信號的頻率值，但是它并非信號的功率譜，一般將PMUSIC(ω)稱為 MUSIC 譜［6］.

3 改進(jìn)的MUSIC譜估計(jì)方法

傳統(tǒng)MUSIC方法是通過對數(shù)據(jù)序列X(n)的自相關(guān)矩陣R進(jìn)行特征分解，然后由譜函數(shù)的峰值位置搜索來進(jìn)行譜估計(jì).但是當(dāng)X(n)是含有α穩(wěn)定分布噪聲的復(fù)正弦信號時(shí)，數(shù)據(jù)序列X(n)的自相關(guān)矩陣不存在，這種基于自相關(guān)矩陣分解的MUSIC譜估計(jì)方法在這種情況下則不再適用.

3.1 FLOC-MUSIC譜估計(jì)方法

由于α穩(wěn)定分布隨機(jī)變量不存在二階及二階以上的統(tǒng)計(jì)量，因此只能采用分?jǐn)?shù)低階統(tǒng)計(jì)量對含有α穩(wěn)定分布的信號進(jìn)行分析和處理.對于SαS分布隨機(jī)變量X和Y，滿足0＜α≤2，則X和Y的分?jǐn)?shù)低階協(xié)方差(FLOC)定義［7］為

可以看出，分?jǐn)?shù)低階協(xié)方差相當(dāng)于對 x(n)和y(n)先進(jìn)行非線性預(yù)處理，然后再計(jì)算互相關(guān).

文獻(xiàn)［8］對上述非線性處理過程進(jìn)行了數(shù)學(xué)證明，結(jié)論為:如果x(n)和y(n)為μ=0的SαS分布過程，則(x(n))＜A＞和(y(n))＜A＞是滿足零均值概率密度函數(shù)的對稱分布，且當(dāng)式(13)中的約束條件滿足時(shí)具有有限的二階矩.

本文依據(jù)α穩(wěn)定分布信號的特性，利用分?jǐn)?shù)低階協(xié)方差代替MUSIC譜估計(jì)方法中的相關(guān)，可以得到一種改進(jìn)的MUSIC譜估計(jì)方法(FLOC－MUSIC).下面給出FLOC－MUSIC譜估計(jì)方法步驟:

步驟1 由X(n)的N個(gè)觀測樣本，先對其進(jìn)行如式(14)的非線性變換得到C(n)，然后構(gòu)造數(shù)據(jù)序列C(n)的自相關(guān)矩陣RCC.

步驟2 對RCC進(jìn)行特征值分解，得到M－K個(gè)最小特征值對應(yīng)的特征向量，得到噪聲子空間的向量，根據(jù)式(10)構(gòu)造矩陣G.

步驟3 定義信號向量a(ω)=［1e－jωi…e－j(K－1)ωi］T，i=1，2，…，M，利用式(12)進(jìn)行譜估計(jì)，峰值位置就是待求復(fù)正弦信號的頻率.

3.2 FLOC-ROOT-MUSIC譜估計(jì)方法

對于長序列信號源，傳統(tǒng)的MUSIC譜估計(jì)方法能很好進(jìn)行頻率估計(jì)，但當(dāng)信號源是短序列時(shí)，傳統(tǒng)的MUSIC譜估計(jì)方法性能退化.求根MUSIC譜估計(jì)方法將MUSIC方法轉(zhuǎn)化為一種多項(xiàng)式求根形式，能對短序列信號源進(jìn)行很好地估計(jì)［6］.

定義向量a(z)為

當(dāng)z=ejω時(shí)，向量a(z)是頻率為ω的信號頻率向量.定義多項(xiàng)式

利用噪聲子空間的向量構(gòu)成矩陣G，式(17)可以表示為

由(17)可 知 zk=ejωk，k=1，…，K 是 方 程PROOT－MUSIC(z)=0的根.由于與復(fù)正弦信號頻率有關(guān)的K個(gè)根zk=ejωk，k=1，…，K都位于單位圓|z|=1上，單位圓上的復(fù)數(shù)z=ejω應(yīng)滿足z*=z－1，所以有

將式(19)代入式(18)，得到修正方程

于是信號頻率估計(jì)問題轉(zhuǎn)化為一元高次方程求根問題，這種方法稱為ROOT－MUSIC譜估計(jì)方法.

若信號源含有α穩(wěn)定分布噪聲，由于信號的自相關(guān)矩陣不存在，上述ROOT－MUSIC譜估計(jì)方法同樣會(huì)出現(xiàn)退化現(xiàn)象.利用分?jǐn)?shù)低階協(xié)方差代替ROOT－MUSIC譜估計(jì)方法中的自相關(guān)，可以得到另一種改進(jìn)的譜估計(jì)方法(FLOC－ROOTMUSIC).下面給出 FLOC－ROOT－MUSIC 譜估計(jì)方法步驟:

步驟1 由X(n)的N個(gè)觀測樣本，先對其進(jìn)行如式(14)的非線性變換得到C(n)，然后如式(15)構(gòu)造數(shù)據(jù)序列C(n)的自相關(guān)矩陣RCC.

步驟2 對RCC進(jìn)行特征值分解，得到M－K個(gè)最小特征值對應(yīng)的特征向量，得到噪聲子空間的向量，根據(jù)式(10)構(gòu)造矩陣G.

步驟3 求解式(20)，找出最接近單位圓的K個(gè)根，這些根的相位就是信號頻率的估計(jì).

4 仿真分析

設(shè)v(n)為服從α穩(wěn)定分布噪聲序列，信號

L為歸一化頻率.觀測序列為

X(n)=A*S(n)+V(n)，其中A為M×K維矩陣，對帶α穩(wěn)定分布噪聲的信號分別用傳統(tǒng)MUSIC譜估計(jì)方法和本文提出的改進(jìn)MUSIC譜估計(jì)方法進(jìn)行頻率估計(jì)，并比較三種方法的性能.信號歸一化頻率fi=，其中f1=0.1，f2=0.2，M=10.

4.1 高斯噪聲和α穩(wěn)定分布噪聲環(huán)境下三種算法的性能比較

采樣信號長度L=2 000，設(shè)定混合信噪比MSNR=0 dB，α分別設(shè)定為2、1.2，獨(dú)立運(yùn)行20次結(jié)果.仿真結(jié)果分別如圖1～圖3所示.

圖1 傳統(tǒng)MUSIC譜估計(jì)結(jié)果

圖2 FLOC-MUSIC譜估計(jì)結(jié)果

圖3 FLOC-ROOT-MUSIC譜估計(jì)結(jié)果

從圖1～圖3中可以看出，當(dāng)α=2時(shí)傳統(tǒng)MUSIC譜估計(jì)方法、FLOC－MUSIC譜估計(jì)方法及FLOC－ROOT－MUSIC譜估計(jì)方法都能很好進(jìn)行頻率估計(jì);但當(dāng)α=1.2時(shí)，由于不滿足高斯分布，傳統(tǒng)MUSIC譜估計(jì)方法不再適用估計(jì)性能嚴(yán)重退化，而 FLOC－MUSIC譜估計(jì)方法及 FLOCROOT－MUSIC譜估計(jì)方法仍能很好進(jìn)行頻率估計(jì).

4.2 不同混合信噪比下三種算法的性能比較

采樣信號長度L=2000，設(shè)定α =1.2，混合信噪比MSNR 分別設(shè)定為 －10、－5、0、5 和10 dB，獨(dú)立運(yùn)行20次，計(jì)算每種方法的均值和標(biāo)準(zhǔn)差.仿真結(jié)果如表1所示.

表1 不同信噪比下三種算法的均值和標(biāo)準(zhǔn)差比較

由表1可以看出，由于α=1.2，不滿足高斯分布，所以基于相關(guān)矩陣分解的傳統(tǒng)MUSIC譜估計(jì)方法的估計(jì)效果很差，而FLOC－MUSIC譜估計(jì)方法及FLOC－ROOT－MUSIC譜估計(jì)方法能對于不同混合信噪比下的帶噪信號進(jìn)行有效的頻率估計(jì)，混合信噪比越大，估計(jì)的均值越接近于真實(shí)值.

4.3 不同α值下三種算法的性能比較

采樣信號長度L=2 000，設(shè)定MSNR=0 dB，分別設(shè)定 α 為0.8、1.1、1.4、1.7、和 2，獨(dú)立運(yùn)行20次結(jié)果，計(jì)算每種方法的均值和標(biāo)準(zhǔn)差.仿真結(jié)果如表2所示.

由表2可以看出，當(dāng)α=2時(shí)，傳統(tǒng)MUSIC譜估計(jì)方法和FLOC－MUSIC譜估計(jì)方法及 FLOCROOT－MUSIC譜估計(jì)方法均能很好的進(jìn)行頻率估計(jì).但是當(dāng)α＜2時(shí)，傳統(tǒng)MUSIC譜估計(jì)方法估計(jì)性能變差，而 FLOC－MUSIC譜估計(jì)方法及FLOC－ROOT－MUSIC譜估計(jì)方法仍能很好的進(jìn)行頻率估計(jì).

表2 不同α值時(shí)三種算法的均值和標(biāo)準(zhǔn)差比較

4.4 不同信號長度時(shí) FLOC-MUSIC算法和FLOC-ROOT-MUSIC算法性能比較

設(shè)定MSNR=0 dB、α =1.1時(shí)，分別取采樣信號長度L=2 000和L=200，采用FLOC－MUSIC譜估計(jì)方法及FLOC－ROOT－MUSIC譜估計(jì)方法進(jìn)行頻率估計(jì)，獨(dú)立運(yùn)行20次結(jié)果，計(jì)算每種方法的均值和標(biāo)準(zhǔn)差.仿真結(jié)果如表3、表4所示.

由表3、表4可以看出，對于長序列的信號源，F(xiàn)LOC－MUSIC 譜估計(jì)方法及 FLOC－ROOT－MUSIC譜估計(jì)方法均能很好的進(jìn)行信號頻率估計(jì);對于短序列，當(dāng)混合信噪比較低時(shí)，F(xiàn)LOC－ROOTMUSIC譜估計(jì)方法的性能明顯好于FLOC－MUSIC譜估計(jì)方法.

表3 α=1.1、L=2 000時(shí)不同MSNR的仿真結(jié)果

表4 α=1.1、L=200時(shí)不同MSNR的仿真結(jié)果

5 結(jié)論

由于不存在二階統(tǒng)計(jì)量，所以在脈沖噪聲環(huán)境下，傳統(tǒng)MUSIC譜估計(jì)方法效果變差，本文利用分?jǐn)?shù)低階統(tǒng)計(jì)量的處理方法，對帶噪信號進(jìn)行了非線性處理，提出了改進(jìn)的MUSIC譜估計(jì)方法.仿真結(jié)果表明，F(xiàn)LOC－MUSIC譜估計(jì)方法和FLOC－ROOT－MUSIC譜估計(jì)方法在高斯噪聲和脈沖分布噪聲環(huán)境下均具有良好的性能.對于帶有α穩(wěn)定分布噪聲的短序列，當(dāng)混合信噪比較低時(shí)，F(xiàn)LOC－ROOT－MUSIC譜估計(jì)方法的性能明顯優(yōu)于FLOC－MUSIC譜估計(jì)方法.

［1］NIKIAS C L，SHAO M.Signal Processing with Alpha－Stable Distribution and Applications［M］.New York:John Wiley＆Sons，Inc，1995.

［2］SHAO M，NIKIAS C L.Signal processing with fractional lower order moments:stable processes and their applications［J］.Proceedings of the IEEE，1993，81(7):987－1010.

［3］MOUNIR D，MESSAOUD B.Robust Polynomial Wigner－Ville Distribution for the analysis of polynomial phase signals in α －stable noise［J］.IEEE international conference，2004(2):613－616.

［4］MA X，NIKIAS C L.Parameter estimation and bland channel identification in impulsive signal environments［J］.IEEE Trans on signal processing，1995，43(11):2884－2897.

［5］孫永梅.穩(wěn)定分布參數(shù)估計(jì)與譜分析理論及應(yīng)用研究［D］.大連:大連理工大學(xué)，2006.

［6］何子述.現(xiàn)代數(shù)字信號處理及其應(yīng)用［M］.北京:清華大學(xué)出版社，2009:110－116.

［7］MA X，NIKIAS C L.Joint estimation of time delay and frequency delay in impulsive noise Using fractional lower order statistics［J］.IEEE Transon Signal Processing，1996，44(11):2669－2687.

［8］QIU T，LI X，WANG H.An AFLC algorithm for latency change estimation of EP under alpha－stable noise condition［C］.Proce of the 25th Annual International Conference of the IEEE on Engineering in Medicine and Biology Society，Cancun，Mexico，2003:2635－2638.