找回密码
 立即注册

QQ登录

只需一步,快速开始

搜索
查看: 7023|回复: 5
收起左侧

STM32单片机卡尔曼滤波程序分享

  [复制链接]
ID:582568 发表于 2019-8-6 17:42 | 显示全部楼层 |阅读模式
在网上找的卡尔曼滤波在自己理解的基础上更改了一下代码!实测有用

单片机源程序如下:
  1. /* Includes ------------------------------------------------------------------*/
  2. #include "stm32f10x.h"       
  3. #include "kalman.h"       
  4. #include "matrix.h"                                                                                 

  5. /*==============================================================================
  6. 1.预估计
  7.    X(k|k-1) = F(k,k-1)*X(k-1|k-1)        //控制量为0


  8. 2.计算预估计协方差矩阵
  9.    P(k|k-1) = F(k,k-1)*P(k-1|k-1)*F(k,k-1)'+Q(k)
  10.    Q(k) = U(k)×U(k)'


  11. 3.计算卡尔曼增益矩阵
  12.    Kg(k) = P(k|k-1)*H' / (H*P(k|k-1)*H' + R(k))
  13.    R(k) = N(k)×N(k)'


  14. 4.更新估计
  15.    X(k|k) = X(k|k-1)+Kg(k)*(Z(k)-H*X(k|k-1))


  16. 5.计算更新后估计协防差矩阵
  17.    P(k|k) =(I-Kg(k)*H)*P(k|k-1)


  18. 6. 更新最优值


  19. F(k,k-1):     状态转移矩阵
  20. X(k|k-1):     根据k-1时刻的最优值估计k时刻的值
  21. X(k-1|k-1):   k-1时刻的最优值
  22. P(k|k-1):     X(k|k-1)对应的covariance
  23. P(k-1|k-1):   X(k-1|k-1)对应的covariance
  24. Q(k):         系统过程的covariance
  25. R(k):         测量过程的协方差
  26. H(k):         观测矩阵转移矩阵
  27. Z(k):         k时刻的测量值


  28. 基本思路: 首先根据上一次(如果是第一次则根据预赋值计算)的数据计算出本次的估计值,
  29.           同理,根据上一次的数据计算出本次估计值的协方差;  接着,由本次估计值的协
  30.           方差计算出卡尔曼增益;  最后,根据估测值和测量值计算当前最优值及其协方差
  31. ==============================================================================*/



  32. //================================================//
  33. //==             最优值方差结构体               ==//
  34. //================================================//
  35. typedef struct  _tCovariance
  36. {
  37.   float PNowOpt[LENGTH];
  38.   float PPreOpt[LENGTH];
  39. }tCovariance;







  40. static tOptimal      tOpt;
  41. static tCovariance   tCov;
  42. //float         Z[LENGTH]  = {4000};           //  测量值(每次测量的数据需要存入该数组)
  43. static float         I[LENGTH]  = {1};              //  单位矩阵
  44. static float         X[LENGTH]  = {9.8};              //  当前状态的预测值
  45. static float         P[LENGTH]  = {0};              //  当前状态的预测值的协方差
  46. static float         K[LENGTH]  = {0};              //  卡尔曼增益
  47. static float         Temp3[LENGTH] = {0};           //  辅助变量
  48. //============================================================================//
  49. //==                    卡尔曼滤波需要配置的变量                            ==//
  50. //============================================================================//
  51. static float         F[LENGTH]  = {1};              //  状态转移矩阵   即:坚信当前的状态与上一次状态的关系,如果坚信是一样的,则状态转移矩阵就为单位矩阵
  52. static float         Q[LENGTH]  = {0.0001f};//0.0001f              //  系统过程的协方差        协方差的定义:真实值与期望值之差的平方的期望值
  53. static float         R[LENGTH]  = {2};              //  测量过程的协方差        协方差的定义:真实值与期望值之差的平方的期望值   
  54. //如果你需要滤波结果更依赖于观测量,那就调小R,增大Q;反之,调大R,调小Q,这样估计值就取决于系统。
  55. //如果R大Q小,就是说,状态估计值比测量值要可靠,这时,所得出的结果就是更接近估计值;
  56. //如果R小Q大,这时,计算出来的结果就会更接近测量值。
  57. static float         H[LENGTH]  = {1};              //  观测矩阵转移矩阵        测量值与状态预测值之间的单位换算关系,即把预测值单位换算成测量值单位
  58. static float         Temp1[LENGTH] = {1};           //  辅助变量, 同时保存tOpt.XPreOpt[]的初始化值
  59. static float         Temp2[LENGTH] = {10000};       //  辅助变量, 同时保存tCov.PPreOpt[]的初始化值


  60. void KalMan_PramInit(void)
  61. {
  62.   unsigned char   i;
  63.   
  64.   for (i=0; i<LENGTH; i++)
  65.   {
  66.     tOpt.XPreOpt[i] = Temp1[i];           //零值初始化
  67.   }
  68.   for (i=0; i<LENGTH; i++)
  69.   {
  70.     tCov.PPreOpt[i] = Temp2[i];           //零值初始化
  71.   }
  72. }


  73. //============================================================================//
  74. //==                          卡尔曼滤波                                    ==//
  75. //============================================================================//
  76. //==入口参数: 当前时刻的测量值                                                            ==//
  77. //==出口参数: 当前时刻的最优值                                                            ==//
  78. //==返回值:   当前时刻的最优值                                                            ==//
  79. //============================================================================//
  80. float KalMan_Update(float *Z)
  81. {
  82.         u8 i;  
  83.         MatrixMul(F, tOpt.XPreOpt, X, ORDER, ORDER, ORDER);       //  基于系统的上一状态而预测现在状态; X(k|k-1) = F(k,k-1)*X(k-1|k-1)

  84.         MatrixCal(F, tCov.PPreOpt, Temp1, ORDER);
  85.         MatrixAdd(Temp1, Q, P, ORDER, ORDER);                     //  预测数据的协方差矩阵; P(k|k-1) = F(k,k-1)*P(k-1|k-1)*F(k,k-1)'+Q

  86.         MatrixCal(H, P, Temp1, ORDER);
  87.         MatrixAdd(Temp1, R, Temp1, ORDER, ORDER);
  88.         Gauss_Jordan(Temp1, ORDER);
  89.         MatrixTrans(H, Temp2, ORDER, ORDER);
  90.         MatrixMul(P, Temp2, Temp3, ORDER, ORDER, ORDER);
  91.         MatrixMul(Temp1, Temp3, K, ORDER, ORDER, ORDER);          //  计算卡尔曼增益; Kg(k) = P(k|k-1)*H' / (H*P(k|k-1)*H' + R)

  92.         MatrixMul(H, X, Temp1, ORDER, ORDER, ORDER);
  93.         MatrixMinus(Z, Temp1, Temp1, ORDER, ORDER);
  94.         MatrixMul(K, Temp1, Temp2, ORDER, ORDER, ORDER);
  95.         MatrixAdd(X, Temp2, tOpt.XNowOpt, ORDER, ORDER);          //  根据估测值和测量值计算当前最优值; X(k|k) = X(k|k-1)+Kg(k)*(Z(k)-H*X(k|k-1))

  96.         MatrixMul(K, H, Temp1, ORDER, ORDER, ORDER);
  97.         MatrixMinus(I, Temp1, Temp1, ORDER, ORDER);
  98.         MatrixMul(Temp1, P, tCov.PNowOpt, ORDER, ORDER, ORDER);   //  计算更新后估计协防差矩阵; P(k|k) =(I-Kg(k)*H)*P(k|k-1)

  99.         for (i=0; i<LENGTH; i++)
  100.         {
  101.           tOpt.XPreOpt[i] = tOpt.XNowOpt[i];
  102.           tCov.PPreOpt[i] = tCov.PNowOpt[i];
  103.         }
  104.        
  105.         return tOpt.XNowOpt[0];
  106. }
复制代码

所有资料51hei提供下载:
kalman.7z (3.71 KB, 下载次数: 252)
回复

使用道具 举报

ID:261849 发表于 2020-3-23 15:24 | 显示全部楼层
谢谢分享,比较容易理解.
回复

使用道具 举报

ID:401936 发表于 2020-4-2 20:14 | 显示全部楼层
首先感谢楼主分享!然后小白想问下卡尔曼滤波返回值不应该是下一时刻的最优值吗
回复

使用道具 举报

ID:920454 发表于 2021-5-13 09:52 | 显示全部楼层
tOptimal 这个结构体的定义在哪儿?
回复

使用道具 举报

ID:319585 发表于 2021-5-13 17:38 来自手机 | 显示全部楼层
这个不错,感觉论坛代码的质量越来越高了。
回复

使用道具 举报

ID:829464 发表于 2021-7-21 16:07 | 显示全部楼层
6啊 楼主
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

手机版|小黑屋|51黑电子论坛 |51黑电子论坛6群 QQ 管理员QQ:125739409;技术交流QQ群281945664

Powered by 单片机教程网

快速回复 返回顶部 返回列表