esprit.m

谱估计源程序：是谱估计中的经典算法esprit,MUSIC等算法的源程序

clear all; close all;clc;
J=sqrt(-1);
source_number=3;
source_doa=[40 60 80];
sensor_number=10;
snapshort_number=100;
snr=10;
lamda=0.6;%波长
array_distance=0.1;%子阵的距离

A=exp(-J*(0:sensor_number-1)'*pi*sin(source_doa*pi/180));
%theta=[20,40,60];
%fei=[10,30,50];
%for i=1:3
%A(:,i)=exp(-j*pi*cos(theta(i)*pi/180-2*pi*(0:9)/10)*sin(fei(i)*pi/180));
%end

s=(randn(source_number,snapshort_number)+J*randn(source_number,snapshort_number));
x=A*s;
X=awgn(x,snr);
M=[exp(j*pi/3*sin(source_doa(1)*pi/180)),exp(j*pi/3*sin(source_doa(2)*pi/180)),exp(j*pi/3*sin(source_doa(3)*pi/180))];
M=diag(M);
y=A*M*s ;
Y=awgn(y,snr);

% % %TLS-ESPRIT
 Z=[X',Y']';
 Rzz=Z*Z'/snapshort_number;
 [Uz,Sz,Vz]=svd(Rzz);
% % [Uz,Sz]=eig(Rzz);
% 
% 
 E=Uz(:,1:source_number);
 Ex=E(1:sensor_number,:);
 Ey=E(sensor_number+1:2*sensor_number,:);
 fy=pinv(Ex)*Ey;
%PRO-ESPRIT
%[Ux,Sx,Vx]=svd(X);
%Uxs=Ux(:,1:q);
%Sxs=Sx(1:q,1:q);
%Vxs=Vx(:,1:q);
%fy=inv(Sxs)*Uxs'*Y*Vxs;
% compute thetae
[Ufy,Sfy]=eig(fy);
for i=1:3
    thetae(1,i)=asin(3*angle(Sfy(i,i))/pi)*180/pi;
end
disp(thetae)


 
 %LS-ESPRIT
%Rx=X*X'/snapshort_number;
%Ry=Y*Y'/snapshort_number;
[%Ux,Sx,Vx]=svd(Rx);
%[Ux,Sx]=eig(Rx);
%Uxs=Ux(:,1:source_number);
%Uxn=Ux(:,source_number+1:sensor_number);

%[Uy,Sy,Vy]=svd(Ry);
 %% [Uy,Sy]=eig(Rx); 
 %Uys=Uy(:,1:source_number);
 %Uyn=Uy(:,source_number+1:sensor_number);
% %fy=pinv(Uys)*Uxs;
 %fy=pinv(Uxs)*Uys;
%[Ufy,Sfy]=eig(fy)
%for i=1:3
   % doa(1,i)=asin(3*angle(Sfy(i,i))/pi)*180/pi;
%end
%disp(doa)

% % %TLS-ESPRIT
 Z=[X',Y']';
 Rzz=Z*Z'/nn;
 [Uz,Sz,Vz]=svd(Rzz);
% % [Uz,Sz]=eig(Rzz);
% 
% 
 E=Uz(:,1:source_number)
 Ex=E(1:source_number,:);
 Ey=E(source_number+1:2*source_number,:);
 fy=pinv(Ex)*Ey;

%[Ux,Sx,Vx]=svd(X);
%Uxs=Ux(:,1:source_number);
%Sxs=Sx(1:source_number,1:source_number);
%Vxs=Vx(:,1:source_number);
%fy=inv(Sxs)*Uxs'*Y*Vxs;
% compute thetae
[Ufy,Sfy]=eig(fy);
for i=1:2
    thetae(1,i)=asin(3*angle(Sfy(i,i))/pi)*180/pi;
end
disp(thetae)

%R=X*X'/snapshort_number;
%[U,S,V]=svd(R);
%Un=U(:,source_number+1:sensor_number);
%Gn=Un*Un';
%searching_doa=-90:0.1:90;
%for i=1:length(searching_doa)
   % a_theta=exp(-J*(0:sensor_number-1)'*pi*sin(pi*searching_doa(i)/180));
  %  P_music(i)=1./abs((a_theta'*Gn*a_theta));
%end
%plot(searching_doa,P_music);
%axis([-90 90 -90 90]);
%xlabel('入射角');
%ylabel('谱峰');
%legend('Music spectrum');