music 方位估計 實驗報告三 MATLAB 代碼

1、實驗報告三 實驗目的： 實現(xiàn)常規(guī)波束形成及基于MUSIC方法的方位估計。實驗內(nèi)容： 1）若干陣元的接收陣，信號頻率為10KHz，波束主軸12度，仿真給出常規(guī)波束形成的波束圖。 2）16個陣元的均勻線列陣，信號頻率為10KHz，信號方位為12度，用MUSIC方法完成目標定向，信噪比-5dB，0dB，5dB。 i) 波束形成時的陣型設計為兩種，一種是均勻線列陣，陣元16個；一種是均勻圓陣，陣元數(shù)為16個，比較這兩種陣型的波束圖。ii）比較不同信噪比下MUSIC方法估計的性能（統(tǒng)計100次）。實驗原理： i)常規(guī)波束形成： 圖 1 波束形成器基本原理圖 如圖所示，基陣的輸出。 采用向量符號則有， 式

2、中，x(t)和w(q )分別為觀測數(shù)據(jù)向量和加權(quán)系數(shù)向量，基陣輸出端的空間功率譜表示為：式中，R為觀測數(shù)據(jù)的協(xié)方差矩陣。 ii）基于MUSIC方法的方位估計： 假設: (1 ) 信號源的數(shù)目d 是已知的, 且d < M; (2 ) 各信號的方向矢量是相互獨立的, 即是一個列滿秩矩陣; (3 ) 噪聲是空間平穩(wěn)隨機過程, 為具有各態(tài)歷經(jīng)性的均值為零、方差為2n 的高斯過程; (4 ) 噪聲各取樣間是統(tǒng)計獨立的。 在上述假設條件下, 基陣輸出的協(xié)方差矩陣可表示為： 其中, Rs 為信號的協(xié)方差矩陣；I 為單位矩陣。對R 進行特征分解, 并以特征值降值排列可得 信號子空間與噪聲子空間正交。 若

3、噪聲子空間記為EN , 即 則有：因此，對矩陣R 進行特征分解后, 取噪聲特征矢量, 就可得到信號源的到達角。 定義多重信號分類法的空間譜函數(shù)為：實驗過程及結(jié)論：1.常規(guī)波束形成，陣列為均勻線列陣： 2.常規(guī)波束形成，陣列為均勻圓陣： 結(jié)論：由上面兩個仿真結(jié)果可以看出，采用均勻線列陣和均勻圓陣都可以完成波束形成，都可以獲得陣處理增益。3. MUSIC算法方位估計： 結(jié)論：從上圖得出，信噪比越高，MUSIC方位估計越準確。附錄：MATLAB仿真程序代碼1.常規(guī)波束形成%clear all;close all;clc%常規(guī)波束形成 %信號參數(shù)f0=10e3;fs=10*f0;T=2e-3;t=1:

