MUSIC算法頻率估計(jì)

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

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

3、，是空間譜估計(jì)技術(shù)的典型代表。其測(cè)向原理是根據(jù)矩陣特征分解的理論，對(duì)陣列輸出協(xié)方差矩陣進(jìn)行特征分解，將信號(hào)空間分解為噪聲子空間G和信號(hào)子空間S,利用噪聲子空間G與陣列的方向矩陣A的列矢量正交的性質(zhì)，構(gòu)造空間譜函數(shù)P(w)并進(jìn)行譜峰搜索，從而估計(jì)出波達(dá)方向信息。設(shè)空間有p個(gè)互不相關(guān)的信號(hào)以方位角9,9,9入射到具有m個(gè)接收陣元12p的接收陣元陣列中，入射信號(hào)的數(shù)目p小于陣列的陣元數(shù)m。則此陣列系統(tǒng)的信號(hào)模型為：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對(duì)由上式描述的陣列信號(hào)觀測(cè)模型做以下假設(shè)：假設(shè)1：對(duì)于不同的w值，向量a(w)相互線性獨(dú)立；ii假設(shè)2：加性噪聲向量u(n)的每個(gè)原色都是零均值的復(fù)白噪聲，它們不相關(guān)，并且具有相同的方差a2；u假設(shè)3：矩陣非奇異P=Es(n)sh(n)，即rank(P)=p。上述三個(gè)假設(shè)都只是一般的假設(shè)，在實(shí)際中容易得到滿足?？梢?，R=Ex(n)xH(n)xx=A(w)Es(n)sH(n)AH(w)+a2Iu=APAH+a2Iu于是R是個(gè)對(duì)稱矩陣，令其特征值分解為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是無加性噪聲時(shí)的觀測(cè)信號(hào)Ax(n)的自相關(guān)矩陣APAh的特征值。1p同時(shí)用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ù)信號(hào)特征值和噪聲特征值，將特征矩陣U的列向量分成兩部分，即U=S#G其中S=s,s=u,u1p1pG=g,g=u,u1m-pp+1m其中分別由信號(hào)特征向量和噪聲特征向量組成?？疾霷G=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時(shí)，ah(w)G豐0t。1p將上式改成標(biāo)量形式，可以定義一種類似于功率譜的函數(shù)1ah(w)GGha(w)上式取峰值的p個(gè)w值給出p個(gè)信號(hào)的波達(dá)方向9,9,9。12p這樣定義的函數(shù)描述了空間參數(shù)(即波達(dá)方向)的分布，因此成為空間譜。由于它能夠?qū)Χ鄠€(gè)信號(hào)進(jìn)行識(shí)別，故以此得名MUSIC方法(multiplesignalclassification)。白

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

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

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

10、個(gè)峰。研究特征向量的正交性可知，特征矢量是兩兩正交的。下面的程序可以證實(shí)：對(duì)特征分解得到的特征矢量組成的矩陣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。這說明特征矢量矩陣的行列矢量(即特征矢量)之間，行矢量之間都是兩兩正交的。特征矢量?jī)蓛烧?，也就是說eHVi應(yīng)該有N-1個(gè)零點(diǎn)，因此應(yīng)該有N-1個(gè)峰值，這是虛假峰產(chǎn)

11、生的真正原因！為此特作了多次實(shí)驗(yàn)加以驗(yàn)證，模型仍用前面1個(gè)復(fù)正弦信號(hào)，31=n/4，N=6，結(jié)果如下圖所示：mymusics.m白噪聲中單個(gè)正弦信號(hào)的頻率檢測(cè)與估計(jì)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白噪聲中多個(gè)正弦信號(hào)的頻率檢測(cè)與估計(jì)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