声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1244|回复: 0

[综合] 随机信号的Levinson处理的程序,附上原程序

[复制链接]
发表于 2009-3-18 09:32 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

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

因为误差功率肯定不可能为负的,请大侠们帮我看看,程序错在哪里!谢谢了
回复
分享到:

使用道具 举报

您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-26 09:36 , Processed in 0.060370 second(s), 19 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表