clear;
format long;
c=3*10.^8;
L=10;
sam=128*2;
N=128*2;
w1=1;
w2=1;
phase=0;

snr=10;
Amp=sqrt(2*10^(snr/10));
lamta=c/w1;
len=lamta/2;

P1=42*pi/180;
P2=47*pi/180;
%P;

for t=1:sam,
    s1(t)=Amp*(exp(j*2*pi*w1*0.1*t));
    s2(t)=Amp*(exp(j*2*pi*w2*0.1*t));
end
  s=[s1(1:N);s2(1:N)];
%t=1:sam ; 
%plot(t,s1)

  %n=1:1800;
  %aa=(0.1*n)*pi/180;
  
  i=1:L;
  x1=exp(j*2*pi*w1*(len*(i-1)*sin(P1))/c);
  x2=exp(j*2*pi*w2*(len*(i-1)*sin(P2))/c);

  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;
 J = zhihuan(L)
 Rb=J*(Rz').'*J;
 Rfb=0.50*(Rz+Rb);
 [e,v]=eig(Rfb);

 es=e(:,1:2);
 en=e(:,(3:L));
   
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*len*(k-1)*sin(0.1*h*pi/180)/c);
     end
 end
b=zeros(L);
for i=1:L;
 b(i,i)=1;
end
 c=b-es*es';
 for m=1:9000;
     aac=aaaa(:,m);
     pp=aac'*aac/((aac)'*c*aac); 
     %pp=1/((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;
plot(u,p);