jiaodu.rar_DOA_DOA估计_matlab

% -5 10 40 clc clear all close all; warning off L=16; %阵元数 d=0.5; M=3; sigma=1;%噪声功率 SNR=-10:0.1:30; for ii=1:length(SNR) % SNR=0; for jj=1:500 SNR1=10; %干噪比 SNR2=10; sigma0=sigma*10^(SNR(ii)/10); %信号 sigma1=sigma*10^(SNR1/10); sigma2=sigma*10^(SNR2/10); %% 保护信号和干扰信号的波达方向角 theta0=0; %保护信号 theta0_0=2; %期望信号波达方向标称值为2度，存在指向误差 theta1=45; %干扰信号1 theta2=-45; %干扰信号2 %% 保护信号和干扰信号相对于观测阵列的导向矢量 a0=exp(1i*2*pi*d*sin(theta0*pi/180)*(0:L-1)'); a0_0=exp(1i*2*pi*d*sin(theta0_0*pi/180)*(0:L-1)'); %标称 a1=exp(1i*2*pi*d*sin(theta1*pi/180)*(0:L-1)'); a2=exp(1i*2*pi*d*sin(theta2*pi/180)*(0:L-1)'); %% 生成仿真信号 sig_len=20; s0=(randn(1,sig_len)+1i*randn(1,sig_len))*sqrt(sigma0/2); s1=(randn(1,sig_len)+1i*randn(1,sig_len))*sqrt(sigma1/2); s2=(randn(1,sig_len)+1i*randn(1,sig_len))*sqrt(sigma2/2); n=(randn(1,sig_len)+1i*randn(1,sig_len))*sqrt(sigma/2); x=a0_0*s0+a1*s1+a2*s2+n; % Rxx=sigma0^2*(a0_0*a0_0')+sigma1^2*(a1*a1')+sigma2^2*(a2*a2')+sigma^2*eye(L); %协方差矩阵 R_in=sigma1^2*(a1*a1')+sigma2^2*(a2*a2')+sigma^2*eye(L); Rxx_est=x*x'/length(x); %样本协方差矩阵估计 %SMI法 w_smi_est=Rxx*inv(Rxx_est)*a0_0; %SSIP法 [V,D]=eig(Rxx_est); d1=diag(D); [B,I]=sort(d1); %升序排列,I为索引 dd=B(L-M+1:end); %大特征值 vv=V(:,I(L-M+1:end)); %大特征值对应的特征矢量 vvv=vv(:,1)*vv(:,1)'/dd(1)+vv(:,2)*vv(:,2)'/dd(2)+vv(:,3)*vv(:,3)'/dd(3); w_ssip_est=Rxx*inv(Rxx_est)*vvv*a0_0; %DL法 lamda=0.01; w_dl_est=Rxx*inv(Rxx_est+lamda*eye(L))*a0_0; %计算三种方法的失配广义角 P_smi(jj)=acos(sqrt((abs(w_smi_est'*inv(R_in)*a0))^2/(w_smi_est'*inv(R_in)*w_smi_est)/(a0'*inv(R_in)*a0)))*180/pi; P_ssip(jj)=acos(sqrt((abs(w_ssip_est'*inv(R_in)*a0))^2/(w_ssip_est'*inv(R_in)*w_ssip_est)/(a0'*inv(R_in)*a0)))*180/pi; P_dl(jj)=acos(sqrt((abs(w_dl_est'*inv(R_in)*a0))^2/(w_dl_est'*inv(R_in)*w_dl_est)/(a0'*inv(R_in)*a0)))*180/pi; end smi(ii)=mean(P_smi); ssip(ii)=mean(P_ssip); dl(ii)=mean(P_dl); end figure; plot(SNR,smi,'r-','linewidth',1.5); hold on plot(SNR,ssip,'b:','linewidth',1.5); hold on plot(SNR,dl,'g--','linewidth',1.5); legend('SMI','SSIP','DL'); xlabel('输入信噪比/dB');ylabel('失配广义夹角/°'); grid on