Rev 973 | Details | Compare with Previous | Last modification | View Log | RSS feed
Rev | Author | Line No. | Line |
---|---|---|---|
966 | MikeW | 1 | /* |
2 | Copyright 2008, by Michael Walter |
||
3 | |||
4 | All functions written by Michael Walter are free software and can be redistributed and/or modified under the terms of the GNU Lesser |
||
5 | General Public License as published by the Free Software Foundation. This program is distributed in the hope that it will be useful, but |
||
6 | WITHOUT ANY WARRANTY, without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. |
||
7 | See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public |
||
8 | License along with this program. If not, see <http://www.gnu.org/licenses/>. |
||
9 | |||
10 | Please note: The software is based on the framework provided by H. Buss and I. Busker in their Mikrokopter projekt. All functions that |
||
11 | are not written by Michael Walter are under the license by H. Buss and I. Busker (license_buss.txt) published by www.mikrokopter.de |
||
12 | unless it is stated otherwise. |
||
13 | */ |
||
14 | |||
838 | MikeW | 15 | /***************************************************************************** |
16 | INCLUDES |
||
17 | **************************************************************************** */ |
||
18 | #include "kafi.h" |
||
966 | MikeW | 19 | #include "analog.h" |
20 | #include "main.h" |
||
21 | #include "FlightControl.h" |
||
22 | #include "mm3.h" |
||
838 | MikeW | 23 | |
24 | /***************************************************************************** |
||
25 | (SYMBOLIC) CONSTANTS |
||
26 | *****************************************************************************/ |
||
27 | |||
966 | MikeW | 28 | /***************************************************************************** |
29 | VARIABLES |
||
30 | *****************************************************************************/ |
||
31 | volatile int KompassValue = 0; |
||
838 | MikeW | 32 | |
966 | MikeW | 33 | /* Kalman filter */ |
34 | f32_t Phi, Theta, sinPhi, cosPhi, tanTheta, secTheta; |
||
35 | f32_t sinPhi, cosPhi, sinTheta, cosTheta; |
||
36 | f32_t sinPhi_P, cosPhi_P, sinTheta_P, cosTheta_P, sinPsi_P, cosPsi_P; |
||
838 | MikeW | 37 | |
966 | MikeW | 38 | f32_t **matPs; /* Predicted system Noise */ |
39 | f32_t **matP; /* Predicted system Noise */ |
||
40 | f32_t **matP_tmp; |
||
41 | f32_t **matQ; /* system noise covariance */ |
||
42 | f32_t **matR; /* measurement noise covariance */ |
||
43 | f32_t **matXd; /* estimated state */ |
||
44 | f32_t **matXs; /* predicted state */ |
||
45 | f32_t **matX_tmp1; |
||
46 | f32_t **matX_tmp2; |
||
47 | |||
48 | f32_t **matPhi; /* transition matrix to predict new state from current state */ |
||
49 | f32_t **matB ; /* control matrix to predict new state from ext. control vector U */ |
||
50 | f32_t **matU; /* control vector */ |
||
51 | f32_t **matC; /* Jacobi matrix */ |
||
52 | f32_t **matY; /* real measurements */ |
||
53 | f32_t **matYs; /* predicted measurements */ |
||
54 | f32_t **matdY; /* innovation */ |
||
55 | f32_t **temp_meas_meas; |
||
56 | f32_t **temp_meas_state; |
||
57 | f32_t **temp_state_meas; |
||
58 | f32_t **temp_meas_2; |
||
59 | f32_t **temp_state_state; |
||
60 | f32_t **matK; /* the Kalman gain matrix */ |
||
61 | |||
62 | char fYValid[KAFI_DIM_Y]; /* measurement validity flag */ |
||
63 | |||
64 | |||
65 | |||
66 | /* **************************************************************************** |
||
67 | Functionname: AttitudeEstimation */ /*! |
||
68 | Description: Main Kalman Filter Loop |
||
69 | |||
70 | @param[in] |
||
71 | |||
72 | @return void |
||
73 | @pre - |
||
74 | @post - |
||
75 | @author Michael Walter |
||
76 | **************************************************************************** */ |
||
77 | void FlightAttitudeEstimation() |
||
78 | { |
||
79 | /* Predict the Measurent */ |
||
80 | trPredictMeasurement(); |
||
81 | |||
82 | /* Update the Jacobi Matrix */ |
||
83 | trUpdateJacobiMatrix(); |
||
84 | |||
85 | /* Extract measurements and store them in vector matY. |
||
86 | Measurement validity is indicated by fYValid[] */ |
||
87 | trMeasure(); |
||
88 | |||
89 | /* Innovation: calculate Xd, and Pd */ |
||
90 | trInnovate(); |
||
91 | |||
92 | /* Update transition matrix Phi */ |
||
93 | //trUpdatePhi(); // Is an Identity matrix which will not change |
||
94 | |||
95 | /* Update control matrix B */ |
||
96 | trUpdateBU(); |
||
97 | |||
98 | /* Predict new state Xs and Ps */ |
||
99 | KAFIPrediction(); |
||
100 | |||
101 | /* store innovated state vector in current state */ |
||
102 | matGetFull(matXs, _Phi , 0, &status.Phi); |
||
103 | matGetFull(matXs, _Theta , 0, &status.Theta); |
||
104 | matGetFull(matXs, _Psi , 0, &status.Psi); |
||
105 | |||
106 | /* Limit Angles to [-2Pi;...;+2Pi] */ |
||
107 | trLimitAngles(); |
||
108 | } |
||
109 | |||
110 | /* **************************************************************************** |
||
111 | Functionname: Kafi_Init */ /*! |
||
112 | Description: Generate Kalman Matrizes |
||
113 | |||
114 | @param[in] |
||
115 | |||
116 | @return void |
||
117 | @pre - |
||
118 | @post - |
||
119 | @author Michael Walter |
||
120 | **************************************************************************** */ |
||
121 | void Kafi_Init() |
||
122 | { |
||
123 | /* create kalman filter and associated matrizes */ |
||
124 | ui32_t i, j; |
||
125 | matPs = matrix( 1, KAFI_DIM_X, 1, KAFI_DIM_X ); |
||
126 | matQ = matrix( 1, KAFI_DIM_X, 1,KAFI_DIM_X); |
||
127 | matR = matrix( 1, KAFI_DIM_Y, 1,KAFI_DIM_Y); |
||
128 | matXd = matrix( 1, KAFI_DIM_X, 1,1); /* estimated state */ |
||
129 | matXs = matrix( 1, KAFI_DIM_X, 1,1); /* predicted state */ |
||
130 | matX_tmp1 = matrix( 1, KAFI_DIM_X,1, 1); |
||
131 | matX_tmp2 = matrix( 1, KAFI_DIM_X, 1,1); |
||
132 | matPhi = matrix( 1, KAFI_DIM_X, 1,KAFI_DIM_X);/* transition matrix to predict new state from current state */ |
||
133 | matB = matrix( 1, KAFI_DIM_X, 1,KAFI_DIM_U); /* control matrix to predict new state from ext. control vector U */ |
||
134 | matU = matrix( 1, KAFI_DIM_U, 1,1); /* control vector */ |
||
135 | matC = matrix( 1, KAFI_DIM_Y, 1,KAFI_DIM_X); /* Jacobi matrix */ |
||
136 | matY = matrix( 1, KAFI_DIM_Y, 1,1); /* real measurements */ |
||
137 | matdY = matrix( 1, KAFI_DIM_Y, 1,1); /* innovation */ |
||
138 | matYs = matrix( 1, KAFI_DIM_Y, 1,1); /* predicted measurements */ |
||
139 | matP = matrix( 1, KAFI_DIM_X, 1,KAFI_DIM_X); |
||
140 | matK = matrix(1, KAFI_DIM_X, 1, KAFI_DIM_Y); |
||
141 | temp_state_state = matrix( 1, KAFI_DIM_X, 1, KAFI_DIM_X ); |
||
142 | temp_meas_state = matrix( 1, KAFI_DIM_Y, 1, KAFI_DIM_X ); |
||
143 | temp_state_meas = matrix( 1, KAFI_DIM_X, 1, KAFI_DIM_Y ); |
||
144 | temp_meas_meas = matrix( 1, KAFI_DIM_Y, 1, KAFI_DIM_Y ); |
||
145 | temp_meas_2 = matrix( 1, KAFI_DIM_Y, 1, KAFI_DIM_Y ); |
||
146 | matP_tmp = matrix( 1, KAFI_DIM_X, 1, KAFI_DIM_X); |
||
147 | |||
148 | /* set phi to identity matrix */ |
||
149 | for ( j = 0; j < KAFI_DIM_X; j++ ) |
||
150 | { |
||
151 | for ( i = 0; i < KAFI_DIM_X; i++ ) |
||
152 | { |
||
153 | matSetFull(matPhi, j, i, (j == i) ? 1.F : 0.F); |
||
154 | } |
||
155 | } |
||
156 | |||
157 | /* set diagonal of measurement noise covariance R */ |
||
158 | matSetDiagonal(matR, _ax, _ax, 10.0F); |
||
159 | matSetDiagonal(matR, _ay, _ay, 10.0F); |
||
160 | matSetDiagonal(matR, _az, _az, 10.0F); |
||
161 | matSetDiagonal(matR, _compass, _compass, 5.0F); |
||
162 | |||
163 | #ifdef MELEXIS_GYRO |
||
164 | /* set diagonal of system noise covariance Q, dim(X,X) */ |
||
165 | matSetDiagonal(matQ, _Phi, _Phi, 0.000006F); |
||
166 | matSetDiagonal(matQ, _Theta, _Theta, 0.000006F); |
||
167 | matSetDiagonal(matQ, _Psi, _Psi, 0.000006F); |
||
168 | #else |
||
169 | /* set diagonal of system noise covariance Q */ |
||
170 | matSetDiagonal(matQ, _Phi, _Phi, 0.00001F); |
||
171 | matSetDiagonal(matQ, _Theta, _Theta, 0.00001F); |
||
172 | matSetDiagonal(matQ, _Psi, _Psi, 0.00005F); |
||
173 | #endif |
||
174 | |||
982 | MikeW | 175 | matSetDiagonal(matQ, _Phi, _Phi, 0.000006F); |
176 | matSetDiagonal(matQ, _Theta, _Theta, 0.000006F); |
||
177 | matSetDiagonal(matQ, _Psi, _Psi, 0.000006F); |
||
178 | |||
966 | MikeW | 179 | /* set diagonal init. estimation error (covariance) P0 */ |
180 | matSetDiagonal(matPs, _Phi, _Phi, 0.00001F); |
||
181 | matSetDiagonal(matPs, _Theta, _Theta, 0.00001F); |
||
182 | matSetDiagonal(matPs, _Psi, _Psi, 0.00001F); |
||
183 | |||
184 | matSetFull(matXd, _Phi, 0, 0.0F); |
||
185 | matSetFull(matXd, _Theta, 0, 0.0F); |
||
186 | matSetFull(matXd, _Psi, 0, 0.0F); |
||
187 | |||
188 | matSetFull(matXs, _Phi, 0, 0.0F); |
||
189 | matSetFull(matXs, _Theta, 0, 0.0F); |
||
190 | matSetFull(matXs, _Psi, 0, 0.0F); |
||
191 | } |
||
192 | |||
193 | /* **************************************************************************** |
||
194 | Functionname: KAFIPrediction */ /*! |
||
195 | Description: Predict new state Xs and Ps |
||
196 | |||
197 | X_k = A * X_(k-1) + B * u_k |
||
198 | P_k = A * P_(k-1) * A^T + Q |
||
199 | |||
200 | @param[in] |
||
201 | |||
202 | @return void |
||
203 | @pre - |
||
204 | @post - |
||
205 | @author Michael Walter |
||
206 | **************************************************************************** */ |
||
207 | void KAFIPrediction() |
||
208 | { |
||
209 | /*--- (SYMBOLIC) CONSTANTS ---*/ |
||
210 | |||
211 | /*--- VARIABLES ---*/ |
||
212 | |||
213 | /* Innovation: calculate Xd, and Pd */ |
||
214 | |||
215 | /* X_k = A * X_(k-1) + B * u_k */ |
||
216 | |||
217 | /* matX_tmp1 = matPhi * matXd */ |
||
218 | /* |
||
219 | matrix_mult(matPhi,matXd, matX_tmp1, |
||
220 | KAFI_DIM_X, KAFI_DIM_X, 1); |
||
221 | */ |
||
222 | /* matX_tmp2 = matB * matU */ |
||
223 | matrix_mult( matB,matU, matX_tmp2, |
||
224 | KAFI_DIM_X, KAFI_DIM_U, 1); |
||
225 | |||
226 | /* matXs = matPhi * matXd + matB * matU */ |
||
227 | /* |
||
228 | matrix_add( matX_tmp1, matX_tmp2, matXs, KAFI_DIM_X, 1 ); |
||
229 | matrix_add( matXd, matX_tmp2, matXs, KAFI_DIM_X, 1 ); |
||
230 | */ |
||
231 | |||
232 | matXs[1][1] = matXd[1][1] + matX_tmp2[1][1]; |
||
233 | matXs[2][1] = matXd[2][1] + matX_tmp2[2][1]; |
||
234 | matXs[3][1] = matXd[3][1] + matX_tmp2[3][1]; |
||
235 | |||
236 | /*P_k = A * P_(k-1)*A^T + Q */ |
||
237 | /* |
||
238 | matrix_mult_transpose( matP, matPhi, matP_tmp, // Can be skiped PHI is a ident matrix |
||
239 | KAFI_DIM_X, KAFI_DIM_X, KAFI_DIM_X ); |
||
240 | matrix_mult( matPhi, matP_tmp, matP, |
||
241 | KAFI_DIM_X, KAFI_DIM_X, KAFI_DIM_X ); |
||
242 | matrix_add( matP, matQ, matPs, KAFI_DIM_X, KAFI_DIM_X ); |
||
243 | */ |
||
244 | |||
245 | matPs[1][1] = matP[1][1] + matQ[1][1]; |
||
246 | matPs[2][2] = matP[2][2] + matQ[2][2]; |
||
247 | matPs[3][3] = matP[3][3] + matQ[3][3]; |
||
248 | } |
||
249 | |||
250 | /* **************************************************************************** |
||
251 | Functionname: trInnovate */ /*! |
||
252 | Description: Innovation: calculate Xd, and Pd |
||
253 | |||
254 | K = Ps * C^T * (C * Ps * C^T + R)^-1 |
||
255 | X = Xs + K (Y - Ys) |
||
256 | |||
257 | @param[in] |
||
258 | |||
259 | @return void |
||
260 | @pre - |
||
261 | @post - |
||
262 | @author Michael Walter |
||
263 | **************************************************************************** */ |
||
264 | void trInnovate() |
||
265 | { |
||
266 | /*--- (SYMBOLIC) CONSTANTS ---*/ |
||
267 | |||
268 | /*--- VARIABLES ---*/ |
||
269 | #if 0 /* Runtime optimised. To be verified */ |
||
270 | /* temp_meas_state = matC * matPs */ |
||
271 | matrix_mult( matC, matPs, temp_meas_state, |
||
272 | KAFI_DIM_Y, KAFI_DIM_X, KAFI_DIM_X ); |
||
273 | |||
274 | /* temp_meas_meas = matC * matPs^T * matC^T */ |
||
275 | matrix_mult_transpose( matC, temp_meas_state, temp_meas_meas, |
||
276 | KAFI_DIM_Y, KAFI_DIM_X, KAFI_DIM_Y ); |
||
277 | |||
278 | /* temp_meas_meas = matC * matPs^T * matC^T + matR */ |
||
279 | /* |
||
280 | matrix_add( temp_meas_meas, matR, temp_meas_meas, |
||
281 | KAFI_DIM_Y, KAFI_DIM_Y ); |
||
282 | */ |
||
283 | temp_meas_meas[1][1] += matR[1][1]; |
||
284 | temp_meas_meas[2][2] += matR[2][2]; |
||
285 | temp_meas_meas[3][3] += matR[3][3]; |
||
286 | temp_meas_meas[4][4] += matR[4][4]; |
||
287 | |||
288 | /* temp_meas_2 = (matC * matPs^T * matC^T + matR)^-1 */ |
||
289 | take_inverse( temp_meas_meas, temp_meas_2, KAFI_DIM_Y ); |
||
290 | |||
291 | /* matK = matPs^T * matC^T * (matC * matPs^T * matC^T + matR)^-1 */ |
||
292 | matrix_transpose_mult( temp_meas_state, temp_meas_2, matK, |
||
293 | KAFI_DIM_X, KAFI_DIM_Y, KAFI_DIM_Y ); |
||
294 | |||
295 | /* temp_state_state = matK * matC * matPs */ |
||
296 | matrix_mult( matK, temp_meas_state, temp_state_state, |
||
297 | KAFI_DIM_X, KAFI_DIM_Y, KAFI_DIM_X ); |
||
298 | |||
299 | /* matP = matPs - matK * matC * matPs */ |
||
300 | matrix_sub( matPs, temp_state_state, matP, KAFI_DIM_X, KAFI_DIM_X ); |
||
301 | #else |
||
302 | /* temp_state_meas = matPs * matC^T */ |
||
303 | matrix_mult_transpose( matPs, matC, temp_state_meas, |
||
304 | KAFI_DIM_X, KAFI_DIM_X, KAFI_DIM_Y ); |
||
305 | |||
306 | /* temp_meas_meas = matC * matPs * matC^T */ |
||
307 | matrix_mult( matC, temp_state_meas, temp_meas_meas, |
||
308 | KAFI_DIM_Y, KAFI_DIM_X, KAFI_DIM_Y ); |
||
309 | |||
310 | /* temp_meas_meas = matC * matPs * matC^T + matR */ |
||
311 | /* |
||
312 | matrix_add( temp_meas_meas, matR, temp_meas_meas, |
||
313 | KAFI_DIM_Y, KAFI_DIM_Y ); |
||
314 | */ |
||
315 | temp_meas_meas[1][1] += matR[1][1]; |
||
316 | temp_meas_meas[2][2] += matR[2][2]; |
||
317 | temp_meas_meas[3][3] += matR[3][3]; |
||
318 | temp_meas_meas[4][4] += matR[4][4]; |
||
319 | |||
320 | /* temp_meas_2 = (matC * matPs * matC^T + matR)^-1 */ |
||
321 | take_inverse( temp_meas_meas, temp_meas_2, KAFI_DIM_Y ); |
||
322 | |||
323 | /* matK = matPs * matC^T * (matC * matPs^T * matC^T + matR)^-1 */ |
||
324 | matrix_mult( temp_state_meas, temp_meas_2, matK, |
||
325 | KAFI_DIM_X, KAFI_DIM_Y, KAFI_DIM_Y ); |
||
326 | |||
327 | /* temp_meas_state = matC * matPs*/ |
||
328 | matrix_mult( matC, matPs, temp_meas_state, |
||
329 | KAFI_DIM_Y, KAFI_DIM_X, KAFI_DIM_X ); |
||
330 | |||
331 | /* temp_state_state = matK * matC * matPs */ |
||
332 | matrix_mult( matK, temp_meas_state, temp_state_state, |
||
333 | KAFI_DIM_X, KAFI_DIM_Y, KAFI_DIM_X ); |
||
334 | |||
335 | /* matP = matPs - matK * matC * matPs */ |
||
336 | matrix_sub( matPs, temp_state_state, matP, KAFI_DIM_X, KAFI_DIM_X ); |
||
337 | #endif |
||
338 | |||
339 | /* matdY = matY - matYs */ |
||
340 | /* |
||
341 | matrix_sub( matY, matYs, matdY, |
||
342 | KAFI_DIM_Y, 1 ); |
||
343 | */ |
||
344 | matdY[1][1] = matY[1][1] - matYs[1][1]; |
||
345 | matdY[2][1] = matY[2][1] - matYs[2][1]; |
||
346 | matdY[3][1] = matY[3][1] - matYs[3][1]; |
||
347 | |||
348 | /* We will only update the Heading if we got a valid compass signal */ |
||
349 | if(fYValid[_compass] == 1) |
||
350 | { |
||
351 | matdY[4][1] = matY[4][1] - matYs[4][1]; |
||
352 | } |
||
353 | else |
||
354 | { |
||
355 | matdY[4][1] = 0.0F; |
||
356 | } |
||
357 | |||
358 | /* matX_tmp1 = matK * (Y -Ys) */ |
||
359 | matrix_mult( matK, matdY, matX_tmp1, |
||
360 | KAFI_DIM_X, KAFI_DIM_Y, 1 ); |
||
361 | |||
362 | /* matXd = matXs + matK * (Y -Ys) */ |
||
363 | /* |
||
364 | matrix_add( matXs, matX_tmp1, matXd, |
||
365 | KAFI_DIM_X, 1 ); |
||
366 | */ |
||
367 | matXd[1][1] = matXs[1][1] + matX_tmp1[1][1]; |
||
368 | matXd[2][1] = matXs[2][1] + matX_tmp1[2][1]; |
||
369 | matXd[3][1] = matXs[3][1] + matX_tmp1[3][1]; |
||
370 | matXd[4][1] = matXs[4][1] + matX_tmp1[4][1]; |
||
371 | } |
||
372 | |||
373 | |||
374 | /* **************************************************************************** |
||
375 | Functionname: trUpdatePhi */ /*! |
||
376 | Description: Setup matPhi, used to calculate the predicted state |
||
377 | from the current state |
||
378 | |||
379 | @param[in] |
||
380 | @param[in] fCycleTime |
||
381 | |||
382 | @return void |
||
383 | @pre - |
||
384 | @post - |
||
385 | @author Michael Walter |
||
386 | **************************************************************************** */ |
||
387 | void trUpdatePhi() |
||
388 | { |
||
389 | /*--- (SYMBOLIC) CONSTANTS ---*/ |
||
390 | |||
391 | /*--- VARIABLES ---*/ |
||
392 | matSetFull(matPhi, _Phi, _Phi, 1.0F); |
||
393 | matSetFull(matPhi, _Theta, _Theta, 1.0F); |
||
394 | matSetFull(matPhi, _Psi, _Psi, 1.0F); |
||
395 | } |
||
396 | |||
397 | /* **************************************************************************** |
||
398 | Functionname: trUpdateJacobiMatrix */ /*! |
||
399 | Description: |
||
400 | |||
401 | @param[in] void |
||
402 | |||
403 | @return void |
||
404 | @pre - |
||
405 | @post - |
||
406 | @author Michael Walter |
||
407 | *****************************************************************************/ |
||
408 | void trUpdateJacobiMatrix(void) |
||
409 | { |
||
410 | /*--- (SYMBOLIC) CONSTANTS ---*/ |
||
411 | /* |
||
412 | Simplified Version d/dt |
||
413 | Pre_ax = du - sinTheta * g; |
||
414 | Pre_ay = dv + cosTheta_sinPhi * g; |
||
415 | Pre_az = dw + cosTheta_cosPhi * g; |
||
416 | */ |
||
417 | |||
418 | /*--- VARIABLES ---*/ |
||
419 | f32_t Phi, Theta, Psi; |
||
420 | f32_t g = 9.81F; |
||
421 | |||
422 | Phi = status.Phi; |
||
423 | Theta = status.Theta; |
||
424 | Psi = status.Psi; |
||
425 | |||
426 | sinPhi = sin(Phi); |
||
427 | cosPhi = cos(Phi); |
||
428 | sinTheta = sin(Theta); |
||
429 | cosTheta = cos(Theta); |
||
430 | |||
431 | //Pre_ax = du - sinTheta * g; |
||
432 | //matSetFull(matC, _ax, _Phi , 0); |
||
433 | matSetFull(matC, _ax, _Theta, -cosTheta * g); |
||
434 | //matSetFull(matC, _ax, _Psi , 0); |
||
435 | |||
436 | //Pre_ay = dv + cosTheta_sinPhi * g; |
||
437 | matSetFull(matC, _ay, _Phi, cosTheta * cosPhi * g); |
||
438 | matSetFull(matC, _ay, _Theta , -sinTheta * sinPhi * g); |
||
439 | //matSetFull(matC, _ay, _Psi , 0); |
||
440 | |||
441 | //Pre_az = dw + cosTheta_cosPhi * g; |
||
442 | matSetFull(matC, _az, _Phi , -cosTheta * sinPhi * g); |
||
443 | matSetFull(matC, _az, _Theta , -sinTheta * cosPhi * g); |
||
444 | //matSetFull(matC, _az, _Psi , 0); |
||
445 | |||
446 | matSetFull(matC, _compass, _Psi , 1.0F); |
||
447 | } |
||
448 | |||
838 | MikeW | 449 | /***************************************************************************** |
966 | MikeW | 450 | Functionname: trPredictMeasurement */ /*! |
451 | Description:Predict Measurements: AX, AY, AZ, Compass |
||
452 | |||
453 | |||
454 | @param[in] void |
||
455 | |||
456 | @return void |
||
457 | @pre - |
||
458 | @post - |
||
459 | @author Michael Walter |
||
838 | MikeW | 460 | *****************************************************************************/ |
966 | MikeW | 461 | void trPredictMeasurement(void) |
462 | { |
||
463 | /*--- (SYMBOLIC) CONSTANTS ---*/ |
||
464 | |||
465 | /* |
||
466 | Simplified Version |
||
467 | Pre_ax = du - sinTheta * g; |
||
468 | Pre_ay = dv + cosTheta_sinPhi * g; |
||
469 | Pre_az = dw + cosTheta_cosPhi * g; |
||
470 | */ |
||
471 | |||
472 | /*--- VARIABLES ---*/ |
||
473 | f32_t g = 9.81F; |
||
474 | f32_t Pre_ax, Pre_ay, Pre_az; |
||
475 | f32_t Phi, Theta, Psi; |
||
476 | |||
477 | matGetFull(matXs, _Phi, 0, &Phi); |
||
478 | matGetFull(matXs, _Theta, 0, &Theta); |
||
479 | matGetFull(matXs, _Psi, 0, &Psi); |
||
480 | |||
481 | sinPhi_P = sin(Phi); |
||
482 | cosPhi_P = cos(Phi); |
||
483 | sinTheta_P = sin(Theta); |
||
484 | cosTheta_P = cos(Theta); |
||
485 | |||
486 | /* Simplified Version */ |
||
487 | Pre_ax = -sinTheta_P * g; |
||
488 | Pre_ay = cosTheta_P * sinPhi_P * g; |
||
489 | Pre_az = cosTheta_P * cosPhi_P * g; |
||
490 | |||
491 | matSetFull(matYs, _ax, 0, Pre_ax); |
||
492 | matSetFull(matYs, _ay, 0, Pre_ay); |
||
493 | matSetFull(matYs, _az, 0, Pre_az); |
||
494 | matSetFull(matYs, _compass, 0, Psi); |
||
495 | } |
||
496 | |||
497 | /* **************************************************************************** |
||
498 | Functionname: trMeasure */ /*! |
||
499 | Description: |
||
500 | |||
501 | @param[in] |
||
502 | |||
503 | @return void |
||
504 | @pre - |
||
505 | @post - |
||
506 | @author Michael Walter |
||
507 | **************************************************************************** */ |
||
508 | void trMeasure() |
||
509 | { |
||
510 | /*--- (SYMBOLIC) CONSTANTS ---*/ |
||
511 | |||
512 | /*--- VARIABLES ---*/ |
||
513 | #ifdef USE_COMPASS |
||
514 | static int updCompass = 10; |
||
515 | f32_t Psi, Compass; |
||
516 | #endif |
||
517 | |||
518 | f32_t ADC_ax = 0.0F; |
||
519 | f32_t ADC_ay = 0.0F; |
||
520 | f32_t ADC_az = 0.0F; |
||
521 | |||
522 | fYValid[_ax] = 1; |
||
523 | fYValid[_ay] = 1; |
||
524 | fYValid[_az] = 1; |
||
525 | |||
526 | #ifdef INTERNAL_REFERENCE |
||
527 | ADC_ax = AdWertAccNick * 0.044F; |
||
528 | ADC_ay = -AdWertAccRoll * 0.044F; |
||
529 | ADC_az = AdWertAccHoch * 0.044F; |
||
530 | #else |
||
531 | ADC_ax = AdWertAccNick * 0.047F; |
||
532 | ADC_ay = -AdWertAccRoll * 0.047F; |
||
533 | ADC_az = AdWertAccHoch * 0.047F; |
||
534 | #endif |
||
535 | |||
536 | matSetFull(matY, _ax, 0, ADC_ax); |
||
537 | matSetFull(matY, _ay, 0, ADC_ay); |
||
538 | matSetFull(matY, _az, 0, ADC_az); |
||
539 | |||
540 | #ifdef USE_COMPASS |
||
541 | if (!updCompass--) |
||
542 | { |
||
543 | updCompass = 10; // update only at 2ms*50 = 100ms (10Hz) |
||
544 | KompassValue = MM3_Heading(); |
||
545 | } |
||
546 | Compass = KompassValue / 180.F * PI; |
||
547 | |||
548 | matGetFull(matXs, _Psi, 0, &Psi); |
||
549 | |||
550 | if (fabs(Compass + 2*PI - Psi) < fabs(Compass - Psi)) |
||
551 | { |
||
552 | Compass += 2 * PI; |
||
553 | } |
||
554 | else if (fabs(Compass - 2*PI - Psi) < fabs(Compass - Psi)) |
||
555 | { |
||
556 | Compass -= 2 * PI; |
||
557 | } |
||
558 | |||
559 | matSetFull(matY, _compass, 0, Compass); |
||
560 | |||
561 | /* Roll and Yaw angle are smaller 8 Deg*/ |
||
562 | if ((abs(status.iPhi10) < 80 ) && (abs(status.iTheta10) < 80)) |
||
563 | { |
||
564 | fYValid[_compass] = 1; |
||
565 | } |
||
566 | else |
||
567 | { |
||
568 | fYValid[_compass] = 0; |
||
569 | } |
||
570 | #else |
||
571 | fYValid[_compass] = 0; |
||
572 | #endif |
||
573 | } |
||
574 | |||
575 | /* **************************************************************************** |
||
576 | Functionname: trUpdateBU */ /*! |
||
577 | Description: Update control matrix B and vector U |
||
578 | |||
579 | |||
580 | @param[in] |
||
581 | |||
582 | @return void |
||
583 | @pre - |
||
584 | @post - |
||
585 | @author Michael Walter |
||
586 | **************************************************************************** */ |
||
587 | void trUpdateBU() |
||
588 | { |
||
589 | /*--- (SYMBOLIC) CONSTANTS ---*/ |
||
590 | |||
591 | /*--- VARIABLES ---*/ |
||
592 | /* |
||
593 | dPhi | 1 sinPhi_tanTheta cosPhi_tanTheta | p |
||
594 | dTheta = | 0 cosPhi -sinPhi | * q |
||
595 | dPsi | 0 sinPhi_secTheta cosPhi_secTheta | w |
||
596 | */ |
||
597 | |||
598 | static f32_t fSumGier = 0.0F; |
||
599 | static f32_t fSumNick = 0.0F; |
||
600 | static f32_t fSumRoll = 0.0F; |
||
601 | |||
602 | f32_t ADC_p = 0.0F; |
||
603 | f32_t ADC_q = 0.0F; |
||
604 | f32_t ADC_r = 0.0F; |
||
605 | |||
606 | /* store predicted state vector in current state */ |
||
607 | matGetFull(matXd, _Phi , 0, &Phi); |
||
608 | matGetFull(matXd, _Theta , 0, &Theta); |
||
609 | |||
610 | sinPhi = sin(Phi); |
||
611 | cosPhi = cos(Phi); |
||
612 | |||
613 | tanTheta = tan(Theta); |
||
614 | secTheta = 1.0F / cos(Theta); |
||
615 | |||
616 | matSetFull(matB, _Phi, _p, 1.0F); |
||
617 | matSetFull(matB, _Phi, _q, sinPhi * tanTheta ); |
||
618 | matSetFull(matB, _Phi, _r, cosPhi * tanTheta ); |
||
619 | |||
620 | //matSetFull(matB, _Theta, _p, 0.0F); |
||
621 | matSetFull(matB, _Theta, _q, cosPhi ); |
||
622 | matSetFull(matB, _Theta, _r, -sinPhi ); |
||
623 | |||
624 | //matSetFull(matB, _Psi, _p, 0.0F); |
||
625 | matSetFull(matB, _Psi, _q, sinPhi * secTheta ); |
||
626 | matSetFull(matB, _Psi, _r, cosPhi * secTheta ); |
||
627 | |||
973 | MikeW | 628 | ADC_p = -(float) AverageRoll * 0.00001F * fCycleTime; |
629 | ADC_q = -(float) AverageNick * 0.00001F * fCycleTime; |
||
630 | ADC_r = -(float) AverageGier * 0.00001F * fCycleTime; |
||
966 | MikeW | 631 | |
632 | fSumGier += ADC_r * 180.F / 3.14F; |
||
633 | fSumNick += ADC_q * 180.F / 3.14F; |
||
634 | fSumRoll += ADC_p * 180.F / 3.14F; |
||
635 | |||
636 | DebugOut.Analog[3] = fSumNick; |
||
637 | DebugOut.Analog[4] = fSumRoll; |
||
638 | DebugOut.Analog[5] = fSumGier; |
||
639 | |||
640 | matSetFull(matU, _p, 0, ADC_p ); |
||
641 | matSetFull(matU, _q, 0, ADC_q); |
||
642 | matSetFull(matU, _r, 0, ADC_r); |
||
643 | } |
||
644 | |||
645 | |||
646 | |||
647 | /* **************************************************************************** |
||
648 | Functionname: trLimitAngles */ /*! |
||
649 | Description: |
||
650 | |||
651 | @param[in] |
||
652 | |||
653 | @return void |
||
654 | @pre - |
||
655 | @post - |
||
656 | @author |
||
657 | **************************************************************************** */ |
||
658 | void trLimitAngles() |
||
659 | { |
||
660 | /*--- (SYMBOLIC) CONSTANTS ---*/ |
||
661 | |||
662 | /*--- VARIABLES ---*/ |
||
663 | f32_t Psi; |
||
664 | |||
665 | if (status.Phi > 2.0F * PI) |
||
666 | { |
||
667 | status.Phi -= 2.0F * PI; |
||
668 | matSetFull(matXs, _Phi, 0, status.Phi); |
||
669 | } |
||
670 | if (status.Phi < -2.0F * PI) |
||
671 | { |
||
672 | status.Phi += 2.0F * PI; |
||
673 | matSetFull(matXs, _Phi, 0, status.Phi); |
||
674 | } |
||
675 | if (status.Theta > 2.0F * PI) |
||
676 | { |
||
677 | status.Theta -= 2.0F * PI; |
||
678 | matSetFull(matXs, _Theta, 0, status.Theta); |
||
679 | } |
||
680 | if (status.Theta < -2.0F * PI) |
||
681 | { |
||
682 | status.Theta += 2.0F * PI; |
||
683 | matSetFull(matXs, _Theta, 0, status.Theta); |
||
684 | } |
||
685 | if (status.Psi > 2.0F * PI) |
||
686 | { |
||
687 | status.Psi -= 2.0F * PI; |
||
688 | sollGier -= 3600; |
||
689 | |||
690 | matGetFull(matXs, _Psi, 0, &Psi); |
||
691 | matSetFull(matXs, _Psi , 0, Psi - 2.0F * PI); |
||
692 | matGetFull(matXd, _Psi, 0, &Psi); |
||
693 | matSetFull(matXd, _Psi , 0, Psi - 2.0F * PI); |
||
694 | } |
||
695 | if (status.Psi < 0.0F * PI) |
||
696 | { |
||
697 | status.Psi += 2.0F * PI; |
||
698 | sollGier += 3600; |
||
699 | |||
700 | matGetFull(matXs, _Psi, 0, &Psi); |
||
701 | matSetFull(matXs, _Psi , 0, Psi + 2.0F * PI); |
||
702 | matGetFull(matXd, _Psi, 0, &Psi); |
||
703 | matSetFull(matXd, _Psi , 0, Psi + 2.0F * PI); |
||
704 | } |
||
705 | |||
706 | status.iPhi10 = (int) (status.Phi * 573.0F); |
||
707 | status.iTheta10 = (int) (status.Theta * 573.0F); |
||
708 | status.iPsi10 = (int) (status.Psi * 573.0F); |
||
709 | |||
710 | if ((sollGier - status.iPsi10) > 3600) |
||
711 | { |
||
712 | sollGier -= 3600; |
||
713 | } |
||
714 | |||
715 | if ((sollGier - status.iPsi10) < -3600) |
||
716 | { |
||
717 | sollGier += 3600; |
||
718 | } |
||
719 | } |
||
720 | |||
721 | |||
722 | /* **************************************************************************** |
||
723 | Functionname: trEstimateVelocity */ /*! |
||
724 | Description: |
||
725 | |||
726 | @param[in] |
||
727 | |||
728 | @return void |
||
729 | @pre - |
||
730 | @post - |
||
731 | @author |
||
732 | **************************************************************************** */ |
||
733 | void trEstimateVelocity() |
||
734 | { |
||
735 | /*--- (SYMBOLIC) CONSTANTS ---*/ |
||
736 | |||
737 | /*--- VARIABLES ---*/ |
||
738 | } |