找回密码
 注册
查看: 721|回复: 1
打印 上一主题 下一主题

基于9轴惯性运动传感器的三阶卡尔曼滤波器算法

[复制链接]

该用户从未签到

跳转到指定楼层
1#
发表于 2019-2-21 12:30 | 只看该作者 回帖奖励 |倒序浏览 |阅读模式

EDA365欢迎您登录!

您需要 登录 才可以下载或查看,没有帐号?注册

x
基于9轴惯性运动传感器的三阶卡尔曼滤波器算法
* @: D+ ^$ ~* x, a/ ]+ l- t
! R3 k* W( D; {2 f5 S2 h7 a
最近在玩九轴的惯性传感器,很是有挑战性.九轴说的是三轴的加速度计、三轴的陀螺仪以及三轴的磁场传感器。但是只是单纯的测出九个轴的数据没什么用,关键是要能够融合这九轴数据得出我们想要的结果。这里就运用三阶卡尔曼滤波算法来融合这九轴运动数据为三轴的角度。运用这三个角度可以用来做自平衡车或者四轴飞行器.
一、卡尔曼算法理解
其实如果不去考虑kalman算法是怎么来的,我们只需要知道有下面几个式子就可以了,具体意思可以看上面的wikipedia链接
二 卡尔曼滤波算法的实现
这里我的算法是运行在avr单片机上的,所以采用的是c语言写的。下面的代码是要放到avr的定时器中断测试刷新的。用示波器测试了一下,这个算法在16M晶振下的运行时间需要0.35ms,而数据采集需要3ms左右,所以选定定时器时间为8ms.之前也写过一阶的kalman算法,运用在自平衡车上,这边是三阶的,主要是矩阵运算.
//kalman.c

$ C4 R8 X8 \$ K, V
float dtTimer   = 0.008;

" H$ D/ M3 C0 ~4 y
float xk[9] = {0,0,0,0,0,0,0,0,0};

3 s9 x* K/ _1 e& E4 ?! K, ]
float pk[9] = {1,0,0,0,1,0,0,0,0};
0 I+ C) P- c( t
float I[9]  = {1,0,0,0,1,0,0,0,1};
; J4 v1 E) x! T  x$ ]
float R[9]  = {0.5,0,0,0,0.5,0,0,0,0.01};
# w8 j" U" n+ s' k" D; J
float Q[9] = {0.005,0,0,0,0.005,0,0,0,0.001};
  H/ }+ f  z1 s) n; v
9 }2 z! S4 }( ?! Y/ }
void matrix_add(float* mata,float* matb,float* matc){

  b5 _" }. {5 w+ x; o9 o/ W: G( }
    uint8_t i,j;

! p* E  E5 V' `- j$ B5 M
    for (i=0; i<3; i++){

9 d  ]5 D9 A7 e
       for (j=0; j<3; j++){
5 G- I1 b, L" C8 H# b3 s8 L
          matc[i*3+j] = mata[i*3+j] + matb[i*3+j];
( b9 Z+ R# L+ P" R% [" ~+ ^7 {% C+ V9 H
       }
7 y& t$ r6 {& K: ]( G
    }

7 b7 j; C2 l- }- U) |
}

) ?4 Z! O& P% N8 l0 g! t# N& n5 g( N# a' ?$ R' |/ ~: Y4 `
void matrix_sub(float* mata,float* matb,float* matc){
) c& o6 t9 ~3 v; _8 Z+ L% ?
    uint8_t i,j;

" Z- o; h6 _8 z
    for (i=0; i<3; i++){
7 N$ g, Y( u3 i1 O; }
       for (j=0; j<3; j++){

4 A* ?9 x, ^$ {9 c: e5 r% ^
          matc[i*3+j] = mata[i*3+j] - matb[i*3+j];

. O! r# i5 o7 F# D( k$ B6 o& U- y
       }
$ U' r7 k! o# U! x( D2 r
    }
0 l; }8 Y7 b( Z5 ?8 v7 C
}
- ?+ |0 w. p6 f1 @6 I1 D

, o- `1 P: O  m5 D: x
void matrix_multi(float* mata,float* matb,float* matc){
* r# a1 p( C$ X* F
  uint8_t i,j,m;
4 P, |7 X* O0 a2 G! F* ~
  for (i=0; i<3; i++)

6 s# Q6 z9 b7 @5 E
  {
/ Q- y0 S) ~. h( ]1 I8 p
    for (j=0; j<3; j++)

. n. K2 H# A5 T7 b
   {
& S  `# R" x1 F' U( o1 C) J( J0 T
      matc[i*3+j]=0.0;

* E. B! F4 _4 ]0 Y; U, m
      for (m=0; m<3; m++)
