%  the  data of signal
clear;
format long;
c=3*10.^8;%光速
L=9;%陣元數
sam=1024;%取樣點
N=1024;
w1=1;%第一個信號的角頻率
w2=1.1;%第二個信號的角頻率
phase1=0;%初始相位
phase2=0;
  
snr1=0;
snr2=0;
Amp1=sqrt(10^(snr1/10));%信號幅度
Amp2=sqrt(10^(snr2/10));

%lamta=c/w1;
%len=lamta/2;
len=1/2;%陣元間距半波長
P1=30*pi/180;
P2=60*pi/180;
%P;

for t=1:sam,
    s1(t)=Amp1*(exp(-j*(2*pi*w1*0.1*t+phase1)));
    s2(t)=Amp2*(exp(-j*(2*pi*w2*0.1*t+phase2)));
end
tt=1:sam;
%plot(tt,s1,'r--',tt,s2,'b--');

  s=[s1(1:N);s2(1:N)];

  %n=1:1800;
  %aa=(0.1*n)*pi/180;
  
  i=1:L;
  x1=exp(-j*2*pi*(len*(i-1)*sin(P1)));%陣列上關于第一個信號的方向向量
  x2=exp(-j*2*pi*(len*(i-1)*sin(P2)));%陣列上關于第二個信號的方向向量

  a1=x1.';
  a2=x2.';
  
  a=[a1,a2];
  
  %a=[a1.';a2.'];
   
 noise=randn(L,N)+j*randn(L,N);%零均值、方差為 的白噪聲，且與信號源不相關
 
 z=a*s+noise;%總的陣列輸出向量
 
 Rz=(z*z')/N;%取z的自相關矩陣
 
 [e,v]=eig(Rz);%對Rz進行特征值分解

 es=e(:,(L-1:L));%信號子空間
 en=e(:,(1:L-2));%噪聲子空間
  
aaaa=zeros(L,9000);
 for k=1:L;
     for h=1:9000;
         aaaa(k,h)=exp(-j*pi*(k-1)*sin(0.01*h*pi/180));
         %aaaa(k,h)=exp(j*2*pi*w1*ln*(k-1)*sein(0.1*h*pi/180)/c);
     end
 end
b=eye(L);
 c=en*en';
 for m=1:9000;
     aac=aaaa(:,m);
     pp=aac'*aac/(aac'*c*aac); 
     p(m)=real(pp);
 end
 
 p1=p(1:4500);
 p2=p(4501:9000);
[m1,n1]=max(p1);
doa1=n1/100
[m2,n2]=max(p2);
doa2=n2/100+45
u=0:0.01:90-0.01;
p=10*log10(p);
plot(u,p,'r-');
grid on;