HAPSRG
Eigen
Dependencies: Eigen
Dependents: optWingforHAPS_Eigen hexaTest_Eigen
solaESKF.cpp@77:780ce6556041, 2021-11-30 (annotated)
- Committer:
- NaotoMorita
- Date:
- Tue Nov 30 11:26:56 2021 +0000
- Revision:
- 77:780ce6556041
- Parent:
- 75:e2c825cdc511
- Child:
- 78:e36a7b844fb5
teleme usapack
Who changed what in which revision?
| User | Revision | Line number | New contents of line |
|---|---|---|---|
| NaotoMorita | 44:7d82e63b6a86 | 1 | #include "solaESKF.hpp" |
| NaotoMorita | 54:cd514d9d4b19 | 2 | |
| NaotoMorita | 19:3fae66745363 | 3 | |
| NaotoMorita | 19:3fae66745363 | 4 | |
| NaotoMorita | 44:7d82e63b6a86 | 5 | solaESKF::solaESKF() |
| NaotoMorita | 74:f5fe7fecbd3c | 6 | :pihat(3,1),vihat(3,1),qhat(4,1),accBias(3,1),gyroBias(3,1),gravity(3,1),errState(18,1),Phat(18,18),Q(18,18) |
| NaotoMorita | 62:5519d34eb6e8 | 7 | { |
| NaotoMorita | 46:15988dc41923 | 8 | pihat << 0.0f << 0.0f << 0.0f; |
| NaotoMorita | 46:15988dc41923 | 9 | vihat << 0.0f << 0.0f << 0.0f; |
| NaotoMorita | 46:15988dc41923 | 10 | qhat << 1.0f << 0.0f << 0.0f << 0.0f; |
| NaotoMorita | 46:15988dc41923 | 11 | accBias << 0.0f << 0.0f << 0.0f; |
| NaotoMorita | 46:15988dc41923 | 12 | gyroBias << 0.0f << 0.0f << 0.0f; |
| NaotoMorita | 46:15988dc41923 | 13 | gravity << 0.0f << 0.0f << 9.8f; |
| NaotoMorita | 58:93ba28cf5cb3 | 14 | |
| NaotoMorita | 62:5519d34eb6e8 | 15 | nState = errState.getRows(); |
| NaotoMorita | 62:5519d34eb6e8 | 16 | |
| NaotoMorita | 58:93ba28cf5cb3 | 17 | for (int i = 1; i < nState+1; i++){ |
| NaotoMorita | 46:15988dc41923 | 18 | errState(i,1) = 0.0f; |
| NaotoMorita | 58:93ba28cf5cb3 | 19 | for (int j = 1; j < nState+1; j++){ |
| NaotoMorita | 46:15988dc41923 | 20 | Phat(i,j) = 0.0f; |
| NaotoMorita | 46:15988dc41923 | 21 | Q(i,j) = 0.0f; |
| NaotoMorita | 46:15988dc41923 | 22 | } |
| NaotoMorita | 46:15988dc41923 | 23 | } |
| NaotoMorita | 62:5519d34eb6e8 | 24 | |
| NaotoMorita | 46:15988dc41923 | 25 | setBlockDiag(Phat,0.1f,1,3);//position |
| NaotoMorita | 46:15988dc41923 | 26 | setBlockDiag(Phat,0.1f,4,6);//velocity |
| NaotoMorita | 46:15988dc41923 | 27 | setBlockDiag(Phat,0.1f,7,9);//angle error |
| NaotoMorita | 47:2467de40951f | 28 | setBlockDiag(Phat,0.1f,10,12);//acc bias |
| NaotoMorita | 47:2467de40951f | 29 | setBlockDiag(Phat,0.1f,13,15);//gyro bias |
| NaotoMorita | 74:f5fe7fecbd3c | 30 | setBlockDiag(Phat,0.00000001f,16,18);//gravity |
| NaotoMorita | 47:2467de40951f | 31 | setBlockDiag(Q,0.00025f,4,6);//velocity |
| NaotoMorita | 55:21611d4cf7e8 | 32 | setBlockDiag(Q,0.005f/57.0f,7,9);//angle error |
| NaotoMorita | 46:15988dc41923 | 33 | setBlockDiag(Q,0.001f,10,12);//acc bias |
| NaotoMorita | 70:d12e46fdc2f0 | 34 | setBlockDiag(Q,0.001f,13,15);//gyro bias//positionとgravityはQ項なし |
| NaotoMorita | 19:3fae66745363 | 35 | |
| NaotoMorita | 19:3fae66745363 | 36 | } |
| NaotoMorita | 19:3fae66745363 | 37 | |
| NaotoMorita | 47:2467de40951f | 38 | |
| NaotoMorita | 44:7d82e63b6a86 | 39 | void solaESKF::updateNominal(Matrix acc, Matrix gyro,float att_dt) |
| NaotoMorita | 19:3fae66745363 | 40 | { |
| NaotoMorita | 44:7d82e63b6a86 | 41 | Matrix gyrom = gyro - gyroBias; |
| NaotoMorita | 44:7d82e63b6a86 | 42 | Matrix accm = acc - accBias; |
| NaotoMorita | 44:7d82e63b6a86 | 43 | |
| NaotoMorita | 44:7d82e63b6a86 | 44 | Matrix qint(4,1); |
| NaotoMorita | 46:15988dc41923 | 45 | qint << 1.0f << 0.5f*gyrom(1,1)*att_dt << 0.5f*gyrom(2,1)*att_dt << 0.5f*gyrom(3,1)*att_dt; |
| NaotoMorita | 44:7d82e63b6a86 | 46 | qhat = quatmultiply(qhat,qint); |
| NaotoMorita | 19:3fae66745363 | 47 | float qnorm = sqrt(MatrixMath::dot(MatrixMath::Transpose(qhat),qhat)); |
| NaotoMorita | 19:3fae66745363 | 48 | qhat *= (1.0f/ qnorm); |
| NaotoMorita | 23:1509648c2318 | 49 | |
| NaotoMorita | 23:1509648c2318 | 50 | Matrix dcm(3,3); |
| NaotoMorita | 23:1509648c2318 | 51 | computeDcm(dcm, qhat); |
| NaotoMorita | 46:15988dc41923 | 52 | Matrix accned = dcm*accm+gravity; |
| NaotoMorita | 44:7d82e63b6a86 | 53 | vihat += accned*att_dt; |
| NaotoMorita | 44:7d82e63b6a86 | 54 | |
| NaotoMorita | 44:7d82e63b6a86 | 55 | pihat += vihat*att_dt+0.5f*accned*att_dt*att_dt; |
| NaotoMorita | 19:3fae66745363 | 56 | } |
| NaotoMorita | 19:3fae66745363 | 57 | |
| NaotoMorita | 44:7d82e63b6a86 | 58 | void solaESKF::updateErrState(Matrix acc,Matrix gyro,float att_dt) |
| NaotoMorita | 23:1509648c2318 | 59 | { |
| NaotoMorita | 44:7d82e63b6a86 | 60 | Matrix gyrom = gyro - gyroBias; |
| NaotoMorita | 44:7d82e63b6a86 | 61 | Matrix accm = acc - accBias; |
| NaotoMorita | 75:e2c825cdc511 | 62 | |
| NaotoMorita | 23:1509648c2318 | 63 | Matrix dcm(3,3); |
| NaotoMorita | 23:1509648c2318 | 64 | computeDcm(dcm, qhat); |
| NaotoMorita | 75:e2c825cdc511 | 65 | Matrix a2v = -dcm*MatrixMath::Matrixcross(accm(1,1),accm(2,1),accm(3,1))*att_dt; |
| NaotoMorita | 75:e2c825cdc511 | 66 | Matrix a2v2 = 0.5f*a2v*att_dt; |
| NaotoMorita | 75:e2c825cdc511 | 67 | |
| NaotoMorita | 75:e2c825cdc511 | 68 | Matrix Fx(nState,nState); |
| NaotoMorita | 44:7d82e63b6a86 | 69 | //position |
| NaotoMorita | 75:e2c825cdc511 | 70 | Fx(1,1) = 1.0f; |
| NaotoMorita | 75:e2c825cdc511 | 71 | Fx(2,2) = 1.0f; |
| NaotoMorita | 75:e2c825cdc511 | 72 | Fx(3,3) = 1.0f; |
| NaotoMorita | 75:e2c825cdc511 | 73 | Fx(1,4) = 1.0f*att_dt; |
| NaotoMorita | 75:e2c825cdc511 | 74 | Fx(2,5) = 1.0f*att_dt; |
| NaotoMorita | 75:e2c825cdc511 | 75 | Fx(3,6) = 1.0f*att_dt; |
| NaotoMorita | 46:15988dc41923 | 76 | for (int i = 1; i < 4; i++){ |
| NaotoMorita | 46:15988dc41923 | 77 | for (int j = 1; j < 4; j++){ |
| NaotoMorita | 75:e2c825cdc511 | 78 | Fx(i,j+6) = a2v2(i,j); |
| NaotoMorita | 75:e2c825cdc511 | 79 | Fx(i,j+9) = -0.5f*dcm(i,j)*att_dt*att_dt; |
| NaotoMorita | 75:e2c825cdc511 | 80 | } |
| NaotoMorita | 75:e2c825cdc511 | 81 | Fx(i,i+15) = 0.5f*att_dt*att_dt; |
| NaotoMorita | 75:e2c825cdc511 | 82 | } |
| NaotoMorita | 75:e2c825cdc511 | 83 | |
| NaotoMorita | 75:e2c825cdc511 | 84 | |
| NaotoMorita | 75:e2c825cdc511 | 85 | //velocity |
| NaotoMorita | 75:e2c825cdc511 | 86 | Fx(4,4) = 1.0f; |
| NaotoMorita | 75:e2c825cdc511 | 87 | Fx(5,5) = 1.0f; |
| NaotoMorita | 75:e2c825cdc511 | 88 | Fx(6,6) = 1.0f; |
| NaotoMorita | 75:e2c825cdc511 | 89 | for (int i = 1; i < 4; i++){ |
| NaotoMorita | 75:e2c825cdc511 | 90 | for (int j = 1; j < 4; j++){ |
| NaotoMorita | 75:e2c825cdc511 | 91 | Fx(i+3,j+6) = a2v(i,j); |
| NaotoMorita | 75:e2c825cdc511 | 92 | Fx(i+3,j+9) = -dcm(i,j)*att_dt; |
| NaotoMorita | 75:e2c825cdc511 | 93 | Fx(i+3,j+12) = -a2v2(i,j); |
| NaotoMorita | 23:1509648c2318 | 94 | } |
| NaotoMorita | 23:1509648c2318 | 95 | } |
| NaotoMorita | 75:e2c825cdc511 | 96 | Fx(4,16) = 1.0f*att_dt; |
| NaotoMorita | 75:e2c825cdc511 | 97 | Fx(5,17) = 1.0f*att_dt; |
| NaotoMorita | 75:e2c825cdc511 | 98 | Fx(6,18) = 1.0f*att_dt; |
| NaotoMorita | 75:e2c825cdc511 | 99 | |
| NaotoMorita | 44:7d82e63b6a86 | 100 | //angulat error |
| NaotoMorita | 75:e2c825cdc511 | 101 | Fx(7,7) = 1.0f; |
| NaotoMorita | 75:e2c825cdc511 | 102 | Fx(8,8) = 1.0f; |
| NaotoMorita | 75:e2c825cdc511 | 103 | Fx(9,9) = 1.0f; |
| NaotoMorita | 75:e2c825cdc511 | 104 | Fx(7,8) = gyrom(3,1)*att_dt; |
| NaotoMorita | 75:e2c825cdc511 | 105 | Fx(7,9) = -gyrom(2,1)*att_dt; |
| NaotoMorita | 75:e2c825cdc511 | 106 | Fx(8,7) = -gyrom(3,1)*att_dt; |
| NaotoMorita | 75:e2c825cdc511 | 107 | Fx(8,9) = gyrom(1,1)*att_dt; |
| NaotoMorita | 75:e2c825cdc511 | 108 | Fx(9,7) = gyrom(2,1)*att_dt; |
| NaotoMorita | 75:e2c825cdc511 | 109 | Fx(9,8) = -gyrom(1,1)*att_dt; |
| NaotoMorita | 75:e2c825cdc511 | 110 | Fx(7,13) = -1.0f*att_dt; |
| NaotoMorita | 75:e2c825cdc511 | 111 | Fx(8,14) = -1.0f*att_dt; |
| NaotoMorita | 75:e2c825cdc511 | 112 | Fx(9,15) = -1.0f*att_dt; |
| NaotoMorita | 47:2467de40951f | 113 | |
| NaotoMorita | 75:e2c825cdc511 | 114 | //acc bias |
| NaotoMorita | 75:e2c825cdc511 | 115 | Fx(10,10) = 1.0f; |
| NaotoMorita | 75:e2c825cdc511 | 116 | Fx(11,11) = 1.0f; |
| NaotoMorita | 75:e2c825cdc511 | 117 | Fx(12,12) = 1.0f; |
| NaotoMorita | 75:e2c825cdc511 | 118 | |
| NaotoMorita | 75:e2c825cdc511 | 119 | //gyro bias |
| NaotoMorita | 75:e2c825cdc511 | 120 | Fx(13,13) = 1.0f; |
| NaotoMorita | 75:e2c825cdc511 | 121 | Fx(14,14) = 1.0f; |
| NaotoMorita | 75:e2c825cdc511 | 122 | Fx(15,15) = 1.0f; |
| NaotoMorita | 75:e2c825cdc511 | 123 | |
| NaotoMorita | 75:e2c825cdc511 | 124 | //gravity bias |
| NaotoMorita | 75:e2c825cdc511 | 125 | Fx(16,16) = 1.0f; |
| NaotoMorita | 75:e2c825cdc511 | 126 | Fx(17,17) = 1.0f; |
| NaotoMorita | 75:e2c825cdc511 | 127 | Fx(18,18) = 1.0f; |
| NaotoMorita | 75:e2c825cdc511 | 128 | |
| NaotoMorita | 75:e2c825cdc511 | 129 | errState = Fx * errState; |
| NaotoMorita | 75:e2c825cdc511 | 130 | Phat = Fx*Phat*MatrixMath::Transpose(Fx); |
| NaotoMorita | 59:03fe5e16a33c | 131 | for (int i = 1; i < nState+1; i++){ |
| NaotoMorita | 59:03fe5e16a33c | 132 | if(i>3 && i<10){ |
| NaotoMorita | 75:e2c825cdc511 | 133 | Phat(i,i) += Q(i,i)*att_dt; |
| NaotoMorita | 59:03fe5e16a33c | 134 | }else if(i>9 && i<16){ |
| NaotoMorita | 75:e2c825cdc511 | 135 | Phat(i,i) += Q(i,i)* att_dt*att_dt; |
| NaotoMorita | 59:03fe5e16a33c | 136 | } |
| NaotoMorita | 59:03fe5e16a33c | 137 | } |
| NaotoMorita | 25:07ac5c6cd61c | 138 | } |
| NaotoMorita | 25:07ac5c6cd61c | 139 | |
| NaotoMorita | 61:5e5c4fe12440 | 140 | void solaESKF::updateAcc(Matrix acc,Matrix R) |
| NaotoMorita | 46:15988dc41923 | 141 | { |
| NaotoMorita | 46:15988dc41923 | 142 | Matrix accm = acc - accBias; |
| NaotoMorita | 56:c10f1168bd4a | 143 | Matrix dcm(3,3); |
| NaotoMorita | 56:c10f1168bd4a | 144 | computeDcm(dcm, qhat); |
| NaotoMorita | 56:c10f1168bd4a | 145 | Matrix tdcm = MatrixMath::Transpose(dcm); |
| NaotoMorita | 46:15988dc41923 | 146 | Matrix tdcm_g = tdcm*gravity; |
| NaotoMorita | 59:03fe5e16a33c | 147 | Matrix rotgrav = MatrixMath::Matrixcross(tdcm_g(1,1),tdcm_g(2,1),tdcm_g(3,1)); |
| NaotoMorita | 46:15988dc41923 | 148 | |
| NaotoMorita | 61:5e5c4fe12440 | 149 | Matrix H(3,nState); |
| NaotoMorita | 59:03fe5e16a33c | 150 | for (int i = 1; i < 4; i++){ |
| NaotoMorita | 59:03fe5e16a33c | 151 | for (int j = 1; j < 4; j++){ |
| NaotoMorita | 59:03fe5e16a33c | 152 | H(i,j+6) = rotgrav(i,j); |
| NaotoMorita | 74:f5fe7fecbd3c | 153 | H(i,j+15) = tdcm(i,j); |
| NaotoMorita | 46:15988dc41923 | 154 | } |
| NaotoMorita | 58:93ba28cf5cb3 | 155 | } |
| NaotoMorita | 56:c10f1168bd4a | 156 | |
| NaotoMorita | 46:15988dc41923 | 157 | H(1,10) = -1.0f; |
| NaotoMorita | 46:15988dc41923 | 158 | H(2,11) = -1.0f; |
| NaotoMorita | 46:15988dc41923 | 159 | H(3,12) = -1.0f; |
| NaotoMorita | 59:03fe5e16a33c | 160 | |
| NaotoMorita | 61:5e5c4fe12440 | 161 | Matrix K = (Phat*MatrixMath::Transpose(H))*MatrixMath::Inv(H*Phat*MatrixMath::Transpose(H)+R); |
| NaotoMorita | 61:5e5c4fe12440 | 162 | Matrix zacc = -accm-tdcm*gravity; |
| NaotoMorita | 61:5e5c4fe12440 | 163 | Matrix z(3,1); |
| NaotoMorita | 61:5e5c4fe12440 | 164 | z << zacc(1,1) << zacc(2,1) << zacc(3,1); |
| NaotoMorita | 61:5e5c4fe12440 | 165 | errState = K * z; |
| NaotoMorita | 61:5e5c4fe12440 | 166 | Phat = (MatrixMath::Eye(nState)-K*H)*Phat; |
| NaotoMorita | 61:5e5c4fe12440 | 167 | |
| NaotoMorita | 61:5e5c4fe12440 | 168 | fuseErr2Nominal(); |
| NaotoMorita | 61:5e5c4fe12440 | 169 | } |
| NaotoMorita | 70:d12e46fdc2f0 | 170 | |
| NaotoMorita | 73:5770a0d470c0 | 171 | void solaESKF::updateHeading(float a,Matrix R) |
| NaotoMorita | 70:d12e46fdc2f0 | 172 | { |
| NaotoMorita | 70:d12e46fdc2f0 | 173 | float q0 = qhat(1,1); |
| NaotoMorita | 70:d12e46fdc2f0 | 174 | float q1 = qhat(2,1); |
| NaotoMorita | 70:d12e46fdc2f0 | 175 | float q2 = qhat(3,1); |
| NaotoMorita | 70:d12e46fdc2f0 | 176 | float q3 = qhat(4,1); |
| NaotoMorita | 73:5770a0d470c0 | 177 | |
| NaotoMorita | 73:5770a0d470c0 | 178 | bool canUseA = false; |
| NaotoMorita | 73:5770a0d470c0 | 179 | const float SA0 = 2*q3; |
| NaotoMorita | 73:5770a0d470c0 | 180 | const float SA1 = 2*q2; |
| NaotoMorita | 73:5770a0d470c0 | 181 | const float SA2 = SA0*q0 + SA1*q1; |
| NaotoMorita | 73:5770a0d470c0 | 182 | const float SA3 = q0*q0 + q1*q1 - q2*q2 - q3*q3; |
| NaotoMorita | 73:5770a0d470c0 | 183 | float SA4, SA5_inv; |
| NaotoMorita | 73:5770a0d470c0 | 184 | if ((SA3*SA3) > 1e-6f) { |
| NaotoMorita | 73:5770a0d470c0 | 185 | SA4 = 1.0F/(SA3*SA3); |
| NaotoMorita | 73:5770a0d470c0 | 186 | SA5_inv = SA2*SA2*SA4 + 1; |
| NaotoMorita | 73:5770a0d470c0 | 187 | canUseA = fabsf(SA5_inv) > 1e-6f; |
| NaotoMorita | 73:5770a0d470c0 | 188 | } |
| NaotoMorita | 73:5770a0d470c0 | 189 | |
| NaotoMorita | 73:5770a0d470c0 | 190 | bool canUseB = false; |
| NaotoMorita | 73:5770a0d470c0 | 191 | const float SB0 = 2*q0; |
| NaotoMorita | 73:5770a0d470c0 | 192 | const float SB1 = 2*q1; |
| NaotoMorita | 73:5770a0d470c0 | 193 | const float SB2 = SB0*q3 + SB1*q2; |
| NaotoMorita | 73:5770a0d470c0 | 194 | const float SB4 = q0*q0 + q1*q1 - q2*q2 - q3*q3; |
| NaotoMorita | 73:5770a0d470c0 | 195 | float SB3, SB5_inv; |
| NaotoMorita | 73:5770a0d470c0 | 196 | if ((SB2*SB2) > 1e-6f) { |
| NaotoMorita | 73:5770a0d470c0 | 197 | SB3 = 1.0F/(SB2*SB2); |
| NaotoMorita | 73:5770a0d470c0 | 198 | SB5_inv = SB3*SB4*SB4 + 1; |
| NaotoMorita | 73:5770a0d470c0 | 199 | canUseB = fabsf(SB5_inv) > 1e-6f; |
| NaotoMorita | 73:5770a0d470c0 | 200 | } |
| NaotoMorita | 73:5770a0d470c0 | 201 | |
| NaotoMorita | 73:5770a0d470c0 | 202 | Matrix Hh(1,4); |
| NaotoMorita | 73:5770a0d470c0 | 203 | |
| NaotoMorita | 73:5770a0d470c0 | 204 | if (canUseA && (!canUseB || fabsf(SA5_inv) >= fabsf(SB5_inv))) { |
| NaotoMorita | 73:5770a0d470c0 | 205 | const float SA5 = 1.0F/SA5_inv; |
| NaotoMorita | 73:5770a0d470c0 | 206 | const float SA6 = 1.0F/SA3; |
| NaotoMorita | 73:5770a0d470c0 | 207 | const float SA7 = SA2*SA4; |
| NaotoMorita | 73:5770a0d470c0 | 208 | const float SA8 = 2*SA7; |
| NaotoMorita | 73:5770a0d470c0 | 209 | const float SA9 = 2*SA6; |
| NaotoMorita | 73:5770a0d470c0 | 210 | |
| NaotoMorita | 73:5770a0d470c0 | 211 | Hh(1,1) = SA5*(SA0*SA6 - SA8*q0); |
| NaotoMorita | 73:5770a0d470c0 | 212 | Hh(1,2) = SA5*(SA1*SA6 - SA8*q1); |
| NaotoMorita | 73:5770a0d470c0 | 213 | Hh(1,3) = SA5*(SA1*SA7 + SA9*q1); |
| NaotoMorita | 73:5770a0d470c0 | 214 | Hh(1,4) = SA5*(SA0*SA7 + SA9*q0); |
| NaotoMorita | 73:5770a0d470c0 | 215 | } else if (canUseB && (!canUseA || fabsf(SB5_inv) > fabsf(SA5_inv))) { |
| NaotoMorita | 73:5770a0d470c0 | 216 | const float SB5 = 1.0F/SB5_inv; |
| NaotoMorita | 73:5770a0d470c0 | 217 | const float SB6 = 1.0F/SB2; |
| NaotoMorita | 73:5770a0d470c0 | 218 | const float SB7 = SB3*SB4; |
| NaotoMorita | 73:5770a0d470c0 | 219 | const float SB8 = 2*SB7; |
| NaotoMorita | 73:5770a0d470c0 | 220 | const float SB9 = 2*SB6; |
| NaotoMorita | 73:5770a0d470c0 | 221 | |
| NaotoMorita | 73:5770a0d470c0 | 222 | Hh(1,1) = -SB5*(SB0*SB6 - SB8*q3); |
| NaotoMorita | 73:5770a0d470c0 | 223 | Hh(1,2) = -SB5*(SB1*SB6 - SB8*q2); |
| NaotoMorita | 73:5770a0d470c0 | 224 | Hh(1,3) = -SB5*(-SB1*SB7 - SB9*q2); |
| NaotoMorita | 73:5770a0d470c0 | 225 | Hh(1,4) = -SB5*(-SB0*SB7 - SB9*q3); |
| NaotoMorita | 73:5770a0d470c0 | 226 | } else { |
| NaotoMorita | 73:5770a0d470c0 | 227 | return; |
| NaotoMorita | 73:5770a0d470c0 | 228 | } |
| NaotoMorita | 70:d12e46fdc2f0 | 229 | |
| NaotoMorita | 70:d12e46fdc2f0 | 230 | Matrix Hdq(4,3); |
| NaotoMorita | 70:d12e46fdc2f0 | 231 | Hdq << -0.5f*q1 << -0.5f*q2 << -0.5f*q3 |
| NaotoMorita | 72:8ebae608ae12 | 232 | << 0.5f*q0 << -0.5f*q3 << 0.5f*q2 |
| NaotoMorita | 72:8ebae608ae12 | 233 | << 0.5f*q3 << 0.5f*q0 << -0.5f*q1 |
| NaotoMorita | 73:5770a0d470c0 | 234 | << -0.5f*q2 << 0.5f*q1 << 0.5f*q0; |
| NaotoMorita | 70:d12e46fdc2f0 | 235 | |
| NaotoMorita | 70:d12e46fdc2f0 | 236 | Matrix Hpart = Hh*Hdq; |
| NaotoMorita | 73:5770a0d470c0 | 237 | Matrix H(1,nState); |
| NaotoMorita | 70:d12e46fdc2f0 | 238 | for(int j=1;j<4;j++){ |
| NaotoMorita | 70:d12e46fdc2f0 | 239 | H(1,j+6) = Hpart(1,j); |
| NaotoMorita | 70:d12e46fdc2f0 | 240 | } |
| NaotoMorita | 70:d12e46fdc2f0 | 241 | |
| NaotoMorita | 74:f5fe7fecbd3c | 242 | const float psi = atan2f(qhat(2,1)*qhat(3,1) + qhat(1,1)*qhat(4,1), 0.5f - qhat(3,1)*qhat(3,1) - qhat(4,1)*qhat(4,1)); |
| NaotoMorita | 73:5770a0d470c0 | 243 | Matrix z(1,1); |
| NaotoMorita | 74:f5fe7fecbd3c | 244 | z << atan2f(sin(a-psi),cos(a-psi)); |
| NaotoMorita | 70:d12e46fdc2f0 | 245 | Matrix K = (Phat*MatrixMath::Transpose(H))*MatrixMath::Inv(H*Phat*MatrixMath::Transpose(H)+R); |
| NaotoMorita | 70:d12e46fdc2f0 | 246 | errState = K * z; |
| NaotoMorita | 70:d12e46fdc2f0 | 247 | Phat = (MatrixMath::Eye(nState)-K*H)*Phat; |
| NaotoMorita | 70:d12e46fdc2f0 | 248 | |
| NaotoMorita | 70:d12e46fdc2f0 | 249 | fuseErr2Nominal(); |
| NaotoMorita | 70:d12e46fdc2f0 | 250 | } |
| NaotoMorita | 74:f5fe7fecbd3c | 251 | void solaESKF::updateGPSPosition(Matrix posgps,float palt,Matrix R) |
| NaotoMorita | 65:c25d7810de44 | 252 | { |
| NaotoMorita | 74:f5fe7fecbd3c | 253 | Matrix H(3,nState); |
| NaotoMorita | 74:f5fe7fecbd3c | 254 | H(1,1) = 1.0f; |
| NaotoMorita | 74:f5fe7fecbd3c | 255 | H(2,2) = 1.0f; |
| NaotoMorita | 74:f5fe7fecbd3c | 256 | H(3,3) = 1.0f; |
| NaotoMorita | 65:c25d7810de44 | 257 | |
| NaotoMorita | 74:f5fe7fecbd3c | 258 | //Matrix A = H*Phat*MatrixMath::Transpose(H)+R; |
| NaotoMorita | 74:f5fe7fecbd3c | 259 | //Matrix K = (Phat*MatrixMath::Transpose(H))*(MatrixMath::Inv(MatrixMath::Transpose(A)*A+1000.0f*MatrixMath::Eye(3))*MatrixMath::Transpose(A)); |
| NaotoMorita | 70:d12e46fdc2f0 | 260 | Matrix K = (Phat*MatrixMath::Transpose(H))*MatrixMath::Inv(H*Phat*MatrixMath::Transpose(H)+R); |
| NaotoMorita | 74:f5fe7fecbd3c | 261 | Matrix z(3,1); |
| NaotoMorita | 74:f5fe7fecbd3c | 262 | z << posgps(1,1) - pihat(1,1) << posgps(2,1)-pihat(2,1) << palt - pihat(3,1); |
| NaotoMorita | 65:c25d7810de44 | 263 | errState = K * z; |
| NaotoMorita | 65:c25d7810de44 | 264 | Phat = (MatrixMath::Eye(nState)-K*H)*Phat; |
| NaotoMorita | 74:f5fe7fecbd3c | 265 | //Phat = Phat-K*(H*Phat*MatrixMath::Transpose(H)+R)*MatrixMath::Transpose(K); |
| NaotoMorita | 74:f5fe7fecbd3c | 266 | //Phat = (MatrixMath::Eye(nState)-K*H)*Phat*MatrixMath::Transpose(MatrixMath::Eye(nState)-K*H)+K*R*MatrixMath::Transpose(K); |
| NaotoMorita | 65:c25d7810de44 | 267 | fuseErr2Nominal(); |
| NaotoMorita | 65:c25d7810de44 | 268 | } |
| NaotoMorita | 75:e2c825cdc511 | 269 | void solaESKF::updateGPSVelocity(Matrix velgps,Matrix R) |
| NaotoMorita | 75:e2c825cdc511 | 270 | { |
| NaotoMorita | 75:e2c825cdc511 | 271 | Matrix H(3,nState); |
| NaotoMorita | 75:e2c825cdc511 | 272 | H(1,4) = 1.0f; |
| NaotoMorita | 75:e2c825cdc511 | 273 | H(2,5) = 1.0f; |
| NaotoMorita | 75:e2c825cdc511 | 274 | H(3,6) = 1.0f; |
| NaotoMorita | 75:e2c825cdc511 | 275 | |
| NaotoMorita | 75:e2c825cdc511 | 276 | //Matrix A = H*Phat*MatrixMath::Transpose(H)+R; |
| NaotoMorita | 75:e2c825cdc511 | 277 | //Matrix K = (Phat*MatrixMath::Transpose(H))*(MatrixMath::Inv(MatrixMath::Transpose(A)*A+1000.0f*MatrixMath::Eye(3))*MatrixMath::Transpose(A)); |
| NaotoMorita | 75:e2c825cdc511 | 278 | Matrix K = (Phat*MatrixMath::Transpose(H))*MatrixMath::Inv(H*Phat*MatrixMath::Transpose(H)+R); |
| NaotoMorita | 75:e2c825cdc511 | 279 | Matrix z(3,1); |
| NaotoMorita | 75:e2c825cdc511 | 280 | z << velgps(1,1) - vihat(1,1) << velgps(2,1)-vihat(2,1)<< velgps(3,1)-vihat(3,1); |
| NaotoMorita | 75:e2c825cdc511 | 281 | errState = K * z; |
| NaotoMorita | 75:e2c825cdc511 | 282 | Phat = (MatrixMath::Eye(nState)-K*H)*Phat; |
| NaotoMorita | 75:e2c825cdc511 | 283 | //Phat = Phat-K*(H*Phat*MatrixMath::Transpose(H)+R)*MatrixMath::Transpose(K); |
| NaotoMorita | 75:e2c825cdc511 | 284 | //Phat = (MatrixMath::Eye(nState)-K*H)*Phat*MatrixMath::Transpose(MatrixMath::Eye(nState)-K*H)+K*R*MatrixMath::Transpose(K); |
| NaotoMorita | 75:e2c825cdc511 | 285 | fuseErr2Nominal(); |
| NaotoMorita | 75:e2c825cdc511 | 286 | } |
| NaotoMorita | 65:c25d7810de44 | 287 | |
| NaotoMorita | 59:03fe5e16a33c | 288 | void solaESKF::updateGPS(Matrix posgps,float palt,Matrix velgps,Matrix R) |
| NaotoMorita | 58:93ba28cf5cb3 | 289 | { |
| NaotoMorita | 77:780ce6556041 | 290 | Matrix H(5,nState); |
| NaotoMorita | 46:15988dc41923 | 291 | H(1,1) = 1.0f; |
| NaotoMorita | 46:15988dc41923 | 292 | H(2,2) = 1.0f; |
| NaotoMorita | 56:c10f1168bd4a | 293 | H(3,3) = 1.0f; |
| NaotoMorita | 56:c10f1168bd4a | 294 | H(4,4) = 1.0f; |
| NaotoMorita | 56:c10f1168bd4a | 295 | H(5,5) = 1.0f; |
| NaotoMorita | 53:8b551358a7e3 | 296 | Matrix K = (Phat*MatrixMath::Transpose(H))*MatrixMath::Inv(H*Phat*MatrixMath::Transpose(H)+R); |
| NaotoMorita | 77:780ce6556041 | 297 | Matrix z(5,1); |
| NaotoMorita | 77:780ce6556041 | 298 | z << posgps(1,1) - pihat(1,1) << posgps(2,1)-pihat(2,1) << palt - pihat(3,1) << velgps(1,1) - vihat(1,1) << velgps(2,1)-vihat(2,1); |
| NaotoMorita | 58:93ba28cf5cb3 | 299 | errState = K * z; |
| NaotoMorita | 58:93ba28cf5cb3 | 300 | Phat = (MatrixMath::Eye(nState)-K*H)*Phat; |
| NaotoMorita | 58:93ba28cf5cb3 | 301 | fuseErr2Nominal(); |
| NaotoMorita | 58:93ba28cf5cb3 | 302 | } |
| NaotoMorita | 58:93ba28cf5cb3 | 303 | |
| NaotoMorita | 46:15988dc41923 | 304 | |
| NaotoMorita | 44:7d82e63b6a86 | 305 | Matrix solaESKF::computeAngles() |
| NaotoMorita | 25:07ac5c6cd61c | 306 | { |
| NaotoMorita | 25:07ac5c6cd61c | 307 | |
| NaotoMorita | 44:7d82e63b6a86 | 308 | Matrix euler(3,1); |
| NaotoMorita | 44:7d82e63b6a86 | 309 | euler(1,1) = atan2f(qhat(1,1)*qhat(2,1) + qhat(3,1)*qhat(4,1), 0.5f - qhat(2,1)*qhat(2,1) - qhat(3,1)*qhat(3,1)); |
| NaotoMorita | 44:7d82e63b6a86 | 310 | euler(2,1) = asinf(-2.0f * (qhat(2,1)*qhat(4,1) - qhat(1,1)*qhat(3,1))); |
| NaotoMorita | 44:7d82e63b6a86 | 311 | euler(3,1) = atan2f(qhat(2,1)*qhat(3,1) + qhat(1,1)*qhat(4,1), 0.5f - qhat(3,1)*qhat(3,1) - qhat(4,1)*qhat(4,1)); |
| NaotoMorita | 44:7d82e63b6a86 | 312 | return euler; |
| NaotoMorita | 19:3fae66745363 | 313 | } |
| NaotoMorita | 19:3fae66745363 | 314 | |
| NaotoMorita | 21:d6079def0473 | 315 | |
| NaotoMorita | 44:7d82e63b6a86 | 316 | void solaESKF::fuseErr2Nominal() |
| NaotoMorita | 19:3fae66745363 | 317 | { |
| NaotoMorita | 44:7d82e63b6a86 | 318 | //position |
| NaotoMorita | 44:7d82e63b6a86 | 319 | pihat(1,1) += errState(1,1); |
| NaotoMorita | 44:7d82e63b6a86 | 320 | pihat(2,1) += errState(2,1); |
| NaotoMorita | 44:7d82e63b6a86 | 321 | pihat(3,1) += errState(3,1); |
| NaotoMorita | 44:7d82e63b6a86 | 322 | |
| NaotoMorita | 44:7d82e63b6a86 | 323 | //velocity |
| NaotoMorita | 44:7d82e63b6a86 | 324 | vihat(1,1) += errState(4,1); |
| NaotoMorita | 44:7d82e63b6a86 | 325 | vihat(2,1) += errState(5,1); |
| NaotoMorita | 44:7d82e63b6a86 | 326 | vihat(3,1) += errState(6,1); |
| NaotoMorita | 44:7d82e63b6a86 | 327 | |
| NaotoMorita | 44:7d82e63b6a86 | 328 | //angle error |
| NaotoMorita | 19:3fae66745363 | 329 | Matrix qerr(4,1); |
| NaotoMorita | 46:15988dc41923 | 330 | qerr << 1.0f << 0.5f*errState(7,1) << 0.5f*errState(8,1) << 0.5f*errState(9,1); |
| NaotoMorita | 19:3fae66745363 | 331 | qhat = quatmultiply(qhat, qerr); |
| NaotoMorita | 19:3fae66745363 | 332 | float qnorm = sqrt(MatrixMath::dot(MatrixMath::Transpose(qhat),qhat)); |
| NaotoMorita | 19:3fae66745363 | 333 | qhat *= (1.0f/ qnorm); |
| NaotoMorita | 23:1509648c2318 | 334 | |
| NaotoMorita | 44:7d82e63b6a86 | 335 | //acc bias |
| NaotoMorita | 44:7d82e63b6a86 | 336 | accBias(1,1) += errState(10,1); |
| NaotoMorita | 44:7d82e63b6a86 | 337 | accBias(2,1) += errState(11,1); |
| NaotoMorita | 44:7d82e63b6a86 | 338 | accBias(3,1) += errState(12,1); |
| NaotoMorita | 44:7d82e63b6a86 | 339 | |
| NaotoMorita | 44:7d82e63b6a86 | 340 | //gyro bias |
| NaotoMorita | 44:7d82e63b6a86 | 341 | gyroBias(1,1) += errState(13,1); |
| NaotoMorita | 44:7d82e63b6a86 | 342 | gyroBias(2,1) += errState(14,1); |
| NaotoMorita | 44:7d82e63b6a86 | 343 | gyroBias(3,1) += errState(15,1); |
| NaotoMorita | 44:7d82e63b6a86 | 344 | |
| NaotoMorita | 44:7d82e63b6a86 | 345 | //gravity bias |
| NaotoMorita | 74:f5fe7fecbd3c | 346 | gravity(1,1) += errState(16,1); |
| NaotoMorita | 74:f5fe7fecbd3c | 347 | gravity(2,1) += errState(17,1); |
| NaotoMorita | 74:f5fe7fecbd3c | 348 | gravity(3,1) += errState(18,1); |
| NaotoMorita | 70:d12e46fdc2f0 | 349 | |
| NaotoMorita | 62:5519d34eb6e8 | 350 | |
| NaotoMorita | 56:c10f1168bd4a | 351 | |
| NaotoMorita | 58:93ba28cf5cb3 | 352 | for (int i = 1; i < nState+1; i++){ |
| NaotoMorita | 47:2467de40951f | 353 | errState(i,1) = 0.0f; |
| NaotoMorita | 47:2467de40951f | 354 | } |
| NaotoMorita | 44:7d82e63b6a86 | 355 | |
| NaotoMorita | 19:3fae66745363 | 356 | } |
| NaotoMorita | 19:3fae66745363 | 357 | |
| NaotoMorita | 44:7d82e63b6a86 | 358 | Matrix solaESKF::quatmultiply(Matrix p, Matrix q) |
| NaotoMorita | 19:3fae66745363 | 359 | { |
| NaotoMorita | 44:7d82e63b6a86 | 360 | |
| NaotoMorita | 19:3fae66745363 | 361 | Matrix qout(4,1); |
| NaotoMorita | 44:7d82e63b6a86 | 362 | qout(1,1) = p(1,1)*q(1,1) - p(2,1)*q(2,1) - p(3,1)*q(3,1) - p(4,1)*q(4,1); |
| NaotoMorita | 44:7d82e63b6a86 | 363 | qout(2,1) = p(1,1)*q(2,1) + p(2,1)*q(1,1) + p(3,1)*q(4,1) - p(4,1)*q(3,1); |
| NaotoMorita | 44:7d82e63b6a86 | 364 | qout(3,1) = p(1,1)*q(3,1) - p(2,1)*q(4,1) + p(3,1)*q(1,1) + p(4,1)*q(2,1); |
| NaotoMorita | 44:7d82e63b6a86 | 365 | qout(4,1) = p(1,1)*q(4,1) + p(2,1)*q(3,1) - p(3,1)*q(2,1) + p(4,1)*q(1,1); |
| NaotoMorita | 19:3fae66745363 | 366 | float qnorm = sqrt(MatrixMath::dot(MatrixMath::Transpose(qout),qout)); |
| NaotoMorita | 19:3fae66745363 | 367 | qout *= (1.0f/ qnorm); |
| NaotoMorita | 19:3fae66745363 | 368 | return qout; |
| NaotoMorita | 19:3fae66745363 | 369 | } |
| NaotoMorita | 19:3fae66745363 | 370 | |
| NaotoMorita | 44:7d82e63b6a86 | 371 | void solaESKF::computeDcm(Matrix& dcm, Matrix quat) |
| NaotoMorita | 19:3fae66745363 | 372 | { |
| NaotoMorita | 19:3fae66745363 | 373 | |
| NaotoMorita | 44:7d82e63b6a86 | 374 | dcm(1,1) = quat(1,1)*quat(1,1) + quat(2,1)*quat(2,1) - quat(3,1)*quat(3,1) - quat(4,1)*quat(4,1); |
| NaotoMorita | 44:7d82e63b6a86 | 375 | dcm(1,2) = 2.0f*(quat(2,1)*quat(3,1) - quat(1,1)*quat(4,1)); |
| NaotoMorita | 44:7d82e63b6a86 | 376 | dcm(1,3) = 2.0f*(quat(2,1)*quat(4,1) + quat(1,1)*quat(3,1)); |
| NaotoMorita | 44:7d82e63b6a86 | 377 | dcm(2,1) = 2.0f*(quat(2,1)*quat(3,1) + quat(1,1)*quat(4,1)); |
| NaotoMorita | 44:7d82e63b6a86 | 378 | dcm(2,2) = quat(1,1)*quat(1,1) - quat(2,1)*quat(2,1) + quat(3,1)*quat(3,1) - quat(4,1)*quat(4,1); |
| NaotoMorita | 44:7d82e63b6a86 | 379 | dcm(2,3) = 2.0f*(quat(3,1)*quat(4,1) - quat(1,1)*quat(2,1)); |
| NaotoMorita | 44:7d82e63b6a86 | 380 | dcm(3,1) = 2.0f*(quat(2,1)*quat(4,1) - quat(1,1)*quat(3,1)); |
| NaotoMorita | 44:7d82e63b6a86 | 381 | dcm(3,2) = 2.0f*(quat(3,1)*quat(4,1) + quat(1,1)*quat(2,1)); |
| NaotoMorita | 44:7d82e63b6a86 | 382 | dcm(3,3) = quat(1,1)*quat(1,1) - quat(2,1)*quat(2,1) - quat(3,1)*quat(3,1) + quat(4,1)*quat(4,1); |
| NaotoMorita | 19:3fae66745363 | 383 | |
| NaotoMorita | 19:3fae66745363 | 384 | } |
| NaotoMorita | 19:3fae66745363 | 385 | |
| NaotoMorita | 44:7d82e63b6a86 | 386 | void solaESKF::setQhat(float ex,float ey,float ez) |
| NaotoMorita | 44:7d82e63b6a86 | 387 | { |
| NaotoMorita | 44:7d82e63b6a86 | 388 | float cos_z_2 = cosf(0.5f*ez); |
| NaotoMorita | 44:7d82e63b6a86 | 389 | float cos_y_2 = cosf(0.5f*ey); |
| NaotoMorita | 44:7d82e63b6a86 | 390 | float cos_x_2 = cosf(0.5f*ex); |
| NaotoMorita | 19:3fae66745363 | 391 | |
| NaotoMorita | 44:7d82e63b6a86 | 392 | float sin_z_2 = sinf(0.5f*ez); |
| NaotoMorita | 44:7d82e63b6a86 | 393 | float sin_y_2 = sinf(0.5f*ey); |
| NaotoMorita | 44:7d82e63b6a86 | 394 | float sin_x_2 = sinf(0.5f*ex); |
| NaotoMorita | 19:3fae66745363 | 395 | |
| NaotoMorita | 19:3fae66745363 | 396 | // and now compute quaternion |
| NaotoMorita | 19:3fae66745363 | 397 | qhat(1,1) = cos_z_2*cos_y_2*cos_x_2 + sin_z_2*sin_y_2*sin_x_2; |
| NaotoMorita | 19:3fae66745363 | 398 | qhat(2,1) = cos_z_2*cos_y_2*sin_x_2 - sin_z_2*sin_y_2*cos_x_2; |
| NaotoMorita | 19:3fae66745363 | 399 | qhat(3,1) = cos_z_2*sin_y_2*cos_x_2 + sin_z_2*cos_y_2*sin_x_2; |
| NaotoMorita | 19:3fae66745363 | 400 | qhat(4,1) = sin_z_2*cos_y_2*cos_x_2 - cos_z_2*sin_y_2*sin_x_2; |
| NaotoMorita | 19:3fae66745363 | 401 | } |
| NaotoMorita | 19:3fae66745363 | 402 | |
| NaotoMorita | 75:e2c825cdc511 | 403 | Matrix solaESKF::calcDynAcc(Matrix acc) |
| NaotoMorita | 75:e2c825cdc511 | 404 | { |
| NaotoMorita | 75:e2c825cdc511 | 405 | Matrix accm = acc - accBias; |
| NaotoMorita | 75:e2c825cdc511 | 406 | Matrix tdcm(3,3); |
| NaotoMorita | 75:e2c825cdc511 | 407 | computeDcm(tdcm, qhat); |
| NaotoMorita | 75:e2c825cdc511 | 408 | tdcm = MatrixMath::Transpose(tdcm); |
| NaotoMorita | 75:e2c825cdc511 | 409 | |
| NaotoMorita | 75:e2c825cdc511 | 410 | Matrix dynAcc = accm+tdcm*gravity; |
| NaotoMorita | 75:e2c825cdc511 | 411 | return dynAcc; |
| NaotoMorita | 75:e2c825cdc511 | 412 | } |
| NaotoMorita | 47:2467de40951f | 413 | |
| NaotoMorita | 47:2467de40951f | 414 | |
| NaotoMorita | 44:7d82e63b6a86 | 415 | void solaESKF::setGravity(float gx,float gy,float gz) |
| NaotoMorita | 23:1509648c2318 | 416 | { |
| NaotoMorita | 44:7d82e63b6a86 | 417 | gravity(1,1) = gx; |
| NaotoMorita | 44:7d82e63b6a86 | 418 | gravity(2,1) = gy; |
| NaotoMorita | 44:7d82e63b6a86 | 419 | gravity(3,1) = gz; |
| NaotoMorita | 23:1509648c2318 | 420 | } |
| NaotoMorita | 62:5519d34eb6e8 | 421 | |
| NaotoMorita | 47:2467de40951f | 422 | Matrix solaESKF::getPihat() |
| NaotoMorita | 47:2467de40951f | 423 | { |
| NaotoMorita | 47:2467de40951f | 424 | return pihat; |
| NaotoMorita | 47:2467de40951f | 425 | } |
| NaotoMorita | 47:2467de40951f | 426 | Matrix solaESKF::getVihat() |
| NaotoMorita | 47:2467de40951f | 427 | { |
| NaotoMorita | 47:2467de40951f | 428 | return vihat; |
| NaotoMorita | 47:2467de40951f | 429 | } |
| NaotoMorita | 47:2467de40951f | 430 | Matrix solaESKF::getQhat() |
| NaotoMorita | 47:2467de40951f | 431 | { |
| NaotoMorita | 47:2467de40951f | 432 | return qhat; |
| NaotoMorita | 47:2467de40951f | 433 | } |
| NaotoMorita | 47:2467de40951f | 434 | Matrix solaESKF::getAccBias() |
| NaotoMorita | 47:2467de40951f | 435 | { |
| NaotoMorita | 47:2467de40951f | 436 | return accBias; |
| NaotoMorita | 47:2467de40951f | 437 | } |
| NaotoMorita | 47:2467de40951f | 438 | Matrix solaESKF::getGyroBias() |
| NaotoMorita | 47:2467de40951f | 439 | { |
| NaotoMorita | 47:2467de40951f | 440 | return gyroBias; |
| NaotoMorita | 47:2467de40951f | 441 | } |
| NaotoMorita | 47:2467de40951f | 442 | Matrix solaESKF::getGravity() |
| NaotoMorita | 47:2467de40951f | 443 | { |
| NaotoMorita | 47:2467de40951f | 444 | return gravity; |
| NaotoMorita | 47:2467de40951f | 445 | } |
| NaotoMorita | 74:f5fe7fecbd3c | 446 | |
| NaotoMorita | 47:2467de40951f | 447 | Matrix solaESKF::getErrState() |
| NaotoMorita | 47:2467de40951f | 448 | { |
| NaotoMorita | 47:2467de40951f | 449 | return errState; |
| NaotoMorita | 47:2467de40951f | 450 | } |
| NaotoMorita | 23:1509648c2318 | 451 | |
| NaotoMorita | 47:2467de40951f | 452 | void solaESKF::setPhatPosition(float val) |
| NaotoMorita | 47:2467de40951f | 453 | { |
| NaotoMorita | 47:2467de40951f | 454 | setBlockDiag(Phat,val,1,3); |
| NaotoMorita | 47:2467de40951f | 455 | } |
| NaotoMorita | 47:2467de40951f | 456 | void solaESKF::setPhatVelocity(float val) |
| NaotoMorita | 47:2467de40951f | 457 | { |
| NaotoMorita | 47:2467de40951f | 458 | setBlockDiag(Phat,val,4,6); |
| NaotoMorita | 47:2467de40951f | 459 | } |
| NaotoMorita | 47:2467de40951f | 460 | void solaESKF::setPhatAngleError(float val) |
| NaotoMorita | 47:2467de40951f | 461 | { |
| NaotoMorita | 47:2467de40951f | 462 | setBlockDiag(Phat,val,7,9); |
| NaotoMorita | 47:2467de40951f | 463 | } |
| NaotoMorita | 47:2467de40951f | 464 | void solaESKF::setPhatAccBias(float val) |
| NaotoMorita | 47:2467de40951f | 465 | { |
| NaotoMorita | 47:2467de40951f | 466 | setBlockDiag(Phat,val,10,12); |
| NaotoMorita | 47:2467de40951f | 467 | } |
| NaotoMorita | 47:2467de40951f | 468 | void solaESKF::setPhatGyroBias(float val) |
| NaotoMorita | 47:2467de40951f | 469 | { |
| NaotoMorita | 47:2467de40951f | 470 | setBlockDiag(Phat,val,13,15); |
| NaotoMorita | 47:2467de40951f | 471 | } |
| NaotoMorita | 47:2467de40951f | 472 | void solaESKF::setPhatGravity(float val) |
| NaotoMorita | 47:2467de40951f | 473 | { |
| NaotoMorita | 74:f5fe7fecbd3c | 474 | setBlockDiag(Phat,val,16,18); |
| NaotoMorita | 62:5519d34eb6e8 | 475 | } |
| NaotoMorita | 47:2467de40951f | 476 | |
| NaotoMorita | 47:2467de40951f | 477 | |
| NaotoMorita | 47:2467de40951f | 478 | void solaESKF::setQVelocity(float val) |
| NaotoMorita | 47:2467de40951f | 479 | { |
| NaotoMorita | 47:2467de40951f | 480 | setBlockDiag(Q,val,4,6); |
| NaotoMorita | 47:2467de40951f | 481 | } |
| NaotoMorita | 47:2467de40951f | 482 | void solaESKF::setQAngleError(float val) |
| NaotoMorita | 47:2467de40951f | 483 | { |
| NaotoMorita | 47:2467de40951f | 484 | setBlockDiag(Q,val,7,9); |
| NaotoMorita | 47:2467de40951f | 485 | } |
| NaotoMorita | 47:2467de40951f | 486 | void solaESKF::setQAccBias(float val) |
| NaotoMorita | 47:2467de40951f | 487 | { |
| NaotoMorita | 47:2467de40951f | 488 | setBlockDiag(Q,val,10,12); |
| NaotoMorita | 47:2467de40951f | 489 | } |
| NaotoMorita | 47:2467de40951f | 490 | void solaESKF::setQGyroBias(float val) |
| NaotoMorita | 47:2467de40951f | 491 | { |
| NaotoMorita | 47:2467de40951f | 492 | setBlockDiag(Q,val,13,15); |
| NaotoMorita | 47:2467de40951f | 493 | } |
| NaotoMorita | 19:3fae66745363 | 494 | |
| NaotoMorita | 19:3fae66745363 | 495 | |
| NaotoMorita | 44:7d82e63b6a86 | 496 | void solaESKF::setDiag(Matrix& mat, float val){ |
| NaotoMorita | 44:7d82e63b6a86 | 497 | for (int i = 1; i < mat.getCols()+1; i++){ |
| NaotoMorita | 44:7d82e63b6a86 | 498 | mat(i,i) = val; |
| NaotoMorita | 22:7d84b8bc20b4 | 499 | } |
| NaotoMorita | 22:7d84b8bc20b4 | 500 | } |
| NaotoMorita | 22:7d84b8bc20b4 | 501 | |
| NaotoMorita | 45:df4618814803 | 502 | void solaESKF::setBlockDiag(Matrix& mat, float val,int startIndex, int endIndex){ |
| NaotoMorita | 44:7d82e63b6a86 | 503 | for (int i = startIndex; i < endIndex+1; i++){ |
| NaotoMorita | 44:7d82e63b6a86 | 504 | mat(i,i) = val; |
| NaotoMorita | 22:7d84b8bc20b4 | 505 | } |
| NaotoMorita | 22:7d84b8bc20b4 | 506 | } |
| osaka | 51:a5af3b280d23 | 507 | |
| NaotoMorita | 58:93ba28cf5cb3 | 508 | void solaESKF::setPihat(float pi_x, float pi_y) |
| osaka | 51:a5af3b280d23 | 509 | { |
| NaotoMorita | 56:c10f1168bd4a | 510 | pihat(1,1) = pi_x; |
| NaotoMorita | 56:c10f1168bd4a | 511 | pihat(2,1) = pi_y; |
| NaotoMorita | 58:93ba28cf5cb3 | 512 | //pihat(3,1) = pi_z; |
| NaotoMorita | 56:c10f1168bd4a | 513 | } |
| NaotoMorita | 73:5770a0d470c0 | 514 | |
| NaotoMorita | 56:c10f1168bd4a | 515 | /* |
| NaotoMorita | 56:c10f1168bd4a | 516 | void solaESKF::updateAccConstraints(Matrix acc,float palt,Matrix R) |
| NaotoMorita | 56:c10f1168bd4a | 517 | { |
| NaotoMorita | 56:c10f1168bd4a | 518 | Matrix accm = acc - accBias; |
| NaotoMorita | 56:c10f1168bd4a | 519 | Matrix tdcm(3,3); |
| NaotoMorita | 56:c10f1168bd4a | 520 | computeDcm(tdcm, qhat); |
| NaotoMorita | 56:c10f1168bd4a | 521 | tdcm = MatrixMath::Transpose(tdcm); |
| NaotoMorita | 56:c10f1168bd4a | 522 | Matrix tdcm_g = tdcm*gravity; |
| NaotoMorita | 56:c10f1168bd4a | 523 | Matrix a2v = MatrixMath::Matrixcross(tdcm_g(1,1),tdcm_g(2,1),tdcm_g(3,1)); |
| NaotoMorita | 56:c10f1168bd4a | 524 | |
| NaotoMorita | 56:c10f1168bd4a | 525 | Matrix H(4,nState); |
| NaotoMorita | 56:c10f1168bd4a | 526 | for (int i = 1; i < 4; i++){ |
| NaotoMorita | 56:c10f1168bd4a | 527 | for (int j = 1; j < 4; j++){ |
| NaotoMorita | 56:c10f1168bd4a | 528 | H(i,j+6) = a2v(i,j); |
| NaotoMorita | 56:c10f1168bd4a | 529 | H(i,j+15) = tdcm(i,j); |
| NaotoMorita | 56:c10f1168bd4a | 530 | } |
| NaotoMorita | 56:c10f1168bd4a | 531 | } |
| NaotoMorita | 56:c10f1168bd4a | 532 | H(1,10) = -1.0f; |
| NaotoMorita | 56:c10f1168bd4a | 533 | H(2,11) = -1.0f; |
| NaotoMorita | 56:c10f1168bd4a | 534 | H(3,12) = -1.0f; |
| NaotoMorita | 56:c10f1168bd4a | 535 | H(4,3) = 1.0f; |
| NaotoMorita | 56:c10f1168bd4a | 536 | |
| NaotoMorita | 56:c10f1168bd4a | 537 | //Matrix A = H*Phat*MatrixMath::Transpose(H)+R; |
| NaotoMorita | 56:c10f1168bd4a | 538 | //Matrix K = (Phat*MatrixMath::Transpose(H))*(MatrixMath::Inv(MatrixMath::Transpose(A)*A+10.0f*MatrixMath::Eye(3))*MatrixMath::Transpose(A)); |
| NaotoMorita | 56:c10f1168bd4a | 539 | Matrix K = (Phat*MatrixMath::Transpose(H))*MatrixMath::Inv(H*Phat*MatrixMath::Transpose(H)+R); |
| NaotoMorita | 56:c10f1168bd4a | 540 | Matrix zacc = -accm-tdcm*gravity; |
| NaotoMorita | 56:c10f1168bd4a | 541 | Matrix z(4,1); |
| NaotoMorita | 56:c10f1168bd4a | 542 | z << zacc(1,1) << zacc(2,1) << zacc(3,1) << palt - pihat(3,1); |
| NaotoMorita | 56:c10f1168bd4a | 543 | errState = K * z; |
| NaotoMorita | 56:c10f1168bd4a | 544 | Phat = (MatrixMath::Eye(nState)-K*H)*Phat; |
| NaotoMorita | 56:c10f1168bd4a | 545 | //Phat = Phat-K*(H*Phat*MatrixMath::Transpose(H)+R)*MatrixMath::Transpose(K); |
| NaotoMorita | 56:c10f1168bd4a | 546 | //Phat = (MatrixMath::Eye(nState)-K*H)*Phat*MatrixMath::Transpose(MatrixMath::Eye(nState)-K*H)+K*R*MatrixMath::Transpose(K); |
| NaotoMorita | 56:c10f1168bd4a | 547 | fuseErr2Nominal(); |
| NaotoMorita | 56:c10f1168bd4a | 548 | } |
| NaotoMorita | 56:c10f1168bd4a | 549 | |
| NaotoMorita | 56:c10f1168bd4a | 550 | void solaESKF::updateGyroConstraints(Matrix gyro,Matrix R) |
| NaotoMorita | 56:c10f1168bd4a | 551 | { |
| NaotoMorita | 56:c10f1168bd4a | 552 | Matrix gyrom = gyro - gyroBias; |
| NaotoMorita | 56:c10f1168bd4a | 553 | Matrix dcm(3,3); |
| NaotoMorita | 56:c10f1168bd4a | 554 | computeDcm(dcm, qhat); |
| NaotoMorita | 56:c10f1168bd4a | 555 | Matrix a2v = dcm*MatrixMath::Matrixcross(gyrom(1,1),gyrom(2,1),gyrom(3,1)); |
| NaotoMorita | 56:c10f1168bd4a | 556 | |
| NaotoMorita | 56:c10f1168bd4a | 557 | Matrix H(2,nState); |
| NaotoMorita | 56:c10f1168bd4a | 558 | for (int i = 1; i < 3; i++){ |
| NaotoMorita | 56:c10f1168bd4a | 559 | for (int j = 1; j < 4; j++){ |
| NaotoMorita | 56:c10f1168bd4a | 560 | H(i,j+6) = a2v(i,j); |
| NaotoMorita | 56:c10f1168bd4a | 561 | H(i,j+12) = -dcm(i,j); |
| NaotoMorita | 56:c10f1168bd4a | 562 | } |
| NaotoMorita | 56:c10f1168bd4a | 563 | } |
| NaotoMorita | 56:c10f1168bd4a | 564 | //Matrix A = H*Phat*MatrixMath::Transpose(H)+R; |
| NaotoMorita | 56:c10f1168bd4a | 565 | //Matrix K = (Phat*MatrixMath::Transpose(H))*(MatrixMath::Inv(MatrixMath::Transpose(A)*A+10.0f*MatrixMath::Eye(2))*MatrixMath::Transpose(A)); |
| NaotoMorita | 56:c10f1168bd4a | 566 | Matrix K = (Phat*MatrixMath::Transpose(H))*MatrixMath::Inv(H*Phat*MatrixMath::Transpose(H)+R); |
| NaotoMorita | 56:c10f1168bd4a | 567 | |
| NaotoMorita | 56:c10f1168bd4a | 568 | Matrix z3 = -dcm*gyrom; |
| NaotoMorita | 56:c10f1168bd4a | 569 | Matrix z(2,1); |
| NaotoMorita | 56:c10f1168bd4a | 570 | z << z3(1,1) << z3(2,1); |
| NaotoMorita | 56:c10f1168bd4a | 571 | errState = K * z; |
| NaotoMorita | 56:c10f1168bd4a | 572 | //Phat = (MatrixMath::Eye(nState)-K*H)*Phat; |
| NaotoMorita | 56:c10f1168bd4a | 573 | //Phat = Phat-K*(H*Phat*MatrixMath::Transpose(H)+R)*MatrixMath::Transpose(K); |
| NaotoMorita | 56:c10f1168bd4a | 574 | Phat = (MatrixMath::Eye(nState)-K*H)*Phat*MatrixMath::Transpose(MatrixMath::Eye(nState)-K*H)+K*R*MatrixMath::Transpose(K); |
| NaotoMorita | 56:c10f1168bd4a | 575 | fuseErr2Nominal(); |
| NaotoMorita | 56:c10f1168bd4a | 576 | } |
| NaotoMorita | 56:c10f1168bd4a | 577 | void solaESKF::updateMag(Matrix mag, Matrix R) |
| NaotoMorita | 56:c10f1168bd4a | 578 | { |
| NaotoMorita | 56:c10f1168bd4a | 579 | Matrix dcm(3,3); |
| NaotoMorita | 56:c10f1168bd4a | 580 | computeDcm(dcm, qhat); |
| NaotoMorita | 56:c10f1168bd4a | 581 | Matrix a2v = -dcm*MatrixMath::Matrixcross(mag(1,1),mag(2,1),mag(3,1)); |
| NaotoMorita | 56:c10f1168bd4a | 582 | |
| NaotoMorita | 56:c10f1168bd4a | 583 | Matrix H(2,nState); |
| NaotoMorita | 56:c10f1168bd4a | 584 | for (int i = 1; i < 3; i++){ |
| NaotoMorita | 56:c10f1168bd4a | 585 | for (int j = 1; j < 4; j++){ |
| NaotoMorita | 56:c10f1168bd4a | 586 | H(i,j+6) = a2v(i,j); |
| NaotoMorita | 56:c10f1168bd4a | 587 | } |
| NaotoMorita | 56:c10f1168bd4a | 588 | } |
| NaotoMorita | 56:c10f1168bd4a | 589 | H(1,19) = 1.0f; |
| NaotoMorita | 56:c10f1168bd4a | 590 | //H(3,20) = 1.0f; |
| NaotoMorita | 56:c10f1168bd4a | 591 | |
| NaotoMorita | 56:c10f1168bd4a | 592 | //Matrix A = H*Phat*MatrixMath::Transpose(H)+R; |
| NaotoMorita | 56:c10f1168bd4a | 593 | //Matrix K = (Phat*MatrixMath::Transpose(H))*(MatrixMath::Inv(MatrixMath::Transpose(A)*A+10.0f*MatrixMath::Eye(3))*MatrixMath::Transpose(A)); |
| NaotoMorita | 56:c10f1168bd4a | 594 | Matrix K = (Phat*MatrixMath::Transpose(H))*MatrixMath::Inv(H*Phat*MatrixMath::Transpose(H)+R); |
| NaotoMorita | 56:c10f1168bd4a | 595 | Matrix zmag = -dcm*mag-magField; |
| NaotoMorita | 56:c10f1168bd4a | 596 | Matrix z(2,1); |
| NaotoMorita | 56:c10f1168bd4a | 597 | z << zmag(1,1) << zmag(2,1); |
| NaotoMorita | 56:c10f1168bd4a | 598 | errState = K * z; |
| NaotoMorita | 56:c10f1168bd4a | 599 | Phat = (MatrixMath::Eye(nState)-K*H)*Phat; |
| NaotoMorita | 56:c10f1168bd4a | 600 | //Phat = Phat-K*(H*Phat*MatrixMath::Transpose(H)+R)*MatrixMath::Transpose(K); |
| NaotoMorita | 56:c10f1168bd4a | 601 | //Phat = (MatrixMath::Eye(nState)-K*H)*Phat*MatrixMath::Transpose(MatrixMath::Eye(nState)-K*H)+K*R*MatrixMath::Transpose(K); |
| NaotoMorita | 56:c10f1168bd4a | 602 | fuseErr2Nominal(); |
| NaotoMorita | 56:c10f1168bd4a | 603 | } |
| NaotoMorita | 59:03fe5e16a33c | 604 | void solaESKF::updateMag(Matrix mag,float palt, Matrix R) |
| NaotoMorita | 59:03fe5e16a33c | 605 | { |
| NaotoMorita | 59:03fe5e16a33c | 606 | Matrix dcm(3,3); |
| NaotoMorita | 59:03fe5e16a33c | 607 | computeDcm(dcm, qhat); |
| NaotoMorita | 59:03fe5e16a33c | 608 | Matrix a2v = -dcm*MatrixMath::Matrixcross(mag(1,1),mag(2,1),mag(3,1)); |
| NaotoMorita | 59:03fe5e16a33c | 609 | |
| NaotoMorita | 59:03fe5e16a33c | 610 | Matrix H(3,nState); |
| NaotoMorita | 59:03fe5e16a33c | 611 | for (int i = 1; i < 3; i++){ |
| NaotoMorita | 59:03fe5e16a33c | 612 | for (int j = 1; j < 4; j++){ |
| NaotoMorita | 59:03fe5e16a33c | 613 | H(i,j+6) = a2v(i,j); |
| NaotoMorita | 59:03fe5e16a33c | 614 | } |
| NaotoMorita | 59:03fe5e16a33c | 615 | } |
| NaotoMorita | 59:03fe5e16a33c | 616 | H(1,19) = 1.0f; |
| NaotoMorita | 59:03fe5e16a33c | 617 | H(3,3) = 1.0f; |
| NaotoMorita | 59:03fe5e16a33c | 618 | //Matrix A = H*Phat*MatrixMath::Transpose(H)+R; |
| NaotoMorita | 59:03fe5e16a33c | 619 | //Matrix K = (Phat*MatrixMath::Transpose(H))*(MatrixMath::Inv(MatrixMath::Transpose(A)*A+10.0f*MatrixMath::Eye(3))*MatrixMath::Transpose(A)); |
| NaotoMorita | 59:03fe5e16a33c | 620 | Matrix K = (Phat*MatrixMath::Transpose(H))*MatrixMath::Inv(H*Phat*MatrixMath::Transpose(H)+R); |
| NaotoMorita | 59:03fe5e16a33c | 621 | Matrix zmag = -dcm*mag-magField; |
| NaotoMorita | 59:03fe5e16a33c | 622 | Matrix z(3,1); |
| NaotoMorita | 59:03fe5e16a33c | 623 | z << zmag(1,1) << zmag(2,1) << palt - pihat(3,1); |
| NaotoMorita | 59:03fe5e16a33c | 624 | errState = K * z; |
| NaotoMorita | 59:03fe5e16a33c | 625 | Phat = (MatrixMath::Eye(nState)-K*H)*Phat; |
| NaotoMorita | 59:03fe5e16a33c | 626 | //Phat = Phat-K*(H*Phat*MatrixMath::Transpose(H)+R)*MatrixMath::Transpose(K); |
| NaotoMorita | 59:03fe5e16a33c | 627 | //Phat = (MatrixMath::Eye(nState)-K*H)*Phat*MatrixMath::Transpose(MatrixMath::Eye(nState)-K*H)+K*R*MatrixMath::Transpose(K); |
| NaotoMorita | 59:03fe5e16a33c | 628 | |
| NaotoMorita | 59:03fe5e16a33c | 629 | fuseErr2Nominal(); |
| NaotoMorita | 59:03fe5e16a33c | 630 | } |
| NaotoMorita | 59:03fe5e16a33c | 631 | void solaESKF::updateGPSVelocity(Matrix velgps,Matrix mag,Matrix R) |
| NaotoMorita | 59:03fe5e16a33c | 632 | { |
| NaotoMorita | 59:03fe5e16a33c | 633 | |
| NaotoMorita | 59:03fe5e16a33c | 634 | Matrix dcm(3,3); |
| NaotoMorita | 59:03fe5e16a33c | 635 | computeDcm(dcm, qhat); |
| NaotoMorita | 59:03fe5e16a33c | 636 | Matrix a2v = -dcm*MatrixMath::Matrixcross(mag(1,1),mag(2,1),mag(3,1)); |
| NaotoMorita | 59:03fe5e16a33c | 637 | |
| NaotoMorita | 59:03fe5e16a33c | 638 | |
| NaotoMorita | 59:03fe5e16a33c | 639 | Matrix H(3,nState); |
| NaotoMorita | 59:03fe5e16a33c | 640 | H(1,4) = 1.0f; |
| NaotoMorita | 59:03fe5e16a33c | 641 | H(2,5) = 1.0f; |
| NaotoMorita | 59:03fe5e16a33c | 642 | |
| NaotoMorita | 59:03fe5e16a33c | 643 | for (int j = 1; j < 4; j++){ |
| NaotoMorita | 59:03fe5e16a33c | 644 | H(3,j+6) = a2v(2,j); |
| NaotoMorita | 59:03fe5e16a33c | 645 | } |
| NaotoMorita | 59:03fe5e16a33c | 646 | //H(3,19) = 1.0f; |
| NaotoMorita | 59:03fe5e16a33c | 647 | |
| NaotoMorita | 59:03fe5e16a33c | 648 | Matrix zmag = -dcm*mag-magField; |
| NaotoMorita | 59:03fe5e16a33c | 649 | |
| NaotoMorita | 59:03fe5e16a33c | 650 | //Matrix A = H*Phat*MatrixMath::Transpose(H)+R; |
| NaotoMorita | 59:03fe5e16a33c | 651 | //Matrix K = (Phat*MatrixMath::Transpose(H))*(MatrixMath::Inv(MatrixMath::Transpose(A)*A+10.0f*MatrixMath::Eye(3))*MatrixMath::Transpose(A)); |
| NaotoMorita | 59:03fe5e16a33c | 652 | Matrix K = (Phat*MatrixMath::Transpose(H))*MatrixMath::Inv(H*Phat*MatrixMath::Transpose(H)+R); |
| NaotoMorita | 59:03fe5e16a33c | 653 | Matrix z(3,1); |
| NaotoMorita | 59:03fe5e16a33c | 654 | z << velgps(1,1) - vihat(1,1) << velgps(2,1)-vihat(2,1) << zmag(2,1); |
| NaotoMorita | 59:03fe5e16a33c | 655 | errState = K * z; |
| NaotoMorita | 59:03fe5e16a33c | 656 | Phat = (MatrixMath::Eye(nState)-K*H)*Phat; |
| NaotoMorita | 59:03fe5e16a33c | 657 | //Phat = Phat-K*(H*Phat*MatrixMath::Transpose(H)+R)*MatrixMath::Transpose(K); |
| NaotoMorita | 59:03fe5e16a33c | 658 | //Phat = (MatrixMath::Eye(nState)-K*H)*Phat*MatrixMath::Transpose(MatrixMath::Eye(nState)-K*H)+K*R*MatrixMath::Transpose(K); |
| NaotoMorita | 59:03fe5e16a33c | 659 | fuseErr2Nominal(); |
| NaotoMorita | 59:03fe5e16a33c | 660 | } |
| NaotoMorita | 56:c10f1168bd4a | 661 | |
| NaotoMorita | 59:03fe5e16a33c | 662 | void solaESKF::updateGPSPosition(Matrix posgps,float palt,Matrix R) |
| NaotoMorita | 59:03fe5e16a33c | 663 | { |
| NaotoMorita | 59:03fe5e16a33c | 664 | Matrix H(3,nState); |
| NaotoMorita | 59:03fe5e16a33c | 665 | H(1,1) = 1.0f; |
| NaotoMorita | 59:03fe5e16a33c | 666 | H(2,2) = 1.0f; |
| NaotoMorita | 59:03fe5e16a33c | 667 | H(3,3) = 1.0f; |
| NaotoMorita | 59:03fe5e16a33c | 668 | |
| NaotoMorita | 59:03fe5e16a33c | 669 | //Matrix A = H*Phat*MatrixMath::Transpose(H)+R; |
| NaotoMorita | 59:03fe5e16a33c | 670 | //Matrix K = (Phat*MatrixMath::Transpose(H))*(MatrixMath::Inv(MatrixMath::Transpose(A)*A+1000.0f*MatrixMath::Eye(3))*MatrixMath::Transpose(A)); |
| NaotoMorita | 59:03fe5e16a33c | 671 | Matrix K = (Phat*MatrixMath::Transpose(H))*MatrixMath::Inv(H*Phat*MatrixMath::Transpose(H)+R); |
| NaotoMorita | 59:03fe5e16a33c | 672 | Matrix z(3,1); |
| NaotoMorita | 59:03fe5e16a33c | 673 | z << posgps(1,1) - pihat(1,1) << posgps(2,1)-pihat(2,1) << palt - pihat(3,1); |
| NaotoMorita | 59:03fe5e16a33c | 674 | errState = K * z; |
| NaotoMorita | 59:03fe5e16a33c | 675 | Phat = (MatrixMath::Eye(nState)-K*H)*Phat; |
| NaotoMorita | 59:03fe5e16a33c | 676 | //Phat = Phat-K*(H*Phat*MatrixMath::Transpose(H)+R)*MatrixMath::Transpose(K); |
| NaotoMorita | 59:03fe5e16a33c | 677 | //Phat = (MatrixMath::Eye(nState)-K*H)*Phat*MatrixMath::Transpose(MatrixMath::Eye(nState)-K*H)+K*R*MatrixMath::Transpose(K); |
| NaotoMorita | 59:03fe5e16a33c | 678 | fuseErr2Nominal(); |
| NaotoMorita | 59:03fe5e16a33c | 679 | } |
| NaotoMorita | 61:5e5c4fe12440 | 680 | |
| NaotoMorita | 61:5e5c4fe12440 | 681 | void solaESKF::updateImuConstraints(Matrix acc,Matrix mag,Matrix R) |
| NaotoMorita | 61:5e5c4fe12440 | 682 | { |
| NaotoMorita | 61:5e5c4fe12440 | 683 | Matrix accm = acc - accBias; |
| NaotoMorita | 61:5e5c4fe12440 | 684 | Matrix magm = mag - magBias; |
| NaotoMorita | 61:5e5c4fe12440 | 685 | Matrix dcm(3,3); |
| NaotoMorita | 61:5e5c4fe12440 | 686 | computeDcm(dcm, qhat); |
| NaotoMorita | 61:5e5c4fe12440 | 687 | Matrix tdcm = MatrixMath::Transpose(dcm); |
| NaotoMorita | 61:5e5c4fe12440 | 688 | Matrix tdcm_g = tdcm*gravity; |
| NaotoMorita | 61:5e5c4fe12440 | 689 | Matrix rotgrav = MatrixMath::Matrixcross(tdcm_g(1,1),tdcm_g(2,1),tdcm_g(3,1)); |
| NaotoMorita | 61:5e5c4fe12440 | 690 | Matrix rotmag = -dcm*MatrixMath::Matrixcross(magm(1,1),magm(2,1),magm(3,1)); |
| NaotoMorita | 61:5e5c4fe12440 | 691 | |
| NaotoMorita | 61:5e5c4fe12440 | 692 | Matrix H(5,nState); |
| NaotoMorita | 61:5e5c4fe12440 | 693 | for (int i = 1; i < 4; i++){ |
| NaotoMorita | 61:5e5c4fe12440 | 694 | for (int j = 1; j < 4; j++){ |
| NaotoMorita | 61:5e5c4fe12440 | 695 | H(i,j+6) = rotgrav(i,j); |
| NaotoMorita | 61:5e5c4fe12440 | 696 | } |
| NaotoMorita | 61:5e5c4fe12440 | 697 | H(i,16) = tdcm(i,3); |
| NaotoMorita | 61:5e5c4fe12440 | 698 | } |
| NaotoMorita | 61:5e5c4fe12440 | 699 | |
| NaotoMorita | 61:5e5c4fe12440 | 700 | H(1,10) = -1.0f; |
| NaotoMorita | 61:5e5c4fe12440 | 701 | H(2,11) = -1.0f; |
| NaotoMorita | 61:5e5c4fe12440 | 702 | H(3,12) = -1.0f; |
| NaotoMorita | 61:5e5c4fe12440 | 703 | |
| NaotoMorita | 61:5e5c4fe12440 | 704 | Matrix magned = dcm*magm; |
| NaotoMorita | 61:5e5c4fe12440 | 705 | float hx = sqrt(magned(1,1)*magned(1,1)+magned(2,1)*magned(2,1)); |
| NaotoMorita | 61:5e5c4fe12440 | 706 | |
| NaotoMorita | 61:5e5c4fe12440 | 707 | for(int j = 1; j < 4; j++){ |
| NaotoMorita | 61:5e5c4fe12440 | 708 | H(4,j+6) = rotmag(1,j)-(rotmag(1,j)+rotmag(2,j))/hx; |
| NaotoMorita | 61:5e5c4fe12440 | 709 | H(4,j+16) = -dcm(1,j)+(dcm(1,j)+dcm(2,j))/hx; |
| NaotoMorita | 61:5e5c4fe12440 | 710 | H(5,j+6) = rotmag(2,j); |
| NaotoMorita | 61:5e5c4fe12440 | 711 | H(5,j+16) = -dcm(2,j); |
| NaotoMorita | 61:5e5c4fe12440 | 712 | } |
| NaotoMorita | 61:5e5c4fe12440 | 713 | |
| NaotoMorita | 61:5e5c4fe12440 | 714 | Matrix K = (Phat*MatrixMath::Transpose(H))*MatrixMath::Inv(H*Phat*MatrixMath::Transpose(H)+R); |
| NaotoMorita | 61:5e5c4fe12440 | 715 | Matrix zacc = -accm-tdcm*gravity; |
| NaotoMorita | 61:5e5c4fe12440 | 716 | Matrix zmag = dcm*magm; |
| NaotoMorita | 61:5e5c4fe12440 | 717 | Matrix z(5,1); |
| NaotoMorita | 61:5e5c4fe12440 | 718 | z << zacc(1,1) << zacc(2,1) << zacc(3,1) << -(zmag(1,1) - hx) << -zmag(2,1); |
| NaotoMorita | 61:5e5c4fe12440 | 719 | twelite.printf("%f %f\r\n",hx,(zmag(1,1) - hx)); |
| NaotoMorita | 61:5e5c4fe12440 | 720 | errState = K * z; |
| NaotoMorita | 61:5e5c4fe12440 | 721 | Phat = (MatrixMath::Eye(nState)-K*H)*Phat; |
| NaotoMorita | 61:5e5c4fe12440 | 722 | |
| NaotoMorita | 61:5e5c4fe12440 | 723 | fuseErr2Nominal(); |
| NaotoMorita | 61:5e5c4fe12440 | 724 | } |
| NaotoMorita | 73:5770a0d470c0 | 725 | float q0 = qhat(1,1); |
| NaotoMorita | 73:5770a0d470c0 | 726 | float q1 = qhat(2,1); |
| NaotoMorita | 73:5770a0d470c0 | 727 | float q2 = qhat(3,1); |
| NaotoMorita | 73:5770a0d470c0 | 728 | float q3 = qhat(4,1); |
| NaotoMorita | 73:5770a0d470c0 | 729 | |
| NaotoMorita | 73:5770a0d470c0 | 730 | float d0 = (-q3*q3-q2*q2+q1*q1+q0*q0); |
| NaotoMorita | 73:5770a0d470c0 | 731 | float q0q3q1q2 = (q0*q3+q1*q2); |
| NaotoMorita | 73:5770a0d470c0 | 732 | float h1lower = 2.0f*(4.0f*q0q3q1q2*q0q3q1q2/d0/d0+1.0f)*sqrt(4.0f*q0q3q1q2*q0q3q1q2/d0/d0+1.0f); |
| NaotoMorita | 73:5770a0d470c0 | 733 | |
| NaotoMorita | 73:5770a0d470c0 | 734 | float d1 = d0*sqrt(4.0f*q0q3q1q2*q0q3q1q2/d0/d0+1.0f); |
| NaotoMorita | 73:5770a0d470c0 | 735 | float d2 = d0*d0*sqrt(4.0f*q0q3q1q2*q0q3q1q2/d0/d0+1.0f); |
| NaotoMorita | 73:5770a0d470c0 | 736 | float d3 = d0*sqrt(4.0f*q0q3q1q2*q0q3q1q2/d0/d0+1.0f)*(4.0f*q0q3q1q2*q0q3q1q2/d0/d0+1.0f); |
| NaotoMorita | 73:5770a0d470c0 | 737 | |
| NaotoMorita | 73:5770a0d470c0 | 738 | |
| NaotoMorita | 73:5770a0d470c0 | 739 | |
| NaotoMorita | 73:5770a0d470c0 | 740 | Matrix Hh(2,4); |
| NaotoMorita | 73:5770a0d470c0 | 741 | Hh(1,1) = -(8.0f*q3*q0q3q1q2/d0/d0-16.0f*q0*q0q3q1q2*q0q3q1q2/d0/d0/d0)/h1lower; |
| NaotoMorita | 73:5770a0d470c0 | 742 | Hh(1,2) = -(8.0f*q2*q0q3q1q2/d0/d0-16.0f*q1*q0q3q1q2*q0q3q1q2/d0/d0/d0)/h1lower; |
| NaotoMorita | 73:5770a0d470c0 | 743 | Hh(1,3) = -(8.0f*q1*q0q3q1q2/d0/d0-16.0f*q2*q0q3q1q2*q0q3q1q2/d0/d0/d0)/h1lower; |
| NaotoMorita | 73:5770a0d470c0 | 744 | Hh(1,4) = -(8.0f*q0*q0q3q1q2/d0/d0-16.0f*q3*q0q3q1q2*q0q3q1q2/d0/d0/d0)/h1lower; |
| NaotoMorita | 73:5770a0d470c0 | 745 | |
| NaotoMorita | 73:5770a0d470c0 | 746 | Hh(2,1) = 2.0f*q3/d1-4.0f*q0*q0q3q1q2/d2-q0q3q1q2*(8.0f*q3*q0q3q1q2/d0/d0-16.0f*q0*q0q3q1q2*q0q3q1q2/d0/d0/d0)/d3; |
| NaotoMorita | 73:5770a0d470c0 | 747 | Hh(2,2) = 2.0f*q2/d1-4.0f*q1*q0q3q1q2/d2-q0q3q1q2*(8.0f*q2*q0q3q1q2/d0/d0-16.0f*q1*q0q3q1q2*q0q3q1q2/d0/d0/d0)/d3; |
| NaotoMorita | 73:5770a0d470c0 | 748 | Hh(2,3) = 2.0f*q1/d1+4.0f*q2*q0q3q1q2/d2-q0q3q1q2*(8.0f*q1*q0q3q1q2/d0/d0-16.0f*q2*q0q3q1q2*q0q3q1q2/d0/d0/d0)/d3; |
| NaotoMorita | 73:5770a0d470c0 | 749 | Hh(2,4) = 2.0f*q0/d1+4.0f*q3*q0q3q1q2/d2-q0q3q1q2*(8.0f*q0*q0q3q1q2/d0/d0-16.0f*q3*q0q3q1q2*q0q3q1q2/d0/d0/d0)/d3; |
| NaotoMorita | 56:c10f1168bd4a | 750 | */ |