for q=1:50
N=1024;
n=1:N;
w1=0.2*pi;w2=0.3*pi;w3=0.8*pi;w4=1.2*pi;
A1=1;A2=2;A3=3;A4=4;
x=A1*exp(j*(n-1)*w1)+A2*exp(j*(n-1)*w2)+A3*exp(j*(n-1)*w3)+A4*exp(j*(n-1)*w4)+randn(1,N);
N1=64;
for i=1:N1
    r(i)=1/N*x(:,i:N-1)*x(:,1:N-i)';
end
D=zeros(N1,N1);
for k=1:N1
    for l=1:N1
        if k-l==0,
            D(k,l)=r(1);
        elseif  k-l>0,
            D(k,l)=r(k-l+1);
        else
            D(k,l)=conj(r(l-k+1));
        end
    end
end
[V,L]=eig(D);
lam=diag(L);
for m=1:N1-1
    max=m;
    for n=m+1:N1
        if lam(max)<lam(n),
            max=n;
            temp=lam(m);
            lam(m)=lam(n);
            lam(n)=temp;
            Temp=V(:,m);
            V(:,m)=V(:,n);
            V(:,n)=Temp;
        end
    end
end
Vn=V(:,5:N1);
Pn=Vn*(Vn)';
U1=zeros(1,N1-1);
U1=cat(2,1,U1);
gama=U1*Pn*U1';
a=Pn*U1'/gama;
K=2000;
k=0:K;
W=2*pi*k/(K+1);
m=1:length(a);
fenmu=(conj(a))'*exp(-j*(m-1)'*W);
P=1./(abs(fenmu)).^2;
subplot(2,2,1);
xlabel('W/pi');
ylabel('10*log10(P) 0.1');

plot(W/pi,10*log10(P));
grid on;
hold on;
end