ESPRIT算法MATLAB仿真程序

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,c2c3];

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([050-8080]);

xlabel('蒙特卡洛次数');

ylabel('DOA/(°)');

legend('DOA1','DOA2','DOA3',4);

title('ESPRTI Algorithm');

grid on;