diff options
| author | Dawsyn Schraiber <[email protected]> | 2024-05-09 02:05:35 -0400 |
|---|---|---|
| committer | GitHub <[email protected]> | 2024-05-09 02:05:35 -0400 |
| commit | 93acde052369568beaefb0d99629d8797f5c191f (patch) | |
| tree | a3fb96ddad2d289aa7f8bf410c60cf6289bca7a1 /src/AltEst/filters.cpp | |
| parent | 5f68c7a1b5c8dec82d4a2e1e12443a41b5196b1d (diff) | |
| download | active-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.cpp | 202 |
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); +} |
