summaryrefslogtreecommitdiff
path: root/src/AltEst/filters.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'src/AltEst/filters.cpp')
-rw-r--r--src/AltEst/filters.cpp202
1 files changed, 0 insertions, 202 deletions
diff --git a/src/AltEst/filters.cpp b/src/AltEst/filters.cpp
deleted file mode 100644
index 7902065..0000000
--- a/src/AltEst/filters.cpp
+++ /dev/null
@@ -1,202 +0,0 @@
-/*
- filters.cpp: Filter class implementations
- */
-
-//#include <cmath>
-#include <stdlib.h> // XXX eventually use fabs() instead of abs() ?
-
-#include "filters.h"
-
-void KalmanFilter::getPredictionCovariance(float covariance[3][3], float previousState[3], float deltat)
-{
- // required matrices for the operations
- float sigma[3][3];
- float identity[3][3];
- identityMatrix3x3(identity);
- float skewMatrix[3][3];
- skew(skewMatrix, previousState);
- float tmp[3][3];
- // Compute the prediction covariance matrix
- scaleMatrix3x3(sigma, pow(sigmaGyro, 2), identity);
- matrixProduct3x3(tmp, skewMatrix, sigma);
- matrixProduct3x3(covariance, tmp, skewMatrix);
- scaleMatrix3x3(covariance, -pow(deltat, 2), covariance);
-}
-
-void KalmanFilter::getMeasurementCovariance(float covariance[3][3])
-{
- // required matrices for the operations
- float sigma[3][3];
- float identity[3][3];
- identityMatrix3x3(identity);
- float norm;
- // Compute measurement covariance
- scaleMatrix3x3(sigma, pow(sigmaAccel, 2), identity);
- vectorLength(& norm, previousAccelSensor);
- scaleAndAccumulateMatrix3x3(sigma, (1.0/3.0)*pow(ca, 2)*norm, identity);
- copyMatrix3x3(covariance, sigma);
-}
-
-void KalmanFilter::predictState(float predictedState[3], float gyro[3], float deltat)
-{
- // helper matrices
- float identity[3][3];
- identityMatrix3x3(identity);
- float skewFromGyro[3][3];
- skew(skewFromGyro, gyro);
- // Predict state
- scaleAndAccumulateMatrix3x3(identity, -deltat, skewFromGyro);
- matrixDotVector3x3(predictedState, identity, currentState);
- normalizeVector(predictedState);
-}
-
-void KalmanFilter::predictErrorCovariance(float covariance[3][3], float gyro[3], float deltat)
-{
- // required matrices
- float Q[3][3];
- float identity[3][3];
- identityMatrix3x3(identity);
- float skewFromGyro[3][3];
- skew(skewFromGyro, gyro);
- float tmp[3][3];
- float tmpTransposed[3][3];
- float tmp2[3][3];
- // predict error covariance
- getPredictionCovariance(Q, currentState, deltat);
- scaleAndAccumulateMatrix3x3(identity, -deltat, skewFromGyro);
- copyMatrix3x3(tmp, identity);
- transposeMatrix3x3(tmpTransposed, tmp);
- matrixProduct3x3(tmp2, tmp, currErrorCovariance);
- matrixProduct3x3(covariance, tmp2, tmpTransposed);
- scaleAndAccumulateMatrix3x3(covariance, 1.0, Q);
-}
-
-void KalmanFilter::updateGain(float gain[3][3], float errorCovariance[3][3])
-{
- // required matrices
- float R[3][3];
- float HTransposed[3][3];
- transposeMatrix3x3(HTransposed, H);
- float tmp[3][3];
- float tmp2[3][3];
- float tmp2Inverse[3][3];
- // update kalman gain
- // P.dot(H.T).dot(inv(H.dot(P).dot(H.T) + R))
- getMeasurementCovariance(R);
- matrixProduct3x3(tmp, errorCovariance, HTransposed);
- matrixProduct3x3(tmp2, H, tmp);
- scaleAndAccumulateMatrix3x3(tmp2, 1.0, R);
- invert3x3(tmp2Inverse, tmp2);
- matrixProduct3x3(gain, tmp, tmp2Inverse);
-}
-
-void KalmanFilter::updateState(float updatedState[3], float predictedState[3], float gain[3][3], float accel[3])
-{
- // required matrices
- float tmp[3];
- float tmp2[3];
- float measurement[3];
- scaleVector(tmp, ca, previousAccelSensor);
- subtractVectors(measurement, accel, tmp);
- // update state with measurement
- // predicted_state + K.dot(measurement - H.dot(predicted_state))
- matrixDotVector3x3(tmp, H, predictedState);
- subtractVectors(tmp, measurement, tmp);
- matrixDotVector3x3(tmp2, gain, tmp);
- sumVectors(updatedState, predictedState, tmp2);
- normalizeVector(updatedState);
-}
-
-void KalmanFilter::updateErrorCovariance(float covariance[3][3], float errorCovariance[3][3], float gain[3][3])
-{
- // required matrices
- float identity[3][3];
- identityMatrix3x3(identity);
- float tmp[3][3];
- float tmp2[3][3];
- // update error covariance with measurement
- matrixProduct3x3(tmp, gain, H);
- matrixProduct3x3(tmp2, tmp, errorCovariance);
- scaleAndAccumulateMatrix3x3(identity, -1.0, tmp2);
- copyMatrix3x3(covariance, tmp2);
-}
-
-
-KalmanFilter::KalmanFilter(float ca, float sigmaGyro, float sigmaAccel)
-{
- this->ca = ca;
- this->sigmaGyro = sigmaGyro;
- this->sigmaAccel = sigmaAccel;
-}
-
-float KalmanFilter::estimate(float gyro[3], float accel[3], float deltat)
-{
- float predictedState[3];
- float updatedState[3];
- float errorCovariance[3][3];
- float updatedErrorCovariance[3][3];
- float gain[3][3];
- float accelSensor[3];
- float tmp[3];
- float accelEarth;
- scaleVector(accel, 9.81, accel); // Scale accel readings since they are measured in gs
- // perform estimation
- // predictions
- predictState(predictedState, gyro, deltat);
- predictErrorCovariance(errorCovariance, gyro, deltat);
- // updates
- updateGain(gain, errorCovariance);
- updateState(updatedState, predictedState, gain, accel);
- updateErrorCovariance(updatedErrorCovariance, errorCovariance, gain);
- // Store required values for next iteration
- copyVector(currentState, updatedState);
- copyMatrix3x3(currErrorCovariance, updatedErrorCovariance);
- // return vertical acceleration estimate
- scaleVector(tmp, 9.81, updatedState);
- subtractVectors(accelSensor, accel, tmp);
- copyVector(previousAccelSensor, accelSensor);
- dotProductVectors(& accelEarth, accelSensor, updatedState);
- return accelEarth;
-}
-
-
-float ComplementaryFilter::ApplyZUPT(float accel, float vel)
-{
- // first update ZUPT array with latest estimation
- ZUPT[ZUPTIdx] = accel;
- // and move index to next slot
- uint8_t nextIndex = (ZUPTIdx + 1) % ZUPT_SIZE;
- ZUPTIdx = nextIndex;
- // Apply Zero-velocity update
- for (uint8_t k = 0; k < ZUPT_SIZE; ++k) {
- if (abs(ZUPT[k]) > accelThreshold) return vel;
- }
- return 0.0;
-}
-
-
-ComplementaryFilter::ComplementaryFilter(float sigmaAccel, float sigmaBaro, float accelThreshold)
-{
- // Compute the filter gain
- gain[0] = sqrt(2 * sigmaAccel / sigmaBaro);
- gain[1] = sigmaAccel / sigmaBaro;
- // If acceleration is below the threshold the ZUPT counter
- // will be increased
- this->accelThreshold = accelThreshold;
- // initialize zero-velocity update
- ZUPTIdx = 0;
- for (uint8_t k = 0; k < ZUPT_SIZE; ++k) {
- ZUPT[k] = 0;
- }
-}
-
-void ComplementaryFilter::estimate(float * velocity, float * altitude, float baroAltitude,
- float pastAltitude, float pastVelocity, float accel, float deltat)
-{
- // Apply complementary filter
- *altitude = pastAltitude + deltat*(pastVelocity + (gain[0] + gain[1]*deltat/2)*(baroAltitude-pastAltitude))+
- accel*pow(deltat, 2)/2;
- *velocity = pastVelocity + deltat*(gain[1]*(baroAltitude-pastAltitude) + accel);
- // Compute zero-velocity update
- *velocity = ApplyZUPT(accel, *velocity);
-}