卡尔曼滤波算法实现代码

1、卡尔曼滤波算法实现代码其c语言实现代码如下：#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,*

2、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

3、+)      for (j=0; j<=n-1; j+)         ii=i*l+j; aii=0.0;          for (kk=0; kk<=n-1; kk+)       &#

4、160;    aii=aii+pi*n+kk*fj*n+kk;            for (i=0; i<=n-1; i+)      for (j=0; j<=n-1; j+)         ii=i*n+j

5、; pii=qii;          for (kk=0; kk<=n-1; kk+)            pii=pii+fi*n+kk*akk*l+j;            for (ii=2;

6、 ii<=k; ii+)       for (i=0; i<=n-1; i+)        for (j=0; j<=m-1; j+)           jj=i*l+j; ajj=0.0;   

7、         for (kk=0; kk<=n-1; kk+)              ajj=ajj+pi*n+kk*hj*n+kk;                

8、  for (i=0; i<=m-1; i+)        for (j=0; j<=m-1; j+)           jj=i*m+j; ejj=rjj;            f

9、or (kk=0; kk<=n-1; kk+)              ejj=ejj+hi*n+kk*akk*l+j;                  js=rinv(e,m);     

10、   if (js=0)            free(e); free(a); free(b); return(js);        for (i=0; i<=n-1; i+)        for 

11、;(j=0; j<=m-1; j+)           jj=i*m+j; gjj=0.0;            for (kk=0; kk<=m-1; kk+)           

12、   gjj=gjj+ai*l+kk*ej*m+kk;                  for (i=0; i<=n-1; i+)           jj=(ii-1)*n+i; xjj=0.0;   

13、;         for (j=0; j<=n-1; j+)              xjj=xjj+fi*n+j*x(ii-2)*n+j;                

14、;  for (i=0; i<=m-1; i+)           jj=i*l; bjj=y(ii-1)*m+i;            for (j=0; j<=n-1; j+)       

15、       bjj=bjj-hi*n+j*x(ii-1)*n+j;                  for (i=0; i<=n-1; i+)           jj=(ii-1)*n+i;

16、0;           for (j=0; j<=m-1; j+)              xjj=xjj+gi*m+j*bj*l;               

17、;   if (ii<k)           for (i=0; i<=n-1; i+)            for (j=0; j<=n-1; j+)        

18、       jj=i*l+j; ajj=0.0;                for (kk=0; kk<=m-1; kk+)                &

19、#160; ajj=ajj-gi*m+kk*hkk*n+j;                if (i=j) ajj=1.0+ajj;                       

20、   for (i=0; i<=n-1; i+)            for (j=0; j<=n-1; j+)               jj=i*l+j; bjj=0.0;   &

21、#160;            for (kk=0; kk<=n-1; kk+)                  bjj=bjj+ai*l+kk*pkk*n+j;        &

22、#160;                 for (i=0; i<=n-1; i+)            for (j=0; j<=n-1; j+)       &

23、#160;       jj=i*l+j; ajj=0.0;                for (kk=0; kk<=n-1; kk+)               &#

24、160;  ajj=ajj+bi*l+kk*fj*n+kk;                          for (i=0; i<=n-1; i+)           

25、; for (j=0; j<=n-1; j+)               jj=i*n+j; pjj=qjj;                for (kk=0; kk<=n-1; kk+)

26、60;                 pjj=pjj+fi*n+kk*aj*l+kk;                            

27、60;     free(e); free(a); free(b);    return(js);  C+实现代码如下：=kalman.h=/ kalman.h: interface for the kalman class./#if !defined(AFX_KALMAN_H_ED3D740F_01D2_4616_8B74_8BF57636F2C0_INCLUDED_)#define 

28、AFX_KALMAN_H_ED3D740F_01D2_4616_8B74_8BF57636F2C0_INCLUDED_#if _MSC_VER > 1000#pragma once#endif / _MSC_VER > 1000#include <math.h>#include "cv.h" class kalman  public: void init_kalman(int x,i

