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]>

6 files changed, 768 insertions, 0 deletions

diff --git a/src/AltEst/algebra.cpp b/src/AltEst/algebra.cpp
new file mode 100644
index 0000000..653c3b9
--- /dev/null
+++ b/src/AltEst/algebra.cpp
@@ -0,0 +1,292 @@
+/*
+  algebra.cpp: This file contains a number of utilities useful for handling
+  3D vectors
+
+  This work is an adaptation from vvector.h, written by Linas Vepstras. The
+  original code can be found at:
+
+  https://github.com/markkilgard/glut/blob/master/lib/gle/vvector.h
+
+  HISTORY:
+  Written by Linas Vepstas, August 1991
+  Added 2D code, March 1993
+  Added Outer products, C++ proofed, Linas Vepstas October 1993
+  Adapted for altitude estimation tasks by Juan Gallostra June 2018
+*/
+
+//#include <cmath>
+#include <stdio.h>
+
+#include "algebra.h"
+
+// Copy 3D vector
+void copyVector(float b[3],float a[3])
+{
+   b[0] = a[0];
+   b[1] = a[1];
+   b[2] = a[2];
+}
+
+
+// Vector difference
+void subtractVectors(float v21[3], float v2[3], float v1[3])
+{
+   v21[0] = v2[0] - v1[0];
+   v21[1] = v2[1] - v1[1];
+   v21[2] = v2[2] - v1[2];
+}
+
+// Vector sum
+void sumVectors(float v21[3], float v2[3], float v1[3])
+{
+   v21[0] = v2[0] + v1[0];
+   v21[1] = v2[1] + v1[1];
+   v21[2] = v2[2] + v1[2];
+}
+
+// scalar times vector
+void scaleVector(float c[3],float a, float b[3])
+{
+   (c)[0] = a*b[0];
+   (c)[1] = a*b[1];
+   (c)[2] = a*b[2];
+}
+
+// accumulate scaled vector
+void accumulateScaledVector(float c[3], float a, float b[3])
+{
+   (c)[0] += a*b[0];
+   (c)[1] += a*b[1];
+   (c)[2] += a*b[2];
+}
+
+// Vector dot product
+void dotProductVectors(float * c, float a[3], float b[3])
+{
+   *c = a[0]*b[0] + a[1]*b[1] + a[2]*b[2];
+}
+
+// Vector length
+void vectorLength(float * len, float a[3])
+{
+   float tmp;
+   tmp = a[0]*a[0] + a[1]*a[1]+a[2]*a[2];
+   *len = sqrt(tmp);
+}
+
+// Normalize vector
+void normalizeVector(float a[3])
+{
+   float len;
+   vectorLength(& len,a);
+   if (len != 0.0) {
+      len = 1.0 / len;
+      a[0] *= len;
+      a[1] *= len;
+      a[2] *= len;
+   }
+}
+
+// 3D Vector cross product yeilding vector
+void crossProductVectors(float c[3], float a[3], float b[3])
+{
+   c[0] = a[1] * b[2] - a[2] * b[1];
+   c[1] = a[2] * b[0] - a[0] * b[2];
+   c[2] = a[0] * b[1] - a[1] * b[0];
+}
+
+// initialize matrix
+void identityMatrix3x3(float m[3][3])
+{
+   m[0][0] = 1.0;
+   m[0][1] = 0.0;
+   m[0][2] = 0.0;
+
+   m[1][0] = 0.0;
+   m[1][1] = 1.0;
+   m[1][2] = 0.0;
+
+   m[2][0] = 0.0;
+   m[2][1] = 0.0;
+   m[2][2] = 1.0;
+}
+
+// matrix copy
+void copyMatrix3x3(float b[3][3], float a[3][3])
+{
+   b[0][0] = a[0][0];
+   b[0][1] = a[0][1];
+   b[0][2] = a[0][2];
+
+   b[1][0] = a[1][0];
+   b[1][1] = a[1][1];
+   b[1][2] = a[1][2];
+
+   b[2][0] = a[2][0];
+   b[2][1] = a[2][1];
+   b[2][2] = a[2][2];
+}
+
+// matrix transpose
+void transposeMatrix3x3(float b[3][3], float a[3][3])
+{
+   b[0][0] = a[0][0];
+   b[0][1] = a[1][0];
+   b[0][2] = a[2][0];
+
+   b[1][0] = a[0][1];
+   b[1][1] = a[1][1];
+   b[1][2] = a[2][1];
+
+   b[2][0] = a[0][2];
+   b[2][1] = a[1][2];
+   b[2][2] = a[2][2];
+}
+
+// multiply matrix by scalar
+void scaleMatrix3x3(float b[3][3], float s, float a[3][3])
+{
+   b[0][0] = (s) * a[0][0];
+   b[0][1] = (s) * a[0][1];
+   b[0][2] = (s) * a[0][2];
+
+   b[1][0] = (s) * a[1][0];
+   b[1][1] = (s) * a[1][1];
+   b[1][2] = (s) * a[1][2];
+
+   b[2][0] = (s) * a[2][0];
+   b[2][1] = (s) * a[2][1];
+   b[2][2] = (s) * a[2][2];
+}
+
+// multiply matrix by scalar and add result to another matrix
+void scaleAndAccumulateMatrix3x3(float b[3][3], float s, float a[3][3])
+{
+   b[0][0] += s * a[0][0];
+   b[0][1] += s * a[0][1];
+   b[0][2] += s * a[0][2];
+
+   b[1][0] += s * a[1][0];
+   b[1][1] += s * a[1][1];
+   b[1][2] += s * a[1][2];
+
+   b[2][0] += s * a[2][0];
+   b[2][1] += s * a[2][1];
+   b[2][2] += s * a[2][2];
+}
+
+// matrix product
+// c[x][y] = a[x][0]*b[0][y]+a[x][1]*b[1][y]+a[x][2]*b[2][y]+a[x][3]*b[3][y]
+void matrixProduct3x3(float c[3][3], float a[3][3], float b[3][3])
+{
+   c[0][0] = a[0][0]*b[0][0]+a[0][1]*b[1][0]+a[0][2]*b[2][0];
+   c[0][1] = a[0][0]*b[0][1]+a[0][1]*b[1][1]+a[0][2]*b[2][1];
+   c[0][2] = a[0][0]*b[0][2]+a[0][1]*b[1][2]+a[0][2]*b[2][2];
+
+   c[1][0] = a[1][0]*b[0][0]+a[1][1]*b[1][0]+a[1][2]*b[2][0];
+   c[1][1] = a[1][0]*b[0][1]+a[1][1]*b[1][1]+a[1][2]*b[2][1];
+   c[1][2] = a[1][0]*b[0][2]+a[1][1]*b[1][2]+a[1][2]*b[2][2];
+
+   c[2][0] = a[2][0]*b[0][0]+a[2][1]*b[1][0]+a[2][2]*b[2][0];
+   c[2][1] = a[2][0]*b[0][1]+a[2][1]*b[1][1]+a[2][2]*b[2][1];
+   c[2][2] = a[2][0]*b[0][2]+a[2][1]*b[1][2]+a[2][2]*b[2][2];
+}
+
+// matrix times vector
+void matrixDotVector3x3(float p[3], float m[3][3], float v[3])
+{
+   p[0] = m[0][0]*v[0] + m[0][1]*v[1] + m[0][2]*v[2];
+   p[1] = m[1][0]*v[0] + m[1][1]*v[1] + m[1][2]*v[2];
+   p[2] = m[2][0]*v[0] + m[2][1]*v[1] + m[2][2]*v[2];
+}
+
+// determinant of matrix
+// Computes determinant of matrix m, returning d
+void determinant3x3(float * d, float m[3][3])
+{
+   *d = m[0][0] * (m[1][1]*m[2][2] - m[1][2] * m[2][1]);
+   *d -= m[0][1] * (m[1][0]*m[2][2] - m[1][2] * m[2][0]);
+   *d += m[0][2] * (m[1][0]*m[2][1] - m[1][1] * m[2][0]);
+}
+
+// adjoint of matrix
+// Computes adjoint of matrix m, returning a
+// (Note that adjoint is just the transpose of the cofactor matrix)
+void adjoint3x3(float a[3][3], float m[3][3])
+{
+   a[0][0] = m[1][1]*m[2][2] - m[1][2]*m[2][1];
+   a[1][0] = - (m[1][0]*m[2][2] - m[2][0]*m[1][2]);
+   a[2][0] = m[1][0]*m[2][1] - m[1][1]*m[2][0];
+   a[0][1] = - (m[0][1]*m[2][2] - m[0][2]*m[2][1]);
+   a[1][1] = m[0][0]*m[2][2] - m[0][2]*m[2][0];
+   a[2][1] = - (m[0][0]*m[2][1] - m[0][1]*m[2][0]);
+   a[0][2] = m[0][1]*m[1][2] - m[0][2]*m[1][1];
+   a[1][2] = - (m[0][0]*m[1][2] - m[0][2]*m[1][0]);
+   a[2][2] = m[0][0]*m[1][1] - m[0][1]*m[1][0];
+}
+
+// compute adjoint of matrix and scale
+// Computes adjoint of matrix m, scales it by s, returning a
+void scaleAdjoint3x3(float a[3][3], float s, float m[3][3])
+{
+   a[0][0] = (s) * (m[1][1] * m[2][2] - m[1][2] * m[2][1]);
+   a[1][0] = (s) * (m[1][2] * m[2][0] - m[1][0] * m[2][2]);
+   a[2][0] = (s) * (m[1][0] * m[2][1] - m[1][1] * m[2][0]);
+
+   a[0][1] = (s) * (m[0][2] * m[2][1] - m[0][1] * m[2][2]);
+   a[1][1] = (s) * (m[0][0] * m[2][2] - m[0][2] * m[2][0]);
+   a[2][1] = (s) * (m[0][1] * m[2][0] - m[0][0] * m[2][1]);
+
+   a[0][2] = (s) * (m[0][1] * m[1][2] - m[0][2] * m[1][1]);
+   a[1][2] = (s) * (m[0][2] * m[1][0] - m[0][0] * m[1][2]);
+   a[2][2] = (s) * (m[0][0] * m[1][1] - m[0][1] * m[1][0]);
+}
+
+// inverse of matrix
+// Compute inverse of matrix a, returning determinant m and
+// inverse b
+void invert3x3(float b[3][3], float a[3][3])
+{
+   float tmp;
+   determinant3x3(& tmp, a);
+   tmp = 1.0 / (tmp);
+   scaleAdjoint3x3(b, tmp, a);
+}
+
+// skew matrix from vector
+void skew(float a[3][3], float v[3])
+{
+  a[0][1] = -v[2];
+  a[0][2] = v[1];
+  a[1][2] = -v[0];
+  a[1][0] = v[2];
+  a[2][0] = -v[1];
+  a[2][1] = v[0];
+  // set diagonal to 0
+  a[0][0] = 0.0;
+  a[1][1] = 0.0;
+  a[2][2] = 0.0;
+}
+
+void printMatrix3X3(float mmm[3][3])
+{
+   int i,j;
+   printf ("matrix mmm is \n");
+   if (mmm == NULL) {
+      printf (" Null \n");
+   } else {
+      for (i=0; i<3; i++) {
+         for (j=0; j<3; j++) {
+            printf ("%f ", mmm[i][j]);
+         }
+         printf (" \n");
+      }
+   }
+}
+
+void vecPrint(float a[3])
+{
+   float len;
+   vectorLength(& len, a);
+   printf(" a is %f %f %f length of a is %f \n", a[0], a[1], a[2], len);
+}
diff --git a/src/AltEst/algebra.h b/src/AltEst/algebra.h
new file mode 100644
index 0000000..382103e
--- /dev/null
+++ b/src/AltEst/algebra.h
@@ -0,0 +1,96 @@
+/*
+  algebra.h: This file contains a number of utilities useful for handling
+  3D vectors
+
+  This work is an adaptation from vvector.h, written by Linas Vepstras. The
+  original code can be found at:
+
+  https://github.com/markkilgard/glut/blob/master/lib/gle/vvector.h
+
+  HISTORY:
+  Written by Linas Vepstas, August 1991
+  Added 2D code, March 1993
+  Added Outer products, C++ proofed, Linas Vepstas October 1993
+  Adapted for altitude estimation tasks by Juan Gallostra June 2018
+  Separated .h, .cpp by Simon D. Levy July 2018
+*/
+
+#pragma once
+
+//#include <cmath>
+#include <math.h>
+
+// Copy 3D vector
+void copyVector(float b[3],float a[3]);
+
+
+// Vector difference
+void subtractVectors(float v21[3], float v2[3], float v1[3]);
+
+// Vector sum
+void sumVectors(float v21[3], float v2[3], float v1[3]);
+
+// scalar times vector
+void scaleVector(float c[3],float a, float b[3]);
+
+// accumulate scaled vector
+void accumulateScaledVector(float c[3], float a, float b[3]);
+
+// Vector dot product
+void dotProductVectors(float * c, float a[3], float b[3]);
+
+// Vector length
+void vectorLength(float * len, float a[3]);
+
+// Normalize vector
+void normalizeVector(float a[3]);
+
+// 3D Vector cross product yeilding vector
+void crossProductVectors(float c[3], float a[3], float b[3]);
+
+// initialize matrix
+void identityMatrix3x3(float m[3][3]);
+
+// matrix copy
+void copyMatrix3x3(float b[3][3], float a[3][3]);
+
+// matrix transpose
+void transposeMatrix3x3(float b[3][3], float a[3][3]);
+
+// multiply matrix by scalar
+void scaleMatrix3x3(float b[3][3], float s, float a[3][3]);
+
+// multiply matrix by scalar and add result to another matrix
+void scaleAndAccumulateMatrix3x3(float b[3][3], float s, float a[3][3]);
+
+// matrix product
+// c[x][y] = a[x][0]*b[0][y]+a[x][1]*b[1][y]+a[x][2]*b[2][y]+a[x][3]*b[3][y]
+void matrixProduct3x3(float c[3][3], float a[3][3], float b[3][3]);
+
+// matrix times vector
+void matrixDotVector3x3(float p[3], float m[3][3], float v[3]);
+
+// determinant of matrix
+// Computes determinant of matrix m, returning d
+void determinant3x3(float * d, float m[3][3]);
+
+// adjoint of matrix
+// Computes adjoint of matrix m, returning a
+// (Note that adjoint is just the transpose of the cofactor matrix);
+void adjoint3x3(float a[3][3], float m[3][3]);
+
+// compute adjoint of matrix and scale
+// Computes adjoint of matrix m, scales it by s, returning a
+void scaleAdjoint3x3(float a[3][3], float s, float m[3][3]);
+
+// inverse of matrix
+// Compute inverse of matrix a, returning determinant m and
+// inverse b
+void invert3x3(float b[3][3], float a[3][3]);
+
+// skew matrix from vector
+void skew(float a[3][3], float v[3]);
+
+void printMatrix3X3(float mmm[3][3]);
+
+void vecPrint(float a[3]);
diff --git a/src/AltEst/altitude.cpp b/src/AltEst/altitude.cpp
new file mode 100644
index 0000000..8838b36
--- /dev/null
+++ b/src/AltEst/altitude.cpp
@@ -0,0 +1,58 @@
+/*
+    altitude.cpp: Altitude estimation via barometer/accelerometer fusion
+*/
+
+#include "filters.h"
+#include "algebra.h"
+#include "altitude.h"
+
+AltitudeEstimator::AltitudeEstimator(float sigmaAccel, float sigmaGyro, float sigmaBaro,
+                                     float ca, float accelThreshold)
+:kalman(ca, sigmaGyro, sigmaAccel), complementary(sigmaAccel, sigmaBaro, accelThreshold)
+{
+      this->sigmaAccel = sigmaAccel;
+      this->sigmaGyro = sigmaGyro;
+      this->sigmaBaro = sigmaBaro;
+      this->ca = ca;
+      this->accelThreshold = accelThreshold;
+}
+
+void AltitudeEstimator::estimate(float accel[3], float gyro[3], float baroHeight, uint32_t timestamp)
+{
+        float deltat = (float)(timestamp-previousTime)/1000000.0f;
+        float verticalAccel = kalman.estimate(pastGyro,
+                                              pastAccel,
+                                              deltat);
+        complementary.estimate(& estimatedVelocity,
+                               & estimatedAltitude,
+                               baroHeight,
+                               pastAltitude,
+                               pastVerticalVelocity,
+                               pastVerticalAccel,
+                               deltat);
+        // update values for next iteration
+        copyVector(pastGyro, gyro);
+        copyVector(pastAccel, accel);
+        pastAltitude = estimatedAltitude;
+        pastVerticalVelocity = estimatedVelocity;
+        pastVerticalAccel = verticalAccel;
+        previousTime = timestamp;
+}
+
+float AltitudeEstimator::getAltitude()
+{
+        // return the last estimated altitude
+        return estimatedAltitude;
+}
+
+float AltitudeEstimator::getVerticalVelocity()
+{
+        // return the last estimated vertical velocity
+        return estimatedVelocity;
+}
+
+float AltitudeEstimator::getVerticalAcceleration()
+{
+        // return the last estimated vertical acceleration
+        return pastVerticalAccel;
+}
diff --git a/src/AltEst/altitude.h b/src/AltEst/altitude.h
new file mode 100644
index 0000000..1ca6cb0
--- /dev/null
+++ b/src/AltEst/altitude.h
@@ -0,0 +1,55 @@
+/*
+    altitude.h: Altitude estimation via barometer/accelerometer fusion
+*/
+
+# pragma once
+
+#include "filters.h"
+#include "algebra.h"
+#include "pico/time.h"
+#include "pico/types.h"
+
+class AltitudeEstimator {
+
+  private:
+    // required parameters for the filters used for the estimations
+    // sensor's standard deviations
+    float sigmaAccel;
+    float sigmaGyro;
+    float sigmaBaro;
+    // Acceleration markov chain model state transition constant
+    float ca;
+    // Zero-velocity update acceleration threshold
+    float accelThreshold;
+    // gravity
+    float g = 9.81;
+    // For computing the sampling period
+    absolute_time_t prevTime = get_absolute_time();
+    uint32_t previousTime = to_us_since_boot(prevTime);
+    // required filters for altitude and vertical velocity estimation
+    KalmanFilter kalman;
+    ComplementaryFilter complementary;
+    // Estimated past vertical acceleration
+    float pastVerticalAccel = 0;
+    float pastVerticalVelocity = 0;
+    float pastAltitude = 0;
+    float pastGyro[3] = {0, 0, 0};
+    float pastAccel[3] = {0, 0, 0};
+    // estimated altitude and vertical velocity
+    float estimatedAltitude = 0;
+    float estimatedVelocity = 0;
+
+  public:
+
+    AltitudeEstimator(float sigmaAccel, float sigmaGyro, float sigmaBaro,
+                      float ca, float accelThreshold);
+
+    void estimate(float accel[3], float gyro[3], float baroHeight, uint32_t timestamp);
+
+    float getAltitude();
+
+    float getVerticalVelocity();
+
+    float getVerticalAcceleration();
+
+}; // class AltitudeEstimator
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);
+}
diff --git a/src/AltEst/filters.h b/src/AltEst/filters.h
new file mode 100644
index 0000000..2e316a3
--- /dev/null
+++ b/src/AltEst/filters.h
@@ -0,0 +1,65 @@
+/*
+   filters.h: Filter class declarations
+ */
+
+#pragma once
+
+//#include <cmath>
+#include <math.h>
+#include <stdint.h>
+
+#include "algebra.h"
+
+class KalmanFilter {
+  private:
+    float currentState[3] = {0, 0, 1};
+    float currErrorCovariance[3][3] = {{100, 0, 0},{0, 100, 0},{0, 0, 100}};
+    float H[3][3] = {{9.81, 0, 0}, {0, 9.81, 0}, {0, 0, 9.81}};
+    float previousAccelSensor[3] = {0, 0, 0};
+    float ca;
+    float sigmaGyro;
+    float sigmaAccel;
+
+    void getPredictionCovariance(float covariance[3][3], float previousState[3], float deltat);
+
+    void getMeasurementCovariance(float covariance[3][3]);
+
+    void predictState(float predictedState[3], float gyro[3], float deltat);
+
+    void predictErrorCovariance(float covariance[3][3], float gyro[3], float deltat);
+
+    void updateGain(float gain[3][3], float errorCovariance[3][3]);
+
+    void updateState(float updatedState[3], float predictedState[3], float gain[3][3], float accel[3]);
+
+    void updateErrorCovariance(float covariance[3][3], float errorCovariance[3][3], float gain[3][3]);
+
+  public:
+
+    KalmanFilter(float ca, float sigmaGyro, float sigmaAccel);
+
+    float estimate(float gyro[3], float accel[3], float deltat);
+
+}; // Class KalmanFilter
+
+class ComplementaryFilter {
+
+  private:
+
+    // filter gain
+    float gain[2];
+    // Zero-velocity update
+    float accelThreshold;
+    static const uint8_t ZUPT_SIZE = 12;
+    uint8_t ZUPTIdx;
+    float   ZUPT[ZUPT_SIZE];
+
+    float ApplyZUPT(float accel, float vel);
+
+  public:
+
+    ComplementaryFilter(float sigmaAccel, float sigmaBaro, float accelThreshold);
+
+    void estimate(float * velocity, float * altitude, float baroAltitude,
+                  float pastAltitude, float pastVelocity, float accel, float deltat);
+}; // Class ComplementaryFilter