( `1 I& o" E) g- q
     {

. x$ C* E) Q: L7 D
        matc[i*3+j] += mata[i*3+m] *  matb[m*3+j];
) d# |% I+ ]" Y% w8 O$ p. S' h
      }
/ |+ B$ {* ~8 O2 I. g6 A  {
    }
7 T5 F( w/ I& ?- Y0 A
}
% L6 E& u8 P* K* {
}
- i$ ]2 o* A/ L. K/ M
% ?6 L4 R& d# s# R4 a5 X
void KalmanFilter(float* am_angle_mat,float* gyro_angle_mat){
7 E: r* Y! h  K5 ?1 q1 j0 `$ }" a
uint8_t i,j;

. M0 ^' ]0 v$ y) @# P
float yk[9];
8 s$ Y9 r& H4 V5 t
float pk_new[9];

! ]( p; S0 d: h
float K[9];

1 H0 E6 z  }$ h% D3 G0 j& i7 M
float KxYk[9];

& L9 o6 v" i$ c7 }% y: f# Q, @
float I_K[9];

2 `( {- A9 a+ a5 s
float S[9];
- G1 T5 }# L1 Y$ ^+ ?
float S_invert[9];

3 H5 u8 t7 S- L( h
float sdet;
4 j  g' T3 k! j" T% N

( F1 S* S* ?. ?
//xk = xk + uk
- v$ m- T+ ~, D# l+ d; a
matrix_add(xk,gyro_angle_mat,xk);

+ Q/ X( \- `" n
//pk = pk + Q
5 K6 c; S& A7 y/ \, D
matrix_add(pk,Q,pk);

. w: ^! {4 R- c9 N& h# Z3 W: Q
//yk =  xnew - xk
* e, ?" b1 h: R9 b- d/ ?% x
matrix_sub(am_angle_mat,xk,yk);
! }  ~: j) l) o' `) ]& W3 o1 G
//S=Pk + R
) Z  Y3 r( B" L, `6 B) p9 q$ N
matrix_add(pk,R,S);
9 |( x& y/ T/ M! H
//S求逆invert

. t" @9 b3 S) t4 [  r) p/ p( i
sdet = S[0] * S[4] * S[8]
! L: ?8 e2 ?2 M# p* Z
          + S[1] * S[5] * S[6]
8 _  x! |  s, C
          + S[2] * S[3] * S[7]

4 P4 Q5 z1 t. H0 H& k  x! p
          - S[2] * S[4] * S[6]
5 B2 {& U( Z* u( O
          - S[5] * S[7] * S[0]

- P& _; a; M8 V  B
          - S[8] * S[1] * S[3];
8 ]9 I. l) S/ f# Q5 Z7 g1 a% n5 U* Z

& S0 K# v! N# R1 J! e6 j
S_invert[0] = (S[4] * S[8] - S[5] *  S[7])/sdet;
1 e5 F. g- T6 ?: F) A
S_invert[1] = (S[2] * S[7] - S[1] *  S[8])/sdet;

& P4 |0 W' C+ M# \% _; r6 J
S_invert[2] = (S[1] * S[7] - S[4] *  S[6])/sdet;

0 i% E; @9 r, l' R' A, d! v$ p0 U2 f
S_invert[3] = (S[5] * S[6] - S[3] *  S[8])/sdet;
7 L# O" p4 X# M, o+ x
S_invert[4] = (S[0] * S[8] - S[2] *  S[6])/sdet;

  B5 y8 L3 B( l* A
S_invert[5] = (S[2] * S[3] - S[0] *  S[5])/sdet;

2 ?' k, V1 q. G0 S: S7 a5 k3 M* s* X; B. D/ v
S_invert[6] = (S[3] * S[7] - S[4] *  S[6])/sdet;
" y  A% j9 Z. ~1 i2 K
S_invert[7] = (S[1] * S[6] - S[0] *  S[7])/sdet;
1 T3 u& h4 q; z1 C& B( Q% }
S_invert[8] = (S[0] * S[4] - S[1] *  S[3])/sdet;

3 V/ w! d2 w) C+ q* o! y. V: p/ @
//K = Pk * S_invert
2 r2 K/ q0 j8 T( H+ \9 |
matrix_multi(pk,S_invert,K);
- p& S$ a/ ]. A* E! T
//KxYk = K * Yk

; ?7 V9 s2 @  i& w3 ?- W3 P$ s
matrix_multi(K,yk,KxYk);

5 f& [% T2 t7 D; m+ Z1 p
//xk = xk + K * yk

. c4 W, P1 |- v1 I1 G4 w
matrix_add(xk,KxYk,xk);
$ _$ N; J+ K2 e
//pk = (I - K)*(pk)

7 f3 R6 H) |. y+ F- h! ^
matrix_sub(I,K,I_K);
# Q# E2 X$ C, D$ M
matrix_multi(I_K,pk,pk_new);

0 ?% b8 x* V+ t% m9 r
//update pk
$ o$ d  ]7 L* J- `8 f% q. Z; Q+ J) p! t
//pk = pk_new;

+ x5 A! L, X' k1 m) _, l
for (i=0; i<3; i++){
5 E  G) O- G6 T& G/ e# r5 N% i
    for (j=0; j<3; j++){
# f6 X2 t$ P: z
        pk[i*3+j] = pk_new[i*3+j];
/ W" k. y. `9 {" V  O
    }
1 l2 }* w0 Z$ M
  }
