esprit.m

esprit算法DOA估计esprit算法DOA估计esprit算法DOA估计

clear all;
j=sqrt(-1);
N=8;%天线阵元数
P=3;%信号源个数
angle1=-40;angle2=20;angle3=60;%信号达到方向
th=[angle1,angle2,angle3];
th=sort(th);
L=1024;%采样点数,即快拍数
rat=0.5;%rat 代表阵元距离是信号波长的一半
degrad=pi/180;
AM=[(10),(10),(10)];%每个信号源的幅度
th=sin(th*degrad);
% sn=[10 9 8];
randn('state',sum(100*clock));%Resets it to a different state each time.
pn=randn(1,P*L);
count=1;
for k=1:P
   for m=1:L
      s(k,m)=AM(1,k)*pn(1,count);%s transmited signal matrix
      count=count+1;
   end
end
randn('state',sum(100*clock));
pn1=randn(1,(N-1)*L);
count=1;
for k=1:N-1
   for m=1:L
      n1(k,m)=pn1(1,count);%s transmited signal matrix
      count=count+1;
   end
end
randn('state',sum(100*clock));
pn2=randn(1,(N-1)*L);
count=1;
for k=1:N-1
   for m=1:L
      n2(k,m)=pn2(1,count);%s transmited signal matrix
      count=count+1;
   end
end
sn=[-4:2:45];
suma=zeros(1,25);
for l=1:25
Ps=s*s'/L;%求S的均值
ps=diag(Ps);%得Ps的主对角线元素
refp=2*10.^(sn(l)/10);
tmp=sqrt(refp./ps');
s=diag(tmp)*s;
for k=1:N-1
   for m=1:P
      A1(k,m)=exp(-j*2*pi*rat*(k-1)*th(m));
   end
end
for k=1:P
      E(k)=exp(j*2*pi*rat*th(k));
  end
E1=diag(E);
A2=A1*E1;
X1=A1*s+n1;
X2=A2*s+n2;
X=[X1;X2];
 Rxx=X*X'/L;
[V,D]=eig(Rxx);
D=diag(D);
minD=min(D);
I=eye(size(Rxx));
Rxx=Rxx-minD*I;
[V1,D1]=eig(Rxx);
    Vs=V1(:,2*(N-1)-P+1:2*(N-1));%信号子空间 
    Vs1=Vs(1:(N-1),:);
    Vs2=Vs(N:2*(N-1),:);
J=inv(Vs1'*Vs1)*Vs1'*Vs2;
d1=eig(J);
d1=angle(d1);
doa=asin(d1/pi)/pi*180;
doa=sort(doa);
% stem(doa,'r');
% title('LS-ESPRIT算法');
% ylabel('角度');
% grid on;
t=doa'-th;
t=abs(t);
suma(l)=sum(t);
end
suma=log(suma);
plot(sn,suma);
ylabel('dB');
grid on;