-
Notifications
You must be signed in to change notification settings - Fork 26
/
kalman_filter.cpp
71 lines (62 loc) · 1.49 KB
/
kalman_filter.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
#include "kalman_filter.h"
using Eigen::MatrixXd;
using Eigen::VectorXd;
KalmanFilter::KalmanFilter() {}
KalmanFilter::~KalmanFilter() {}
void KalmanFilter::Init(VectorXd &x_in, MatrixXd &P_in, MatrixXd &F_in,
MatrixXd &H_in, MatrixXd &R_in, MatrixXd &Q_in) {
x_ = x_in;
P_ = P_in;
F_ = F_in;
H_ = H_in;
R_ = R_in;
Q_ = Q_in;
}
void KalmanFilter::Predict() {
x_ = F_ * x_ ;
MatrixXd Ft = F_.transpose();
P_ = F_ * P_ * Ft + Q_;
}
void KalmanFilter::Update(const VectorXd &z) {
/**
TODO:
* update the state by using Kalman Filter equations
*/
VectorXd y = z - H_ * x_;
UpdateWithY(y);
}
void KalmanFilter::UpdateEKF(const VectorXd &z) {
/**
TODO:
* update the state by using Extended Kalman Filter equations
*/
double px = x_(0);
double py = x_(1);
double vx = x_(2);
double vy = x_(3);
double rho = sqrt(px*px + py*py);
double theta = atan2(py, px);
double rho_dot = (px*vx + py*vy) / rho;
VectorXd h = VectorXd(3);
h << rho, theta, rho_dot;
VectorXd y = z - h;
while ( y(1) > M_PI || y(1) < -M_PI ) {
if ( y(1) > M_PI ) {
y(1) -= M_PI;
} else {
y(1) += M_PI;
}
}
UpdateWithY(y);
}
void KalmanFilter::UpdateWithY(const VectorXd &y){
MatrixXd Ht = H_.transpose();
MatrixXd S = H_ * P_ * Ht + R_;
MatrixXd Si = S.inverse();
MatrixXd K = P_ * Ht * Si;
// New state
x_ = x_ + (K * y);
int x_size = x_.size();
MatrixXd I = MatrixXd::Identity(x_size, x_size);
P_ = (I - K * H_) * P_;
}