%esprit算法子程序
function DOA_guji=esprit(y,snr,M,fc,c,q,d,N)
%y阵列接收数据，snr信噪比，M阵元个数，fc中心频率，c光速
%q信号源个数，d阵元间距，N采样次数
y=awgn(y,snr,0);%按信噪比snr给个阵元加上高斯白噪声
R=y*y'/N;%协方差矩阵
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[v,dd]=eig(R);%特征值分解
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%求取信号子空间
if dd(1,1)>dd(2,2)
u=v(:,1:q);
else
u=v(:,M-q+1:M);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
u1=u(1:M-1,:);
u2=u(2:M,:);
phi=pinv(u1)*u2;%估计旋转矩阵
dg=eig(phi);
DOA_guji=zeros(q,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%利用旋转矩阵估计DOA
for k=1:q
DOA_guji(k,1)=-asin((angle(dg(k))/(2*pi*d*fc)*c))*180/pi;
end
% disp('esprit 测向结果');
% DOA_guji
