esprit.m

阵列信号处理中经典的的ESPRIT算法

%%%%%%%%%%%%% TLS-ESPRIT algorithm %%%%%%%%%%%%%

clear;
clc;

%---first value---

M=6;
K=2;
d=0.5;
N=1024;
theta1=30*pi/180;  theta111=30;
theta2=50*pi/180;  theta222=50;


%----s---
for SSS=1:7, 
   snr=0+(SSS-1)*5;
   Amp=sqrt(2*10^(snr/10));

for MM=1:50,
    SSS
    MM

   s1=(exp(j*2*pi*rand(1,N)))*Amp;

   s2=(exp(j*2*pi*rand(1,N)))*Amp;

s=[s1;s2];


%----a---

a1=exp(-j*2*pi*d*[0:M-1].'*sin(theta1));

a2=exp(-j*2*pi*d*[0:M-1].'*sin(theta2));

a=[a1,a2];
%----noise---

noise=randn(M,N)+j*randn(M,N);

%-----z-------

z=a*s+noise;

%----caculate R----
R=z*z'/N;

% ----do the eigendecomposition---
[U,D,V]=svd(R);

%----get the signal subspace-----
Ep1=U(1:(M-1),1:K);
Ep2=U(2:M,1:K);

temp1=[Ep1,Ep2];
temp2=[Ep1';Ep2'];
temp=temp2*temp1;

[u_p,lpp]=eig(temp);
E12=u_p(1:2,3:4);
E22=u_p(3:4,3:4);
psaip=-E12/(E22);
[u,faip]=eig(psaip);


%----estimate  DOA ----

theta11(MM)=-asin(angle(faip(1,1))/pi/2/d)*180/pi;
theta22(MM)=-asin(angle(faip(2,2))/pi/2/d)*180/pi;

if abs(theta11(MM)-theta111)>abs(theta22(MM)-theta111)
  temp=theta11(MM);
  theta11(MM)=theta22(MM);
  theta22(MM)=temp;
  
end

end

%  mtheta1(SSS)=mean(theta11);
%  mtheta2(SSS)=mean(theta22);
%  
%  vart1(SSS)=std(theta11);
%  vart2(SSS)=std(theta22);
rmset1(SSS)=sqrt(sum((theta11-theta111).^2)/MM)
rmset2(SSS)=sqrt(sum((theta22-theta222).^2)/MM)

end

%---plot command---
i=1:SSS;
snr=0+(i-1)*5;

% figure(1)
% subplot(1,2,1),plot(snr,abs(mtheta1(i)-theta111),'k-s');
% legend({'VESPA'});
% xlabel('信噪比/dB'); ylabel('\theta_{1}的偏差/\circ');
% subplot(1,2,2),plot(snr,vart1(i),'k-s');
% legend({'VESPA'});
% xlabel('信噪比/dB'); ylabel('\theta_{1}的标准差/\circ');
% 
% figure(2)
% subplot(1,2,1),plot(snr,abs(mtheta2(i)-theta222),'k-s');
% legend({'VESPA'});
% xlabel('信噪比/dB'); ylabel('\theta_{2}的偏差/\circ');
% subplot(1,2,2),plot(snr,vart2(i),'k-s');
% legend({'VESPA'});
% xlabel('信噪比/dB'); ylabel('\theta_{2}的标准差/\circ');

figure(1)
plot(snr,rmset1,'k-o');%,snr,rmsea2,'k-s');
xlabel('信噪比/dB'); ylabel('\theta_{1}估计的RMSE/\circ');
legend({'VESPA'});

figure(2)
plot(snr,rmset2,'k-o');%,snr,rmsea2,'k-s');
xlabel('信噪比/dB'); ylabel('\theta_{2}估计的RMSE/\circ');
legend({'VESPA'});