4、1/fs:T-1/fs; %均勻線列陣陣列參數(shù)c=1500;lamda=c/f0;d=lamda/2;M=16;direct=12;%波束指向 %陣列流行theta=-90:90;for i_theta=1:length(theta) fai=2*pi*d*sin(theta(i_theta)*pi/180)/lamda; for i_M=1:M A(i_M,i_theta)=exp(-sqrt(-1)*(i_M-1)*fai); endend %加權(quán)向量wfai_direct=2*pi*f0*d*sin(direct*pi/180)/c;for i_M=1:M w(i_M)=exp(sqrt

5、(-1)*(i_M-1)*fai_direct); endw=w' %無噪聲輸出P=w'*A;P=20*log10(abs(P)/max(abs(P);figure(1)plot(theta,P);title('均勻線列陣 12度主軸');xlabel('方位角/度');ylabel('歸一化功率/dB');grid hold on clear all;clc%常規(guī)波束形成 %信號參數(shù)f0=10e3;fs=10*f0;T=2e-3;t=1:1/fs:T-1/fs; %元陣陣列參數(shù)c=1500;lamda=c/f0;d=lamda/

6、2;M=16;i_M=0:M-1;fai=2*pi/M.*i_M;B=d*cos(fai);d*sin(fai);d*zeros(1,M);direct=12;%波束指向 %陣列流行theta=-180:180;theta=theta*pi/180;for i_theta=1:length(theta) e=cos(theta(i_theta);sin(theta(i_theta);0; tao=-e.'*B./c; for i_M=1:M A(i_M,i_theta)=exp(sqrt(-1)*2*pi*f0*tao(i_M); endend %加權(quán)向量wfai_direct=2*p

7、i*f0*d*sin(direct*pi*180)/c;for i_M=1:M w(i_M)=exp(sprt(-1)*2*pi*f0*d*cos(fai(i_M)-direct/180*pi)/c);endw=w'%無噪聲輸出P=w'*A;P=20*log10(abs(P)/max(abs(P);figure(2)plot(theta./pi*180,P);title('均勻圓陣 12度主軸');xlabel('方位角/度');ylabel('歸一化功率/dB');grid 2.MUSIC方位估計% clc; clear; cl

8、ose all;M=16;c=1500;f0=10e3;lamda=c/f0;d=lamda/2;fs=10*f0;T=2;t=0:1/fs:T;N=fs*T;signalnum=1;%信噪比SNR=-5; theta1=12;num=100;theta1=theta1*pi/180;tau0=d*sin(theta1)/c;for i_M=1:M S(i_M,:)=exp(sqrt(-1)*2*pi*f0*t)*exp(-sqrt(-1)*2*pi*f0*(i_M-1)*tau0);endP=zeros(1,181); %length(theta)for i=1:num X=awgn(S,S

9、NR);R=X*X'/N; %協(xié)方差矩陣V,I=eig(R); Order Index=sort(diag(I),'descend'); %特征值降值排列VS=; %向量子空間P_music=;for i_M=1:M VS=VS,V(:,Index(i_M); %排序后的特征向量endE=VS(:,signalnum+1:M); %后M-signalnum向量構(gòu)成向量子空間En=E*E' %噪聲子空間theta=-90:1:90; fai=2*pi*f0*d*sin(theta*pi/180)/c;for i_theta=1:length(theta) for

10、i_M=1:M w(i_theta,i_M)=exp(sqrt(-1)*(i_M-1)*fai(i_theta);%方向矢量 end P_music(i,i_theta)=1./(w(i_theta,:)*En*En'*w(i_theta,:)'); endendP=abs(P)/num;P=P+P_music(i,:);figure(3)plot(theta,10*log10(P/max(P),'r');xlabel('方位角 / 度');ylabel('歸一化功率 / dB');grid title('不同信噪比下MU

11、SIC方法估計的性能（統(tǒng)計100次）')hold on %信噪比SNR=0; theta1=12;num=100;theta1=theta1*pi/180;tau0=d*sin(theta1)/c;for i_M=1:M S(i_M,:)=exp(sqrt(-1)*2*pi*f0*t)*exp(-sqrt(-1)*2*pi*f0*(i_M-1)*tau0);endP=zeros(1,181); %length(theta)for i=1:num X=awgn(S,SNR);R=X*X'/N; %協(xié)方差矩陣V,I=eig(R); Order Index=sort(diag(I),

12、'descend'); %特征值降值排列VS=; %向量子空間P_music=;for i_M=1:M VS=VS,V(:,Index(i_M); %排序后的特征向量endE=VS(:,signalnum+1:M); %后M-signalnum向量構(gòu)成向量子空間En=E*E' %噪聲子空間theta=-90:1:90; fai=2*pi*f0*d*sin(theta*pi/180)/c;for i_theta=1:length(theta) for i_M=1:M w(i_theta,i_M)=exp(sqrt(-1)*(i_M-1)*fai(i_theta);%方向矢

13、量 end P_music(i,i_theta)=1./(w(i_theta,:)*En*En'*w(i_theta,:)'); endendP=abs(P)/num;P=P+P_music(i,:);figure(3)plot(theta,10*log10(P/max(P),'*');xlabel('方位角 / 度');ylabel('歸一化功率 / dB');grid title('不同信噪比下MUSIC方法估計的性能（統(tǒng)計100次）')hold on %信噪比SNR=5; theta1=12;num=100;

14、theta1=theta1*pi/180;tau0=d*sin(theta1)/c;for i_M=1:M S(i_M,:)=exp(sqrt(-1)*2*pi*f0*t)*exp(-sqrt(-1)*2*pi*f0*(i_M-1)*tau0);endP=zeros(1,181); %length(theta)for i=1:num X=awgn(S,SNR);R=X*X'/N; %協(xié)方差矩陣V,I=eig(R); Order Index=sort(diag(I),'descend'); %特征值降值排列VS=; %向量子空間P_music=;for i_M=1:M V

15、S=VS,V(:,Index(i_M); %排序后的特征向量endE=VS(:,signalnum+1:M); %后M-signalnum向量構(gòu)成向量子空間En=E*E' %噪聲子空間theta=-90:1:90; fai=2*pi*f0*d*sin(theta*pi/180)/c;for i_theta=1:length(theta) for i_M=1:M w(i_theta,i_M)=exp(sqrt(-1)*(i_M-1)*fai(i_theta);%方向矢量 end P_music(i,i_theta)=1./(w(i_theta,:)*En*En'*w(i_theta,:)'); endendP=ab