贝叶斯辨识.m

一个用MATLAB实现的贝叶斯辨识的代码

clear
L=60; %四位移位寄存器产生的M序列的长度
y1=1; y2=1; y3=1; y4=0; %四个移位寄存器的输出初始值
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; %M序列的值为"1"时,辨识的输入信号取“-1”
    else u(i)=1; %M序列的值为"0"时,辨识的输入信号取“1”
    end
    y1=x1; y2=x2; y3=x3; y4=x4; %为下一次的输入信号作准备
end
figure(1); subplot(2,1,1); stem(u), grid on %画第一个图形的第一个子图——M序列输入信号
v=randn(1,60); subplot(2,1,2); %产生一组长度为60的正态分布随机噪声
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); %直接给出初始状态P0，即一个充分大的实数单位矩阵
E=0.00000000005; %取相对误差E<=0.000000005为停机准则
c=[c0,zeros(7,14)]; %被辨识参数矩阵的初始值及大小
e=zeros(7,15); %相对误差的初始值及大小
for k=3:20; %开始求K 
    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); %系统在M序列输入下的输出采样信号    
    h1=[-z(k-1),-z(k-2),u(k-1),u(k-2),v(k),v(k-1),v(k-2)]'; %为求K(k)作准备
    x=h1'*p0*h1+rv;
    x1=inv(x);
k1=p0*h1*x1; %K
    d1=z(k)-h1'*c0; %开始求被辨识参数c
    c1=c0+k1*d1; %辨识参数c
    zs(k)=-1.5*z(k-1)+0.7*z(k-2)+u(k-1)+0.5*u(k-2); %系统在M序列输入下的输出响应
    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; %把辨识参数c 列向量加入辨识参数矩阵 
    p1=p0-k1*k1'*[h1'*p0*h1+1];%find p(k)
    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,:); %分离出a1、a2、b1、b2
d1=c(5,:); d2=c(6,:); d3=c(7,:); %分离出d1、 d2、 d3
ea1=e(1,:); ea2=e(2,:); eb1=e(3,:); eb2=e(4,:); %分离出a1、 a2、b1、b2的收敛情况
ed1=e(5,:); ed2=e(6,:); ed3=e(7,:); %分离出d1、 d2、 d3的收敛情况
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('Parameter Identification with Recursive Least Squares Method') %标题
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,ed2,'r+') %画出各个参数收敛情况
title('Identification Error') %标题
%Response %响应
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') %画出模型去除噪声后的输出响应