卡尔曼滤波算法(C--C++两种实现代码)
无迹卡尔曼滤波c语言

无迹卡尔曼滤波c语言卡尔曼滤波是一种高效的递归滤波器,它使用状态方程和测量方程来估计系统的状态。
在C语言中实现卡尔曼滤波器需要定义状态变量、测量变量、状态转移矩阵、测量矩阵和过程噪声协方差矩阵等。
以下是一个简单的无迹卡尔曼滤波器的C语言实现示例:```cinclude <>include <>define Q // 过程噪声协方差矩阵define R // 测量噪声协方差矩阵typedef struct {double x; // 状态变量double P; // 状态协方差矩阵double K; // 卡尔曼增益double Q; // 过程噪声协方差矩阵double R; // 测量噪声协方差矩阵} KalmanFilter;void KalmanInit(KalmanFilter kf, double init_x, double init_P) {kf->x = init_x;kf->P = init_P;kf->K = 0;}void KalmanUpdate(KalmanFilter kf, double measurement) {// 预测步骤double P_predict = kf->P + Q; // P的预测值double K = P_predict / (P_predict + R); // 卡尔曼增益double x_predict = kf->x + K (measurement - kf->x); // x的预测值 double P_update = (1 - K) P_predict; // P的更新值// 更新状态变量和协方差矩阵kf->x = x_predict;kf->P = P_update;}int main() {KalmanFilter kf;double measurements[10] = {, , , , , , , , , }; // 测量数据double true_x = ; // 真实状态变量double init_x = true_x; // 初始状态变量值double init_P = ; // 初始协方差矩阵值KalmanInit(&kf, init_x, init_P); // 初始化卡尔曼滤波器for (int i = 0; i < 10; i++) {KalmanUpdate(&kf, measurements[i]); // 进行卡尔曼滤波更新步骤printf("State variable: %f\n", ); // 输出估计的状态变量值true_x += ; // 真实状态变量进行微小变化,模拟真实情况下的状态变化}return 0;}```。
卡尔曼滤波器实现代码

VectorXd k = Cal_kalmen_gain_1_2(P_predict_i, H, Rnoise);
VectorXd x_predict_calibration_i = Correctpredcit_x_2_1(y, x_predict_i, k);//F is the relation of this time and next time
printf("hello world\n"); VectorXd my_vector(2); my_vector << 10, 20; MatrixXd my_matrix(2, 2); my_matrix << 1, 2, 3, 4; cout << my_matrix << endl; cout << my_vector << endl;
R = MatrixXd(1, 1); R << 1;
I = MatrixXd::Identity(2, 2);
Q = MatrixXd(2, 2); Q << 0, 0, 0, 0;
//create a list of measurements VectorXd single_meas(1); single_meas << 1; measurements.push_back(single_meas); single_meas << 5; measurements.push_back(single_meas); single_meas << 9; measurements.push_back(single_meas); single_meas << 13; measurements.push_back(single_meas); single_meas << 17; measurements.push_back(single_meas); single_meas << 21; measurements.push_back(single_meas); single_meas << 25; measurements.push_back(single_meas); single_meas << 29; measurements.push_back(single_meas); filter(x, P,u,F,H,R,I, Q,8); ; }
自适应卡尔曼滤波c语言实现

