MUSIC算法頻率估計

1、采用MUSIC方法的白噪聲頻率檢測仿真本試驗提供了一種使用MUSIC方法的白噪聲中一個正弦信號和M個正弦信號的特征分解頻率估計的仿真試驗，并討論了虛假峰的成因并給出了實驗證明。問題描述假定仿真的觀測數(shù)據(jù)分別由單個正弦信號檢測的情況vnn、x(n)=4ej(4n+3)+u(n)多個正弦信號檢測的情況nnnn5nn、x(n)=4ej(4n+3)+2ej(3n+6)+3ej(4n+5)+u(n)產(chǎn)生，其中u(n)是一高斯白噪聲，其均值為0,方差為1。用MUSIC方法估計觀測數(shù)據(jù)中正弦波的頻率，并給出白噪聲方差Q2與復(fù)正弦波的振幅A的估計值。u多重信號分類的MUSIC方法實際應(yīng)用中常常需要對空間中存在

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

3、，是空間譜估計技術(shù)的典型代表。其測向原理是根據(jù)矩陣特征分解的理論，對陣列輸出協(xié)方差矩陣進(jìn)行特征分解，將信號空間分解為噪聲子空間G和信號子空間S,利用噪聲子空間G與陣列的方向矩陣A的列矢量正交的性質(zhì)，構(gòu)造空間譜函數(shù)P(w)并進(jìn)行譜峰搜索，從而估計出波達(dá)方向信息。設(shè)空間有p個互不相關(guān)的信號以方位角9,9,9入射到具有m個接收陣元12p的接收陣元陣列中，入射信號的數(shù)目p小于陣列的陣元數(shù)m。則此陣列系統(tǒng)的信號模型為：x(n)=區(qū)a(w)s(n)+u(n)=A(w)s(n)+u(n)iii=1其中A(w)=a(w),a(w)1p_111_e-jWe-丿叫e-jwpe-j(m-i)We-j(m-i)w?e

4、j(m-1)wps(n)=s(n),s(n)1p對由上式描述的陣列信號觀測模型做以下假設(shè)：假設(shè)1：對于不同的w值，向量a(w)相互線性獨立；ii假設(shè)2：加性噪聲向量u(n)的每個原色都是零均值的復(fù)白噪聲，它們不相關(guān)，并且具有相同的方差a2；u假設(shè)3：矩陣非奇異P=Es(n)sh(n)，即rank(P)=p。上述三個假設(shè)都只是一般的假設(shè)，在實際中容易得到滿足?？梢?，R=Ex(n)xH(n)xx=A(w)Es(n)sH(n)AH(w)+a2Iu=APAH+a2Iu于是R是個對稱矩陣，令其特征值分解為xxR=UZUhxx式中Z=diag(a2,.,a2)。1m由于A滿列秩，故rank(APAh)=r

5、ank(P)=p，這里假定pm。于是有UhAPAhU=diag(a2,a2,0,0)1p其中a2,a2是無加性噪聲時的觀測信號Ax(n)的自相關(guān)矩陣APAh的特征值。1p同時用Uh左乘和U右乘上式，得UhRU=UhAPAhU+a2UhUxxu=diag(a2,a2,0，,0)+a211pu這表明，自相關(guān)矩陣R的特征值為xxra2+a2,i=1,p九=O2=iuio2i=p+1,mu根據(jù)信號特征值和噪聲特征值，將特征矩陣U的列向量分成兩部分，即U=S#G其中S=s,s=u,u1p1pG=g,g=u,u1m-pp+1m其中分別由信號特征向量和噪聲特征向量組成。考察RG=S#GZxxSHGH=O2G

6、u又由R=APAh+o21有RG=APAhG+o2G,利用上式的結(jié)果，得到xxuxxuAPAHG=O進(jìn)而有GHAPAHG=O上式成立的充分必要條件是AHG=O將A=a(w),a(w)要代入上式得1paH(w)G=0T,w=w,w1p顯然，當(dāng)w豐w，,w時，ah(w)G豐0t。1p將上式改成標(biāo)量形式，可以定義一種類似于功率譜的函數(shù)1ah(w)GGha(w)上式取峰值的p個w值給出p個信號的波達(dá)方向9,9,9。12p這樣定義的函數(shù)描述了空間參數(shù)(即波達(dá)方向)的分布，因此成為空間譜。由于它能夠?qū)Χ鄠€信號進(jìn)行識別，故以此得名MUSIC方法(multiplesignalclassification)。白

7、噪聲中單個正弦信號的頻率檢測與估計1實驗結(jié)果(1)短數(shù)據(jù)時的估計結(jié)果，N_x=64.S=1.4058A=3.9835(2)長數(shù)據(jù)時的估計結(jié)果，N_x=1024。噪聲功率)信號振幅估計)wl_estimate=0.2500(n)(估計頻率)0.9646(噪聲功率)00pi/42對結(jié)果的討論pi/23pi/4(信號振幅估計)5pi/43pi/27pr/42pi(1)由給出的試驗結(jié)果可知，在n/4處有明顯的峰值，其它地方都是平坦的，頻率估計效果非常好。(2)在信號長度很短(64)的情況下，a2和A的估計誤差比較大。造成誤差大的原因是：由于信號短，加窗效應(yīng)較明顯，自相關(guān)矩陣的估計誤差比較大。增加信號長

