%MUSIC; WSF; etc;
clear all; close all;clc;
J=sqrt(-1);
source_number=2;
source_doa=[40 70];
sensor_number=13;
snapshot_number=100 ;
snr=10;
 
A=exp(-J*(0:sensor_number-1)'*pi*sin(source_doa*pi/180));

        
b=[1e6 ,5e6];
Fs=10e6;%采樣周期
Ts=1/Fs;
 
 
t=linspace(Ts,snapshot_number*Ts,snapshot_number);
for n=1:snapshot_number 
   S1(1,n)=exp(j*2*pi*1e6*t(1,n));
   S2(1,n)=exp(j*2*pi*5e6*t(1,n));
end
 
s=[S1;S2];%%s=[S1.',S2.']'
 
 

%s=(randn(source_number,snapshot_number)+J*randn(source_number,snapshot_number))/sqrt(2);

x=A*s;
y=awgn(x,snr);

R=y*y'/snapshot_number;
 
[V,D]=eig(R);
Un=V(:,1:sensor_number-source_number);
Un1=Un(1:source_number,:);
Un2=Un(source_number+1:sensor_number,:);

T=[1,0,0,0,0,0,0,0,0,0,0]';
c=Un1*inv(Un2)*T;


c=[1,c(2,1),c(1,1)];
f=roots(c);
esource_doa=asin(angle(f)/pi)*180/pi;
esource_doa=sort(esource_doa);
disp('esource_doa');
disp(esource_doa);