git.gag.com Git - fw/altos/blob - src/kernel/ao_kalman.c

   1 /*
   2  * Copyright © 2011 Keith Packard <keithp@keithp.com>
   3  *
   4  * This program is free software; you can redistribute it and/or modify
   5  * it under the terms of the GNU General Public License as published by
   6  * the Free Software Foundation; either version 2 of the License, or
   7  * (at your option) any later version.
   8  *
   9  * This program is distributed in the hope that it will be useful, but
  10  * WITHOUT ANY WARRANTY; without even the implied warranty of
  11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12  * General Public License for more details.
  13  *
  14  * You should have received a copy of the GNU General Public License along
  15  * with this program; if not, write to the Free Software Foundation, Inc.,
  16  * 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA.
  17  */
  18
  19 #ifndef AO_FLIGHT_TEST
  20 #include "ao.h"
  21 #include "ao_flight.h"
  22 #endif
  23
  24 #include "ao_sample.h"
  25 #include "ao_kalman.h"
  26
  27 static ao_k_t           ao_k_height;
  28 static ao_k_t           ao_k_speed;
  29 static ao_k_t           ao_k_accel;
  30
  31 #define AO_K_STEP_100           to_fix_v(0.01)
  32 #define AO_K_STEP_2_2_100       to_fix_v(0.00005)
  33
  34 #define AO_K_STEP_10            to_fix_v(0.1)
  35 #define AO_K_STEP_2_2_10        to_fix_v(0.005)
  36
  37 #define AO_K_STEP_1             to_fix_v(1)
  38 #define AO_K_STEP_2_2_1         to_fix_v(0.5)
  39
  40 ao_v_t                  ao_height;
  41 ao_v_t                  ao_speed;
  42 ao_v_t                  ao_accel;
  43 ao_v_t                  ao_max_height;
  44 static ao_k_t           ao_avg_height_scaled;
  45 ao_v_t                  ao_avg_height;
  46
  47 ao_v_t                  ao_error_h;
  48 #if !HAS_ACCEL || AO_FLIGHT_TEST
  49 #define AO_ERROR_H_SQ_AVG       1
  50 #endif
  51
  52 #if AO_ERROR_H_SQ_AVG
  53 ao_v_t                  ao_error_h_sq_avg;
  54 #endif
  55
  56 #if HAS_ACCEL
  57 ao_v_t                  ao_error_a;
  58 #endif
  59
  60 static void
  61 ao_kalman_predict(void)
  62 {
  63 #ifdef AO_FLIGHT_TEST
  64         if ((int16_t) (ao_sample_tick - ao_sample_prev_tick) > 50) {
  65                 ao_k_height += ((ao_k_t) ao_speed * AO_K_STEP_1 +
  66                                 (ao_k_t) ao_accel * AO_K_STEP_2_2_1) >> 4;
  67                 ao_k_speed += (ao_k_t) ao_accel * AO_K_STEP_1;
  68
  69                 return;
  70         }
  71         if ((int16_t) (ao_sample_tick - ao_sample_prev_tick) > 5) {
  72                 ao_k_height += ((ao_k_t) ao_speed * AO_K_STEP_10 +
  73                                 (ao_k_t) ao_accel * AO_K_STEP_2_2_10) >> 4;
  74                 ao_k_speed += (ao_k_t) ao_accel * AO_K_STEP_10;
  75
  76                 return;
  77         }
  78         if (ao_flight_debug) {
  79                 printf ("predict speed %g + (%g * %g) = %g\n",
  80                         ao_k_speed / (65536.0 * 16.0), ao_accel / 16.0, AO_K_STEP_100 / 65536.0,
  81                         (ao_k_speed + (ao_k_t) ao_accel * AO_K_STEP_100) / (65536.0 * 16.0));
  82         }
  83 #endif
  84         ao_k_height += ((ao_k_t) ao_speed * AO_K_STEP_100 +
  85                         (ao_k_t) ao_accel * AO_K_STEP_2_2_100) >> 4;
  86         ao_k_speed += (ao_k_t) ao_accel * AO_K_STEP_100;
  87 }
  88
  89 #if HAS_BARO
  90 static void
  91 ao_kalman_err_height(void)
  92 {
  93 #if AO_ERROR_H_SQ_AVG
  94         ao_v_t  e;
  95 #endif
  96         ao_v_t height_distrust;
  97 #if HAS_ACCEL
  98         ao_v_t  speed_distrust;
  99 #endif
 100
 101         ao_error_h = ao_sample_height - (ao_v_t) (ao_k_height >> 16);
 102
 103 #if AO_ERROR_H_SQ_AVG
 104         e = ao_error_h;
 105         if (e < 0)
 106                 e = -e;
 107         if (e > 127)
 108                 e = 127;
 109         ao_error_h_sq_avg -= ao_error_h_sq_avg >> 4;
 110         ao_error_h_sq_avg += (e * e) >> 4;
 111 #endif
 112
 113         if (ao_flight_state >= ao_flight_drogue)
 114                 return;
 115         height_distrust = ao_sample_alt - AO_MAX_BARO_HEIGHT;
 116 #if HAS_ACCEL
 117         /* speed is stored * 16, but we need to ramp between 248 and 328, so
 118          * we want to multiply by 2. The result is a shift by 3.
 119          */
 120         speed_distrust = (ao_speed - AO_MS_TO_SPEED(AO_MAX_BARO_SPEED)) >> (4 - 1);
 121         if (speed_distrust > AO_MAX_SPEED_DISTRUST)
 122                 speed_distrust = AO_MAX_SPEED_DISTRUST;
 123         if (speed_distrust > height_distrust)
 124                 height_distrust = speed_distrust;
 125 #endif
 126         if (height_distrust > 0) {
 127 #ifdef AO_FLIGHT_TEST
 128                 int     old_ao_error_h = ao_error_h;
 129 #endif
 130                 if (height_distrust > 0x100)
 131                         height_distrust = 0x100;
 132                 ao_error_h = (ao_v_t) (((ao_k_t) ao_error_h * (0x100 - height_distrust)) >> 8);
 133 #ifdef AO_FLIGHT_TEST
 134                 if (ao_flight_debug) {
 135                         printf("over height %g over speed %g distrust: %g height: error %d -> %d\n",
 136                                (double) (ao_sample_alt - AO_MAX_BARO_HEIGHT),
 137                                (ao_speed - AO_MS_TO_SPEED(AO_MAX_BARO_SPEED)) / 16.0,
 138                                height_distrust / 256.0,
 139                                old_ao_error_h, ao_error_h);
 140                 }
 141 #endif
 142         }
 143 }
 144 #endif
 145
 146 #if HAS_BARO
 147 static void
 148 ao_kalman_correct_baro(void)
 149 {
 150         ao_kalman_err_height();
 151 #ifdef AO_FLIGHT_TEST
 152         if ((int16_t) (ao_sample_tick - ao_sample_prev_tick) > 50) {
 153                 ao_k_height += (ao_k_t) AO_BARO_K0_1 * ao_error_h;
 154                 ao_k_speed  += (ao_k_t) AO_BARO_K1_1 * ao_error_h;
 155                 ao_k_accel  += (ao_k_t) AO_BARO_K2_1 * ao_error_h;
 156                 return;
 157         }
 158         if ((int16_t) (ao_sample_tick - ao_sample_prev_tick) > 5) {
 159                 ao_k_height += (ao_k_t) AO_BARO_K0_10 * ao_error_h;
 160                 ao_k_speed  += (ao_k_t) AO_BARO_K1_10 * ao_error_h;
 161                 ao_k_accel  += (ao_k_t) AO_BARO_K2_10 * ao_error_h;
 162                 return;
 163         }
 164 #endif
 165         ao_k_height += (ao_k_t) AO_BARO_K0_100 * ao_error_h;
 166         ao_k_speed  += (ao_k_t) AO_BARO_K1_100 * ao_error_h;
 167         ao_k_accel  += (ao_k_t) AO_BARO_K2_100 * ao_error_h;
 168 }
 169 #endif
 170
 171 #if HAS_ACCEL
 172
 173 static void
 174 ao_kalman_err_accel(void)
 175 {
 176         ao_k_t  accel;
 177
 178         accel = (ao_config.accel_plus_g - ao_sample_accel) * ao_accel_scale;
 179
 180         /* Can't use ao_accel here as it is the pre-prediction value still */
 181         ao_error_a = (accel - ao_k_accel) >> 16;
 182 }
 183
 184 #if !defined(FORCE_ACCEL) && HAS_BARO
 185 static void
 186 ao_kalman_correct_both(void)
 187 {
 188         ao_kalman_err_height();
 189         ao_kalman_err_accel();
 190
 191 #ifdef AO_FLIGHT_TEST
 192         if ((int16_t) (ao_sample_tick - ao_sample_prev_tick) > 50) {
 193                 if (ao_flight_debug) {
 194                         printf ("correct speed %g + (%g * %g) + (%g * %g) = %g\n",
 195                                 ao_k_speed / (65536.0 * 16.0),
 196                                 (double) ao_error_h, AO_BOTH_K10_1 / 65536.0,
 197                                 (double) ao_error_a, AO_BOTH_K11_1 / 65536.0,
 198                                 (ao_k_speed +
 199                                  (ao_k_t) AO_BOTH_K10_1 * ao_error_h +
 200                                  (ao_k_t) AO_BOTH_K11_1 * ao_error_a) / (65536.0 * 16.0));
 201                 }
 202                 ao_k_height +=
 203                         (ao_k_t) AO_BOTH_K00_1 * ao_error_h +
 204                         (ao_k_t) AO_BOTH_K01_1 * ao_error_a;
 205                 ao_k_speed +=
 206                         (ao_k_t) AO_BOTH_K10_1 * ao_error_h +
 207                         (ao_k_t) AO_BOTH_K11_1 * ao_error_a;
 208                 ao_k_accel +=
 209                         (ao_k_t) AO_BOTH_K20_1 * ao_error_h +
 210                         (ao_k_t) AO_BOTH_K21_1 * ao_error_a;
 211                 return;
 212         }
 213         if ((int16_t) (ao_sample_tick - ao_sample_prev_tick) > 5) {
 214                 if (ao_flight_debug) {
 215                         printf ("correct speed %g + (%g * %g) + (%g * %g) = %g\n",
 216                                 ao_k_speed / (65536.0 * 16.0),
 217                                 (double) ao_error_h, AO_BOTH_K10_10 / 65536.0,
 218                                 (double) ao_error_a, AO_BOTH_K11_10 / 65536.0,
 219                                 (ao_k_speed +
 220                                  (ao_k_t) AO_BOTH_K10_10 * ao_error_h +
 221                                  (ao_k_t) AO_BOTH_K11_10 * ao_error_a) / (65536.0 * 16.0));
 222                 }
 223                 ao_k_height +=
 224                         (ao_k_t) AO_BOTH_K00_10 * ao_error_h +
 225                         (ao_k_t) AO_BOTH_K01_10 * ao_error_a;
 226                 ao_k_speed +=
 227                         (ao_k_t) AO_BOTH_K10_10 * ao_error_h +
 228                         (ao_k_t) AO_BOTH_K11_10 * ao_error_a;
 229                 ao_k_accel +=
 230                         (ao_k_t) AO_BOTH_K20_10 * ao_error_h +
 231                         (ao_k_t) AO_BOTH_K21_10 * ao_error_a;
 232                 return;
 233         }
 234         if (ao_flight_debug) {
 235                 printf ("correct speed %g + (%g * %g) + (%g * %g) = %g\n",
 236                         ao_k_speed / (65536.0 * 16.0),
 237                         (double) ao_error_h, AO_BOTH_K10_100 / 65536.0,
 238                         (double) ao_error_a, AO_BOTH_K11_100 / 65536.0,
 239                         (ao_k_speed +
 240                          (ao_k_t) AO_BOTH_K10_100 * ao_error_h +
 241                          (ao_k_t) AO_BOTH_K11_100 * ao_error_a) / (65536.0 * 16.0));
 242         }
 243 #endif
 244         ao_k_height +=
 245                 (ao_k_t) AO_BOTH_K00_100 * ao_error_h +
 246                 (ao_k_t) AO_BOTH_K01_100 * ao_error_a;
 247         ao_k_speed +=
 248                 (ao_k_t) AO_BOTH_K10_100 * ao_error_h +
 249                 (ao_k_t) AO_BOTH_K11_100 * ao_error_a;
 250         ao_k_accel +=
 251                 (ao_k_t) AO_BOTH_K20_100 * ao_error_h +
 252                 (ao_k_t) AO_BOTH_K21_100 * ao_error_a;
 253 }
 254
 255 #else
 256
 257 static void
 258 ao_kalman_correct_accel(void)
 259 {
 260         ao_kalman_err_accel();
 261
 262 #ifdef AO_FLIGHT_TEST
 263         if ((int16_t) (ao_sample_tick - ao_sample_prev_tick) > 5) {
 264                 ao_k_height +=(ao_k_t) AO_ACCEL_K0_10 * ao_error_a;
 265                 ao_k_speed  += (ao_k_t) AO_ACCEL_K1_10 * ao_error_a;
 266                 ao_k_accel  += (ao_k_t) AO_ACCEL_K2_10 * ao_error_a;
 267                 return;
 268         }
 269 #endif
 270         ao_k_height += (ao_k_t) AO_ACCEL_K0_100 * ao_error_a;
 271         ao_k_speed  += (ao_k_t) AO_ACCEL_K1_100 * ao_error_a;
 272         ao_k_accel  += (ao_k_t) AO_ACCEL_K2_100 * ao_error_a;
 273 }
 274
 275 #endif /* else FORCE_ACCEL */
 276 #endif /* HAS_ACCEL */
 277
 278 #if !HAS_BARO
 279 static ao_k_t   ao_k_height_prev;
 280 static ao_k_t   ao_k_speed_prev;
 281
 282 /*
 283  * While in pad mode without a barometric sensor, remove accumulated
 284  * speed and height values to reduce the effect of systematic sensor
 285  * error
 286  */
 287 void
 288 ao_kalman_reset_accumulate(void)
 289 {
 290         ao_k_height -= ao_k_height_prev;
 291         ao_k_speed -= ao_k_speed_prev;
 292         ao_k_height_prev = ao_k_height;
 293         ao_k_speed_prev = ao_k_speed;
 294 }
 295 #endif
 296
 297 void
 298 ao_kalman(void)
 299 {
 300         ao_kalman_predict();
 301 #if HAS_BARO
 302 #if HAS_ACCEL
 303         if (ao_flight_state <= ao_flight_coast) {
 304 #ifdef FORCE_ACCEL
 305                 ao_kalman_correct_accel();
 306 #else
 307                 ao_kalman_correct_both();
 308 #endif
 309         } else
 310 #endif
 311                 ao_kalman_correct_baro();
 312 #else
 313 #if HAS_ACCEL
 314         ao_kalman_correct_accel();
 315 #endif
 316 #endif
 317         ao_height = from_fix(ao_k_height);
 318         ao_speed = from_fix(ao_k_speed);
 319         ao_accel = from_fix(ao_k_accel);
 320         if (ao_height > ao_max_height)
 321                 ao_max_height = ao_height;
 322 #if HAS_BARO
 323         ao_avg_height_scaled = ao_avg_height_scaled - ao_avg_height + ao_sample_height;
 324 #else
 325         ao_avg_height_scaled = ao_avg_height_scaled - ao_avg_height + ao_height;
 326 #endif
 327 #ifdef AO_FLIGHT_TEST
 328         if ((int16_t) (ao_sample_tick - ao_sample_prev_tick) > 50)
 329                 ao_avg_height = (ao_avg_height_scaled + 1) >> 1;
 330         else if ((int16_t) (ao_sample_tick - ao_sample_prev_tick) > 5)
 331                 ao_avg_height = (ao_avg_height_scaled + 7) >> 4;
 332         else
 333 #endif
 334                 ao_avg_height = (ao_avg_height_scaled + 63) >> 7;
 335 }