8、度可以減小自相關(guān)矩陣的估計誤差，可以得到比較理想的結(jié)果，這可以從(N_x=1024)的結(jié)果中看出。5pi/100白噪聲中多個正弦信號的頻率檢測與估計1實驗結(jié)果（1）短數(shù)據(jù)時的估計結(jié)果，N_x=64.700600500400（噪聲功率）300-S=2.1681(2)A=3.88722.18742.910短數(shù)據(jù)時的估計結(jié)果，N_x=1（24.（信號振幅估計）0pi/4pi/23pi/4Pi5p80007000600050004000S=1.03853000-A=4.0098(噪聲功率)2.02182.9917(信號振幅估計)2對結(jié)果的討論(1)由給出的試驗結(jié)果可知，在n/4,n/3,5*n/4處分

9、別都有明顯的峰,與實際值十分吻合。為了予省篇幅沒有趣的可以分段找最大值的方法將估計值與實際值十分吻合計算估計值，只從圖上觀察是否相符，有興5pi/1005pi/100由信號長度導(dǎo)致的噪聲功率信號振幅估計的差別在上述結(jié)果中也有明顯的體現(xiàn)，產(chǎn)生的原因同11L0pi/4pi/23pi/4pi疑問：從圖中發(fā)現(xiàn)三個峰高度有明顯的差異，開始以為與Ai或有關(guān)，經(jīng)過實驗并未發(fā)現(xiàn)有何聯(lián)系，未找到影響峰值的原因。100虛假峰成因與數(shù)量的討論由于在MUSIC方法中，前面的幾個噪聲特征值對應(yīng)的特征向量都是相同的，于是嘗試用單個特征矢量進(jìn)行頻率估計。在用單個特征矢量估計31的時候會出現(xiàn)虛假峰，即不等于31的位置上多出幾

10、個峰。研究特征向量的正交性可知，特征矢量是兩兩正交的。下面的程序可以證實：對特征分解得到的特征矢量組成的矩陣V進(jìn)行處理，判斷列矢量之間，行矢量之間是否正交：%Test_V.m查看特征矢量是否互相垂直。z=zeros(6);結(jié)果為：z=100000fori=1:6010000forj=i:6z(i,j)=V(:,i)*V(:,j);z(j,i)=V(j,:)*V(i,:);001000000100000010end000001end。這說明特征矢量矩陣的行列矢量(即特征矢量)之間，行矢量之間都是兩兩正交的。特征矢量兩兩正交，也就是說eHVi應(yīng)該有N-1個零點，因此應(yīng)該有N-1個峰值，這是虛假峰產(chǎn)

11、生的真正原因！為此特作了多次實驗加以驗證，模型仍用前面1個復(fù)正弦信號，31=n/4，N=6，結(jié)果如下圖所示：mymusics.m白噪聲中單個正弦信號的頻率檢測與估計clear;closeall;%FrequencyEstimationbyEigendecompositionofAutocorrelationMatrixN_x=1024;%LengthofSignalN=10;%SizeofRxMatrixA=423;w=pi/4pi/35*pi/4;phase=pi/3*ones(1,N_x);pi/6*ones(1,N_x);pi/5*ones(1,N_x);M=3;%NumberofSig

12、nalsx=randn(1,N_x)+A*exp(j*(w*0:N_x-1+phase);Cx=xcorr(x,biased);Rxx=Cx(N_x:N_x+N-1);Rx=toeplitz(Rxx);V,D=eig(Rx);%Eigendecomposition特征分解D=sum(D);Nw=128;ww=0:256/128*pi;e=exp(-j*ww*0:N-1);ev=e*V(:,1:N-M);Pw=1./real(diag(ev*ev);figureplot(ww,Pw);xlim(02*pi)set(gca,XTick,0:pi/4:2*pi)set(gca,XTickLabel,

13、0,pi/4,pi/2,3pi/4,pi,5pi/4,3pi/2,7pi/4,2pi)S=mean(D(1:N-M)E=exp(j*0:N-1.*w);P=real(E(Rx-(eye(N).*S)/E);A=sqrt(sum(P)mymusicm.m白噪聲中多個正弦信號的頻率檢測與估計clear;closeall;%FrequencyEstimationbyEigendecompositionofAutocorrelationMatrixN_x=1024;%LengthofSignalN=6;%SizeofRxMatrixw1=pi/4;x=4*exp(j*(w1*0:N_x-1+pi/3)+randn(1,N_x);%Generatethesignal:Cx=xcorr(x,biased);Rxx=Cx(N_x:N_x+N-1);Rx=toeplitz(Rxx);V,D=eig(Rx);%Eigendecomposition特征分解D=sum(D);Nw=128;w=0:256/128*pi;e=exp(-j*w*0:N-1);ev=e*V(:,1:N-1);Pw=1./real(diag(ev*ev);plot(w,Pw);xlim(02*pi)set(gca,XTick,0:pi/4:2*pi)se