%產生信號序列
n=linspace(-pi,pi,1024);
xn=2*sin(n*40)+cos(n*100)+randn(size(n));
N=length(n);
M=round(N/10);
K=4;
%求R矩陣
for i=M:N-1
    x=xn(i:-1:i-M+1);
    a=x';
    y(:,i)=a;
end
r=0;
for i=M:N-1
    r=r+y(:,i)*(y(:,i))';
end
R=r/(N);
[v,d]=eig(R);
G=v(:,1:M-K);
GG=G*G';
%求系數
gg(1:2*M-1)=0;
 for i=1:M
     for j=1:M
         gg(i-j+M)=gg(i-j+M)+GG(i,j);
     end
 end
gg1=fliplr(gg);
r=roots(gg1);
for i=1:2*M-2
    r1(i)=abs(1-abs(r(i)));
end
%搜索排序
[xs,idx]=sort(r1);
r3(1:2*K)=r(idx(1:2*K));
for i=1:2*K
   af(i)=phase(r3(i))/2/pi*(N+1);
   ap(i)=abs(r3(i));   
end
af