|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
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
因为误差功率肯定不可能为负的,请大侠们帮我看看,程序错在哪里!谢谢了 |
|