|
| 1 | +// KalmanFilter.cpp : Uses all available information to make optimal estimation. |
| 2 | +// C++ implementation of Python code. |
| 3 | +// From Udacity CS 373 Artificial Intelliegence for Robots. |
| 4 | +// Taught by Sebastian Thrun. |
| 5 | +// C++ version by Tyler Folsom; March 6, 2012. |
| 6 | +// |
| 7 | + |
| 8 | + |
| 9 | +#include "Matrix.h" |
| 10 | +#include "Kalman.h" |
| 11 | +/* |
| 12 | +Kalman filter |
| 13 | +Matrix x state: longs for position_mm_east, position_mm_north, |
| 14 | + velocity_mm_s_east, velocity_mm_s_north |
| 15 | +*/ |
| 16 | +#define NSTATES (7) |
| 17 | +#define NSEEN (7) |
| 18 | +void Filter(REAL* State, // e.g x, y position, x, y velocity |
| 19 | + REAL* uncertainty, // P uncertainty covariance |
| 20 | + REAL* measurements, // Z: measured x and y position in mm |
| 21 | + REAL deltaT_s, // time from previous state to measurement in seconds |
| 22 | + REAL* variance) // variance of the measurements |
| 23 | +{ |
| 24 | + // State is set for the Sensor Hub; It needs to be modified for Vision. TCF 6/29/20 |
| 25 | + REAL state_transition[NSTATES * NSTATES] = |
| 26 | + { 1., 0 , deltaT_s, 0, 0.5 * deltaT_s * deltaT_s, 0, 0, |
| 27 | + 0, 1., 0, deltaT_s, 0 , 0.5 * deltaT_s * deltaT_s, 0, |
| 28 | + 0, 0 , 1. , 0, deltaT_s, 0, 0, |
| 29 | + 0, 0 , 0, 1, 0, deltaT_s, 0, |
| 30 | + 0, 0 , 0, 0, 1, 0, 0, |
| 31 | + 0, 0 , 0, 0, 0, 1, 0, |
| 32 | + 0, 0 , 0, 0, 0, 0, 1 }; |
| 33 | + //REAL state_transition[] = {1., 0, deltaT_s, 0, |
| 34 | + // 0, 1., 0, deltaT_s, |
| 35 | + // 0, 0, 1, 0, |
| 36 | + // 0, 0, 0, 1}; |
| 37 | + // position' = 1*position + 1*velocity |
| 38 | + // velocity' = 1*velocity |
| 39 | + REAL observable[NSEEN*NSTATES]; // = {1., 0, 0, 0, |
| 40 | + // 0, 1., 0, 0}; |
| 41 | + // only position is measurable; not velocity |
| 42 | + |
| 43 | + /*----------------------------------------------------------------------------- |
| 44 | + // Compiler is confused; it thinks there is just one argument of type matrix TCF 6/11/20 |
| 45 | + // Polymorphism seems to be a C++ featurs not well supported by the Arduino dialect. |
| 46 | + // Aduino code needs to be rewriiten with C style mechanisms. |
| 47 | + ------------------------------------------------------------------------------*/ |
| 48 | + matrix x = matrix( NSTATES, 1, State); // initial state (location and velocity) |
| 49 | + matrix P = matrix(NSTATES,NSTATES, uncertainty); // 4x4 initial uncertainty |
| 50 | + matrix meas = matrix(NSEEN,1,measurements); |
| 51 | + |
| 52 | + matrix u = matrix( NSTATES, 1); // 4x1 external motion; set to 0. |
| 53 | + // u comes from what the vehicle is trying to do, e.g. accelerate. |
| 54 | + matrix F = matrix(NSTATES,NSTATES, state_transition); // next state function |
| 55 | + for (int i = 0; i < NSEEN*NSTATES; i++) |
| 56 | + { |
| 57 | + observable[i] = i%(NSTATES+1) ? 0 : 1; |
| 58 | + } |
| 59 | + matrix H = matrix(NSEEN,NSTATES, observable); // measurement function |
| 60 | + H.show(); |
| 61 | + matrix R = matrix(NSEEN,NSEEN, variance); // measurement variance |
| 62 | + matrix I = matrix(NSTATES); // 4x4 identity matrix |
| 63 | + |
| 64 | + // prediction |
| 65 | + matrix Fx = F * x; |
| 66 | + matrix xNew = Fx + u; |
| 67 | + x = xNew; // x = F*x + u |
| 68 | + matrix Ftrans = F.transpose(); |
| 69 | + matrix PFt = P * Ftrans; |
| 70 | + matrix P2 = F * PFt; |
| 71 | + P = P2; // P = F * P * F.transpose(); |
| 72 | + Show("x= "); |
| 73 | + x.show(); |
| 74 | + Show("dt = "); |
| 75 | + Show(deltaT_s); |
| 76 | + Show("P= "); |
| 77 | + P.show(); |
| 78 | + |
| 79 | + // measurement update |
| 80 | + // Arduino compiler gets confused on combining multiple matrix operations. |
| 81 | + |
| 82 | + matrix Hx = H * x; |
| 83 | + matrix y = meas - Hx; // y = Z- H * x |
| 84 | + matrix transp = H.transpose(); |
| 85 | + matrix PHt = P * transp; |
| 86 | + matrix HPH = H * PHt; |
| 87 | + matrix S = HPH + R; // S = H * P * H.transpose() + R; |
| 88 | + matrix Sinv = S.inverse(); |
| 89 | + matrix K = PHt * Sinv;// K = P * H.transpose()*S.inverse(); |
| 90 | + matrix Ky = K * y; |
| 91 | + matrix x_Ky = x + Ky; |
| 92 | + x = x_Ky; |
| 93 | + matrix KH = K*H; |
| 94 | + matrix IKH = I - KH; |
| 95 | + matrix Pnew = IKH * P; // (I-(K*H)) * P; |
| 96 | + P = Pnew; |
| 97 | + Show("R= "); |
| 98 | + R.show(); |
| 99 | + Show("K= "); |
| 100 | + K.show(); |
| 101 | + Show("x= "); |
| 102 | + x.show(); |
| 103 | + Show("P= "); |
| 104 | + P.show(); |
| 105 | + // Put x into State and P into Uncertainty. |
| 106 | + x.values(State); |
| 107 | + P.values(uncertainty); |
| 108 | + return; |
| 109 | +} |
0 commit comments