clear all;close all;clc;
snapshot_number=500; %快拍数
Ts=0.001; %采样周期
SNR=40; %信噪比
snr=15;%信噪比
tangle1=10; %三个信号的发射角度
tangle2=35; %三个信号的入射角度
tangle3=60;
c1=kron(exp(-1i*pi*(0:9)*sin(tangle1*pi/180))',exp(-1i*pi*(0:9)*sin(tangle1*pi/180))'); %3个信号导向矢量
c2=kron(exp(-1i*pi*(0:9)*sin(tangle2*pi/180))',exp(-1i*pi*(0:9)*sin(tangle2*pi/180))');
c3=kron(exp(-1i*pi*(0:9)*sin(tangle2*pi/180))',exp(-1i*pi*(0:9)*sin(tangle2*pi/180))');
A=[c1,c2 c3];
for r=1:50
s=5*(2*round(rand(3,1*snapshot_number))-1); %仿真信号
z1=A*s;
sigPower=mean(sum(abs(z1).^2)/length(c1));
noisePower=sigPower*10.^(-SNR/10);
x=z1+randn(length(c1),snapshot_number)*sqrt(noisePower); %总输入信号
R=x*x'/snapshot_number;
iR=inv(R);
[U,V]=eig(R);
y=diag(V);
[m,n]=sort(y);
Us=U(:,n(98:100));
Us1=Us(1:12,:);
Us2=Us(4:15,:);
Rr=(Us1'*Us1)\(Us1'*Us2);
[Ur,Vr]=eig(Rr); % Vr为接收角度矩阵
ranglea=asin(1/pi*angle(Vr(1,1)));
rangleb=asin(1/pi*angle(Vr(2,2)));
ranglec=asin(1/pi*angle(Vr(3,3)));
rangle=[ranglea,rangleb,ranglec];
DOA1(r)=min(rangle)*180/pi; % 测得的3个接收角
DOA3(r)=max(rangle)*180/pi;
DOA2(r)=(ranglea+rangleb+ranglec)*180/pi-DOA1(r)-DOA3(r);
end
r=1:1:50;
figure(1)
plot(r,DOA1,'-or',r,DOA2,'-*k',r,DOA3,'-+b');
axis([0 50 -80 80]);
xlabel('蒙特卡洛次数');
ylabel('DOA/(°)');
legend('DOA1','DOA2','DOA3',4);
title('ESPRTI Algorithm');
grid on;
