summaryrefslogtreecommitdiff
path: root/src/AltEst/filters.cpp
diff options
context:
space:
mode:
authorDawsyn Schraiber <[email protected]>2024-05-09 02:05:35 -0400
committerGitHub <[email protected]>2024-05-09 02:05:35 -0400
commit93acde052369568beaefb0d99629d8797f5c191f (patch)
treea3fb96ddad2d289aa7f8bf410c60cf6289bca7a1 /src/AltEst/filters.cpp
parent5f68c7a1b5c8dec82d4a2e1e12443a41b5196b1d (diff)
downloadactive-drag-system-93acde052369568beaefb0d99629d8797f5c191f.tar.gz
active-drag-system-93acde052369568beaefb0d99629d8797f5c191f.tar.bz2
active-drag-system-93acde052369568beaefb0d99629d8797f5c191f.zip
Raspberry Pi Pico (#12)
* Adding a 90% completed, compilable but untested ADS * Made basic changes to actuator & sensor. Also added motor class * Removed unnecessary .cpp files * Updated sensor & actuator classes, finished ads, added variable time step to kalman filter, set up all tests for future assertions * Relocated 'main' to 'active-drag-system.cpp'. Added more info to README * Removed main.cpp * Added more details to README * Changed some function parameters from pass-by-pointer to pass-by-reference. Also removed the std namespace * Started writing the test cases * Updated the .gitignore file * Removed some files that should be gitignored * Up to date with Jazz's pull request * Test Launch Branch Created; PRU Servo Control with Test Program * Added I2C device class and register IDs for MPL [INCOMPLETE SENSOR IMPLEMENTATION] Needs actual data getting function implementation for both sensors and register IDs for BNO, will implement shortly. * Partial implementation of MPL sensor Added startup method, still needs fleshed out data getters and setters and finished I2C implementation. MOST LIKELY WILL HAVE COMPILATION ISSUES. * *Hypothetically* complete MPL implementation NEEDS HARDWARE TESTING * IMU Header and init() method implementation Needs like, all data handling still lol * Hypothetically functional (Definitely won't compile) * We ball? * Conversion to Raspberry Pi Pico Build System; Removed Beaglebone specific code; Simple blinking example in ADS source file; builds for Pico W * Rearranged build so dependent upon cmake file already existing in pico-sdk; current executable prints current altitude, velocity, and time taken to read and calculate said values; ~320 us to do so * Altimeter interrupt callback for Pad to Boost State; dummy templates for other callbacks with comments describing potential implementation details * Altimeter interrupts relatively finished; need to test with vacuum chamber to verify behavior * Established interrupt pins as pullup and active-low; adjusted callback functions to properly use function pointers; still need to verify interrupt system with vacuum chamber * Removed weird artifact in .gitignore, adjust CMakeLists to auto pull pico sdk, added Dockerfile * added Docker dev container file * modified CMakeLists to auto pull sdk if not already downloaded, add build.sh script, fixed Dockerfile * added bno055 support * changed bno055 lin accel struct to use float instead of double * added bno055 support not tested, but compiles, fixed CMakLists to before I messed with it * added absolute quaternion output from bno055 * Added Euler and aboslute linear accelration * Flash implementation for data logging; each log entry is 32 bytes long * added base pwm functions and started on apogee detection * State machine verified functional with logging capabilities; currently on same core * Ooops missed double define, renamed LOOP_HZ to LOOP_PERIOD; State machine functional after merge still * Simple test program to see servo PWM range; logging with semaphores for safe multithreading * Kalman filters generously provided from various sources for temporary replacement; minimum deployment 30 percent; state machine functionality restored; multithreading logging verified; altimeter broke and replaced * Stop logging on END state; provide deployment function with AGL instead of ASL altitude * Various minimal changes; Flash size from 1MB to 8MB; M1939 to M2500T burn time; pin assignments for new PCB; External Status LED to Internal Status LED --------- Co-authored-by: Jazz Jackson <[email protected]> Co-authored-by: Cian Capacci <[email protected]> Co-authored-by: Gregory Wainer <[email protected]>
Diffstat (limited to 'src/AltEst/filters.cpp')
-rw-r--r--src/AltEst/filters.cpp202
1 files changed, 202 insertions, 0 deletions
diff --git a/src/AltEst/filters.cpp b/src/AltEst/filters.cpp
new file mode 100644
index 0000000..7902065
--- /dev/null
+++ b/src/AltEst/filters.cpp
@@ -0,0 +1,202 @@
+/*
+ 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);
+}