From 93acde052369568beaefb0d99629d8797f5c191f Mon Sep 17 00:00:00 2001 From: Dawsyn Schraiber <32221234+dawsynth@users.noreply.github.com> Date: Thu, 9 May 2024 02:05:35 -0400 Subject: 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 Co-authored-by: Cian Capacci <46614295+BeeGuyDude@users.noreply.github.com> Co-authored-by: Gregory Wainer --- src/AltEst/filters.cpp | 202 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 202 insertions(+) create mode 100644 src/AltEst/filters.cpp (limited to 'src/AltEst/filters.cpp') 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 +#include // 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); +} -- cgit v1.2.3