lkk_esprit.m

旋转不变子空间esprit算法对波达方向估计源程序

clc
clear
Pd=1000;

s1=10*exp(j*pi/3*[0:Pd])';
s2=10*exp(j*pi/4*[0:Pd])';
s3=10*exp(j*pi/5*[0:Pd])';

nn=512; 
M=9;%原阵元数
m=8;%子阵元数 
p=3;%信源数 

angle=[30;45;60];
degrad=pi/180;
tt=1:nn;
S=[s1(tt),s2(tt),s3(tt)].';

%引入信噪比
SN1=20;SN2=20;SN3=20;  
SN=[SN1;SN2;SN3];
Ps=S*S'/nn;    
Ps1=diag(Ps);
refp=2*10.^(SN/10);
tmp=sqrt(refp./Ps1);
S1=diag(tmp)*S;

%构造噪声
nr=randn(M,nn);  
ni=randn(M,nn);
E=nr+ni*j;
Ex=E(1:m,:);
Ey=E(2:M,:);

A=zeros(m,p);
k=[0:m-1]';
for t=1:p
    A(:,t)=exp(j*pi*k*sin(angle(t)*degrad));
end

X=A*S1+Ex;
F=diag(exp(j*pi*sin(angle*degrad)));
Y=A*F*S1+Ey;

%Rxx=cov(X');
Rxx=X*X'./nn;%X的自协方差矩阵
Rxy=X*Y'./nn;%XY的互协方差矩阵

[U,s,V_1]=svd(Rxx);
S_1=diag(s);
delta=S_1(m);%为1~2之间的一个数

Cxx=Rxx-delta*eye(m);%%%%错误出在Cxx和Cxy上面，由于Rxx和Rxy远大于1，以至于减号后的东西基本没用...
z=[0 0 0 0 0 0 0 0;
   1 0 0 0 0 0 0 0;
   0 1 0 0 0 0 0 0;
   0 0 1 0 0 0 0 0;
   0 0 0 1 0 0 0 0;
   0 0 0 0 1 0 0 0;
   0 0 0 0 0 1 0 0;
   0 0 0 0 0 0 1 0];
Cxy=Rxy-delta*z;

[V_2,D]=eig(Cxx,Cxy);%做两个矩阵之间的广义特征分解
B=sort(diag(D));
aangle=B(m-p+1:m);%取最大的3个作为信号源对应的特征值
abs1=abs(aangle);
eangle=acos(real(aangle./abs1))/pi*180