music方位估计实验报告三MATLAB代码

实验报告三实验目的：实现常规波束形成及基于MUSIC方法的方位估计。实验内容：1）若干阵元的接收阵，信号频率为10KHz，波束主轴12度，仿真给出常规波束形成的波束图。2）16个阵元的均匀线列阵，信号频率为10KHz，信号方位为12度，用MUSIC方法完成目标定向，信噪比-5dB，0dB，5dB。i)波束形成时的阵型设计为两种，一种是均匀线列阵，阵元16个；一种是均匀圆阵，阵元数为16个，比较这两种阵型的波束图。ii）比较不同信噪比下MUSIC方法估计的性能（统计100次）。实验原理：i)常规波束形成：如图所示，基阵的输出),(tv。MmiiMmiiwtxtxwtv11)()()()(),(采用向量符号则有，)()()()(),(HHwxxwtttv式中，x(t)和w(q)分别为观测数据向量和加权系数向量，)(1w)(2w)(Mw)(1tx)(2tx)(txM),(tv图1波束形成器基本原理图TM21])()()([)(txtxtxtxTM21])()()([)(w基阵输出端的空间功率谱表示为：)()()()]()([)()]()()()([)],(),([]),([)(HHHHH*2RttEttEtvtvEtvEP式中，R为观测数据的协方差矩阵。ii）基于MUSIC方法的方位估计：)()()()(1tntsatxidiTM21])()()([)(txtxtxtx)()()()(tntsAtx假设:(1)信号源的数目d是已知的,且dM;(2)各信号的方向矢量是相互独立的,即)(A是一个列满秩矩阵;(3)噪声)(tn是空间平稳随机过程,为具有各态历经性的均值为零、方差为σ2n的高斯过程;(4)噪声各取样间是统计独立的。在上述假设条件下,基阵输出的协方差矩阵可表示为：IAARtxtxERHsH2])()([其中,Rs为信号的协方差矩阵；I为单位矩阵。对R进行特征分解,并以特征值降值排列可得HmmMdmmHmmdmmeeeeR11信号子空间与噪声子空间正交。若噪声子空间记为EN,即MdmHmmNeeE1则有：0)(iNaE因此，对矩阵R进行特征分解后,取噪声特征矢量,就可得到信号源的到达角。定义多重信号分类法的空间谱函数为：)()(1aEEaPHNNHMUSIC实验过程及结论：1.常规波束形成，阵列为均匀线列阵：-100-80-60-40-20020406080100-70-60-50-40-30-20-100均匀线列阵12度主轴方位角/度归一化功率/dB2.常规波束形成，阵列为均匀圆阵：-200-150-100-50050100150200-80-70-60-50-40-30-20-100均匀圆阵12度主轴方位角/度归一化功率/dB结论：由上面两个仿真结果可以看出，采用均匀线列阵和均匀圆阵都可以完成波束形成，都可以获得阵处理增益。3.MUSIC算法方位估计：-100-80-60-40-20020406080100-60-50-40-30-20-100方位角/度归一化功率/dB不同信噪比下MUSIC方法估计的性能（统计100次）-5dB0dB5dB结论：从上图得出，信噪比越高，MUSIC方位估计越准确。附录：MATLAB仿真程序代码1.常规波束形成%%clearall;closeall;clc%%%%%%常规波束形成%信号参数f0=10e3;fs=10*f0;T=2e-3;t=1:1/fs:T-1/fs;%均匀线列阵阵列参数c=1500;lamda=c/f0;d=lamda/2;M=16;direct=12;%波束指向%阵列流行theta=-90:90;fori_theta=1:length(theta)fai=2*pi*d*sin(theta(i_theta)*pi/180)/lamda;fori_M=1:MA(i_M,i_theta)=exp(-sqrt(-1)*(i_M-1)*fai);endend%加权向量wfai_direct=2*pi*f0*d*sin(direct*pi/180)/c;fori_M=1:Mw(i_M)=exp(sqrt(-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');gridholdonclearall;clc%%%%%%常规波束形成%信号参数f0=10e3;fs=10*f0;T=2e-3;t=1:1/fs:T-1/fs;%元阵阵列参数c=1500;lamda=c/f0;d=lamda/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;fori_theta=1:length(theta)e=[cos(theta(i_theta));sin(theta(i_theta));0];tao=-e.'*B./c;fori_M=1:MA(i_M,i_theta)=exp(sqrt(-1)*2*pi*f0*tao(i_M));endend%加权向量wfai_direct=2*pi*f0*d*sin(direct*pi*180)/c;fori_M=1:Mw(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');grid2.MUSIC方位估计%%clc;clear;closeall;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;fori_M=1:MS(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)fori=1:numX=awgn(S,SNR);R=X*X'/N;%协方差矩阵[V,I]=eig(R);[OrderIndex]=sort(diag(I),'descend');%特征值降值排列VS=[];%向量子空间P_music=[];fori_M=1:MVS=[VS,V(:,Index(i_M))];%排序后的特征向量endE=VS(:,signalnum+1:M);%后M-signalnum向量构成向量子空间En=E*E';%噪声子空间theta=-90:1:90;fai=2*pi*f0*d*sin(theta*pi/180)/c;fori_theta=1:length(theta)fori_M=1:Mw(i_theta,i_M)=exp(sqrt(-1)*(i_M-1)*fai(i_theta));%方向矢量endP_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');gridtitle('不同信噪比下MUSIC方法估计的性能（统计100次）')holdon%%%信噪比SNR=0;theta1=12;num=100;theta1=theta1*pi/180;tau0=d*sin(theta1)/c;fori_M=1:MS(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)fori=1:numX=awgn(S,SNR);R=X*X'/N;%协方差矩阵[V,I]=eig(R);[OrderIndex]=sort(diag(I),'descend');%特征值降值排列VS=[];%向量子空间P_music=[];fori_M=1:MVS=[VS,V(:,Index(i_M))];%排序后的特征向量endE=VS(:,signalnum+1:M);%后M-signalnum向量构成向量子空间En=E*E';%噪声子空间theta=-90:1:90;fai=2*pi*f0*d*sin(theta*pi/180)/c;fori_theta=1:length(theta)fori_M=1:Mw(i_theta,i_M)=exp(sqrt(-1)*(i_M-1)*fai(i_theta));%方向矢量endP_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');gridtitle('不同信噪比下MUSIC方法估计的性能（统计100次）')holdon%%%信噪比SNR=5;theta1=12;num=100;theta1=theta1*pi/180;tau0=d*sin(theta1)/c;fori_M=1:MS(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)fori=1:numX=awgn(S,SNR);R=X*X'/N;%协方差矩阵[V,I]=eig(R);[OrderIndex]=sort(diag(I),'descend');%特征值降值排列VS=[];%向量子空间P_music=[];fori_M=1:MVS=[VS,V(:,Index(i_M))];%排序后的特征向量endE=VS(:,signalnum+1:M);%后M-signalnum向量构成向量子空间En=E*E';%噪声子空间theta=-90:1:90;fai=2*pi*f0*d*sin(theta*pi/180)/c;fori_theta=1:length(theta)fori_M=1:Mw(i_theta,i_M)=exp(sqrt(-1)*(i_M-1)*fai(i_theta));%方向矢量endP_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');gridtitle('不同信噪比下MUSIC方法估