自适应卡尔曼滤波C语言实现1. 什么是卡尔曼滤波?卡尔曼滤波是一种用于估计系统状态的算法,它基于线性动态系统模型和高斯噪声假设。
卡尔曼滤波器通过不断地更新状态估计值,将测量结果和系统动态进行融合,提供更准确的状态估计。
在实际应用中,传感器测量值常常包含噪声或者不完全准确。
卡尔曼滤波通过对测量值进行加权平均,同时考虑系统的动态模型和测量噪声的特性,可以有效地抑制噪声并提高状态估计的精度。
2. 卡尔曼滤波的基本原理卡尔曼滤波器由两个步骤组成:预测步骤和更新步骤。
预测步骤预测步骤用于根据当前时刻的状态估计值和系统动态模型,预测下一时刻的状态估计值。
假设我们有一个状态向量x表示系统的状态,在每个时刻t,我们可以使用状态转移矩阵A将当前时刻的状态向量x(t)预测到下一时刻的状态向量x(t+1):x(t+1) = A * x(t)同时,我们还需要考虑过程噪声的影响。
假设过程噪声服从均值为0、协方差矩阵为Q的高斯分布,我们可以使用协方差矩阵Q来描述过程噪声的特性。
预测步骤可以表示为:P(t+1|t) = A * P(t|t) * A' + Q其中,P(t|t)表示当前时刻的状态估计误差协方差矩阵,P(t+1|t)表示预测步骤中的状态估计误差协方差矩阵。
更新步骤更新步骤用于根据当前时刻的测量值和预测步骤得到的状态估计值,更新系统状态的估计。
假设我们有一个观测向量z表示系统的观测值,在每个时刻t,我们可以使用观测模型C将当前时刻的状态向量x(t)映射到观测空间中:z(t) = C * x(t)同时,我们还需要考虑观测噪声的影响。
假设观测噪声服从均值为0、协方差矩阵为R的高斯分布,我们可以使用协方差矩阵R来描述观测噪声的特性。
更新步骤可以表示为:K(t+1) = P(t+1|t) * C' * (C * P(t+1|t) * C' + R)^(-1)x(t+1|t+1) = x(t+1|t) + K(t+1) * (z(t+1) - C * x(t+1|t))P(t+1|t+1) = (I - K(t+1) * C) * P(t+1|t)其中,K(t+1)表示卡尔曼增益,x(t+1|t+1)表示更新步骤中的状态估计值,P(t+1|t+1)表示更新步骤中的状态估计误差协方差矩阵。
卡尔曼滤波算法及C语言代码

卡尔曼滤波简介及其算法实现代码卡尔曼滤波算法实现代码(C,C++分别实现)卡尔曼滤波器简介近来发现有些问题很多人都很感兴趣。
所以在这里希望能尽自己能力跟大家讨论一些力所能及的算法。
现在先讨论一下卡尔曼滤波器,如果时间和能力允许,我还希望能够写写其他的算法,例如遗传算法,傅立叶变换,数字滤波,神经网络,图像处理等等。
因为这里不能写复杂的数学公式,所以也只能形象的描述。
希望如果哪位是这方面的专家,欢迎讨论更正。
卡尔曼滤波器– Kalman Filter1.什么是卡尔曼滤波器(What is the Kalman Filter?)在学习卡尔曼滤波器之前,首先看看为什么叫“卡尔曼”。
跟其他著名的理论(例如傅立叶变换,泰勒级数等等)一样,卡尔曼也是一个人的名字,而跟他们不同的是,他是个现代人!卡尔曼全名Rudolf Emil Kalman,匈牙利数学家,1930年出生于匈牙利首都布达佩斯。
1953,1954年于麻省理工学院分别获得电机工程学士及硕士学位。
1957年于哥伦比亚大学获得博士学位。
我们现在要学习的卡尔曼滤波器,正是源于他的博士论文和1960年发表的论文《A New Approach to Linear Filtering and Prediction Problems》(线性滤波与预测问题的新方法)。
如果对这编论文有兴趣,可以到这里的地址下载:/~welch/media/pdf/Kalman1960.pdf。
简单来说,卡尔曼滤波器是一个“optimal recursive data processing algorithm(最优化自回归数据处理算法)”。
对于解决很大部分的问题,他是最优,效率最高甚至是最有用的。
他的广泛应用已经超过30年,包括机器人导航,控制,传感器数据融合甚至在军事方面的雷达系统以及导弹追踪等等。
近年来更被应用于计算机图像处理,例如头脸识别,图像分割,图像边缘检测等等。
2.卡尔曼滤波器的介绍(Introduction to the Kalman Filter)为了可以更加容易的理解卡尔曼滤波器,这里会应用形象的描述方法来讲解,而不是像大多数参考书那样罗列一大堆的数学公式和数学符号。
(完整word版)卡尔曼滤波算法(C--C++两种实现代码)

