clear
N=200;
w(1)=0;
x(1)=5;
a=1;
c=1;
Q1 = randn(1,N)*1;%过程噪声
Q2 = randn(1,N);%测量噪声
for k=2:N;x(k)=a*x(k-1)+Q1(k-1);end%状态矩阵
for k=1:N;Y(k)=c*x(k)+Q2(k);end
p(1)=10;
s(1)=1;
for t=2:N;
Rww=cov(Q1(1:t));
Rvv=cov(Q2(1:t));
p1(t)=a.^2*p(t-1)+Rww;
b(t)=c*p1(t)/(c.^2*p1(t)+Rvv);%kalman增益
s(t)=a*s(t-1)+b(t)*(Y(t)-a*c*s(t-1)); p(t)=p1(t)-c*b(t)*p1(t); end
t=1:N;
plot(t,s,'r',t,Y,'g',t,x,'b');%红色卡尔曼,绿色观测值,蓝色状态值 legend('kalman estimate','ovservations','truth');
下面的链接地址就有,希望对您有帮助。
http://wenku.baidu.com/link?url=VExBsM4W1bg9adNs7I6H3EqAST9Jk7GYLpS0XRR0qOyqLucNxa6Jr4240kNsB6CSwVOfLQxPTIJ5I4Q-JM3Uib99uW0gpF8cGvemX7VXPxq