doa.rar_DOA估计_SNR Radar_doa估计算法_music_radar

%classic music clear all clc tic %参数设定 M=8; f=1000;c=1500; lambda=c/f;d=lambda/2; SNR=[5:2:25]; F=length(SNR); N=128;%快拍数。 K=100; doa=[-10 30]; P=length(doa); for g=1:F %信噪比的循环。 for k=1:K %阵列流型A for i=1:P A(:,i)=exp(-j*2*pi*d*[0:M-1]'/lambda*sin(doa(i)*pi/180));%均匀线性阵的阵列流型矢量。 end %信源模型建立 for i=1:P S(i,:)=sqrt(10^(SNR(g)/10))*(randn(1,N)+j*randn(1,N));%sqrt(10.^(snr/10))使信号符合该信噪比。 end %通道幅度相位误差 for m=2:M u1= unifrnd(-pi,pi,1,M-1);%1*7. ph(1,1)=1; ph(:,m)=exp(j*0.2*u1(1,m-1));%1*8相位误差。 u2=unifrnd(-1,1,1,M-1);%1*7。 al(1,1)=1; al(:,m)=1+0.2*u2(1,m-1);%1*8幅度误差。 gamma(:,1)=1; gamma(:,m)=ph(:,m)*al(:,m);%1*8幅相误差。 end %接收信号模型建立 X=A*S+1/sqrt(2)*(randn(M,N)+j*randn(M,N));%1/sqrt(2)使高斯白噪声的能量为1。 %协方差矩阵特征值分解得到噪声子空间 R=X*X'/N; [V,D]=eig(R);% [V,D]=eig(A)：求矩阵R的全部特征值，构成对角阵D，并求A的特征向量构成V的列向量。 [Y,I]=sort(diag(D));%diag是将对角阵的对角元素提取成一个向量.如果A是向量，sort(A) 对A中元素按照升序排列。 % 如果A是矩阵sort(A) 对A按每一列元素按照升序排列。Y为排列后矩阵，I是原矩阵中各元素位置所组成的新矩阵。 Un=V(:,I(1:M-P));%把原矩阵第1到M-P列中的数对应到V矩阵中的列数，建立新的矩阵赋值给Un。 Xe=diag(gamma)*A*S+1/sqrt(2)*(randn(M,N)+j*randn(M,N));%1/sqrt(2)使高斯白噪声的能量为1。 %协方差矩阵特征值分解得到噪声子空间 Re=Xe*Xe'/N; [Ve,De]=eig(Re);% [V,D]=eig(A)：求矩阵R的全部特征值，构成对角阵D，并求A的特征向量构成V的列向量。 [Ye,Ie]=sort(diag(De));%diag是将对角阵的对角元素提取成一个向量.如果A是向量，sort(A) 对A中元素按照升序排列。 % 如果A是矩阵sort(A) 对A按每一列元素按照升序排列。Y为排列后矩阵，I是原矩阵中各元素位置所组成的新矩阵。 Une=Ve(:,Ie(1:M-P));%把原矩阵第1到M-P列中的数对应到V矩阵中的列数，建立新的矩阵赋值给Un。 %谱峰搜索部分 theta=-90:0.1:90;%线阵的搜索范围为-90~90度 for i=1:length(theta) a_theta=exp(-j*(0:M-1)'*2*pi*d*sin(pi/180*theta(i))/lambda); Pmusic(i)=1./abs((a_theta)'*Un*Un'*a_theta); Pmusic2(i)=1./abs((a_theta)'*Une*Une'*a_theta); P1(i)=10*log10(Pmusic(i)/min(Pmusic)); P2(i)=10*log10(Pmusic2(i)/min(Pmusic2)); end [peaks,locals]=findpeaks(Pmusic);%峰值及其位置 b=sortrows( [peaks' locals']);%峰值按从小到大排列 bb=b'; B=bb(:,size(b',2)-1:size(b',2)); BB=B(2,:);%峰值角度对应的数字 BBB=sort(BB);%最大峰值角度对应的数字 M_theta=theta(BBB);%所测得入射角 mistake1(g,2*k-1:2*k)=M_theta-doa;%每次实际入射角与测得入射角得差值，11*200,11对应11个信噪比，200对应200个快拍。 [peaks2,locals2]=findpeaks(Pmusic2);%峰值及其位置 b2=sortrows( [peaks2' locals2']);%峰值按从小到大排列 bb2=b2'; B2=bb2(:,size(b2',2)-1:size(b2',2)); BB2=B2(2,:);%峰值角度对应的数字 BBB2=sort(BB2);%最大峰值角度对应的数字 M_theta2=theta(BBB2);%所测得入射角 mistake2(g,2*k-1:2*k)=M_theta2-doa;%每次实际入射角与测得入射角得差值 end RMSE(g,:)=sqrt(sum(mistake1(g,:).^2,2)/2/K); RMSE2(g,:)=sqrt(sum(mistake2(g,:).^2,2)/2/K); end figure(1) plot(SNR,RMSE,'bd-',SNR,RMSE2,'ro-','MarkerFaceColor','c'); xlabel('SNR/dB'); ylabel('RMSE/°'); legend('无误差均方值RMSE1','有误差均方值RMSE2'); title('存在幅度相位误差与无幅相误差的角度误差均方值比较'); grid on; figure(2) plot(theta,P1,'b--',theta,P2,'r-'); axis([-60 60 0 60]); xlabel('theta/°'); ylabel('功率谱/dB'); %h=legend('g=0','g=0.2',0); %h=legend('g=0','g=0.2',0); h=legend('理想情况下','存在幅相误差',0); set(h,'box','off'); title('存在相位误差与理想情况下的功率谱比较'); %title('仅存在幅度误差与理想情况下的功率谱比较'); %title('存在幅相误差与理想情况下的功率谱比较'); grid on; toc