卡尔曼滤波算法(C--C++两种实现代码)
- 格式:doc
- 大小:176.37 KB
- 文档页数:10
无迹卡尔曼滤波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;}```。
自适应卡尔曼滤波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++分别实现)卡尔曼滤波器简介近来发现有些问题很多人都很感兴趣。
所以在这里希望能尽自己能力跟大家讨论一些力所能及的算法。
现在先讨论一下卡尔曼滤波器,如果时间和能力允许,我还希望能够写写其他的算法,例如遗传算法,傅立叶变换,数字滤波,神经网络,图像处理等等。
因为这里不能写复杂的数学公式,所以也只能形象的描述。
希望如果哪位是这方面的专家,欢迎讨论更正。
卡尔曼滤波器– 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)为了可以更加容易的理解卡尔曼滤波器,这里会应用形象的描述方法来讲解,而不是像大多数参考书那样罗列一大堆的数学公式和数学符号。
卡尔曼滤波算法实现代码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);}。
#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卡尔曼滤波是一种用于估计系统状态的滤波算法。
以下是一个简单的卡尔曼滤波的 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`。
然后,我们打印出当前测量值和估计值。
这只是一个简单的卡尔曼滤波实现,更复杂的应用可能需要更多的参数和计算。
如果你有特定的应用需求,可以根据卡尔曼滤波的原理进行相应的调整和扩展。
卡尔曼滤波原理卡尔曼滤波最早可以追溯到Wiener滤波,不同的是卡尔曼采用状态空间来描述它的滤波器,卡尔曼滤波器同时具有模糊/平滑与预测功能,特别是后者在视频分析与对象跟踪应用场景中被发扬光大,在离散空间(图像或者视频帧)使用卡尔曼滤波器相对简单。
假设我们根据一个处理想知道一个变量值如下:最终卡尔曼滤波完整的评估与空间预测模型工作流程如下:OpenCV APIcv::KalmanFilter::KalmanFilter(int dynamParams,int measureParams,int controlParams = 0,int type = CV_32F)# dynamParams表示state的维度# measureParams表示测量维度# controlParams表示控制向量# type表示创建的matrices代码演示import cv2from math import cos, sin, sqrtimport numpy as npif __name__ == "__main__":img_height = 500img_width = 500kalman = cv2.KalmanFilter(2, 1, 0)dWindow("Kalman", cv2.WINDOW_AUTOSIZE)while True:state = 0.1 * np.random.randn(2, 1)# 初始化kalman.transitionMatrix = np.array([[1., 1.], [0., 1.]])kalman.measurementMatrix = 1. * np.ones((1, 2))kalman.processNoiseCov = 1e-5 * np.eye(2)kalman.measurementNoiseCov = 1e-1 * np.ones((1, 1))kalman.errorCovPost = 1. * np.ones((2, 2))kalman.statePost = 0.1 * np.random.randn(2, 1)while True:def calc_point(angle):return (np.around(img_width/2 + img_width/3*cos(angle), 0).astype(int), np.around(img_height/2 - img_width/3*sin(angle), 1).astype(int)) state_angle = state[0, 0]state_pt = calc_point(state_angle)# 预测prediction = kalman.predict()predict_angle = prediction[0, 0]predict_pt = calc_point(predict_angle)measurement = kalman.measurementNoiseCov * np.random.randn(1, 1)# 生成测量measurement = np.dot(kalman.measurementMatrix, state) + measurementmeasurement_angle = measurement[0, 0]measurement_pt = calc_point(measurement_angle)# plot pointsdef draw_cross(center, color, d):cv2.line(img,(center[0] - d, center[1] - d), (center[0] + d, center[1] + d),color, 1, cv2.LINE_AA, 0)cv2.line(img,(center[0] + d, center[1] - d), (center[0] - d, center[1] + d),color, 1, cv2.LINE_AA, 0)img = np.zeros((img_height, img_width, 3), np.uint8)cv2.line(img, state_pt, measurement_pt, (0, 0, 255), 3, cv2.LINE_AA, 0)cv2.line(img, state_pt, predict_pt, (255, 0, 0), 3, cv2.LINE_AA, 0)# 校正预测与测量值差异kalman.correct(measurement)# 更新noise矩阵与状态process_noise = sqrt(kalman.processNoiseCov[0,0]) * np.random.randn(2, 1) state = np.dot(kalman.transitionMatrix, state) + process_noisecv2.imshow("Kalman", img)code = cv2.waitKey(100)if code != -1:breakif code in [27, ord('q'), ord('Q')]:breakcv2.destroyWindow("Kalman")。
卡尔曼滤波算法实现代码C++实现代码如下:============================kalma n.h=================// kalman.h: interface for the kalman class./////////////////////////////////////////////////////////////////////// /#if !defined(AFX_KALMAN_H__ED3D740F_01D2_4616_8B74_8BF57636F2C0_IN CLUDED」#defi ne AFX_KALMAN_H__ED3D740F_01D2_4616_8B74_8BF57636F2C0__INCLU DED_#if _MSC_VER > 1000#pragma once#en dif // _MSC_VER > 1000#in clude <math.h>#i nclude "cv.h" class kalma n{public :void init_kalman( int x, int xv, int y, int yv);CvKalma n* cvkalma n;CvMat* state;CvMat* process_ no ise;CvMat* measureme nt;const CvMat* prediction;CvPoint2D32f get_predict( float x, float y);kalman( int x=0, int xv=O, int y=0, int yv=O);//virtual ~kalma n();};#en dif // !defi ned(AFX_KALMAN_H_ED3D740F_01D2_4616_8B74_8BF57636F2CO__INCLUDED_)============================kalma n. cpp===============#i nclude "kalma n.h"#i nclude <stdio.h>/* tester de prin ter toutes les valeurs des vecteurs ■■■ *//* tester de cha nger les matrices du no ises *//* replace state by cvkalma n->state_post ??? */CvRa ndState rng;const double T = 0.1;kalman::kalman( int x, int xv, int y, int yv){cvkalma n = cvCreateKalma n( 4, 4, 0 );state = cvCreateMat( 4, 1, CV_32FC1 );process_ noise = cvCreateMat( 4, 1, CV_32FC1 ); measureme nt = cvCreateMat( 4, 1, CV_32FC1 );int code = -1;/* create matrix data *I con st float A[] = {1, T , 0, 0,0,1, 0, 0,0, 0,1, T ,0, 0, 0, 1};con st float H[] = {1, 0, 0, 0,con st 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};con st float R[] = {1, 0, 0, 0,0, 0, 0, 0,0, 0,1,0,0, 0, 0, 0};cvRa ndl nit( &rng, 0, 1, -1, CV_RAND_UNI );cvZero( measureme nt );cvRa ndSetRa nge( &rng, 0, 0.1, 0 );rng.disttype = CV_RAND_NORMAL;cvRa nd( &rng, state );memcpy( cvkalman->transition_matrix->data.fl, A, sizeof (A)); memcpy( cvkalman->measurement_matrix->data.fl, H, sizeof (H)); memcpy( cvkalma n->process_ no ise_cov->data.fl,sizeof (Q));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));//cvSetIde ntity( cvkalma n-> measureme nt_no ise_cov,cvRealScalar(1e-1));/* choose in itial state */state->data.fl[0]=x;state->data.fl[1]=xv;state->data.fl[2]=y;state->data.fl[3]=yv;cvkalma n->state_post->data.fl[0]=x;cvkalma n->state_post->data.fl[1]=xv;cvkalma n->state_post->data.fl[2]=y;cvkalma n->state_post->data.fl[3]=yv;cvRa ndSetRa nge( &rng, 0, sqrt(cvkalma n->process_ no ise_cov->data.fl[0]), 0 );cvRa nd( &rng, process_ no ise );}CvPoint2D32f kalman::get_predict( float x, float y){/* update state with curre nt positi on */state->data.fl[O]=x;state->data.fl[2]=y;/* predict point positi on *//* x'k=A 欽k+B 欽kP'k=A 欽k-1*AT + Q */cvRandSetRange( &rng, 0, sqrt(cvkalman->measurement_noise_cov->data.fl [0]), 0 );cvRa nd( &rng, measureme nt );/* xk=A?xk-1+B?uk+wk */cvMatMulAdd( cvkalma n->tra nsitio n_ matrix, state, process_ no ise, cvkalma n-> state_post );/* zk=H?xk+vk */cvMatMulAdd( cvkalma n->measureme nt_matrix, cvkalma n->state_post, meas ureme nt, measureme nt );cvKalma nCorrect( cvkalma n, measureme nt );float measured_value_x = measurement->data.fl[O];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;cvkalma n->state_post->data.fl[0]=x;cvkalma n->state_post->data.fl[1]=xv;cvkalma n->state_post->data.fl[2]=y;cvkalma n->state_post->data.fl[3]=yv;}c 语言实现代码如下for (j=0; j<=m-1; j++)for (ii=2; ii<=k; ii++) for (ii=2; ii<=k; ii++){ for (i=0; i<=n-1; i++){ jj=i*+ a[jj]=0.0;for (i=0; i<=n-1; i++) for (j=0;j<=m-1; j++){ jj=i*m+j; g[jj]=O.O;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;{ free(e); free(a); free(b); return (js);} for (kk=0; kk<=m-1; kk++)for (j=0; j<=n-1; j++){ jj=i*l; b[jj]=y[(ii-1)*m+i];for (j=0; j<=m-1; j++)I x[jj]=x[jj]+g [i *m+j]*b『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];I if (i==j) a[jj]=1.0+a[jj];■ }I for (i=0; i<=n-1; i++)for (i=0; i<=n-1; i++)for (j=0; jv=n-1; j++){ jj=i*l+j; a[jj]=0.0;for (kk=0; kk<=n-1; kk++)}}}free(e); free(a); free(b);I return (js);}。
卡尔曼滤波简介及其算法实现代码(C++/C/MATLAB)卡尔曼滤波器简介近来发现有些问题很多人都很感兴趣。
所以在这里希望能尽自己能力跟大家讨论一些力所能及的算法。
现在先讨论一下卡尔曼滤波器,如果时间和能力允许,我还希望能够写写其他的算法,例如遗传算法,傅立叶变换,数字滤波,神经网络,图像处理等等。
因为这里不能写复杂的数学公式,所以也只能形象的描述。
希望如果哪位是这方面的专家,欢迎讨论更正。
卡尔曼滤波器– 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)为了可以更加容易的理解卡尔曼滤波器,这里会应用形象的描述方法来讲解,而不是像大多数参考书那样罗列一大堆的数学公式和数学符号。
卡尔曼滤波算法实现代码C++实现代码如下:============================kalman.h================================// kalman.h: interface for the kalman class.////////////////////////////////////////////////////////////////////////#if!defined(AFX_KALMAN_H__ED3D740F_01D2_4616_8B74_8BF57636F2C0__INCLUDED_) #define 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,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_8BF57636F2C0__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, cvkalma n->state_post );/* zk=H?xk+vk */cvMatMulAdd( cvkalman->measurement_matrix, cvkalman->state_post, measure ment, 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);}。