: Y: z% [9 B* B, x6 r! O7 O1 K
}

" [. s/ j2 Y; X/ G% l
8 O; ~3 P7 {7 S" P$ o9 D
三 运用卡尔曼滤波器
这里的kalman滤波器是离散数字滤波器,需要迭代运算。这里把算法放到avr的定时器中断里面执行,进行递归运算.
//isr.c
. x2 p6 [& B4 [; ~; ~. s8 i2 |, M
#include "kalman.h"
- u' P3 V0 f2 C$ U0 E% W; \
float mpu_9dof_values[9]={0};

& t8 l' C' B, Q% ^
float am_angle[3];

& `$ V! h: Z3 l) H1 v1 V1 j
float gyro_angle[3];
) Y7 s) e/ y  y7 `3 @$ p
float am_angle_mat[9]={0,0,0,0,0,0,0,0,0};

1 }9 p6 [" C- U& Z7 j: X% M' d
float gyro_angle_mat[9]={0,0,0,0,0,0,0,0,0};
8 T4 p0 G$ B- C! Y7 h) e! {% M* L2 r
$ r- _( c" J$ J3 J  J! k' h
//8MS

) ]6 m2 o; n3 r5 V; e
ISR(TIMER0_OVF_vect){
) Y# D' T8 t/ n" ~0 R1 c0 z) `
//设置计数开始的初始值

) q! X; M/ A3 G: o* I
TCNT0 = 130 ;  //8ms
, q7 d7 [' _4 `. h: ~/ ?
sei();

4 ~4 Y3 v* }- t7 p
//采集九轴数据
$ v  k! A' q. V# V3 e' m
//mpu_9dof_values 单位为g与度/s

3 J: t/ W4 e" @% i( a7 `4 o
//加速度计和陀螺仪
3 _' o7 z  n8 E, T! n- C9 E3 C. c
mpu_getValue6(&mpu_9dof_values[0],&mpu_9dof_values[1],&mpu_9dof_values[2],&mpu_9dof_values[3],&mpu_hmc_values[4],&mpu_hmc_values[5]);

4 f# U8 B/ g1 c! D
//磁场传感器
6 ?9 q0 X& L- z/ [8 F. i
compass_mgetValues(&mpu_9dof_values[6],&mpu_9dof_values[7],&mpu_9dof_values[8]);

- t) ~* c) T1 s& p$ @' f' F* S
' a! i3 f9 l" v: c7 e+ [" _
accel_compass2angle(&mpu_9dof_values[0],&mpu_9dof_values[6],am_angle);
$ Z0 S% w8 I2 I) v" M+ Q
gyro2angle(&mpu_9dof_values[3],gyro_angle);
. x  \$ m( t6 r$ e; K4 u
  S! k. J6 D- @1 R7 l& i
am_angle_mat[0] = am_angle[0];
, b& n7 u. [0 W: W6 Z$ L% S
am_angle_mat[4] = am_angle[1];

8 R& K1 W3 K( I7 h2 y! ~- t, X
am_angle_mat[8] = am_angle[2];

5 q. f. Q7 C4 |% ~3 r
0 `$ W* O$ Q5 Q" D; f) V. N
gyro_angle_mat[0] = gyro_angle[1];
( S9 x: o  Z3 `. Y$ T
gyro_angle_mat[4] = - gyro_angle[0];
' u# C* Y) p5 i8 E" K+ ~% F. q
gyro_angle_mat[8] = - gyro_angle[2];

+ C- V' W7 p0 _) F5 s
4 \! L8 k$ o# G: o' O. y: J+ e3 D
//卡尔曼
( B, ^; w% \6 Z9 }+ N( a
KalmanFilter(am_angle_mat,gyro_angle_mat);
; A- _, |1 _" N; d6 q7 G1 H
3 c/ u; S1 c! Y! ]8 C1 \  h" {# Y+ _
//输出三轴角度

/ r9 b% B, v3 N2 V& y* `( m1 o7 [8 ]
//xk[0] xk[4] xk[8]

/ J  F) L& r' f* x1 T. F; R# J0 b
}

) e: d/ B( p; X8 q- h4 [! ?
6 Z' T8 a% L/ d% P7 m, q2 W( Y
实测可以准确的输出三轴的角度,为了获得更好的响应速度和跟踪精度还需调整参数.

' o! M" w. g" c+ `4 J) a. ^

该用户从未签到

2#
发表于 2019-2-21 16:18 | 只看该作者
研究一下,谢谢分享
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

推荐内容上一条 /1 下一条

EDA365公众号

关于我们|手机版|EDA365电子论坛网 ( 粤ICP备18020198号-1 )

GMT+8, 2025-5-29 06:33 , Processed in 0.078125 second(s), 23 queries , Gzip On.

深圳市墨知创新科技有限公司

地址:深圳市南山区科技生态园2栋A座805 电话:19926409050

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