%基本ESPRIT算法
clear all;close all;clc;
c=3*10^8;
f=3*10^9;
%%求得信号的波长
lamda=c/f;
%%阵元间距
d=lamda/2;
%%（n-1）为子阵列的个数即阵元数
n=10;
%%信号的数目
signal_number=3;
%%三个信号的角度值
thita1=-25;
thita2=30;
thita3=65;
%%三个信号的中心频率
f1=40;
f2=20;
f3=70;
%%在时域来说是快拍数（一段时间内对阵列数据采样的个数）；在频域来说是DFT的时间子段的个数
snapshot=1:2000;
%%S是信号空间，有三个信号组成
S1=44.67*exp(j*2*pi*f1*snapshot/length(snapshot));
S2=44.67*exp(j*2*pi*f2*snapshot/length(snapshot));
S3=44.67*exp(j*2*pi*f3*snapshot/length(snapshot));
S=[S1;S2;S3];
%%子阵1
A1=exp(-j*2*pi*d*[0:n-1]*sin(thita1*pi/180)/lamda).';%旋转不变量矩阵
A2=exp(-j*2*pi*d*[0:n-1]*sin(thita2*pi/180)/lamda).';
A3=exp(-j*2*pi*d*[0:n-1]*sin(thita3*pi/180)/lamda).';
A=[A1,A2,A3];
%%噪声假设为高斯白噪声，均值为0
N= wgn(10,2000,3);
%%求信噪比的s1、s2、s3 信噪比依次是10、20、30
s_power1=10*log(2.72^2/2);
s_power2=10*log(4.48^2/2);
s_power3=10*log(7.37^2/2);
snr1=s_power1-3;
snr2=s_power2-3;
snr3=s_power3-3;
%%整个阵列接收到的数据0-（n-1）为阵列1，1-n为阵列2；
X=A*S+N;
%%协方差矩阵?
Rxx=X*X'/length(snapshot);
%%对整个数据的协方差矩阵进行特征值分解，从而得到特征值向量D和特征向量V
[V,D]=eig(Rxx);
%[Y,I]=sort(diag(D));
Us=V(:,n-signal_number+1:n);
%%两个方阵张成的两个子空间
U1=Us(1:n-1,:);
U2=Us(2:n,:);
%%利用最小二乘法求得旋转不变关系矩阵，然后进行特征特征分解
[p,q]=eig(inv(U1'*U1)*U1'*U2); %张贤达《矩阵分析与应用》第528页
%%利用上面求得的矩阵来获得角度
for i=1:signal_number;
alpha(i)=real(asin(-j*(log(q(i,i)))*lamda/(-2*pi*d))*180/pi);%q为特征值，对应角度
end;
%%作图
stem(alpha,ones(1,signal_number),'r--');grid;
axis([-90 90 0 2]);
text(alpha(1)-4,1.1,num2str(alpha(1)));text(alpha(1)-15,1.4,'(S1|SNR10)');
text(alpha(2)-4,1.1,num2str(alpha(2)));text(alpha(2)-15,1.4,'(S2|SNR20)');
text(alpha(3)-4,1.1,num2str(alpha(3)));text(alpha(3)-15,1.4,'(S3|SNR30)');
ylabel('DOA估计的角度值');
xlabel('角度');
title('ESPRIT算法DOA估计');
