找到一段卡尔曼滤波的代码 求大神帮看一下

[复制链接]
 楼主| doodle 发表于 2015-5-4 09:46 | 显示全部楼层 |阅读模式
  1. /*
  2. * KalmanFilter.h
  3. * Non-original
  4. * Author: x2d
  5. * Copyright (c) 2012 China
  6. *
  7. */

  8. #ifndef KalmanFilter_h
  9. #define KalmanFilter_h

  10. #include <WProgram.h>

  11. class KalmanFilter
  12. {
  13.   public:
  14.   KalmanFilter();
  15.   
  16.   /*
  17.     卡尔曼融合计算
  18.     angle_m: 加速度计测量并通过atan2(ax,ay)方法计算得到的角度(弧度值)
  19.     gyro_m:陀螺仪测量的角速度值(弧度值)
  20.     dt:采样时间(s)
  21.     outAngle:卡尔曼融合计算出的角度(弧度值)
  22.     outAngleDot:卡尔曼融合计算出的角速度(弧度值)
  23.   */
  24.   void getValue(double angle_m, double gyro_m, double dt, double &outAngle, double &outAngleDot);
  25.   
  26.   private:
  27.   double C_0, Q_angle, Q_gyro, R_angle;
  28.   double q_bias, angle_err, PCt_0, PCt_1, E, K_0, K_1, t_0, t_1;
  29.   double angle, angle_dot;
  30.   double P[2][2];
  31.   double Pdot[4];
  32. };

  33. /*CPP文件  */
  34. /*
  35. * KalmanFilter.cpp
  36. * Non-original
  37. * Author: x2d
  38. * Copyright (c) 2012 China
  39. *
  40. */

  41. #include "KalmanFilter.h"

  42. KalmanFilter::KalmanFilter()
  43. {
  44.    C_0 = 1.0f;
  45.    Q_angle = 0.001f;//加速度计的过程噪声协方差
  46.    Q_gyro = 0.003f;//陀螺仪的过程噪声协方差
  47.    R_angle = 0.5f;//加速度计的测量噪声协方差
  48.    q_bias = angle_err = PCt_0 = PCt_1 = E = K_0 = K_1 = t_0 = t_1 = 0.0f;
  49.    angle = angle_dot = 0.0f;
  50.    P[0][0] = 1.0f;
  51.    P[0][1] = 0.0f;
  52.    P[1][0] = 0.0f;
  53.    P[1][1] = 1.0f;
  54.    Pdot[0] = 0.0f;
  55.    Pdot[1] = 0.0f;
  56.    Pdot[2] = 0.0f;
  57.    Pdot[3] = 0.0f;
  58. }

  59. void  KalmanFilter::getValue(double angle_m, double gyro_m, double dt, double &outAngle, double &outAngleDot)
  60. {
  61. /*
  62.   Serial.print("angle_m = ");
  63.   Serial.print(angle_m);
  64.   Serial.print(";");
  65.   Serial.print("gyro_m = ");
  66.   Serial.print(gyro_m);
  67.   Serial.print(";");
  68. */
  69.   //陀螺仪积分的先验角度
  70.   angle+=(gyro_m-q_bias) * dt;

  71.   //角度偏差(用于计算后验状态估计)
  72.   angle_err = angle_m - angle;

  73.   //先验误差协方差的微分
  74.   Pdot[0] = Q_angle - P[0][1] - P[1][0];
  75.   Pdot[1] = -P[1][1];
  76.   Pdot[2] = -P[1][1];
  77.   Pdot[3] = Q_gyro;

  78.   //先验误差协方差的积分
  79.   P[0][0] += Pdot[0] * dt;//计算加速度计过程噪声协方差
  80.   P[0][1] += Pdot[1] * dt;
  81.   P[1][0] += Pdot[2] * dt;
  82.   P[1][1] += Pdot[3] * dt;//计算陀螺仪过程噪声协方差

  83.   //卡尔曼增益计算
  84.   PCt_0 = C_0 * P[0][0];
  85.   PCt_1 = C_0 * P[1][0];
  86.   E = R_angle + C_0 * PCt_0;//测量余量协方差
  87.   K_0 = PCt_0 / E;//加速度计的卡尔曼增益
  88.   K_1 = PCt_1 / E;//陀螺仪偏差的卡尔曼增益

  89.   //后验误差协方差计算
  90.   t_0 = PCt_0;
  91.   t_1 = C_0 * P[0][1];
  92.   P[0][0] -= K_0 * t_0;
  93.   P[0][1] -= K_0 * t_1;
  94.   P[1][0] -= K_1 * t_0;
  95.   P[1][1] -= K_1 * t_1;

  96.   //更新状态
  97.   angle += K_0 * angle_err;//后验估计最优角度值
  98.   q_bias += K_1 * angle_err;//更新最优估计值偏差
  99.   angle_dot = gyro_m-q_bias;//更新最优角速度(后验估计)
  100.   
  101.   outAngle = angle;
  102.   outAngleDot = angle_dot;
  103. /*
  104.   Serial.print("angle = ");
  105.   Serial.print(angle);
  106.   Serial.print(";");
  107.   Serial.print("angle_dot = ");
  108.   Serial.print(angle_dot);
  109.   Serial.print(";");
  110. */
  111. }
 楼主| doodle 发表于 2015-5-4 09:48 | 显示全部楼层
看了很久 其中的对数组的运算一直无法理解是怎么回事。。。。。求大神解释下
han0097 发表于 2015-5-7 10:53 | 显示全部楼层
什么东西?
您需要登录后才可以回帖 登录 | 注册

本版积分规则

1

主题

16

帖子

0

粉丝
快速回复 在线客服 返回列表 返回顶部

1

主题

16

帖子

0

粉丝
快速回复 在线客服 返回列表 返回顶部