卡尔曼滤波算法实现代码C++实现代码如下: ============================kalman.h================= =============== // kalman.h: interface for the kalman class. // ////////////////////////////////////////////////////////////////////// #if !defined(AFX_KALMAN_H__ED3D740F_01D2_4616_8B74_8BF57636F2C0__IN CLUDED_) #define AFX_KALMAN_H__ED3D740F_01D2_4616_8B74_8BF57636F2C0__INCLU DED_#if _MSC_VER > 1000 #pragma once #endif // _MSC_VER > 1000#include <math.h> #include "cv.h"class kalman { public: void init_kalman(int x,int 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 xv=0,int y=0,int yv=0); //virtual ~kalman();};#endif // !defined(AFX_KALMAN_H__ED3D740F_01D2_4616_8B74_8BF57636F2C 0__INCLUDED_)============================kalman.cpp=============== =================#include "kalman.h" #include <stdio.h>/* tester de printer toutes les valeurs 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 yv) {cvkalman = cvCreateKalman( 4, 4, 0 ); state = cvCreateMat( 4, 1, CV_32FC1 ); process_noise = cvCreateMat( 4, 1, CV_32FC1 ); measurement = cvCreateMat( 4, 1, CV_32FC1 ); int code = -1;/* create matrix data */ const float A[] = { 1, T, 0, 0, 0, 1, 0, 0, 0, 0, 1, T, 0, 0, 0, 1 };const float H[] = { 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0 };const 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)/pow(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, 0, pow(T,2)/2, T };const float R[] = { 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0 }; cvRandInit( &rng, 0, 1, -1, CV_RAND_UNI );cvZero( measurement );cvRandSetRange( &rng, 0, 0.1, 0 ); rng.disttype = CV_RAND_NORMAL;cvRand( &rng, state );memcpy( cvkalman->transition_matrix->data.fl, A, sizeof(A)); memcpy( 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)); memcpy( cvkalman->measurement_noise_cov->data.fl, R, sizeof(R)); //cvSetIdentity( cvkalman->process_noise_cov, cvRealScalar(1e-5) ); //cvSetIdentity( cvkalman->error_cov_post, cvRealScalar(1)); //cvSetIdentity( cvkalman->measurement_noise_cov, cvRealScalar(1e-1) );/* choose initial state */state->data.fl[0]=x; state->data.fl[1]=xv; state->data.fl[2]=y;state->data.fl[3]=yv; cvkalman->state_post->data.fl[0]=x; cvkalman->state_post->data.fl[1]=xv; cvkalman->state_post->data.fl[2]=y; cvkalman->state_post->data.fl[3]=yv;cvRandSetRange( &rng, 0, sqrt(cvkalman->process_noise_cov->data.fl[0]), 0 ); cvRand( &rng, process_noise ); }CvPoint2D32f kalman::get_predict(float x, float y) {/* update state with current position */ state->data.fl[0]=x; state->data.fl[2]=y;/* predict point position */ /* x'k=A 鈥 k+B 鈥 kP'k=A 鈥 k-1*AT + Q */ cvRandSetRange( &rng, 0, sqrt(cvkalman->measurement_noise_cov->data.fl [0]), 0 ); cvRand( &rng, measurement );/* xk=A?xk-1+B?uk+wk */ cvMatMulAdd( cvkalman->transition_matrix, state, process_noise, cvkalman-> state_post );/* zk=H?xk+vk */cvMatMulAdd( cvkalman->measurement_matrix, cvkalman->state_post, meas urement, measurement );cvKalmanCorrect( cvkalman, measurement ); float measured_value_x = measurement->data.fl[0]; float measured_value_y = measurement->data.fl[2];const CvMat* prediction = cvKalmanPredict( cvkalman, 0 ); float predict_value_x = prediction->data.fl[0]; float predict_value_y = prediction->data.fl[2];return(cvPoint2D32f(predict_value_x,predict_value_y)); }void kalman::init_kalman(int x,int xv,int y,int yv) { state->data.fl[0]=x;state->data.fl[1]=xv; state->data.fl[2]=y; state->data.fl[3]=yv; cvkalman->state_post->data.fl[0]=x; cvkalman->state_post->data.fl[1]=xv; cvkalman->state_post->data.fl[2]=y; cvkalman->state_post->data.fl[3]=yv; }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,*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);}。
卡尔曼滤波C程序

