music1snr.m

不同信噪比对MUSIC算法实现波达方向估计的影响

clear all; 
close all;
clc;
theta=-90:0.1:90;%ULA估计角度变化的范围和频率选择
Ltheta=length(theta);
SNR=[-10:30];%信噪比单位是dB.ps已经确定为2.由xxx的计算可知。
ps=2;
Lsnr=length(SNR);
temp=10.^(SNR/10);
sigmaN=sqrt(ps./temp);%由信噪比计算出来的噪声标准差。
DOA=zeros(Ltheta,Lsnr);%用来存储doa结果的矩阵。
N_x=1000; %快拍数 
N=8; % Size of Rx Matrix 
l=1.8;%波长
c=3e8;
w0=2*pi*c/l;
d=0.5*l;%阵元间距为半波长。
M=3; % Number of Signals 
w=[pi/6 pi/10 pi/8]';%信号频率 与2不同之处
xx=pi/30;yy=0;zz=-pi/4; %三个信号的入射角,小于5度就难以辨别了
B=[1 exp(-j*2*pi*d*sin(xx)/l) exp(-j*2*2*pi*d*sin(xx)/l) exp(-j*2*3*pi*d*sin(xx)/l)...
    exp(-j*2*4*pi*d*sin(xx)/l) exp(-j*2*5*pi*d*sin(xx)/l) exp(-j*2*6*pi*d*sin(xx)/l)...
    exp(-j*2*7*pi*d*sin(xx)/l);
   1 exp(-j*2*pi*d*sin(yy)/l) exp(-j*2*2*pi*d*sin(yy)/l) exp(-j*2*3*pi*d*sin(yy)/l)...
   exp(-j*2*4*pi*d*sin(yy)/l) exp(-j*2*5*pi*d*sin(yy)/l) exp(-j*2*6*pi*d*sin(yy)/l)...
   exp(-j*2*7*pi*d*sin(yy)/l);
   1 exp(-j*2*pi*d*sin(zz)/l) exp(-j*2*2*pi*d*sin(zz)/l) exp(-j*2*3*pi*d*sin(zz)/l)...
   exp(-j*2*4*pi*d*sin(zz)/l) exp(-j*2*5*pi*d*sin(zz)/l) exp(-j*2*6*pi*d*sin(zz)/l)...
   exp(-j*2*7*pi*d*sin(zz)/l)]'; %阵列流型,信号源决定行数，阵元数决定列数
   xxx=2*exp(j*(w+w0)*[0:N_x-1]);%仿真信号  
for ii=1:Lsnr
x=B*xxx+sigmaN(ii)*randn(8,N_x)+j*randn(8,N_x)*sigmaN(ii); %加噪声
R=x*x'/N_x;  %%%'就是共轭转置。应该用快拍数来实现循环.
[V D]=eig(R);
[lambda,index] = sort(diag(D));  %%从小到大的排序。
UU=V(:,index(1:5));    
for i = 1:length(theta)
AA=[1 exp(-j*2*pi*d*sin(theta(i)/180*pi)/l) exp(-j*2*2*pi*d*sin(theta(i)/180*pi)/l)...
    exp(-j*2*3*pi*d*sin(theta(i)/180*pi)/l) exp(-j*2*4*pi*d*sin(theta(i)/180*pi)/l)...
    exp(-j*2*5*pi*d*sin(theta(i)/180*pi)/l) exp(-j*2*6*pi*d*sin(theta(i)/180*pi)/l) ...
    exp(-j*2*7*pi*d*sin(theta(i)/180*pi)/l)];%方向矢量
WW=AA*UU*UU'*AA';
Pmusic(i)=abs(1/WW);%角谱
end
Pmusic = 10*log10(Pmusic);
DOA(:,ii)=Pmusic';
end
aa=theta'*ones(1,Lsnr);
bb=ones(Ltheta,1)*SNR;
mesh(aa,bb,DOA);
title('PMUSIC(\theta,SNR)');
xlabel('方位角\theta/degree') 
ylabel('SNR/dB')
zlabel('PMUSIC')