diff options
Diffstat (limited to 'src/AltEst/filters.cpp')
| -rw-r--r-- | src/AltEst/filters.cpp | 202 |
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); -} |
