clear
%取樣點
%定義常數值
N=32;
a(1)=-0.850848;d2=0.101043;
f1=0.05;f2=0.40;f3=0.42;
ur=0.5*d2.*randn(1,N);
ui=0.5*d2.*randn(1,N);
u=ur+ui*i;
%定義32個復數點
z(1)=u(1);
x(1)=6+z(1);
for n=2:N
    z(n)=-a(1)*z(n-1)+u(n);
    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
%定義f范圍
fmin=-0.5;fstep=0.001;fmax=0.5;
f=fmin:fstep:fmax;
nf=(fmax-fmin)/fstep;  
t=sqrt(-1);
%初值
rxx=0;p0=zeros(1,11);ef=zeros(11,N);eb=zeros(11,N);a=zeros(10,10);
for n=1:N
    rxx=rxx+(abs(x(n)))^2;
end
rxx=(1/N)*rxx;p0(1)=rxx;
ef(1,:)=x;eb(1,:)=x;ef(1,1)=0;eb(1,32)=0;
%算法
p=10;kk=zeros(1,10);
for k=1:p
    e1=0;e2=0;
    for n=(k+1):N
        e1=e1+ef(k,n)*(conj(eb(k,n-1)));
        e2=e2+(abs(ef(k,n))^2+abs(eb(k,n-1))^2);
        kk(k)=(-2)*e1/e2;
        ef(k+1,n)=ef(k,n)+kk(k)*eb(k,n-1);
        eb(k+1,n)=eb(k,n-1)+conj(kk(k))*ef(k,n);
    end
    for i=1:(k-1)
        a(k,i)=a(k-1,i)+kk(k)*conj(a(k-1,k-i));
    end
      a(k,k)=kk(k);
       p0(k+1)=(1-abs(kk(k))^2)*p0(k);
end
%功率譜
for j=1:nf+1
    pxx=0;
    for k=1:p
        pxx=pxx+a(10,k)*exp(-t*2*pi*f(j)*k);
    end
    pxxf(j)=d2/[(abs(1+pxx))^2];
    pxxf(j)=10*log10(pxxf(j));
end
plot(f,pxxf)
title('伯格算法');