%  the  data of signal
clear;
format long;
c=3*10.^8;%光速
L=10;%陣元數(shù)
sam=500;%取樣點(diǎn)
N=500;
dx=1/2; %陣元間距
dy=1/2;
w1=0.93;%第一個信號的角頻率
w2=0.98;%第二個信號的角頻率
w3=1.03;
w4=1.08;
phase1=0;%初始相位
phase2=0;
phase3=0;
phase4=0;
singnum=4;
snr1=30;
Amp=sqrt(10^(snr1/10));%信號幅度
cx=0.3527+0.4854j;
cy=0.3527+0.4854j; %互相關(guān)系數(shù)
cxy=0.0927+0.2853j;
theta=[74,35;40,20;54,66;28,41];
%lamta=c/w1;
%len=lamta/2;


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

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

  %n=1:1800;
  %aa=(0.1*n)*pi/180;
  
  i=1:L;
  ax1=exp(-j*2*pi*(dx*(i-1)*cos(theta(1,1)*sin(theta(1,2)))))';%陣列上關(guān)于信號的方向向量
  ax2=exp(-j*2*pi*(dx*(i-1)*cos(theta(2,1)*sin(theta(2,2)))))';
  ax3=exp(-j*2*pi*(dx*(i-1)*cos(theta(3,1)*sin(theta(3,2)))))';
  ax4=exp(-j*2*pi*(dx*(i-1)*cos(theta(4,1)*sin(theta(4,2)))))';
  ay1=exp(-j*2*pi*(dx*(i-1)*cos(theta(1,1)*sin(theta(1,2)))))';
  ay2=exp(-j*2*pi*(dx*(i-1)*cos(theta(2,1)*sin(theta(2,2)))))';
  ay3=exp(-j*2*pi*(dx*(i-1)*cos(theta(3,1)*sin(theta(3,2)))))';
  ay4=exp(-j*2*pi*(dx*(i-1)*cos(theta(4,1)*sin(theta(4,2)))))';
  a=[kron(ay1,ax1),kron(ay2,ax2),kron(ay3,ax3),kron(ay4,ax4)];
  
  c1=toeplitz([1,cx,zeros(1,L-2)]);
  c2=toeplitz([cx,cxy,zeros(1,L-2)]);
  C1=kron(eye(L),c1);
  C2=kron((tril(ones(L),1)-tril(ones(L),-2)-eye(L)),c2);
  C=C1+C2;
  %a=[a1.';a2.'];
  J=[zeros(L-2,1),eye(L-2),zeros(L-2,1)];
  P1=J;
  P=kron(P1,J);
 noise=randn((L-2)*(L-2),N)+j*randn((L-2)*(L-2),N);%零均值、方差為 的白噪聲，且與信號源不相關(guān)
 
 z=P*C*a*s+noise;%總的陣列輸出向量
 
 Rz=(z*z')/N;%取z的自相關(guān)矩陣
 
 [e,v]=eig(Rz);%對Rz進(jìn)行特征值分解

 es=e(:,(L-3:L));%信號子空間
 en=e(:,(1:L-4));%噪聲子空間
 c=en*en';
 b=eye(L);
aaaa=zeros(L,90);
  k=1:L;
     for the1=1:90
         for the2=1:90
         ax=exp(-j*2*pi*(dx*(k-1)*cos(the1)*sin(the2)));
         ay=exp(-j*2*pi*(dx*(k-1)*cos(the1)*sin(the2)));
         %aaaa(k,h)=exp(j*2*pi*w1*ln*(k-1)*sein(0.1*h*pi/180)/c);
         ac=kron(ay',ax');


     pp=ac'*ac/(ac'*c*ac); 
     p(the1,the2)=real(pp);
 
 end
     end
[m1,n1]=max(p(:));
doa1=floor(n1/90);

doa2=n1-floor(n1/90)*90;
u=0:0.01:90-0.01;
p=10*log10(p);
plot(u,p,'r-');
grid on;