bayes.m

用matlab的bayes方法实现系统辨识

%Bayes辨识算法与仿真
clear
L=60;
y1=1;y2=1;y3=1;y4=0;%初始化M序列值
for i=1:L;
    x1=xor(y3,y4);
    x2=y1;
    x3=y2;
    x4=y3;
    y(i)=y4;
    if y(i)>0.5
        u(i)=-1;
    else u(i)=1;
    end
    y1=x1;y2=x2;y3=x3;y4=x4;
end
figure(1);
subplot(2,1,1),stem(u),grid on 
v=randn(1,60);subplot(2,1,2);plot(v),grid on;%噪声的产生
R=corrcoef(u,v);%计算输入信号与随机噪声信号的相关系数
r=R(1,2)
rv=std(v)*std(v)
u,v
%初始化
z(2)=0;z(1)=0;zs(2)=0;zs(1)=0;zm(2)=0;zm(1)=0;
%Bayes辨识方法
c0=[0.001 0.001 0.001 0.001 0.001 0.001 0.001]';
p0=10^6*eye(7,7);
E=0.00000000005;
c=[c0,zeros(7,14)];
e=zeros(7,15);
for k=3:20;
    z(k)=-1.5*z(k-1)+0.7*z(k-2)+u(k-1)+0.5*u(k-2)+v(k)-v(k-1)+0.2*v(k-2);
    hl=[-z(k-1),-z(k-2),u(k-1),u(k-2),v(k),v(k-1),v(k-2)]';
    x=hl'*p0*hl+rv;
    x1=inv(x);
    k1=p0*hl*x1;%K
    d1=z(k)-hl'*c0;
    c1=c0+k1*d1;
    zs(k)=-1.5*z(k-1)+0.7*z(k-2)+u(k-1)+0.5*u(k-2);
    zm(k)=[-z(k-1),-z(k-2),u(k-1),u(k-2)]*[c1(1);c1(2);c1(3);c1(4)];
    e1=c1-c0;
    e2=e1./c0;%相对误差
    e(:,k)=e2;
    c0=c1;
    c(:,k)=c1;
    p1=p0-k1*k1'*[hl'*p0*hl+1];
    p0=p1;
    %停机标准
    if e2<=E break;
    end
end
c,e,z,zs,zm
%输出结果
a1=c(1,:);a2=c(2,:);b1=c(3,:);b2=c(4,:);d1=c(5,:);d2=c(6,:);d3=c(7,:);
ea1=e(1,:);ea2=e(2,:);eb1=e(3,:);eb2=e(4,:);ed1=e(5,:);ed2=e(6,:);ed3=e(7,:);
figure(2);
i=1:20;
plot(i,a1,'r',i,a2,'r:',i,b1,'b',i,b2,'b:',i,d1,'g',i,d2,'g:',i,d3,'g+');
title('Bayes ID')
figure(3);
i=1:20;
plot(i,ea1,'r',i,ea2,'r:',i,eb1,'b',i,eb2,'b:',i,ed1,'g',i,ed2,'g:',i,ed3,'g+');
title('辨识误差')
figure(4);
subplot(4,1,1);i=1:20;plot(i,zs(i),'r');
subplot(4,1,2);i=1:20;plot(i,z(i),'g');
subplot(4,1,3);i=1:20;plot(i,zm(i),'b');
subplot(4,1,4);i=1:20;plot(i,zmy(i),'b');%除去噪声的输出向量