#include <avr/io.h>//#include <avr/delay.h>#include <avr/pgmspace.h>#include <avr/interrupt.h>#include <util/twi.h>#include <math.h>#include <stdlib.h>#include"crc16.c"#include"uart.c"#include"kalman.c"#include"twimaster.c"unsigned char Array[8]={0x01,0x02,0x03,0x04,0x05,0x06,0x07,0x08}; unsigned char Rcvbuf[2];static int i;float y1,y2;int y3=0;//输出端口初始化void PORT_initial(void){DDRB=0B00110000;PINB=0X00;PORTB=0X00;DDRD=0B01100000;PIND=0X00;PORTD=0X00;}//Timer1初始化void Timer1Init(void){TCCR1A=0;TCCR1B=_BV(WGM12)|_BV(CS11);OCR1A=0x9C40; //20ms //0x4E20; //周期是:10毫秒TIMSK1|=_BV(OCIE1A);//Timer1 CTC中断,此中断发生周期是:10毫秒ISR(TIMER1_COMPA_vect){//gyro,acc为外部变量stateUpdate(gyro);//外部函数y1=yyh1+(gyro*0.02);y2=(int)kalmanUpdate(acc);//外部函数Array[0]=acc;Array[1]=acc>>8;Array[2]=(int)y1;Array[3]=(int)y1>>8;Array[4]=y2;Array[5]=y>>8;Array[6]=y3;//Array[7]=y3>>8;CRC16(Array, Rcvbuf,8);//外部函数for(i=0;i<8;i++) USART_Transmit(Array);//外部函数 USART_Transmit(Rcvbuf[1]);//外部函数USART_Transmit(Rcvbuf[0]);//外部函数}int main(void){PORT_initial();Timer1Init();I2C_Init();USART_initial();I2C_Start(); //启动twi的中断模式。
卡尔曼滤波算法C语言实现

0 o
0
2
2
呢, 简单的求平均已经不能满足精度的要求了。 我们可以用他们的协方差 covariance 来判断。 协方差本身就能体现两个信号的相关性, 通过它就能判断到底真值更逼近于预测值还是测量 值。引入 kalman gain( k g ) ,有公式计算 k g ,
2 kalman filter 最优化递归估计
Kalman filter 是一个“optimal recursive data processing algorithm(最优化递归数据处理 方法) ” 。 对于解决很大部分的问题, 他是最优, 效率最高甚至是最有用的方法。 而 kalman filter 最为核心的内容是体现它最优化估计和递归特点的 5 条公式。 举一个例子来详细说明 5 条公 式的物理意义。 假设我们要研究的对象是某一个房间的温度信号。 对于室温来说, 一分钟内或一小段时 间内的值是基本上不变的或者变化范围很小。也就是说 t1 时刻的温度 T1 和 t 2 时刻的温度 T2 基本不变,即 T2 T1 。在这个过程中,因为毕竟温度还是有所改变的,设有几度的偏差。 我们把这几度的偏差看成是高斯白噪声 w(t ) ,也就是说 E[ w(t )] 0 , D[ w(t )] 。除此
T
X (k | k ) X (k | k 1) k g (k ) (Z (k ) H X (k | k 1))
由上面分析可知为了实现递归,每次的 k g 都是实| K 1) H T /( H P(k | k 1) H T R)
4
图 1 真值为 0.56 的运行结果 给定值 z_real=3.32 时的运行结果如图 2
卡尔曼滤波代码c

