kalman滤波程序(matlab实现).txt

kalman filter；matlab implementation

附matlab下面的kalman滤波程序： 

clear 
N=200; 
w(1)=0; 
w=randn(1,N) 
x(1)=0; 
a=1; 
for k=2:N; 
x(k)=a*x(k-1)+w(k-1); 
end 


V=randn(1,N); 
q1=std(V); 
Rvv=q1.^2; 
q2=std(x); 
Rxx=q2.^2; 
q3=std(w); 
Rww=q3.^2; 
c=0.2; 
Y=c*x+V; 

p(1)=0; 
s(1)=0; 
for t=2:N; 
p1(t)=a.^2*p(t-1)+Rww; 
b(t)=c*p1(t)/(c.^2*p1(t)+Rvv); 
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'); 



算法实现
X(k) = F(k,k-1)·X(k-1)+T(k,k-1)·U(k-1) 
Y(k) = H(k)·X(k)+N(k) 
其中 
X(k)和Y(k)分别是k时刻的状态矢量与观测矢量 
F(k,k-1)为状态转移矩阵 
U(k)为k时刻动态噪声 
T(k,k-1)为系统控制矩阵 
H(k)为k时刻观测矩阵 
N(k)为k时刻观测噪声 

则卡尔曼滤波的算法流程为： 

X(k)^= F(k,k-1)·X(k-1) 

计算协方差矩阵 
C(k)^=F(k,k-1)×C(k)×F(k,k-1)'+T(k,k-1)×Q(k)×T(k,k-1)' 
Q(k) = U(k)×U(k)' 

计算卡尔曼增益矩阵 
K(k) = C(k)^×H(k)'×[H(k)×C(k)^×H(k)'+R(k)]^(-1) 
R(k) = N(k)×N(k)' 

更新
X(k)~=X(k)^+K(k)×[Y(k)-H(k)×X(k)^] 

计算更新后协防差矩阵 
C(k)~ = [I-K(k)×H(k)]×C(k)^×[I-K(k)×H(k)]'+K(k)×R(k)×K(k)' 

X(k+1) = X(k)~ 
C(k+1) = C(k)~ 
重复以上步骤 

#include "stdlib.h" 
#include "rinv.c" 
int lman(n,m,k,f,q,r,h,y,x,p,g) 
int n,m,k; 
double f[],q[],r[],h[],y[],x[],p[],g[]; 
{ int i,j,kk,ii,l,jj,js; 
double *e,*a,*b; 
e=malloc(m*m*sizeof(double)); 
l=m; 
if (l<n) l=n; 
a=malloc(l*l*sizeof(double)); 
b=malloc(l*l*sizeof(double)); 
for (i=0; i<=n-1; i++) 
for (j=0; j<=n-1; j++) 
{ ii=i*l+j; a[ii]=0.0; 
for (kk=0; kk<=n-1; kk++) 
a[ii]=a[ii]+p[i*n+kk]*f[j*n+kk]; 
} 
for (i=0; i<=n-1; i++) 
for (j=0; j<=n-1; j++) 
{ ii=i*n+j; p[ii]=q[ii]; 
for (kk=0; kk<=n-1; kk++) 
p[ii]=p[ii]+f[i*n+kk]*a[kk*l+j]; 
} 
for (ii=2; ii<=k; ii++) 
{ for (i=0; i<=n-1; i++) 
for (j=0; j<=m-1; j++) 
{ jj=i*l+j; a[jj]=0.0; 
for (kk=0; kk<=n-1; kk++) 
a[jj]=a[jj]+p[i*n+kk]*h[j*n+kk]; 
} 
for (i=0; i<=m-1; i++) 
for (j=0; j<=m-1; j++) 
{ jj=i*m+j; e[jj]=r[jj]; 
for (kk=0; kk<=n-1; kk++) 
e[jj]=e[jj]+h[i*n+kk]*a[kk*l+j]; 
} 
js=rinv(e,m); 
if (js==0) 
{ free(e); free(a); free(b); return(js);} 
for (i=0; i<=n-1; i++) 
for (j=0; j<=m-1; j++) 
{ jj=i*m+j; g[jj]=0.0; 
for (kk=0; kk<=m-1; kk++) 
g[jj]=g[jj]+a[i*l+kk]*e[j*m+kk]; 
} 
for (i=0; i<=n-1; i++) 
{ jj=(ii-1)*n+i; x[jj]=0.0; 
for (j=0; j<=n-1; j++) 
x[jj]=x[jj]+f[i*n+j]*x[(ii-2)*n+j]; 
} 
for (i=0; i<=m-1; i++) 
{ jj=i*l; b[jj]=y[(ii-1)*m+i]; 
for (j=0; j<=n-1; j++) 
b[jj]=b[jj]-h[i*n+j]*x[(ii-1)*n+j]; 
} 
for (i=0; i<=n-1; i++) 
{ jj=(ii-1)*n+i; 
for (j=0; j<=m-1; j++) 
x[jj]=x[jj]+g[i*m+j]*b[j*l]; 
} 
if (ii<k) 
{ for (i=0; i<=n-1; i++) 
for (j=0; j<=n-1; j++) 
{ jj=i*l+j; a[jj]=0.0; 
for (kk=0; kk<=m-1; kk++) 
a[jj]=a[jj]-g[i*m+kk]*h[kk*n+j]; 
if (i==j) a[jj]=1.0+a[jj]; 
} 
for (i=0; i<=n-1; i++) 
for (j=0; j<=n-1; j++) 
{ jj=i*l+j; b[jj]=0.0; 
for (kk=0; kk<=n-1; kk++) 
b[jj]=b[jj]+a[i*l+kk]*p[kk*n+j]; 
} 
for (i=0; i<=n-1; i++) 
for (j=0; j<=n-1; j++) 
{ jj=i*l+j; a[jj]=0.0; 
for (kk=0; kk<=n-1; kk++) 
a[jj]=a[jj]+b[i*l+kk]*f[j*n+kk]; 
} 
for (i=0; i<=n-1; i++) 
for (j=0; j<=n-1; j++) 
{ jj=i*n+j; p[jj]=q[jj]; 
for (kk=0; kk<=n-1; kk++) 
p[jj]=p[jj]+f[i*n+kk]*a[j*l+kk]; 
} 
} 
} 
free(e); free(a); free(b); 
return(js); 
}