PX4 Firmware
PX4 Autopilot Software http://px4.io
mag_fusion.cpp
Go to the documentation of this file.
1 /****************************************************************************
2  *
3  * Copyright (c) 2015 Estimation and Control Library (ECL). All rights reserved.
4  *
5  * Redistribution and use in source and binary forms, with or without
6  * modification, are permitted provided that the following conditions
7  * are met:
8  *
9  * 1. Redistributions of source code must retain the above copyright
10  * notice, this list of conditions and the following disclaimer.
11  * 2. Redistributions in binary form must reproduce the above copyright
12  * notice, this list of conditions and the following disclaimer in
13  * the documentation and/or other materials provided with the
14  * distribution.
15  * 3. Neither the name ECL nor the names of its contributors may be
16  * used to endorse or promote products derived from this software
17  * without specific prior written permission.
18  *
19  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
20  * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
21  * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
22  * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
23  * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
24  * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
25  * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
26  * OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
27  * AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
28  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
29  * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
30  * POSSIBILITY OF SUCH DAMAGE.
31  *
32  ****************************************************************************/
33 
34 /**
35  * @file heading_fusion.cpp
36  * Magnetometer fusion methods.
37  *
38  * @author Roman Bast <bapstroman@gmail.com>
39  * @author Paul Riseborough <p_riseborough@live.com.au>
40  *
41  */
42 
43 #include "ekf.h"
44 #include <ecl.h>
45 #include <mathlib/mathlib.h>
46 
48 {
49  // assign intermediate variables
50  float q0 = _state.quat_nominal(0);
51  float q1 = _state.quat_nominal(1);
52  float q2 = _state.quat_nominal(2);
53  float q3 = _state.quat_nominal(3);
54 
55  float magN = _state.mag_I(0);
56  float magE = _state.mag_I(1);
57  float magD = _state.mag_I(2);
58 
59  // XYZ Measurement uncertainty. Need to consider timing errors for fast rotations
60  float R_MAG = fmaxf(_params.mag_noise, 0.0f);
61  R_MAG = R_MAG * R_MAG;
62 
63  // intermediate variables from algebraic optimisation
64  float SH_MAG[9];
65  SH_MAG[0] = 2.0f*magD*q3 + 2.0f*magE*q2 + 2.0f*magN*q1;
66  SH_MAG[1] = 2.0f*magD*q0 - 2.0f*magE*q1 + 2.0f*magN*q2;
67  SH_MAG[2] = 2.0f*magD*q1 + 2.0f*magE*q0 - 2.0f*magN*q3;
68  SH_MAG[3] = sq(q3);
69  SH_MAG[4] = sq(q2);
70  SH_MAG[5] = sq(q1);
71  SH_MAG[6] = sq(q0);
72  SH_MAG[7] = 2.0f*magN*q0;
73  SH_MAG[8] = 2.0f*magE*q3;
74 
75  // rotate magnetometer earth field state into body frame
77 
78  Vector3f mag_I_rot = R_to_body * _state.mag_I;
79 
80  // compute magnetometer innovations
81  _mag_innov[0] = (mag_I_rot(0) + _state.mag_B(0)) - _mag_sample_delayed.mag(0);
82  _mag_innov[1] = (mag_I_rot(1) + _state.mag_B(1)) - _mag_sample_delayed.mag(1);
83 
84  // do not use the synthesized measurement for the magnetomter Z component for 3D fusion
86  _mag_innov[2] = 0.0f;
87  } else {
88  _mag_innov[2] = (mag_I_rot(2) + _state.mag_B(2)) - _mag_sample_delayed.mag(2);
89  }
90  // Observation jacobian and Kalman gain vectors
91  float H_MAG[24];
92  float Kfusion[24];
93 
94  // X axis innovation variance
95  _mag_innov_var[0] = (P[19][19] + R_MAG + P[1][19]*SH_MAG[0] - P[2][19]*SH_MAG[1] + P[3][19]*SH_MAG[2] - P[16][19]*(SH_MAG[3] + SH_MAG[4] - SH_MAG[5] - SH_MAG[6]) + (2.0f*q0*q3 + 2.0f*q1*q2)*(P[19][17] + P[1][17]*SH_MAG[0] - P[2][17]*SH_MAG[1] + P[3][17]*SH_MAG[2] - P[16][17]*(SH_MAG[3] + SH_MAG[4] - SH_MAG[5] - SH_MAG[6]) + P[17][17]*(2.0f*q0*q3 + 2.0f*q1*q2) - P[18][17]*(2.0f*q0*q2 - 2.0f*q1*q3) + P[0][17]*(SH_MAG[7] + SH_MAG[8] - 2.0f*magD*q2)) - (2.0f*q0*q2 - 2.0f*q1*q3)*(P[19][18] + P[1][18]*SH_MAG[0] - P[2][18]*SH_MAG[1] + P[3][18]*SH_MAG[2] - P[16][18]*(SH_MAG[3] + SH_MAG[4] - SH_MAG[5] - SH_MAG[6]) + P[17][18]*(2.0f*q0*q3 + 2.0f*q1*q2) - P[18][18]*(2.0f*q0*q2 - 2.0f*q1*q3) + P[0][18]*(SH_MAG[7] + SH_MAG[8] - 2.0f*magD*q2)) + (SH_MAG[7] + SH_MAG[8] - 2.0f*magD*q2)*(P[19][0] + P[1][0]*SH_MAG[0] - P[2][0]*SH_MAG[1] + P[3][0]*SH_MAG[2] - P[16][0]*(SH_MAG[3] + SH_MAG[4] - SH_MAG[5] - SH_MAG[6]) + P[17][0]*(2.0f*q0*q3 + 2.0f*q1*q2) - P[18][0]*(2.0f*q0*q2 - 2.0f*q1*q3) + P[0][0]*(SH_MAG[7] + SH_MAG[8] - 2.0f*magD*q2)) + P[17][19]*(2.0f*q0*q3 + 2.0f*q1*q2) - P[18][19]*(2.0f*q0*q2 - 2.0f*q1*q3) + SH_MAG[0]*(P[19][1] + P[1][1]*SH_MAG[0] - P[2][1]*SH_MAG[1] + P[3][1]*SH_MAG[2] - P[16][1]*(SH_MAG[3] + SH_MAG[4] - SH_MAG[5] - SH_MAG[6]) + P[17][1]*(2.0f*q0*q3 + 2.0f*q1*q2) - P[18][1]*(2.0f*q0*q2 - 2.0f*q1*q3) + P[0][1]*(SH_MAG[7] + SH_MAG[8] - 2.0f*magD*q2)) - SH_MAG[1]*(P[19][2] + P[1][2]*SH_MAG[0] - P[2][2]*SH_MAG[1] + P[3][2]*SH_MAG[2] - P[16][2]*(SH_MAG[3] + SH_MAG[4] - SH_MAG[5] - SH_MAG[6]) + P[17][2]*(2.0f*q0*q3 + 2.0f*q1*q2) - P[18][2]*(2.0f*q0*q2 - 2.0f*q1*q3) + P[0][2]*(SH_MAG[7] + SH_MAG[8] - 2.0f*magD*q2)) + SH_MAG[2]*(P[19][3] + P[1][3]*SH_MAG[0] - P[2][3]*SH_MAG[1] + P[3][3]*SH_MAG[2] - P[16][3]*(SH_MAG[3] + SH_MAG[4] - SH_MAG[5] - SH_MAG[6]) + P[17][3]*(2.0f*q0*q3 + 2.0f*q1*q2) - P[18][3]*(2.0f*q0*q2 - 2.0f*q1*q3) + P[0][3]*(SH_MAG[7] + SH_MAG[8] - 2.0f*magD*q2)) - (SH_MAG[3] + SH_MAG[4] - SH_MAG[5] - SH_MAG[6])*(P[19][16] + P[1][16]*SH_MAG[0] - P[2][16]*SH_MAG[1] + P[3][16]*SH_MAG[2] - P[16][16]*(SH_MAG[3] + SH_MAG[4] - SH_MAG[5] - SH_MAG[6]) + P[17][16]*(2.0f*q0*q3 + 2.0f*q1*q2) - P[18][16]*(2.0f*q0*q2 - 2.0f*q1*q3) + P[0][16]*(SH_MAG[7] + SH_MAG[8] - 2.0f*magD*q2)) + P[0][19]*(SH_MAG[7] + SH_MAG[8] - 2.0f*magD*q2));
96  // check for a badly conditioned covariance matrix
97 
98  if (_mag_innov_var[0] >= R_MAG) {
99  // the innovation variance contribution from the state covariances is non-negative - no fault
100  _fault_status.flags.bad_mag_x = false;
101 
102  } else {
103  // the innovation variance contribution from the state covariances is negative which means the covariance matrix is badly conditioned
105 
106  // we need to re-initialise covariances and abort this fusion step
108  ECL_ERR_TIMESTAMPED("EKF magX fusion numerical error - covariance reset");
109  return;
110 
111  }
112 
113  // Y axis innovation variance
114  _mag_innov_var[1] = (P[20][20] + R_MAG + P[0][20]*SH_MAG[2] + P[1][20]*SH_MAG[1] + P[2][20]*SH_MAG[0] - P[17][20]*(SH_MAG[3] - SH_MAG[4] + SH_MAG[5] - SH_MAG[6]) - (2.0f*q0*q3 - 2.0f*q1*q2)*(P[20][16] + P[0][16]*SH_MAG[2] + P[1][16]*SH_MAG[1] + P[2][16]*SH_MAG[0] - P[17][16]*(SH_MAG[3] - SH_MAG[4] + SH_MAG[5] - SH_MAG[6]) - P[16][16]*(2.0f*q0*q3 - 2.0f*q1*q2) + P[18][16]*(2.0f*q0*q1 + 2.0f*q2*q3) - P[3][16]*(SH_MAG[7] + SH_MAG[8] - 2.0f*magD*q2)) + (2.0f*q0*q1 + 2.0f*q2*q3)*(P[20][18] + P[0][18]*SH_MAG[2] + P[1][18]*SH_MAG[1] + P[2][18]*SH_MAG[0] - P[17][18]*(SH_MAG[3] - SH_MAG[4] + SH_MAG[5] - SH_MAG[6]) - P[16][18]*(2.0f*q0*q3 - 2.0f*q1*q2) + P[18][18]*(2.0f*q0*q1 + 2.0f*q2*q3) - P[3][18]*(SH_MAG[7] + SH_MAG[8] - 2.0f*magD*q2)) - (SH_MAG[7] + SH_MAG[8] - 2.0f*magD*q2)*(P[20][3] + P[0][3]*SH_MAG[2] + P[1][3]*SH_MAG[1] + P[2][3]*SH_MAG[0] - P[17][3]*(SH_MAG[3] - SH_MAG[4] + SH_MAG[5] - SH_MAG[6]) - P[16][3]*(2.0f*q0*q3 - 2.0f*q1*q2) + P[18][3]*(2.0f*q0*q1 + 2.0f*q2*q3) - P[3][3]*(SH_MAG[7] + SH_MAG[8] - 2.0f*magD*q2)) - P[16][20]*(2.0f*q0*q3 - 2.0f*q1*q2) + P[18][20]*(2.0f*q0*q1 + 2.0f*q2*q3) + SH_MAG[2]*(P[20][0] + P[0][0]*SH_MAG[2] + P[1][0]*SH_MAG[1] + P[2][0]*SH_MAG[0] - P[17][0]*(SH_MAG[3] - SH_MAG[4] + SH_MAG[5] - SH_MAG[6]) - P[16][0]*(2.0f*q0*q3 - 2.0f*q1*q2) + P[18][0]*(2.0f*q0*q1 + 2.0f*q2*q3) - P[3][0]*(SH_MAG[7] + SH_MAG[8] - 2.0f*magD*q2)) + SH_MAG[1]*(P[20][1] + P[0][1]*SH_MAG[2] + P[1][1]*SH_MAG[1] + P[2][1]*SH_MAG[0] - P[17][1]*(SH_MAG[3] - SH_MAG[4] + SH_MAG[5] - SH_MAG[6]) - P[16][1]*(2.0f*q0*q3 - 2.0f*q1*q2) + P[18][1]*(2.0f*q0*q1 + 2.0f*q2*q3) - P[3][1]*(SH_MAG[7] + SH_MAG[8] - 2.0f*magD*q2)) + SH_MAG[0]*(P[20][2] + P[0][2]*SH_MAG[2] + P[1][2]*SH_MAG[1] + P[2][2]*SH_MAG[0] - P[17][2]*(SH_MAG[3] - SH_MAG[4] + SH_MAG[5] - SH_MAG[6]) - P[16][2]*(2.0f*q0*q3 - 2.0f*q1*q2) + P[18][2]*(2.0f*q0*q1 + 2.0f*q2*q3) - P[3][2]*(SH_MAG[7] + SH_MAG[8] - 2.0f*magD*q2)) - (SH_MAG[3] - SH_MAG[4] + SH_MAG[5] - SH_MAG[6])*(P[20][17] + P[0][17]*SH_MAG[2] + P[1][17]*SH_MAG[1] + P[2][17]*SH_MAG[0] - P[17][17]*(SH_MAG[3] - SH_MAG[4] + SH_MAG[5] - SH_MAG[6]) - P[16][17]*(2.0f*q0*q3 - 2.0f*q1*q2) + P[18][17]*(2.0f*q0*q1 + 2.0f*q2*q3) - P[3][17]*(SH_MAG[7] + SH_MAG[8] - 2.0f*magD*q2)) - P[3][20]*(SH_MAG[7] + SH_MAG[8] - 2.0f*magD*q2));
115 
116  // check for a badly conditioned covariance matrix
117  if (_mag_innov_var[1] >= R_MAG) {
118  // the innovation variance contribution from the state covariances is non-negative - no fault
119  _fault_status.flags.bad_mag_y = false;
120 
121  } else {
122  // the innovation variance contribution from the state covariances is negtive which means the covariance matrix is badly conditioned
124 
125  // we need to re-initialise covariances and abort this fusion step
127  ECL_ERR_TIMESTAMPED("EKF magY fusion numerical error - covariance reset");
128  return;
129 
130  }
131 
132  // Z axis innovation variance
133  _mag_innov_var[2] = (P[21][21] + R_MAG + P[0][21]*SH_MAG[1] - P[1][21]*SH_MAG[2] + P[3][21]*SH_MAG[0] + P[18][21]*(SH_MAG[3] - SH_MAG[4] - SH_MAG[5] + SH_MAG[6]) + (2.0f*q0*q2 + 2.0f*q1*q3)*(P[21][16] + P[0][16]*SH_MAG[1] - P[1][16]*SH_MAG[2] + P[3][16]*SH_MAG[0] + P[18][16]*(SH_MAG[3] - SH_MAG[4] - SH_MAG[5] + SH_MAG[6]) + P[16][16]*(2.0f*q0*q2 + 2.0f*q1*q3) - P[17][16]*(2.0f*q0*q1 - 2.0f*q2*q3) + P[2][16]*(SH_MAG[7] + SH_MAG[8] - 2.0f*magD*q2)) - (2.0f*q0*q1 - 2.0f*q2*q3)*(P[21][17] + P[0][17]*SH_MAG[1] - P[1][17]*SH_MAG[2] + P[3][17]*SH_MAG[0] + P[18][17]*(SH_MAG[3] - SH_MAG[4] - SH_MAG[5] + SH_MAG[6]) + P[16][17]*(2.0f*q0*q2 + 2.0f*q1*q3) - P[17][17]*(2.0f*q0*q1 - 2.0f*q2*q3) + P[2][17]*(SH_MAG[7] + SH_MAG[8] - 2.0f*magD*q2)) + (SH_MAG[7] + SH_MAG[8] - 2.0f*magD*q2)*(P[21][2] + P[0][2]*SH_MAG[1] - P[1][2]*SH_MAG[2] + P[3][2]*SH_MAG[0] + P[18][2]*(SH_MAG[3] - SH_MAG[4] - SH_MAG[5] + SH_MAG[6]) + P[16][2]*(2.0f*q0*q2 + 2.0f*q1*q3) - P[17][2]*(2.0f*q0*q1 - 2.0f*q2*q3) + P[2][2]*(SH_MAG[7] + SH_MAG[8] - 2.0f*magD*q2)) + P[16][21]*(2.0f*q0*q2 + 2.0f*q1*q3) - P[17][21]*(2.0f*q0*q1 - 2.0f*q2*q3) + SH_MAG[1]*(P[21][0] + P[0][0]*SH_MAG[1] - P[1][0]*SH_MAG[2] + P[3][0]*SH_MAG[0] + P[18][0]*(SH_MAG[3] - SH_MAG[4] - SH_MAG[5] + SH_MAG[6]) + P[16][0]*(2.0f*q0*q2 + 2.0f*q1*q3) - P[17][0]*(2.0f*q0*q1 - 2.0f*q2*q3) + P[2][0]*(SH_MAG[7] + SH_MAG[8] - 2.0f*magD*q2)) - SH_MAG[2]*(P[21][1] + P[0][1]*SH_MAG[1] - P[1][1]*SH_MAG[2] + P[3][1]*SH_MAG[0] + P[18][1]*(SH_MAG[3] - SH_MAG[4] - SH_MAG[5] + SH_MAG[6]) + P[16][1]*(2.0f*q0*q2 + 2.0f*q1*q3) - P[17][1]*(2.0f*q0*q1 - 2.0f*q2*q3) + P[2][1]*(SH_MAG[7] + SH_MAG[8] - 2.0f*magD*q2)) + SH_MAG[0]*(P[21][3] + P[0][3]*SH_MAG[1] - P[1][3]*SH_MAG[2] + P[3][3]*SH_MAG[0] + P[18][3]*(SH_MAG[3] - SH_MAG[4] - SH_MAG[5] + SH_MAG[6]) + P[16][3]*(2.0f*q0*q2 + 2.0f*q1*q3) - P[17][3]*(2.0f*q0*q1 - 2.0f*q2*q3) + P[2][3]*(SH_MAG[7] + SH_MAG[8] - 2.0f*magD*q2)) + (SH_MAG[3] - SH_MAG[4] - SH_MAG[5] + SH_MAG[6])*(P[21][18] + P[0][18]*SH_MAG[1] - P[1][18]*SH_MAG[2] + P[3][18]*SH_MAG[0] + P[18][18]*(SH_MAG[3] - SH_MAG[4] - SH_MAG[5] + SH_MAG[6]) + P[16][18]*(2.0f*q0*q2 + 2.0f*q1*q3) - P[17][18]*(2.0f*q0*q1 - 2.0f*q2*q3) + P[2][18]*(SH_MAG[7] + SH_MAG[8] - 2.0f*magD*q2)) + P[2][21]*(SH_MAG[7] + SH_MAG[8] - 2.0f*magD*q2));
134 
135  // check for a badly conditioned covariance matrix
136  if (_mag_innov_var[2] >= R_MAG) {
137  // the innovation variance contribution from the state covariances is non-negative - no fault
138  _fault_status.flags.bad_mag_z = false;
139 
140  } else if (_mag_innov_var[2] > 0.0f) {
141  // the innovation variance contribution from the state covariances is negative which means the covariance matrix is badly conditioned
143 
144  // we need to re-initialise covariances and abort this fusion step
146  ECL_ERR_TIMESTAMPED("EKF magZ fusion numerical error - covariance reset");
147  return;
148 
149  }
150 
151  // Perform an innovation consistency check and report the result
152  bool healthy = true;
153 
154  for (uint8_t index = 0; index <= 2; index++) {
155  _mag_test_ratio[index] = sq(_mag_innov[index]) / (sq(math::max(_params.mag_innov_gate, 1.0f)) * _mag_innov_var[index]);
156 
157  if (_mag_test_ratio[index] > 1.0f) {
158  healthy = false;
159  _innov_check_fail_status.value |= (1 << (index + 3));
160 
161  } else {
162  _innov_check_fail_status.value &= ~(1 << (index + 3));
163  }
164  }
165 
166  // we are no longer using heading fusion so set the reported test level to zero
167  _yaw_test_ratio = 0.0f;
168 
169  // if any axis fails, abort the mag fusion
170  if (!healthy) {
171  return;
172  }
173 
174  bool update_all_states = !_flt_mag_align_converging;
175 
176  // update the states and covariance using sequential fusion of the magnetometer components
177  for (uint8_t index = 0; index <= 2; index++) {
178 
179  // Calculate Kalman gains and observation jacobians
180  if (index == 0) {
181  // Calculate X axis observation jacobians
182  memset(H_MAG, 0, sizeof(H_MAG));
183  H_MAG[0] = SH_MAG[7] + SH_MAG[8] - 2.0f*magD*q2;
184  H_MAG[1] = SH_MAG[0];
185  H_MAG[2] = -SH_MAG[1];
186  H_MAG[3] = SH_MAG[2];
187  H_MAG[16] = SH_MAG[5] - SH_MAG[4] - SH_MAG[3] + SH_MAG[6];
188  H_MAG[17] = 2.0f*q0*q3 + 2.0f*q1*q2;
189  H_MAG[18] = 2.0f*q1*q3 - 2.0f*q0*q2;
190  H_MAG[19] = 1.0f;
191 
192  // Calculate X axis Kalman gains
193  float SK_MX[5];
194  SK_MX[0] = 1.0f / _mag_innov_var[0];
195  SK_MX[1] = SH_MAG[3] + SH_MAG[4] - SH_MAG[5] - SH_MAG[6];
196  SK_MX[2] = SH_MAG[7] + SH_MAG[8] - 2.0f*magD*q2;
197  SK_MX[3] = 2.0f*q0*q2 - 2.0f*q1*q3;
198  SK_MX[4] = 2.0f*q0*q3 + 2.0f*q1*q2;
199 
200  if (update_all_states) {
201  Kfusion[0] = SK_MX[0]*(P[0][19] + P[0][1]*SH_MAG[0] - P[0][2]*SH_MAG[1] + P[0][3]*SH_MAG[2] + P[0][0]*SK_MX[2] - P[0][16]*SK_MX[1] + P[0][17]*SK_MX[4] - P[0][18]*SK_MX[3]);
202  Kfusion[1] = SK_MX[0]*(P[1][19] + P[1][1]*SH_MAG[0] - P[1][2]*SH_MAG[1] + P[1][3]*SH_MAG[2] + P[1][0]*SK_MX[2] - P[1][16]*SK_MX[1] + P[1][17]*SK_MX[4] - P[1][18]*SK_MX[3]);
203  Kfusion[2] = SK_MX[0]*(P[2][19] + P[2][1]*SH_MAG[0] - P[2][2]*SH_MAG[1] + P[2][3]*SH_MAG[2] + P[2][0]*SK_MX[2] - P[2][16]*SK_MX[1] + P[2][17]*SK_MX[4] - P[2][18]*SK_MX[3]);
204  Kfusion[3] = SK_MX[0]*(P[3][19] + P[3][1]*SH_MAG[0] - P[3][2]*SH_MAG[1] + P[3][3]*SH_MAG[2] + P[3][0]*SK_MX[2] - P[3][16]*SK_MX[1] + P[3][17]*SK_MX[4] - P[3][18]*SK_MX[3]);
205  Kfusion[4] = SK_MX[0]*(P[4][19] + P[4][1]*SH_MAG[0] - P[4][2]*SH_MAG[1] + P[4][3]*SH_MAG[2] + P[4][0]*SK_MX[2] - P[4][16]*SK_MX[1] + P[4][17]*SK_MX[4] - P[4][18]*SK_MX[3]);
206  Kfusion[5] = SK_MX[0]*(P[5][19] + P[5][1]*SH_MAG[0] - P[5][2]*SH_MAG[1] + P[5][3]*SH_MAG[2] + P[5][0]*SK_MX[2] - P[5][16]*SK_MX[1] + P[5][17]*SK_MX[4] - P[5][18]*SK_MX[3]);
207  Kfusion[6] = SK_MX[0]*(P[6][19] + P[6][1]*SH_MAG[0] - P[6][2]*SH_MAG[1] + P[6][3]*SH_MAG[2] + P[6][0]*SK_MX[2] - P[6][16]*SK_MX[1] + P[6][17]*SK_MX[4] - P[6][18]*SK_MX[3]);
208  Kfusion[7] = SK_MX[0]*(P[7][19] + P[7][1]*SH_MAG[0] - P[7][2]*SH_MAG[1] + P[7][3]*SH_MAG[2] + P[7][0]*SK_MX[2] - P[7][16]*SK_MX[1] + P[7][17]*SK_MX[4] - P[7][18]*SK_MX[3]);
209  Kfusion[8] = SK_MX[0]*(P[8][19] + P[8][1]*SH_MAG[0] - P[8][2]*SH_MAG[1] + P[8][3]*SH_MAG[2] + P[8][0]*SK_MX[2] - P[8][16]*SK_MX[1] + P[8][17]*SK_MX[4] - P[8][18]*SK_MX[3]);
210  Kfusion[9] = SK_MX[0]*(P[9][19] + P[9][1]*SH_MAG[0] - P[9][2]*SH_MAG[1] + P[9][3]*SH_MAG[2] + P[9][0]*SK_MX[2] - P[9][16]*SK_MX[1] + P[9][17]*SK_MX[4] - P[9][18]*SK_MX[3]);
211  Kfusion[10] = SK_MX[0]*(P[10][19] + P[10][1]*SH_MAG[0] - P[10][2]*SH_MAG[1] + P[10][3]*SH_MAG[2] + P[10][0]*SK_MX[2] - P[10][16]*SK_MX[1] + P[10][17]*SK_MX[4] - P[10][18]*SK_MX[3]);
212  Kfusion[11] = SK_MX[0]*(P[11][19] + P[11][1]*SH_MAG[0] - P[11][2]*SH_MAG[1] + P[11][3]*SH_MAG[2] + P[11][0]*SK_MX[2] - P[11][16]*SK_MX[1] + P[11][17]*SK_MX[4] - P[11][18]*SK_MX[3]);
213  Kfusion[12] = SK_MX[0]*(P[12][19] + P[12][1]*SH_MAG[0] - P[12][2]*SH_MAG[1] + P[12][3]*SH_MAG[2] + P[12][0]*SK_MX[2] - P[12][16]*SK_MX[1] + P[12][17]*SK_MX[4] - P[12][18]*SK_MX[3]);
214  Kfusion[13] = SK_MX[0]*(P[13][19] + P[13][1]*SH_MAG[0] - P[13][2]*SH_MAG[1] + P[13][3]*SH_MAG[2] + P[13][0]*SK_MX[2] - P[13][16]*SK_MX[1] + P[13][17]*SK_MX[4] - P[13][18]*SK_MX[3]);
215  Kfusion[14] = SK_MX[0]*(P[14][19] + P[14][1]*SH_MAG[0] - P[14][2]*SH_MAG[1] + P[14][3]*SH_MAG[2] + P[14][0]*SK_MX[2] - P[14][16]*SK_MX[1] + P[14][17]*SK_MX[4] - P[14][18]*SK_MX[3]);
216  Kfusion[15] = SK_MX[0]*(P[15][19] + P[15][1]*SH_MAG[0] - P[15][2]*SH_MAG[1] + P[15][3]*SH_MAG[2] + P[15][0]*SK_MX[2] - P[15][16]*SK_MX[1] + P[15][17]*SK_MX[4] - P[15][18]*SK_MX[3]);
217  Kfusion[22] = SK_MX[0]*(P[22][19] + P[22][1]*SH_MAG[0] - P[22][2]*SH_MAG[1] + P[22][3]*SH_MAG[2] + P[22][0]*SK_MX[2] - P[22][16]*SK_MX[1] + P[22][17]*SK_MX[4] - P[22][18]*SK_MX[3]);
218  Kfusion[23] = SK_MX[0]*(P[23][19] + P[23][1]*SH_MAG[0] - P[23][2]*SH_MAG[1] + P[23][3]*SH_MAG[2] + P[23][0]*SK_MX[2] - P[23][16]*SK_MX[1] + P[23][17]*SK_MX[4] - P[23][18]*SK_MX[3]);
219  } else {
220  for (uint8_t i = 0; i < 16; i++) {
221  Kfusion[i] = 0.0f;
222  }
223 
224  Kfusion[22] = 0.0f;
225  Kfusion[23] = 0.0f;
226  }
227 
228  Kfusion[16] = SK_MX[0]*(P[16][19] + P[16][1]*SH_MAG[0] - P[16][2]*SH_MAG[1] + P[16][3]*SH_MAG[2] + P[16][0]*SK_MX[2] - P[16][16]*SK_MX[1] + P[16][17]*SK_MX[4] - P[16][18]*SK_MX[3]);
229  Kfusion[17] = SK_MX[0]*(P[17][19] + P[17][1]*SH_MAG[0] - P[17][2]*SH_MAG[1] + P[17][3]*SH_MAG[2] + P[17][0]*SK_MX[2] - P[17][16]*SK_MX[1] + P[17][17]*SK_MX[4] - P[17][18]*SK_MX[3]);
230  Kfusion[18] = SK_MX[0]*(P[18][19] + P[18][1]*SH_MAG[0] - P[18][2]*SH_MAG[1] + P[18][3]*SH_MAG[2] + P[18][0]*SK_MX[2] - P[18][16]*SK_MX[1] + P[18][17]*SK_MX[4] - P[18][18]*SK_MX[3]);
231  Kfusion[19] = SK_MX[0]*(P[19][19] + P[19][1]*SH_MAG[0] - P[19][2]*SH_MAG[1] + P[19][3]*SH_MAG[2] + P[19][0]*SK_MX[2] - P[19][16]*SK_MX[1] + P[19][17]*SK_MX[4] - P[19][18]*SK_MX[3]);
232  Kfusion[20] = SK_MX[0]*(P[20][19] + P[20][1]*SH_MAG[0] - P[20][2]*SH_MAG[1] + P[20][3]*SH_MAG[2] + P[20][0]*SK_MX[2] - P[20][16]*SK_MX[1] + P[20][17]*SK_MX[4] - P[20][18]*SK_MX[3]);
233  Kfusion[21] = SK_MX[0]*(P[21][19] + P[21][1]*SH_MAG[0] - P[21][2]*SH_MAG[1] + P[21][3]*SH_MAG[2] + P[21][0]*SK_MX[2] - P[21][16]*SK_MX[1] + P[21][17]*SK_MX[4] - P[21][18]*SK_MX[3]);
234 
235  } else if (index == 1) {
236  // Calculate Y axis observation jacobians
237  memset(H_MAG, 0, sizeof(H_MAG));
238  H_MAG[0] = SH_MAG[2];
239  H_MAG[1] = SH_MAG[1];
240  H_MAG[2] = SH_MAG[0];
241  H_MAG[3] = 2.0f*magD*q2 - SH_MAG[8] - SH_MAG[7];
242  H_MAG[16] = 2.0f*q1*q2 - 2.0f*q0*q3;
243  H_MAG[17] = SH_MAG[4] - SH_MAG[3] - SH_MAG[5] + SH_MAG[6];
244  H_MAG[18] = 2.0f*q0*q1 + 2.0f*q2*q3;
245  H_MAG[20] = 1.0f;
246 
247  // Calculate Y axis Kalman gains
248  float SK_MY[5];
249  SK_MY[0] = 1.0f / _mag_innov_var[1];
250  SK_MY[1] = SH_MAG[3] - SH_MAG[4] + SH_MAG[5] - SH_MAG[6];
251  SK_MY[2] = SH_MAG[7] + SH_MAG[8] - 2.0f*magD*q2;
252  SK_MY[3] = 2.0f*q0*q3 - 2.0f*q1*q2;
253  SK_MY[4] = 2.0f*q0*q1 + 2.0f*q2*q3;
254 
255  if (update_all_states) {
256  Kfusion[0] = SK_MY[0]*(P[0][20] + P[0][0]*SH_MAG[2] + P[0][1]*SH_MAG[1] + P[0][2]*SH_MAG[0] - P[0][3]*SK_MY[2] - P[0][17]*SK_MY[1] - P[0][16]*SK_MY[3] + P[0][18]*SK_MY[4]);
257  Kfusion[1] = SK_MY[0]*(P[1][20] + P[1][0]*SH_MAG[2] + P[1][1]*SH_MAG[1] + P[1][2]*SH_MAG[0] - P[1][3]*SK_MY[2] - P[1][17]*SK_MY[1] - P[1][16]*SK_MY[3] + P[1][18]*SK_MY[4]);
258  Kfusion[2] = SK_MY[0]*(P[2][20] + P[2][0]*SH_MAG[2] + P[2][1]*SH_MAG[1] + P[2][2]*SH_MAG[0] - P[2][3]*SK_MY[2] - P[2][17]*SK_MY[1] - P[2][16]*SK_MY[3] + P[2][18]*SK_MY[4]);
259  Kfusion[3] = SK_MY[0]*(P[3][20] + P[3][0]*SH_MAG[2] + P[3][1]*SH_MAG[1] + P[3][2]*SH_MAG[0] - P[3][3]*SK_MY[2] - P[3][17]*SK_MY[1] - P[3][16]*SK_MY[3] + P[3][18]*SK_MY[4]);
260  Kfusion[4] = SK_MY[0]*(P[4][20] + P[4][0]*SH_MAG[2] + P[4][1]*SH_MAG[1] + P[4][2]*SH_MAG[0] - P[4][3]*SK_MY[2] - P[4][17]*SK_MY[1] - P[4][16]*SK_MY[3] + P[4][18]*SK_MY[4]);
261  Kfusion[5] = SK_MY[0]*(P[5][20] + P[5][0]*SH_MAG[2] + P[5][1]*SH_MAG[1] + P[5][2]*SH_MAG[0] - P[5][3]*SK_MY[2] - P[5][17]*SK_MY[1] - P[5][16]*SK_MY[3] + P[5][18]*SK_MY[4]);
262  Kfusion[6] = SK_MY[0]*(P[6][20] + P[6][0]*SH_MAG[2] + P[6][1]*SH_MAG[1] + P[6][2]*SH_MAG[0] - P[6][3]*SK_MY[2] - P[6][17]*SK_MY[1] - P[6][16]*SK_MY[3] + P[6][18]*SK_MY[4]);
263  Kfusion[7] = SK_MY[0]*(P[7][20] + P[7][0]*SH_MAG[2] + P[7][1]*SH_MAG[1] + P[7][2]*SH_MAG[0] - P[7][3]*SK_MY[2] - P[7][17]*SK_MY[1] - P[7][16]*SK_MY[3] + P[7][18]*SK_MY[4]);
264  Kfusion[8] = SK_MY[0]*(P[8][20] + P[8][0]*SH_MAG[2] + P[8][1]*SH_MAG[1] + P[8][2]*SH_MAG[0] - P[8][3]*SK_MY[2] - P[8][17]*SK_MY[1] - P[8][16]*SK_MY[3] + P[8][18]*SK_MY[4]);
265  Kfusion[9] = SK_MY[0]*(P[9][20] + P[9][0]*SH_MAG[2] + P[9][1]*SH_MAG[1] + P[9][2]*SH_MAG[0] - P[9][3]*SK_MY[2] - P[9][17]*SK_MY[1] - P[9][16]*SK_MY[3] + P[9][18]*SK_MY[4]);
266  Kfusion[10] = SK_MY[0]*(P[10][20] + P[10][0]*SH_MAG[2] + P[10][1]*SH_MAG[1] + P[10][2]*SH_MAG[0] - P[10][3]*SK_MY[2] - P[10][17]*SK_MY[1] - P[10][16]*SK_MY[3] + P[10][18]*SK_MY[4]);
267  Kfusion[11] = SK_MY[0]*(P[11][20] + P[11][0]*SH_MAG[2] + P[11][1]*SH_MAG[1] + P[11][2]*SH_MAG[0] - P[11][3]*SK_MY[2] - P[11][17]*SK_MY[1] - P[11][16]*SK_MY[3] + P[11][18]*SK_MY[4]);
268  Kfusion[12] = SK_MY[0]*(P[12][20] + P[12][0]*SH_MAG[2] + P[12][1]*SH_MAG[1] + P[12][2]*SH_MAG[0] - P[12][3]*SK_MY[2] - P[12][17]*SK_MY[1] - P[12][16]*SK_MY[3] + P[12][18]*SK_MY[4]);
269  Kfusion[13] = SK_MY[0]*(P[13][20] + P[13][0]*SH_MAG[2] + P[13][1]*SH_MAG[1] + P[13][2]*SH_MAG[0] - P[13][3]*SK_MY[2] - P[13][17]*SK_MY[1] - P[13][16]*SK_MY[3] + P[13][18]*SK_MY[4]);
270  Kfusion[14] = SK_MY[0]*(P[14][20] + P[14][0]*SH_MAG[2] + P[14][1]*SH_MAG[1] + P[14][2]*SH_MAG[0] - P[14][3]*SK_MY[2] - P[14][17]*SK_MY[1] - P[14][16]*SK_MY[3] + P[14][18]*SK_MY[4]);
271  Kfusion[15] = SK_MY[0]*(P[15][20] + P[15][0]*SH_MAG[2] + P[15][1]*SH_MAG[1] + P[15][2]*SH_MAG[0] - P[15][3]*SK_MY[2] - P[15][17]*SK_MY[1] - P[15][16]*SK_MY[3] + P[15][18]*SK_MY[4]);
272  Kfusion[22] = SK_MY[0]*(P[22][20] + P[22][0]*SH_MAG[2] + P[22][1]*SH_MAG[1] + P[22][2]*SH_MAG[0] - P[22][3]*SK_MY[2] - P[22][17]*SK_MY[1] - P[22][16]*SK_MY[3] + P[22][18]*SK_MY[4]);
273  Kfusion[23] = SK_MY[0]*(P[23][20] + P[23][0]*SH_MAG[2] + P[23][1]*SH_MAG[1] + P[23][2]*SH_MAG[0] - P[23][3]*SK_MY[2] - P[23][17]*SK_MY[1] - P[23][16]*SK_MY[3] + P[23][18]*SK_MY[4]);
274  } else {
275  for (uint8_t i = 0; i < 16; i++) {
276  Kfusion[i] = 0.0f;
277  }
278 
279  Kfusion[22] = 0.0f;
280  Kfusion[23] = 0.0f;
281  }
282 
283  Kfusion[16] = SK_MY[0]*(P[16][20] + P[16][0]*SH_MAG[2] + P[16][1]*SH_MAG[1] + P[16][2]*SH_MAG[0] - P[16][3]*SK_MY[2] - P[16][17]*SK_MY[1] - P[16][16]*SK_MY[3] + P[16][18]*SK_MY[4]);
284  Kfusion[17] = SK_MY[0]*(P[17][20] + P[17][0]*SH_MAG[2] + P[17][1]*SH_MAG[1] + P[17][2]*SH_MAG[0] - P[17][3]*SK_MY[2] - P[17][17]*SK_MY[1] - P[17][16]*SK_MY[3] + P[17][18]*SK_MY[4]);
285  Kfusion[18] = SK_MY[0]*(P[18][20] + P[18][0]*SH_MAG[2] + P[18][1]*SH_MAG[1] + P[18][2]*SH_MAG[0] - P[18][3]*SK_MY[2] - P[18][17]*SK_MY[1] - P[18][16]*SK_MY[3] + P[18][18]*SK_MY[4]);
286  Kfusion[19] = SK_MY[0]*(P[19][20] + P[19][0]*SH_MAG[2] + P[19][1]*SH_MAG[1] + P[19][2]*SH_MAG[0] - P[19][3]*SK_MY[2] - P[19][17]*SK_MY[1] - P[19][16]*SK_MY[3] + P[19][18]*SK_MY[4]);
287  Kfusion[20] = SK_MY[0]*(P[20][20] + P[20][0]*SH_MAG[2] + P[20][1]*SH_MAG[1] + P[20][2]*SH_MAG[0] - P[20][3]*SK_MY[2] - P[20][17]*SK_MY[1] - P[20][16]*SK_MY[3] + P[20][18]*SK_MY[4]);
288  Kfusion[21] = SK_MY[0]*(P[21][20] + P[21][0]*SH_MAG[2] + P[21][1]*SH_MAG[1] + P[21][2]*SH_MAG[0] - P[21][3]*SK_MY[2] - P[21][17]*SK_MY[1] - P[21][16]*SK_MY[3] + P[21][18]*SK_MY[4]);
289 
290  } else if (index == 2) {
291 
292  // we do not fuse synthesized magnetomter measurements when doing 3D fusion
294  continue;
295  }
296 
297  // calculate Z axis observation jacobians
298  memset(H_MAG, 0, sizeof(H_MAG));
299  H_MAG[0] = SH_MAG[1];
300  H_MAG[1] = -SH_MAG[2];
301  H_MAG[2] = SH_MAG[7] + SH_MAG[8] - 2.0f*magD*q2;
302  H_MAG[3] = SH_MAG[0];
303  H_MAG[16] = 2.0f*q0*q2 + 2.0f*q1*q3;
304  H_MAG[17] = 2.0f*q2*q3 - 2.0f*q0*q1;
305  H_MAG[18] = SH_MAG[3] - SH_MAG[4] - SH_MAG[5] + SH_MAG[6];
306  H_MAG[21] = 1.0f;
307 
308  // Calculate Z axis Kalman gains
309  float SK_MZ[5];
310  SK_MZ[0] = 1.0f / _mag_innov_var[2];
311  SK_MZ[1] = SH_MAG[3] - SH_MAG[4] - SH_MAG[5] + SH_MAG[6];
312  SK_MZ[2] = SH_MAG[7] + SH_MAG[8] - 2.0f*magD*q2;
313  SK_MZ[3] = 2.0f*q0*q1 - 2.0f*q2*q3;
314  SK_MZ[4] = 2.0f*q0*q2 + 2.0f*q1*q3;
315 
316  if (update_all_states) {
317  Kfusion[0] = SK_MZ[0]*(P[0][21] + P[0][0]*SH_MAG[1] - P[0][1]*SH_MAG[2] + P[0][3]*SH_MAG[0] + P[0][2]*SK_MZ[2] + P[0][18]*SK_MZ[1] + P[0][16]*SK_MZ[4] - P[0][17]*SK_MZ[3]);
318  Kfusion[1] = SK_MZ[0]*(P[1][21] + P[1][0]*SH_MAG[1] - P[1][1]*SH_MAG[2] + P[1][3]*SH_MAG[0] + P[1][2]*SK_MZ[2] + P[1][18]*SK_MZ[1] + P[1][16]*SK_MZ[4] - P[1][17]*SK_MZ[3]);
319  Kfusion[2] = SK_MZ[0]*(P[2][21] + P[2][0]*SH_MAG[1] - P[2][1]*SH_MAG[2] + P[2][3]*SH_MAG[0] + P[2][2]*SK_MZ[2] + P[2][18]*SK_MZ[1] + P[2][16]*SK_MZ[4] - P[2][17]*SK_MZ[3]);
320  Kfusion[3] = SK_MZ[0]*(P[3][21] + P[3][0]*SH_MAG[1] - P[3][1]*SH_MAG[2] + P[3][3]*SH_MAG[0] + P[3][2]*SK_MZ[2] + P[3][18]*SK_MZ[1] + P[3][16]*SK_MZ[4] - P[3][17]*SK_MZ[3]);
321  Kfusion[4] = SK_MZ[0]*(P[4][21] + P[4][0]*SH_MAG[1] - P[4][1]*SH_MAG[2] + P[4][3]*SH_MAG[0] + P[4][2]*SK_MZ[2] + P[4][18]*SK_MZ[1] + P[4][16]*SK_MZ[4] - P[4][17]*SK_MZ[3]);
322  Kfusion[5] = SK_MZ[0]*(P[5][21] + P[5][0]*SH_MAG[1] - P[5][1]*SH_MAG[2] + P[5][3]*SH_MAG[0] + P[5][2]*SK_MZ[2] + P[5][18]*SK_MZ[1] + P[5][16]*SK_MZ[4] - P[5][17]*SK_MZ[3]);
323  Kfusion[6] = SK_MZ[0]*(P[6][21] + P[6][0]*SH_MAG[1] - P[6][1]*SH_MAG[2] + P[6][3]*SH_MAG[0] + P[6][2]*SK_MZ[2] + P[6][18]*SK_MZ[1] + P[6][16]*SK_MZ[4] - P[6][17]*SK_MZ[3]);
324  Kfusion[7] = SK_MZ[0]*(P[7][21] + P[7][0]*SH_MAG[1] - P[7][1]*SH_MAG[2] + P[7][3]*SH_MAG[0] + P[7][2]*SK_MZ[2] + P[7][18]*SK_MZ[1] + P[7][16]*SK_MZ[4] - P[7][17]*SK_MZ[3]);
325  Kfusion[8] = SK_MZ[0]*(P[8][21] + P[8][0]*SH_MAG[1] - P[8][1]*SH_MAG[2] + P[8][3]*SH_MAG[0] + P[8][2]*SK_MZ[2] + P[8][18]*SK_MZ[1] + P[8][16]*SK_MZ[4] - P[8][17]*SK_MZ[3]);
326  Kfusion[9] = SK_MZ[0]*(P[9][21] + P[9][0]*SH_MAG[1] - P[9][1]*SH_MAG[2] + P[9][3]*SH_MAG[0] + P[9][2]*SK_MZ[2] + P[9][18]*SK_MZ[1] + P[9][16]*SK_MZ[4] - P[9][17]*SK_MZ[3]);
327  Kfusion[10] = SK_MZ[0]*(P[10][21] + P[10][0]*SH_MAG[1] - P[10][1]*SH_MAG[2] + P[10][3]*SH_MAG[0] + P[10][2]*SK_MZ[2] + P[10][18]*SK_MZ[1] + P[10][16]*SK_MZ[4] - P[10][17]*SK_MZ[3]);
328  Kfusion[11] = SK_MZ[0]*(P[11][21] + P[11][0]*SH_MAG[1] - P[11][1]*SH_MAG[2] + P[11][3]*SH_MAG[0] + P[11][2]*SK_MZ[2] + P[11][18]*SK_MZ[1] + P[11][16]*SK_MZ[4] - P[11][17]*SK_MZ[3]);
329  Kfusion[12] = SK_MZ[0]*(P[12][21] + P[12][0]*SH_MAG[1] - P[12][1]*SH_MAG[2] + P[12][3]*SH_MAG[0] + P[12][2]*SK_MZ[2] + P[12][18]*SK_MZ[1] + P[12][16]*SK_MZ[4] - P[12][17]*SK_MZ[3]);
330  Kfusion[13] = SK_MZ[0]*(P[13][21] + P[13][0]*SH_MAG[1] - P[13][1]*SH_MAG[2] + P[13][3]*SH_MAG[0] + P[13][2]*SK_MZ[2] + P[13][18]*SK_MZ[1] + P[13][16]*SK_MZ[4] - P[13][17]*SK_MZ[3]);
331  Kfusion[14] = SK_MZ[0]*(P[14][21] + P[14][0]*SH_MAG[1] - P[14][1]*SH_MAG[2] + P[14][3]*SH_MAG[0] + P[14][2]*SK_MZ[2] + P[14][18]*SK_MZ[1] + P[14][16]*SK_MZ[4] - P[14][17]*SK_MZ[3]);
332  Kfusion[15] = SK_MZ[0]*(P[15][21] + P[15][0]*SH_MAG[1] - P[15][1]*SH_MAG[2] + P[15][3]*SH_MAG[0] + P[15][2]*SK_MZ[2] + P[15][18]*SK_MZ[1] + P[15][16]*SK_MZ[4] - P[15][17]*SK_MZ[3]);
333  Kfusion[22] = SK_MZ[0]*(P[22][21] + P[22][0]*SH_MAG[1] - P[22][1]*SH_MAG[2] + P[22][3]*SH_MAG[0] + P[22][2]*SK_MZ[2] + P[22][18]*SK_MZ[1] + P[22][16]*SK_MZ[4] - P[22][17]*SK_MZ[3]);
334  Kfusion[23] = SK_MZ[0]*(P[23][21] + P[23][0]*SH_MAG[1] - P[23][1]*SH_MAG[2] + P[23][3]*SH_MAG[0] + P[23][2]*SK_MZ[2] + P[23][18]*SK_MZ[1] + P[23][16]*SK_MZ[4] - P[23][17]*SK_MZ[3]);
335  } else {
336  for (uint8_t i = 0; i < 16; i++) {
337  Kfusion[i] = 0.0f;
338  }
339 
340  Kfusion[22] = 0.0f;
341  Kfusion[23] = 0.0f;
342  }
343 
344  Kfusion[16] = SK_MZ[0]*(P[16][21] + P[16][0]*SH_MAG[1] - P[16][1]*SH_MAG[2] + P[16][3]*SH_MAG[0] + P[16][2]*SK_MZ[2] + P[16][18]*SK_MZ[1] + P[16][16]*SK_MZ[4] - P[16][17]*SK_MZ[3]);
345  Kfusion[17] = SK_MZ[0]*(P[17][21] + P[17][0]*SH_MAG[1] - P[17][1]*SH_MAG[2] + P[17][3]*SH_MAG[0] + P[17][2]*SK_MZ[2] + P[17][18]*SK_MZ[1] + P[17][16]*SK_MZ[4] - P[17][17]*SK_MZ[3]);
346  Kfusion[18] = SK_MZ[0]*(P[18][21] + P[18][0]*SH_MAG[1] - P[18][1]*SH_MAG[2] + P[18][3]*SH_MAG[0] + P[18][2]*SK_MZ[2] + P[18][18]*SK_MZ[1] + P[18][16]*SK_MZ[4] - P[18][17]*SK_MZ[3]);
347  Kfusion[19] = SK_MZ[0]*(P[19][21] + P[19][0]*SH_MAG[1] - P[19][1]*SH_MAG[2] + P[19][3]*SH_MAG[0] + P[19][2]*SK_MZ[2] + P[19][18]*SK_MZ[1] + P[19][16]*SK_MZ[4] - P[19][17]*SK_MZ[3]);
348  Kfusion[20] = SK_MZ[0]*(P[20][21] + P[20][0]*SH_MAG[1] - P[20][1]*SH_MAG[2] + P[20][3]*SH_MAG[0] + P[20][2]*SK_MZ[2] + P[20][18]*SK_MZ[1] + P[20][16]*SK_MZ[4] - P[20][17]*SK_MZ[3]);
349  Kfusion[21] = SK_MZ[0]*(P[21][21] + P[21][0]*SH_MAG[1] - P[21][1]*SH_MAG[2] + P[21][3]*SH_MAG[0] + P[21][2]*SK_MZ[2] + P[21][18]*SK_MZ[1] + P[21][16]*SK_MZ[4] - P[21][17]*SK_MZ[3]);
350 
351  } else {
352  return;
353  }
354 
355  // apply covariance correction via P_new = (I -K*H)*P
356  // first calculate expression for KHP
357  // then calculate P - KHP
358  float KHP[_k_num_states][_k_num_states];
359  float KH[10];
360 
361  for (unsigned row = 0; row < _k_num_states; row++) {
362 
363  KH[0] = Kfusion[row] * H_MAG[0];
364  KH[1] = Kfusion[row] * H_MAG[1];
365  KH[2] = Kfusion[row] * H_MAG[2];
366  KH[3] = Kfusion[row] * H_MAG[3];
367  KH[4] = Kfusion[row] * H_MAG[16];
368  KH[5] = Kfusion[row] * H_MAG[17];
369  KH[6] = Kfusion[row] * H_MAG[18];
370  KH[7] = Kfusion[row] * H_MAG[19];
371  KH[8] = Kfusion[row] * H_MAG[20];
372  KH[9] = Kfusion[row] * H_MAG[21];
373 
374  for (unsigned column = 0; column < _k_num_states; column++) {
375  float tmp = KH[0] * P[0][column];
376  tmp += KH[1] * P[1][column];
377  tmp += KH[2] * P[2][column];
378  tmp += KH[3] * P[3][column];
379  tmp += KH[4] * P[16][column];
380  tmp += KH[5] * P[17][column];
381  tmp += KH[6] * P[18][column];
382  tmp += KH[7] * P[19][column];
383  tmp += KH[8] * P[20][column];
384  tmp += KH[9] * P[21][column];
385  KHP[row][column] = tmp;
386  }
387  }
388 
389  // if the covariance correction will result in a negative variance, then
390  // the covariance matrix is unhealthy and must be corrected
391  _fault_status.flags.bad_mag_x = false;
392  _fault_status.flags.bad_mag_y = false;
393  _fault_status.flags.bad_mag_z = false;
394 
395  for (int i = 0; i < _k_num_states; i++) {
396  if (P[i][i] < KHP[i][i]) {
397  // zero rows and columns
398  zeroRows(P, i, i);
399  zeroCols(P, i, i);
400 
401  //flag as unhealthy
402  healthy = false;
403 
404  // update individual measurement health status
405  if (index == 0) {
407 
408  } else if (index == 1) {
410 
411  } else if (index == 2) {
413  }
414  }
415  }
416 
417  // only apply covariance and state corrections if healthy
418  if (healthy) {
419  // apply the covariance corrections
420  for (unsigned row = 0; row < _k_num_states; row++) {
421  for (unsigned column = 0; column < _k_num_states; column++) {
422  P[row][column] = P[row][column] - KHP[row][column];
423  }
424  }
425 
426  // correct the covariance matrix for gross errors
428 
429  // apply the state corrections
430  fuse(Kfusion, _mag_innov[index]);
431 
432  // constrain the declination of the earth field states
434  }
435  }
436 }
437 
439 {
440  // assign intermediate state variables
441  float q0 = _state.quat_nominal(0);
442  float q1 = _state.quat_nominal(1);
443  float q2 = _state.quat_nominal(2);
444  float q3 = _state.quat_nominal(3);
445 
446  float R_YAW = 1.0f;
447  float predicted_hdg;
448  float H_YAW[4];
449  Vector3f mag_earth_pred;
450  float measured_hdg;
451 
452  // determine if a 321 or 312 Euler sequence is best
453  if (fabsf(_R_to_earth(2, 0)) < fabsf(_R_to_earth(2, 1))) {
454  // calculate observation jacobian when we are observing the first rotation in a 321 sequence
455  float t9 = q0*q3;
456  float t10 = q1*q2;
457  float t2 = t9+t10;
458  float t3 = q0*q0;
459  float t4 = q1*q1;
460  float t5 = q2*q2;
461  float t6 = q3*q3;
462  float t7 = t3+t4-t5-t6;
463  float t8 = t7*t7;
464  if (t8 > 1e-6f) {
465  t8 = 1.0f/t8;
466  } else {
467  return;
468  }
469  float t11 = t2*t2;
470  float t12 = t8*t11*4.0f;
471  float t13 = t12+1.0f;
472  float t14;
473  if (fabsf(t13) > 1e-6f) {
474  t14 = 1.0f/t13;
475  } else {
476  return;
477  }
478 
479  H_YAW[0] = t8*t14*(q3*t3-q3*t4+q3*t5+q3*t6+q0*q1*q2*2.0f)*-2.0f;
480  H_YAW[1] = t8*t14*(-q2*t3+q2*t4+q2*t5+q2*t6+q0*q1*q3*2.0f)*-2.0f;
481  H_YAW[2] = t8*t14*(q1*t3+q1*t4+q1*t5-q1*t6+q0*q2*q3*2.0f)*2.0f;
482  H_YAW[3] = t8*t14*(q0*t3+q0*t4-q0*t5+q0*t6+q1*q2*q3*2.0f)*2.0f;
483 
484  // rotate the magnetometer measurement into earth frame
485  Eulerf euler321(_state.quat_nominal);
486  predicted_hdg = euler321(2); // we will need the predicted heading to calculate the innovation
487 
488  // calculate the observed yaw angle
490  // Set the yaw angle to zero and rotate the measurements into earth frame using the zero yaw angle
491  euler321(2) = 0.0f;
492  Dcmf R_to_earth(euler321);
493 
494  // rotate the magnetometer measurements into earth frame using a zero yaw angle
496  // don't apply bias corrections if we are learning them
497  mag_earth_pred = R_to_earth * _mag_sample_delayed.mag;
498 
499  } else {
500  mag_earth_pred = R_to_earth * (_mag_sample_delayed.mag - _state.mag_B);
501  }
502 
503  // the angle of the projection onto the horizontal gives the yaw angle
504  measured_hdg = -atan2f(mag_earth_pred(1), mag_earth_pred(0)) + getMagDeclination();
505 
506  } else if (_control_status.flags.ev_yaw) {
507  // calculate the yaw angle for a 321 sequence
508  // Expressions obtained from yaw_input_321.c produced by https://github.com/PX4/ecl/blob/master/matlab/scripts/Inertial%20Nav%20EKF/quat2yaw321.m
511  measured_hdg = atan2f(Tbn_1_0,Tbn_0_0);
512 
513  } else if (_mag_use_inhibit) {
514  // Special case where we either use the current or last known stationary value
515  // so set to current value as a safe default
516  measured_hdg = predicted_hdg;
517 
518  } else {
519  // Should not be doing yaw fusion
520  return;
521 
522  }
523 
524  } else {
525  // calculate observation jacobian when we are observing a rotation in a 312 sequence
526  float t9 = q0*q3;
527  float t10 = q1*q2;
528  float t2 = t9-t10;
529  float t3 = q0*q0;
530  float t4 = q1*q1;
531  float t5 = q2*q2;
532  float t6 = q3*q3;
533  float t7 = t3-t4+t5-t6;
534  float t8 = t7*t7;
535  if (t8 > 1e-6f) {
536  t8 = 1.0f/t8;
537  } else {
538  return;
539  }
540  float t11 = t2*t2;
541  float t12 = t8*t11*4.0f;
542  float t13 = t12+1.0f;
543  float t14;
544 
545  if (fabsf(t13) > 1e-6f) {
546  t14 = 1.0f/t13;
547  } else {
548  return;
549  }
550 
551  H_YAW[0] = t8*t14*(q3*t3+q3*t4-q3*t5+q3*t6-q0*q1*q2*2.0f)*-2.0f;
552  H_YAW[1] = t8*t14*(q2*t3+q2*t4+q2*t5-q2*t6-q0*q1*q3*2.0f)*-2.0f;
553  H_YAW[2] = t8*t14*(-q1*t3+q1*t4+q1*t5+q1*t6-q0*q2*q3*2.0f)*2.0f;
554  H_YAW[3] = t8*t14*(q0*t3-q0*t4+q0*t5+q0*t6-q1*q2*q3*2.0f)*2.0f;
555 
556  /* Calculate the 312 sequence euler angles that rotate from earth to body frame
557  * Derived from https://github.com/PX4/ecl/blob/master/matlab/scripts/Inertial%20Nav%20EKF/quat2yaw312.m
558  * Body to nav frame transformation using a yaw-roll-pitch rotation sequence is given by:
559  *
560  [ cos(pitch)*cos(yaw) - sin(pitch)*sin(roll)*sin(yaw), -cos(roll)*sin(yaw), cos(yaw)*sin(pitch) + cos(pitch)*sin(roll)*sin(yaw)]
561  [ cos(pitch)*sin(yaw) + cos(yaw)*sin(pitch)*sin(roll), cos(roll)*cos(yaw), sin(pitch)*sin(yaw) - cos(pitch)*cos(yaw)*sin(roll)]
562  [ -cos(roll)*sin(pitch), sin(roll), cos(pitch)*cos(roll)]
563  */
564  float yaw = atan2f(-_R_to_earth(0, 1), _R_to_earth(1, 1)); // first rotation (yaw)
565  float roll = asinf(_R_to_earth(2, 1)); // second rotation (roll)
566  float pitch = atan2f(-_R_to_earth(2, 0), _R_to_earth(2, 2)); // third rotation (pitch)
567 
568  predicted_hdg = yaw; // we will need the predicted heading to calculate the innovation
569 
570  // calculate the observed yaw angle
572  // Set the first rotation (yaw) to zero and rotate the measurements into earth frame
573  yaw = 0.0f;
574 
575  // Calculate the body to earth frame rotation matrix from the euler angles using a 312 rotation sequence
576  // Equations from Tbn_312.c produced by https://github.com/PX4/ecl/blob/master/matlab/scripts/Inertial%20Nav%20EKF/quat2yaw312.m
577  Dcmf R_to_earth;
578  float sy = sinf(yaw);
579  float cy = cosf(yaw);
580  float sp = sinf(pitch);
581  float cp = cosf(pitch);
582  float sr = sinf(roll);
583  float cr = cosf(roll);
584  R_to_earth(0,0) = cy*cp-sy*sp*sr;
585  R_to_earth(0,1) = -sy*cr;
586  R_to_earth(0,2) = cy*sp+sy*cp*sr;
587  R_to_earth(1,0) = sy*cp+cy*sp*sr;
588  R_to_earth(1,1) = cy*cr;
589  R_to_earth(1,2) = sy*sp-cy*cp*sr;
590  R_to_earth(2,0) = -sp*cr;
591  R_to_earth(2,1) = sr;
592  R_to_earth(2,2) = cp*cr;
593 
594  // rotate the magnetometer measurements into earth frame using a zero yaw angle
596  // don't apply bias corrections if we are learning them
597  mag_earth_pred = R_to_earth * _mag_sample_delayed.mag;
598  } else {
599  mag_earth_pred = R_to_earth * (_mag_sample_delayed.mag - _state.mag_B);
600  }
601 
602  // the angle of the projection onto the horizontal gives the yaw angle
603  measured_hdg = -atan2f(mag_earth_pred(1), mag_earth_pred(0)) + getMagDeclination();
604 
605  } else if (_control_status.flags.ev_yaw) {
606  // calculate the yaw angle for a 312 sequence
607  // Values from yaw_input_312.c file produced by https://github.com/PX4/ecl/blob/master/matlab/scripts/Inertial%20Nav%20EKF/quat2yaw312.m
610  measured_hdg = atan2f(Tbn_0_1_neg,Tbn_1_1);
611 
612  } else if (_mag_use_inhibit) {
613  // Special case where we either use the current or last known stationary value
614  // so set to current value as a safe default
615  measured_hdg = predicted_hdg;
616 
617  } else {
618  // Should not be doing yaw fusion
619  return;
620 
621  }
622  }
623 
624  // Calculate the observation variance
626  // using magnetic heading tuning parameter
627  R_YAW = sq(fmaxf(_params.mag_heading_noise, 1.0e-2f));
628 
629  } else if (_control_status.flags.ev_yaw) {
630  // using error estimate from external vision data
631  R_YAW = sq(fmaxf(_ev_sample_delayed.angErr, 1.0e-2f));
632 
633  } else {
634  // default value
635  R_YAW = 0.01f;
636  }
637 
638  // wrap the heading to the interval between +-pi
639  measured_hdg = wrap_pi(measured_hdg);
640 
641  // calculate the innovation and define the innovation gate
642  float innov_gate = math::max(_params.heading_innov_gate, 1.0f);
643  if (_mag_use_inhibit) {
644  // The magnetometer cannot be trusted but we need to fuse a heading to prevent a badly
645  // conditioned covariance matrix developing over time.
646  if (!_vehicle_at_rest) {
647  // Vehicle is not at rest so fuse a zero innovation and record the
648  // predicted heading to use as an observation when movement ceases.
649  _heading_innov = 0.0f;
650  _last_static_yaw = predicted_hdg;
651 
652  } else {
653  // Vehicle is at rest so use the last moving prediction as an observation
654  // to prevent the heading from drifting and to enable yaw gyro bias learning
655  // before takeoff.
656  _heading_innov = predicted_hdg - _last_static_yaw;
657  innov_gate = 5.0f;
658 
659  }
660  } else {
661  _heading_innov = predicted_hdg - measured_hdg;
662  _last_static_yaw = predicted_hdg;
663 
664  }
666 
667  // wrap the innovation to the interval between +-pi
669 
670  // Calculate innovation variance and Kalman gains, taking advantage of the fact that only the first 4 elements in H are non zero
671  // calculate the innovation variance
672  float PH[4];
673  _heading_innov_var = R_YAW;
674 
675  for (unsigned row = 0; row <= 3; row++) {
676  PH[row] = 0.0f;
677 
678  for (uint8_t col = 0; col <= 3; col++) {
679  PH[row] += P[row][col] * H_YAW[col];
680  }
681 
682  _heading_innov_var += H_YAW[row] * PH[row];
683  }
684 
685  float heading_innov_var_inv;
686 
687  // check if the innovation variance calculation is badly conditioned
688  if (_heading_innov_var >= R_YAW) {
689  // the innovation variance contribution from the state covariances is not negative, no fault
690  _fault_status.flags.bad_hdg = false;
691  heading_innov_var_inv = 1.0f / _heading_innov_var;
692 
693  } else {
694  // the innovation variance contribution from the state covariances is negative which means the covariance matrix is badly conditioned
695  _fault_status.flags.bad_hdg = true;
696 
697  // we reinitialise the covariance matrix and abort this fusion step
699  ECL_ERR_TIMESTAMPED("EKF mag yaw fusion numerical error - covariance reset");
700  return;
701  }
702 
703  // calculate the Kalman gains
704  // only calculate gains for states we are using
705  float Kfusion[_k_num_states] = {};
706 
707  for (uint8_t row = 0; row <= 15; row++) {
708  Kfusion[row] = 0.0f;
709 
710  for (uint8_t col = 0; col <= 3; col++) {
711  Kfusion[row] += P[row][col] * H_YAW[col];
712  }
713 
714  Kfusion[row] *= heading_innov_var_inv;
715  }
716 
717  if (_control_status.flags.wind) {
718  for (uint8_t row = 22; row <= 23; row++) {
719  Kfusion[row] = 0.0f;
720 
721  for (uint8_t col = 0; col <= 3; col++) {
722  Kfusion[row] += P[row][col] * H_YAW[col];
723  }
724 
725  Kfusion[row] *= heading_innov_var_inv;
726  }
727  }
728 
729  // innovation test ratio
731 
732  // we are no longer using 3-axis fusion so set the reported test levels to zero
733  memset(_mag_test_ratio, 0, sizeof(_mag_test_ratio));
734 
735  // set the magnetometer unhealthy if the test fails
736  if (_yaw_test_ratio > 1.0f) {
738 
739  // if we are in air we don't want to fuse the measurement
740  // we allow to use it when on the ground because the large innovation could be caused
741  // by interference or a large initial gyro bias
743  return;
744 
745  } else {
746  // constrain the innovation to the maximum set by the gate
747  float gate_limit = sqrtf((sq(innov_gate) * _heading_innov_var));
748  _heading_innov = math::constrain(_heading_innov, -gate_limit, gate_limit);
749  }
750 
751  } else {
753  }
754 
755  // apply covariance correction via P_new = (I -K*H)*P
756  // first calculate expression for KHP
757  // then calculate P - KHP
758  float KHP[_k_num_states][_k_num_states];
759  float KH[4];
760 
761  for (unsigned row = 0; row < _k_num_states; row++) {
762 
763  KH[0] = Kfusion[row] * H_YAW[0];
764  KH[1] = Kfusion[row] * H_YAW[1];
765  KH[2] = Kfusion[row] * H_YAW[2];
766  KH[3] = Kfusion[row] * H_YAW[3];
767 
768  for (unsigned column = 0; column < _k_num_states; column++) {
769  float tmp = KH[0] * P[0][column];
770  tmp += KH[1] * P[1][column];
771  tmp += KH[2] * P[2][column];
772  tmp += KH[3] * P[3][column];
773  KHP[row][column] = tmp;
774  }
775  }
776 
777  // if the covariance correction will result in a negative variance, then
778  // the covariance matrix is unhealthy and must be corrected
779  bool healthy = true;
780  _fault_status.flags.bad_hdg = false;
781 
782  for (int i = 0; i < _k_num_states; i++) {
783  if (P[i][i] < KHP[i][i]) {
784  // zero rows and columns
785  zeroRows(P, i, i);
786  zeroCols(P, i, i);
787 
788  //flag as unhealthy
789  healthy = false;
790 
791  // update individual measurement health status
792  _fault_status.flags.bad_hdg = true;
793 
794  }
795  }
796 
797  // only apply covariance and state corrections if healthy
798  if (healthy) {
799  // apply the covariance corrections
800  for (unsigned row = 0; row < _k_num_states; row++) {
801  for (unsigned column = 0; column < _k_num_states; column++) {
802  P[row][column] = P[row][column] - KHP[row][column];
803  }
804  }
805 
806  // correct the covariance matrix for gross errors
808 
809  // apply the state corrections
810  fuse(Kfusion, _heading_innov);
811 
812  }
813 }
814 
815 void Ekf::fuseDeclination(float decl_sigma)
816 {
817  // assign intermediate state variables
818  float magN = _state.mag_I(0);
819  float magE = _state.mag_I(1);
820 
821  // minimum horizontal field strength before calculation becomes badly conditioned (T)
822  float h_field_min = 0.001f;
823 
824  // observation variance (rad**2)
825  float R_DECL = sq(decl_sigma);
826 
827  // Calculate intermediate variables
828  float t2 = magE*magE;
829  float t3 = magN*magN;
830  float t4 = t2+t3;
831  // if the horizontal magnetic field is too small, this calculation will be badly conditioned
832  if (t4 < h_field_min*h_field_min) {
833  return;
834  }
835  float t5 = P[16][16]*t2;
836  float t6 = P[17][17]*t3;
837  float t7 = t2*t2;
838  float t8 = R_DECL*t7;
839  float t9 = t3*t3;
840  float t10 = R_DECL*t9;
841  float t11 = R_DECL*t2*t3*2.0f;
842  float t14 = P[16][17]*magE*magN;
843  float t15 = P[17][16]*magE*magN;
844  float t12 = t5+t6+t8+t10+t11-t14-t15;
845  float t13;
846  if (fabsf(t12) > 1e-6f) {
847  t13 = 1.0f / t12;
848  } else {
849  return;
850  }
851  float t18 = magE*magE;
852  float t19 = magN*magN;
853  float t20 = t18+t19;
854  float t21;
855  if (fabsf(t20) > 1e-6f) {
856  t21 = 1.0f/t20;
857  } else {
858  return;
859  }
860 
861  // Calculate the observation Jacobian
862  // Note only 2 terms are non-zero which can be used in matrix operations for calculation of Kalman gains and covariance update to significantly reduce cost
863  float H_DECL[24] = {};
864  H_DECL[16] = -magE*t21;
865  H_DECL[17] = magN*t21;
866 
867  // Calculate the Kalman gains
868  float Kfusion[_k_num_states] = {};
869  Kfusion[0] = -t4*t13*(P[0][16]*magE-P[0][17]*magN);
870  Kfusion[1] = -t4*t13*(P[1][16]*magE-P[1][17]*magN);
871  Kfusion[2] = -t4*t13*(P[2][16]*magE-P[2][17]*magN);
872  Kfusion[3] = -t4*t13*(P[3][16]*magE-P[3][17]*magN);
873  Kfusion[4] = -t4*t13*(P[4][16]*magE-P[4][17]*magN);
874  Kfusion[5] = -t4*t13*(P[5][16]*magE-P[5][17]*magN);
875  Kfusion[6] = -t4*t13*(P[6][16]*magE-P[6][17]*magN);
876  Kfusion[7] = -t4*t13*(P[7][16]*magE-P[7][17]*magN);
877  Kfusion[8] = -t4*t13*(P[8][16]*magE-P[8][17]*magN);
878  Kfusion[9] = -t4*t13*(P[9][16]*magE-P[9][17]*magN);
879  Kfusion[10] = -t4*t13*(P[10][16]*magE-P[10][17]*magN);
880  Kfusion[11] = -t4*t13*(P[11][16]*magE-P[11][17]*magN);
881  Kfusion[12] = -t4*t13*(P[12][16]*magE-P[12][17]*magN);
882  Kfusion[13] = -t4*t13*(P[13][16]*magE-P[13][17]*magN);
883  Kfusion[14] = -t4*t13*(P[14][16]*magE-P[14][17]*magN);
884  Kfusion[15] = -t4*t13*(P[15][16]*magE-P[15][17]*magN);
885  Kfusion[16] = -t4*t13*(P[16][16]*magE-P[16][17]*magN);
886  Kfusion[17] = -t4*t13*(P[17][16]*magE-P[17][17]*magN);
887  Kfusion[18] = -t4*t13*(P[18][16]*magE-P[18][17]*magN);
888  Kfusion[19] = -t4*t13*(P[19][16]*magE-P[19][17]*magN);
889  Kfusion[20] = -t4*t13*(P[20][16]*magE-P[20][17]*magN);
890  Kfusion[21] = -t4*t13*(P[21][16]*magE-P[21][17]*magN);
891  Kfusion[22] = -t4*t13*(P[22][16]*magE-P[22][17]*magN);
892  Kfusion[23] = -t4*t13*(P[23][16]*magE-P[23][17]*magN);
893 
894  // calculate innovation and constrain
895  float innovation = atan2f(magE, magN) - getMagDeclination();
896  innovation = math::constrain(innovation, -0.5f, 0.5f);
897 
898  // apply covariance correction via P_new = (I -K*H)*P
899  // first calculate expression for KHP
900  // then calculate P - KHP
901  // take advantage of the empty columns in KH to reduce the number of operations
902  float KHP[_k_num_states][_k_num_states];
903  float KH[2];
904  for (unsigned row = 0; row < _k_num_states; row++) {
905 
906  KH[0] = Kfusion[row] * H_DECL[16];
907  KH[1] = Kfusion[row] * H_DECL[17];
908 
909  for (unsigned column = 0; column < _k_num_states; column++) {
910  float tmp = KH[0] * P[16][column];
911  tmp += KH[1] * P[17][column];
912  KHP[row][column] = tmp;
913  }
914  }
915 
916  // if the covariance correction will result in a negative variance, then
917  // the covariance matrix is unhealthy and must be corrected
918  bool healthy = true;
920  for (int i = 0; i < _k_num_states; i++) {
921  if (P[i][i] < KHP[i][i]) {
922  // zero rows and columns
923  zeroRows(P, i, i);
924  zeroCols(P, i, i);
925 
926  //flag as unhealthy
927  healthy = false;
928 
929  // update individual measurement health status
931 
932  }
933  }
934 
935  // only apply covariance and state corrections if healthy
936  if (healthy) {
937  // apply the covariance corrections
938  for (unsigned row = 0; row < _k_num_states; row++) {
939  for (unsigned column = 0; column < _k_num_states; column++) {
940  P[row][column] = P[row][column] - KHP[row][column];
941  }
942  }
943 
944  // correct the covariance matrix for gross errors
946 
947  // apply the state corrections
948  fuse(Kfusion, innovation);
949 
950  // constrain the declination of the earth field states
952  }
953 }
954 
956 {
957  // get a reference value for the earth field declinaton and minimum plausible horizontal field strength
958  // set to 50% of the horizontal strength from geo tables if location is known
959  float decl_reference;
960  float h_field_min = 0.001f;
962  // use parameter value until GPS is available, then use value returned by geo library
964  decl_reference = _mag_declination_gps;
965  h_field_min = fmaxf(h_field_min , 0.5f * _mag_strength_gps * cosf(_mag_inclination_gps));
966  } else {
967  decl_reference = math::radians(_params.mag_declination_deg);
968  }
969  } else {
970  // always use the parameter value
971  decl_reference = math::radians(_params.mag_declination_deg);
972  }
973 
974  // do not allow the horizontal field length to collapse - this will make the declination fusion badly conditioned
975  // and can result in a reversal of the NE field states which the filter cannot recover from
976  // apply a circular limit
977  float h_field = sqrtf(_state.mag_I(0)*_state.mag_I(0) + _state.mag_I(1)*_state.mag_I(1));
978  if (h_field < h_field_min) {
979  if (h_field > 0.001f * h_field_min) {
980  float h_scaler = h_field_min / h_field;
981  _state.mag_I(0) *= h_scaler;
982  _state.mag_I(1) *= h_scaler;
983  } else {
984  // too small to scale radially so set to expected value
985  float mag_declination = getMagDeclination();
986  _state.mag_I(0) = 2.0f * h_field_min * cosf(mag_declination);
987  _state.mag_I(1) = 2.0f * h_field_min * sinf(mag_declination);
988  }
989  h_field = h_field_min;
990  }
991 
992  // do not allow the declination estimate to vary too much relative to the reference value
993  const float decl_tolerance = 0.5f;
994  const float decl_max = decl_reference + decl_tolerance;
995  const float decl_min = decl_reference - decl_tolerance;
996  const float decl_estimate = atan2f(_state.mag_I(1) , _state.mag_I(0));
997  if (decl_estimate > decl_max) {
998  _state.mag_I(0) = h_field * cosf(decl_max);
999  _state.mag_I(1) = h_field * sinf(decl_max);
1000  } else if (decl_estimate < decl_min) {
1001  _state.mag_I(0) = h_field * cosf(decl_min);
1002  _state.mag_I(1) = h_field * sinf(decl_min);
1003  }
1004 }
1005 
1007 {
1008  // theoretical magnitude of the magnetometer Z component value given X and Y sensor measurement and our knowledge
1009  // of the earth magnetic field vector at the current location
1010  float mag_z_abs = sqrtf(math::max(sq(mag_earth_predicted.length()) - sq(mag_meas(0)) - sq(mag_meas(1)), 0.0f));
1011 
1012  // calculate sign of synthetic magnetomter Z component based on the sign of the predicted magnetomer Z component
1013  float mag_z_body_pred = _R_to_earth(0,2) * mag_earth_predicted(0) + _R_to_earth(1,2) * mag_earth_predicted(1) + _R_to_earth(2,2) * mag_earth_predicted(2);
1014 
1015  return mag_z_body_pred < 0 ? -mag_z_abs : mag_z_abs;
1016 }
Matrix3f quat_to_invrotmat(const Quatf &quat)
void fuseDeclination(float decl_sigma)
Definition: mag_fusion.cpp:815
constexpr _Tp constrain(_Tp val, _Tp min_val, _Tp max_val)
Definition: Limits.hpp:66
struct estimator::filter_control_status_u::@60 flags
t8
Definition: calcH_YAW312.c:9
Adapter / shim layer for system calls needed by ECL.
stateSample _state
state struct of the ekf running at the delayed time horizon
Definition: ekf.h:288
t12
Definition: calcH_YAW312.c:11
#define ECL_ERR_TIMESTAMPED
Definition: ecl.h:95
float mag_heading_noise
measurement noise used for simple heading fusion (rad)
Definition: common.h:269
t9
Definition: calcH_YAW312.c:1
bool bad_mag_x
0 - true if the fusion of the magnetometer X-axis has encountered a numerical error ...
Definition: common.h:380
bool bad_mag_y
1 - true if the fusion of the magnetometer Y-axis has encountered a numerical error ...
Definition: common.h:381
t4
Definition: calcH_YAW312.c:5
Dcm< float > Dcmf
Definition: Dcm.hpp:185
float _last_static_yaw
last yaw angle recorded when on ground motion checks were passing (rad)
Definition: ekf.h:362
bool bad_mag_decl
4 - true if the fusion of the magnetic declination has encountered a numerical error ...
Definition: common.h:384
Quatf quat
quaternion defining rotation from body to earth frame
Definition: common.h:158
float _heading_innov_var
heading measurement innovation variance (rad**2)
Definition: ekf.h:389
uint32_t mag_hdg
4 - true if a simple magnetic yaw heading is being fused
Definition: common.h:445
t6
Definition: calcH_YAW312.c:7
t5
Definition: calcH_YAW312.c:6
t14
Definition: calcH_YAW312.c:13
void initialiseCovariance()
Definition: covariance.cpp:49
Vector3f mag_I
NED earth magnetic field in gauss.
Definition: common.h:373
t19
Definition: calcH_YAWGPS.c:17
struct estimator::fault_status_u::@57 flags
float _mag_innov[3]
earth magnetic field innovations (Gauss)
Definition: ekf.h:376
t21
Definition: calcH_YAWGPS.c:20
float P[_k_num_states][_k_num_states]
state covariance matrix
Definition: ekf.h:367
void zeroCols(float(&cov_mat)[_k_num_states][_k_num_states], uint8_t first, uint8_t last)
float _heading_innov
heading measurement innovation (rad)
Definition: ekf.h:388
#define MASK_USE_GEO_DECL
set to true to use the declination from the geo library when the GPS position becomes available...
Definition: common.h:184
void resetMagRelatedCovariances()
Definition: covariance.cpp:899
float getMagDeclination()
Definition: ekf_helper.cpp:750
uint32_t synthetic_mag_z
25 - true when we are using a synthesized measurement for the magnetometer Z component ...
Definition: common.h:466
void fuse(float *K, float innovation)
t10
Definition: calcH_YAW312.c:2
Dcmf _R_to_earth
transformation matrix from body frame to earth frame from last EKF prediction
Definition: ekf.h:348
bool bad_mag_z
2 - true if the fusion of the magnetometer Z-axis has encountered a numerical error ...
Definition: common.h:382
t15
Definition: calcH_YAWGPS.c:14
float calculate_synthetic_mag_z_measurement(Vector3f mag_meas, Vector3f mag_earth_predicted)
int32_t mag_declination_source
bitmask used to control the handling of declination data
Definition: common.h:274
bool _mag_use_inhibit_prev
true when magnetometer use was being inhibited the previous frame
Definition: ekf.h:360
float mag_noise
measurement noise used for 3-axis magnetoemeter fusion (Gauss)
Definition: common.h:270
bool _flt_mag_align_converging
true when the in-flight mag field post alignment convergence is being performd
Definition: ekf.h:440
t20
Definition: calcH_YAWGPS.c:18
void limitDeclination()
Definition: mag_fusion.cpp:955
constexpr T radians(T degrees)
Definition: Limits.hpp:85
Vector< float, 6 > f(float t, const Matrix< float, 6, 1 > &, const Matrix< float, 3, 1 > &)
Definition: integration.cpp:8
bool bad_hdg
3 - true if the fusion of the heading angle has encountered a numerical error
Definition: common.h:383
t13
Definition: calcH_YAW312.c:12
filter_control_status_u _control_status
Euler< float > Eulerf
Definition: Euler.hpp:156
static constexpr float sq(float var)
Definition: ekf.h:695
float _mag_innov_var[3]
earth magnetic field innovation variance (Gauss**2)
Definition: ekf.h:377
uint32_t ev_yaw
13 - true when yaw data from external vision measurements is being fused
Definition: common.h:454
fault_status_u _fault_status
void fixCovarianceErrors()
Definition: covariance.cpp:728
bool reject_yaw
6 - true if the yaw observation has been rejected
Definition: common.h:410
Type wrap_pi(Type x)
Wrap value in range [-π, π)
Quatf quat_nominal
quaternion defining the rotation from body to earth frame
Definition: common.h:368
void fuseMag()
Definition: mag_fusion.cpp:47
t11
Definition: calcH_YAW312.c:10
float mag_innov_gate
magnetometer fusion innovation consistency gate size (STD)
Definition: common.h:273
t18
Definition: calcH_YAWGPS.c:19
bool _mag_use_inhibit
true when magnetometer use is being inhibited
Definition: ekf.h:359
innovation_fault_status_u _innov_check_fail_status
Vector3< float > Vector3f
Definition: Vector3.hpp:136
extVisionSample _ev_sample_delayed
Kfusion[0]
Definition: K_VELX.c:93
Vector3f mag_B
magnetometer bias estimate in body frame in gauss
Definition: common.h:374
constexpr _Tp max(_Tp a, _Tp b)
Definition: Limits.hpp:60
uint32_t mag_3D
5 - true if 3-axis magnetometer measurement are being fused
Definition: common.h:446
void zeroRows(float(&cov_mat)[_k_num_states][_k_num_states], uint8_t first, uint8_t last)
uint32_t in_air
7 - true when the vehicle is airborne
Definition: common.h:448
float mag_declination_deg
magnetic declination (degrees)
Definition: common.h:271
t2
Definition: calcH_YAW312.c:3
uint32_t wind
8 - true when wind velocity is being estimated
Definition: common.h:449
float angErr
1-Sigma angular error (rad)
Definition: common.h:162
void fuseHeading()
Definition: mag_fusion.cpp:438
struct estimator::innovation_fault_status_u::@58 flags
t7
Definition: calcH_YAW312.c:8
Vector3f mag
NED magnetometer body frame measurements (Gauss)
Definition: common.h:126
static constexpr uint8_t _k_num_states
number of EKF states
Definition: ekf.h:270
float heading_innov_gate
heading fusion innovation consistency gate size (STD)
Definition: common.h:272
Class for core functions for ekf attitude and position estimator.
t3
Definition: calcH_YAW312.c:4