29、nt xv,int y,int yv); CvKalman* cvkalman; CvMat* state;  CvMat* process_noise; CvMat* measurement; const CvMat* prediction; CvPoint2D32f get_predict(float x, float y); kalman(int x=0,int 

30、;xv=0,int y=0,int yv=0); /virtual kalman();#endif / !defined(AFX_KALMAN_H_ED3D740F_01D2_4616_8B74_8BF57636F2C0_INCLUDED_)=kalman.cpp=#include "kalman.h"#include <stdio.h>/* tester de printer toutes les valeurs 

31、;des vecteurs*/* tester de changer les matrices du noises */* replace state by cvkalman->state_post ? */CvRandState rng;const double T = 0.1;kalman:kalman(int x,int xv,int y,int 

32、;yv)         cvkalman = cvCreateKalman( 4, 4, 0 );    state = cvCreateMat( 4, 1, CV_32FC1 );    process_noise = cvCreateMat( 4, 1, 

33、;CV_32FC1 );    measurement = cvCreateMat( 4, 1, CV_32FC1 );    int code = -1;        /* create matrix data */     const 

34、;float A =     1, T, 0, 0,   0, 1, 0, 0,   0, 0, 1, T,   0, 0, 0, 1            const float 

35、;H =      1, 0, 0, 0,    0, 0, 0, 0,   0, 0, 1, 0,   0, 0, 0, 0              const

36、60;float P =     pow(320,2), pow(320,2)/T, 0, 0,   pow(320,2)/T, pow(320,2)/pow(T,2), 0, 0,   0, 0, pow(240,2), pow(240,2)/T,   0, 0, pow(240,2)/T, pow(240,2)/p

37、ow(T,2)         const float Q =    pow(T,3)/3, pow(T,2)/2, 0, 0,   pow(T,2)/2, T, 0, 0,   0, 0, pow(T,3)/3, pow(T,2)/2,   0, 

38、0, pow(T,2)/2, T           const float R =    1, 0, 0, 0,   0, 0, 0, 0,   0, 0, 1, 0,   0, 0, 0,

39、0;0              cvRandInit( &rng, 0, 1, -1, CV_RAND_UNI );    cvZero( measurement );        cvRandSetRange( &rng,

40、60;0, 0.1, 0 );    rng.disttype = CV_RAND_NORMAL;    cvRand( &rng, state );    memcpy( cvkalman->transition_matrix->data.fl, A, sizeof(A);    memcpy( 

41、cvkalman->measurement_matrix->data.fl, H, sizeof(H);    memcpy( cvkalman->process_noise_cov->data.fl, Q, sizeof(Q);    memcpy( cvkalman->error_cov_post->data.fl, P, sizeof(P);    mem

42、cpy( cvkalman->measurement_noise_cov->data.fl, R, sizeof(R);    /cvSetIdentity( cvkalman->process_noise_cov, cvRealScalar(1e-5) );        /cvSetIdentity( cvkalman->error_cov_post, cvRealSc

43、alar(1); /cvSetIdentity( cvkalman->measurement_noise_cov, cvRealScalar(1e-1) );    /* choose initial state */    state->data.fl0=x;    state->data.fl1=xv;    state->dat

44、a.fl2=y;    state->data.fl3=yv;    cvkalman->state_post->data.fl0=x;    cvkalman->state_post->data.fl1=xv;    cvkalman->state_post->data.fl2=y;    cvkalman->state_post->data

45、.fl3=yv; cvRandSetRange( &rng, 0, sqrt(cvkalman->process_noise_cov->data.fl0), 0 );    cvRand( &rng, process_noise );         CvPoint2D32f kalman:get_predict(float x,

46、 float y)        /* update state with current position */    state->data.fl0=x;    state->data.fl2=y;        /* predict poin

47、t position */    /* x'k=A鈥k+B鈥k       P'k=A鈥k-1*AT + Q */    cvRandSetRange( &rng, 0, sqrt(cvkalman->measurement_noise_cov->data.fl0), 0 );  &

48、#160; cvRand( &rng, measurement );         /* xk=A?xk-1+B?uk+wk */    cvMatMulAdd( cvkalman->transition_matrix, state, process_noise, cvkalman->state_post );       &#