clear;
p=26;
pi=3.1415;
f1=0.05;
f2=0.3;
f3=0.42;
a1=-0.850848;
a=0; 
b=sqrt(0.101043/2); 
% 噪聲
u1=randn(1,32); 
u1=u1/std(u1); 
u1=u1-mean(u1); 
u1=a+b*u1; 

u2=randn(1,32); 
u2=u2/std(u2); 
u2=u2-mean(u2); 
u2=a+b*u2; 
u=u1+u2*i;

% 信號產生
z=zeros(1,32);
x=zeros(1,32);
z(1)=u(1);
for n=2:1:32,z(n)=-a1*z(n-1)+u(n);end,z;
for n=1:1:32,x(n)=2*cos(2*pi*f1*(n-1))+2*cos(2*pi*f2*(n-1))+2*cos(2*pi*f3*(n-1))+z(n);end,x;

% 求自相關矩陣Rxx
rxx=zeros(32,1);
for k=0:31
    for n=1:(32-k)
        rxx(k+1)=rxx(k+1)+conj(x(n))*x(n+k);
    end
end
rxx=rxx/32;
rxx=rxx';%????????????
Rxx=toeplitz(rxx);

% 特征分解
[V,D]=eig(Rxx);

f=-0.5:0.001:0.5;
e=zeros(32,1001);
for k=1:32
    e(k,:)=exp(i*2*pi*f*(k-1));
end

Pmu=zeros(1,1001);
for m=1:1001
    for n=1:p
        Pmu(m)=Pmu(m)+abs(e(:,m)'*V(:,n))^2;
    end
end
Pmu=1./Pmu;
Pmu=10*log10(Pmu);
plot(f,Pmu);