%  the  data of signal
clear;
format long;
c=3*10.^8;
L=6;
sam=128*4;
N=128*4;
f0=0.175;
lamta=c/f0;
len=lamta/2;

P1=20*pi/180;
snr=-5;
Amp=sqrt(2*10^(snr/10));
sig=Amp.*fmlin(N,0.1,0.25);
s=sig.';
tt=1:N;
f1=(0.15/N)*tt+0.1;
S=zeros(L,N);
 
  for t=1:N;
      i=1:L;  
      x=exp(j*2*pi*(len*f1(1,t)*(i-1)*sin(P1))/c);
      a=x.';
      S(:,t)=a*s(:,t);
  end
  noise=randn(L,N)+j*randn(L,N);
  z=S+noise;
%[tfrf,t,f]=tfrsp(z(1,:).');

%[m n]=max(tfrf);
%f01=min(n)./N
%f02=max(n)./N
%f0=(f01+f02)./2
  for i=1:L 
      zz=z(i,:).';
      %[tfr,t,f]=tfrwv(zz);
      [tfr,t,f]=tfrwv([z(1,:).',zz]);
      %[WH,rho,theta]=htl(tfr,N,N);
      %[tfr,rtfr,hat]=tfrrsp(zz);
      y=max(tfr);
      yy(i,:)=y;
      %y=max(WH);
      %yy(i,:)=max(y);
  end
%figure(1);
%mesh(t,f,abs(tfr));
%mesh(rho,theta,abs(WH));
%xlabel('極半徑 rho');
%ylabel('角度 theta');

  Rz=(yy*yy')/N;
  [e,v]=eig(Rz);
  J = zhihuan(L);

  if v(1,1)<v(L,L)
     vv=J*v*J;
     ee=J*e*J;
  else  vv=v;
        ee=e;
  end
  es=ee(:,1:1);
  en=ee(:,(2:L));
  
  
  
  %構造多項式
  Un1=ee(1,2:L);
  Un2=ee(2:L,2:L); 
  Li=[1 0 0 0 0];
  c1=Un1*inv(Un2)*Li.';
  ccc=[1 c1];
  R=roots(ccc);
 est1=asin(angle(R)*lamta/(2*pi*len));
 agl1=est1*180/pi
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  aaaa=zeros(L,1800);
  for k=1:L;
      for h=1:1800;
          aaaa(k,h)=exp(j*2*pi*(k-1)*len*f0*sin((0.1*h*pi/180-pi/2))/c);
     end
  end
  b=eye(L);
  cc=b-es*es';
  for m=1:1800;
      aac=aaaa(:,m);
      pp=aac'*aac/((aac)'*cc*aac); 
      p(m)=real(pp);
  end
 
 %p1=p(1:450);
 %p2=p(451:900);
%[m1,n1]=max(p1);
%doa1=n1/10
[m,n]=max(p);
doa=n/10-90
figure(2);
u=-90:0.1:90-0.1;
p=10*log(p/max(p));
plot(u,p,'r--');
grid;