随机信号的Levinson处理的程序，附上原程序

yjvfhp1211

clc
clear all
Gx0=256;                  %参考空间频率n0下的路面功率谱密度
n0=0.1;                   %参考空间频率n0
L=409.6;
l=0.1;
N=L/l;                     %采样点数
n1=0.01;                   %空间频率范围n1--nu
nu=3;
w=2;                       %频率指数
no=1/L;                    %空间频率间隔
Xk=[];
Xm=[];
xm=[];
n=linspace(n1,nu,N/2+1);
GxC=Gx0*(n/n0).^(-w);
k=0:N/2;
fik=unifrnd(0,2*pi,1,N);                            %产生0到2pi的均匀分布的随机序列
pg=GxC(1:N/2+1);
Xk=sqrt((N/2+1)/(2*l)*pg).*exp(j*fik(1:N/2+1));  %调用函数GxC（n）
Xk(1)=sqrt((N/2+1)/(2*l)*pg(1)); Xk(N/2+1)=sqrt((N/2+1)/(2*l)*pg(N/2+1));
Xk=[Xk(1:N/2+1) conj(Xk(N/2:-1:2))];
Xm=ifft(Xk);%逆傅立叶变换后得到复数形式随机序列
x=real(Xm);
N1=150;
for m=1:N1
    R(m)=0;
end
for m=1:N1
    s=0;
    for n=1:N1-(m-1)
        s=s+x(n)*x(n+m-1);
    end
    R(m)=s/N1;
end
%迭代
for M=1:N1-5%阶数M要小于数据长度N，这里取M=N-5
    a(M,M)=0;
    FPE(M)=0;
   P(M)=0;
end
%（一)初始化：由Levinson公式计算一阶模型参数a(1，1)和误差功率P（1）及FPE（1)；
a(1,1)=-R(2)/R(1); %计算一阶模型参数a(1,1)；
P(1)=R(1)*(1-(a(1,1))^2); %计算一阶模型的误差功率P(1)；
FPE(1)=P(1)*(N1+1+1)/(N1-1-1); %计算一阶模型的最终预测误差准则函数FPE(1)；
%（二）递推：设已经计算出M-1阶时a(M-1,k)和P(M-1),计算M阶模型参数a(M,M),a(M,k)和误差功率P（M）及FPE（M)；
sum=0;
for M=2:N1-5
    for k=1:M-1
        sum=sum+a(M-1,k)*R(M-k+1);
    end
    a(M,M)=-(R(M+1)+sum)/P(M-1);%计算反射系数
    for k=1:M-1
        a(M,k)=a(M-1,k)+a(M,M)*a(M-1,M-k);
    end
   P(M)=(1-(abs(a(k,k)))^2)*P(M-1);
    FPE(M)=P(M)*(N1+M+1)/(N1-M-1);
    sum=0;
end %递推计算完毕；
%(三)Model阶次确定：求出使FPE最小时的阶数M1；
min=FPE(1);
for M=2:N-5
    if FPE(M)<min
        min=FPE(M);
        M1=M;
    end
end
disp('阶数M为：')
disp(M1);
disp('模型参数a为：');
for i=1:M1
    disp(a(M1,i));
end
disp('误差功率为：');
disp(P(M1));

我得出的结果：阶数M为：
     7
模型参数a为：
-4.9685e+005
-8.2513e+007
  1.1251e+008
-1.0643e+008
  8.6858e+007
-1.0587e+007
  2.0119e+001
误差功率为：
-3.3531e+012

因为误差功率肯定不可能为负的，请大侠们帮我看看，程序错在哪里！谢谢了