随机信号的Levinson处理的程序，附上原程序
clcclear 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=;
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
因为误差功率肯定不可能为负的，请大侠们帮我看看，程序错在哪里！谢谢了
页:
[1]