MUSIC算法諧波恢復(fù).doc

《MUSIC算法諧波恢復(fù).doc》由會(huì)員上傳分享，免費(fèi)在線閱讀，更多相關(guān)內(nèi)容在工程資料-天天文庫(kù)。

1、MUSIC算法諧波恢復(fù)1.原理設(shè)某個(gè)天線陣列具有M個(gè)陣元，且陣元按照等距直線排列，相鄰兩個(gè)陣元之間的距離為d,天線陣接收到的諧波信號(hào)為，到達(dá)每個(gè)陣元的波達(dá)方向角為，其中P為遠(yuǎn)場(chǎng)信號(hào)的個(gè)數(shù)。則第k個(gè)陣元接收到的第i個(gè)信號(hào)可表示為:，則由于信號(hào)到第k+l個(gè)陣元的傳播距離比到第k個(gè)長(zhǎng)，第k+1個(gè)陣元接收的信號(hào)與第k個(gè)陣元接收的信號(hào)與第k個(gè)陣元接收的信號(hào)的相位差為:其中是相位傳播因子，表示波長(zhǎng).因此，信號(hào)到M個(gè)陣元的相位差所組成的向量為:因此，當(dāng)有p個(gè)信號(hào)位于遠(yuǎn)場(chǎng)時(shí)，且p不大于M，該等距直線陣上各陣元在觀測(cè)加性高斯白噪聲(AWGN,additiveGaussi

2、anwhitenoise)下所接收到的信號(hào)為:用向量形式表示為：（1）其中：在陣列信號(hào)處理中，一次采樣成為一次快拍，假定在每個(gè)陣元上共觀測(cè)到N次快拍的接受信號(hào)Xk(1),….,Xk(n),其中，k=1,…,m?？紤]使N次快拍的輸出能量最小，即（2）設(shè)計(jì)權(quán)向量w的這一準(zhǔn)則稱為最小輸出能量（MOE）準(zhǔn)則，它在通信信號(hào)處理，雷達(dá)信號(hào)處理等中具有重要的應(yīng)用。令（3）經(jīng)過(guò)推導(dǎo)可以得出一種類似于功率譜的函數(shù)：（4）（4）式被稱為信號(hào)子空間算法，相應(yīng)的下式是被稱為噪聲子空間算法：（5）2MUSIC算法步驟1計(jì)算樣本自相關(guān)矩陣Rxx的特征值分解，得到其主特征，….，和

3、次特征值，并儲(chǔ)存主特征值是S1,S2,….Sp。步驟2利用（4）式或（5）式計(jì)算MUSIC譜P(wi)，其中wi=(i-1)w，網(wǎng)格w可取作2∏0.001等。步驟3找出P(w)的p個(gè)峰值，它們就是待求的MUSIC估計(jì)值w1,w2,…wp。3Matlab仿真觀測(cè)信號(hào)為：按照MUSIC算法的步驟，將其代入（3）式計(jì)算出Rxx，對(duì)Rxx進(jìn)行特征值分解，再代入（4）或（5）式即可得出功率譜表達(dá)式。按照上述方法編程進(jìn)行Matlab仿真得出的仿真波形為：m=40,N=80m=30,N=70M=20,N=504功率譜分析通過(guò)以上的理論要求，我們分別設(shè)置m個(gè)等距天線陣

4、元并對(duì)觀測(cè)信號(hào)進(jìn)行N次快拍，對(duì)觀測(cè)信號(hào)進(jìn)行估計(jì)，可以看出，功率譜的分辨率隨著m和N的取值大小而變化，m和N取值越大，分辨率越大；m和N取值偏小時(shí)則得不到想要的分辨率。因此，只要取合適的m和N，我們就可以得到分辨率滿足要求的功率譜波形，進(jìn)而可以恢復(fù)出諧波。附程序：clearalln=1:128;N1=128;m=30;%p=2;N=50;%設(shè)置觀測(cè)數(shù)據(jù)x=sqrt(20)*sin(2*pi*0.2*n)+sqrt(2)*sin(2*pi*0.213*n)+randn(1,N1);%建立自相關(guān)矩陣Rxx=zeros(m);fork=1:N;fori=1:m

5、;xm(i)=x(i+k-1);endrxx=xm'*xm;Rxx=Rxx+rxx;endRx=Rxx/N;%做特征值分解[V,D]=eig(Rx);p=4;%對(duì)特征值和特征向量做處理fori=1:m;d(i)=D(m+1-i,m+1-i);endV1=zeros(m);fori=1:m;V1(1:m,m+1-i)=V(1:m,i);ends=zeros(m,p);fori=1:p,s(1:m,i)=V1(1:m,i);end%用信號(hào)子空間公式求功率譜df=0.001;I=eye(m);f=0.001:0.001:0.5;fori=1:length(f

6、);wi=2*pi*df*i;fort=1:m;aw(t)=exp(-j*(t-1)*wi);endp(i)=1/(aw*(I-s*s')*aw');endpw=10*log10(abs(p));plot(f,pw);