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