%LS_ESPRIT ALGORITHM
%DOA ESTIMATION BY LS_ESPRIT ALOGRITHM
clear all;
%close all;
clc;
source_number=2;%信元数
sensor_number=8;%原阵元数
sub_sensor_number=7;%子阵元数
N_x=1024; %信号长度
snapshot_number=N_x;%快拍数
w=[pi/4 pi/6].';%信号频率
l=sum(2*pi*3e8./w)/length(source_number);%信号波长
d=0.5*l;%阵元间距
snr=0;%信噪比/dB
source_doa=[-10 60];%两个信号的入射角度
A=[exp(-1j*(0:sensor_number-1)*d*2*pi*sin(source_doa(1)*pi/180)/l);exp(-1j*(0:sensor_number-1)*d*2*pi*sin(source_doa(2)*pi/180)/l)].';%阵列流型
for n=1:50
s=sqrt(10.^(snr/10))*exp(1j*w*[0:N_x-1]);%仿真信号
x=A*s+(1/sqrt(2))*(randn(sensor_number,N_x)+1j*randn(sensor_number,N_x));%加了高斯白噪声后的阵列接收信号
x1=x(1:sub_sensor_number,:);%子阵1接受的数据矢量
x2=x(2:(sub_sensor_number+1),:);%子阵2接受的数据矢量
%对两个子阵的模型进行合并
X=[x1;x2];
R=X*X'/snapshot_number;
%对R进行奇异值分解
[U,S,V]=svd(R);
R=R-S(2*sub_sensor_number,2*sub_sensor_number)*eye(2*sub_sensor_number);
[U,S,V]=svd(R);
Us=U(:,1:source_number);
disp(Us);
Us1=Us(1:sub_sensor_number,:);
Us2=Us((sub_sensor_number+1):2*sub_sensor_number,:);
%按照公式得到旋转不变矩阵M
E=pinv(Us1)*Us2;
disp('E');
disp(E);
%对得到的旋转不变矩阵进行特征分解
[V,D]=eig(E);
disp(D);
D=(diag(D)).';
doa=-asin(angle(D)/pi)*180/pi;
doa=sort(doa);
plot(n,doa(1),'*',n,doa(2),'square');
axis([0 51 -100 100]);
grid on;
xlabel('Number of Experiments');
ylabel('DOA Estimation/degree');
title('DOA Estimation by LS ESPRIT Algorithm with Monte Carlo');
hold on;
drawnow;
end