卡尔曼滤波代码c卡尔曼滤波是一种用于估计系统状态的滤波算法。
以下是一个简单的卡尔曼滤波的 C 代码实现:```c#include <stdio.h>#define Q 0.1 // 过程噪声的方差#define R 1.0 // 测量噪声的方差void kalman_filter(float* x, float* P, float z, float* K) {// 预测步骤float x_pred = *x;float P_pred = *P + Q;// 更新步骤float y = z - x_pred;float S = P_pred + R;float K_gain = P_pred / S;*x = x_pred + K_gain * y;*P = (1 - K_gain) * P_pred;*K = K_gain;}int main() {float x = 0; // 初始状态float P = 1; // 初始协方差float K; // 卡尔曼增益float measurements[] = {1, 2, 3, 4, 5}; // 测量值for (int i = 0; i < 5; i++) {kalman_filter(&x, &P, measurements[i], &K);printf("Measurement: %.2f, Estimated: %.2f\n", measurements[i], x);}return 0;}```在这个例子中,我们使用卡尔曼滤波来估计一个连续的测量值序列。
在每次测量后,我们调用`kalman_filter` 函数来更新状态估计值`x` 和协方差 `P`。
然后,我们打印出当前测量值和估计值。
这只是一个简单的卡尔曼滤波实现,更复杂的应用可能需要更多的参数和计算。
如果你有特定的应用需求,可以根据卡尔曼滤波的原理进行相应的调整和扩展。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
卡尔曼滤波算法实现代码C++实现代码如下: ============================kalman.h================= =============== // kalman.h: interface for the kalman class. // ////////////////////////////////////////////////////////////////////// #if !defined(AFX_KALMAN_H__ED3D740F_01D2_4616_8B74_8BF57636F2C0__IN CLUDED_) #define AFX_KALMAN_H__ED3D740F_01D2_4616_8B74_8BF57636F2C0__INCLU DED_#if _MSC_VER > 1000 #pragma once #endif // _MSC_VER > 1000#include <math.h> #include "cv.h"class kalman { public: void init_kalman(int x,int 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 xv=0,int y=0,int yv=0); //virtual ~kalman();};#endif // !defined(AFX_KALMAN_H__ED3D740F_01D2_4616_8B74_8BF57636F2C 0__INCLUDED_)============================kalman.cpp=============== =================#include "kalman.h" #include <stdio.h>/* tester de printer toutes les valeurs 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 yv) {cvkalman = cvCreateKalman( 4, 4, 0 ); state = cvCreateMat( 4, 1, CV_32FC1 ); process_noise = cvCreateMat( 4, 1, CV_32FC1 ); measurement = cvCreateMat( 4, 1, CV_32FC1 ); int code = -1;/* create matrix data */ const float A[] = { 1, T, 0, 0, 0, 1, 0, 0, 0, 0, 1, T, 0, 0, 0, 1 };const float H[] = { 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0 };const 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)/pow(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, 0, pow(T,2)/2, T };const float R[] = { 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0 }; cvRandInit( &rng, 0, 1, -1, CV_RAND_UNI );cvZero( measurement );cvRandSetRange( &rng, 0, 0.1, 0 ); rng.disttype = CV_RAND_NORMAL;cvRand( &rng, state );memcpy( cvkalman->transition_matrix->data.fl, A, sizeof(A)); memcpy( 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)); memcpy( cvkalman->measurement_noise_cov->data.fl, R, sizeof(R)); //cvSetIdentity( cvkalman->process_noise_cov, cvRealScalar(1e-5) ); //cvSetIdentity( cvkalman->error_cov_post, cvRealScalar(1)); //cvSetIdentity( cvkalman->measurement_noise_cov, cvRealScalar(1e-1) );/* choose initial state */state->data.fl[0]=x; state->data.fl[1]=xv; state->data.fl[2]=y;state->data.fl[3]=yv; cvkalman->state_post->data.fl[0]=x; cvkalman->state_post->data.fl[1]=xv; cvkalman->state_post->data.fl[2]=y; cvkalman->state_post->data.fl[3]=yv;cvRandSetRange( &rng, 0, sqrt(cvkalman->process_noise_cov->data.fl[0]), 0 ); cvRand( &rng, process_noise ); }CvPoint2D32f kalman::get_predict(float x, float y) {/* update state with current position */ state->data.fl[0]=x; state->data.fl[2]=y;/* predict point position */ /* x'k=A 鈥 k+B 鈥 kP'k=A 鈥 k-1*AT + Q */ cvRandSetRange( &rng, 0, sqrt(cvkalman->measurement_noise_cov->data.fl [0]), 0 ); cvRand( &rng, measurement );/* xk=A?xk-1+B?uk+wk */ cvMatMulAdd( cvkalman->transition_matrix, state, process_noise, cvkalman-> state_post );/* zk=H?xk+vk */cvMatMulAdd( cvkalman->measurement_matrix, cvkalman->state_post, meas urement, measurement );cvKalmanCorrect( cvkalman, measurement ); float measured_value_x = measurement->data.fl[0]; float measured_value_y = measurement->data.fl[2];const CvMat* prediction = cvKalmanPredict( cvkalman, 0 ); float predict_value_x = prediction->data.fl[0]; float predict_value_y = prediction->data.fl[2];return(cvPoint2D32f(predict_value_x,predict_value_y)); }void kalman::init_kalman(int x,int xv,int y,int yv) { state->data.fl[0]=x;state->data.fl[1]=xv; state->data.fl[2]=y; state->data.fl[3]=yv; cvkalman->state_post->data.fl[0]=x; cvkalman->state_post->data.fl[1]=xv; cvkalman->state_post->data.fl[2]=y; cvkalman->state_post->data.fl[3]=yv; }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,*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);}。