Subversion Repositories FlightCtrl

Rev

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
}