Giter VIP home page Giter VIP logo

inertialnav's Introduction

Inertial Navigation Estimation Library

Build Status

Files for prototype 21, 22, 23 and 24 state Extended Kalman filters designed for APMPlane implementation Author: Paul Riseborough

This is an implementation of a strapdown inertial navigation system with an Extended Kalman Filter algorithm used to provide aiding using the following data sources (depending on filter variant):

3D velocity measurements (taken from UBlox6 or equivalent GPS) 2D position (taken from UBlox6 or equivalent GPS) Height measurement (could be GPS or barometric height or some combination) 3-axis body fixed magnetic flux measurements True airspeed measurement Terrain laser range finder measurements (aligned with vehicle Z axis) Ground optical flow measurements (flow angular rates about vehicle X and Y axes)

The IMU delta angles and delta velocities are assumed to be simple integrals of the corresponding angular rates and axial accelerations, with no coning or sculling compensation applied to them before input to the filter.

The filter estimates the following states: 4 quaternion parameters 3 North,East,Down velocity components 3 North,East,Down position components 3 IMU delta angle bias components 2 IMU delta velocity bias in X and Y (only included in 24 state derivation) 1 IMU delta velocity bias in Z (only included in the 22 and 23 state derivation) 2 North,East wind velocity components 3 North,East,Down earth magnetic flux components 3 X,Y,Z body fixed magnetic flux components (these are opposite sign to the compass offsets used by APM) 1 Offset of terrain along down axis (only included in the 23 state filter derivation)

The filter is designed to run in parallel with the existing APM AHRS complementary filter, firstly to provide a bootstrap for initial alignment, and secondly to provide a watchdog reference to detect filter divergence.

No optimisation of tuning parameters has been performed, other than to roughly set them based on innovation sequence RMS errors obtained from the test data set.

The three test data sets were obtained from a PX4 FMU and digital airspeed sensor, installed in a Skywalker X-8 airframe. The 2nd and 3rd data sets include aerobatic manoeuvres (loops and axial rolls), but no spinning due to airframe limitations.

Some additional data sets incorporating optical flow and range finder measurements are included.

Instructions To Run Simulink Model

Note : Simulink models are only available for 21 and 24 state architecture, and do not include range finder or optical flow measurements.

You will need Matlab + Simulink 2012a or later to run this model

Instructions to run:

  1. Add 'plots', 'scripts' and 'TestData' directories to your Matlab path
  2. Make 'models' your working directory
  3. Run the RunNavFilterTestHarness24.m script file to run the 24-state filter, or RunNavFilterTestHarness21.m to run the 21-state filter
  4. Test data will be loaded and the model will be built, run and plots generated.

You can load other test data by modifying the file load command at the top of the LoadNavFilterTestData.m script file.

Instructions To Run C++ code test harness

  1. Go to the code directory and run make
  2. Make sure the ATT, GPS, IMU, MAG and NTUN.txt files are in the directory with the executable (unzip one of the data ZIP files)
  3. Run: ./estimator_closed_loop_test
  4. The program will put the results into the following space deliminted data files
  5. Run make plots to generate Python plots
  • CovDataOut.txt
  • EulDataOut.txt
  • MagFuse.txt
  • RefVelPosDataOut.txt
  • StateDataOut.txt
  • TasFuse.txt
  • VelPosFuse.txt

inertialnav's People

Contributors

kd0aij avatar lorenzmeier avatar priseborough avatar thomasgubler avatar

Stargazers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

Watchers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

inertialnav's Issues

Coning term is always zero

In the code:

// Save current measurements
    Vector3f  prevDelAng = correctedDelAng;

// Apply corrections for earths rotation rate and coning errors
// * and + operators have been overloaded
    correctedDelAng   = correctedDelAng - Tnb*earthRateNED*dtIMU + 8.333333333333333e-2f*(prevDelAng % correctedDelAng);

the cross product prevDelAng % correctedDelAng will always be zero as they have the same value

Why the magMeas is the sum up of inertial frame mag and body frame mag

Hello, Paul
I am very inteseted in your work and had a little doubt here.

As shown in the matlab code, here you have

%% derive equations for fusion of magnetic field measurement
magMeas = transpose(Tbn)*[magN;magE;magD] + [magX;magY;magZ]; % predicted measurement

I think the magnetometer sensor only measures [magX;magY;magZ], which is earth mag field in body frame.
so I think

magMeas = transpose(Tbn)*[magN;magE;magD]

for the observation function would be correct

and then

innovation = magMeas - [magX;magY;magZ]

Could you explain a little about this?

Help in assembling for platform

Good day @priseborough and @thomasgubler
Sorry for the long text...
I am rewriting the code to assemble it under the stm32f411ceu6 or BlackPill in the Arduino ecosystem on the platform + vscode.
Everything is standard, I added the file main.cpp and includes.
Only instead of

int main() {
 while() {
 }
}

I write in Arduino style

#include <Arduino.h>

void setup() {
 // Your initialization code is here
}

void loop() {
 // Your main program code is here
}

The first compilation showed that the compiler could not determine setup and loop
Added inclusion

#include "Arduino.h"

at the very beginning the code
The second compilation showed that

static inline float sq(float valIn) {return valIn * valIn;}

already used in the file

.platformio\packages\framework-arduinoststm32\cores\arduino\wiring_constants.h

I had to comment it out

// error in inertialNAV
// #define sq(x) ((x)*(x))

The third compilation gave an error about missing millis() definition in the file estimator_22states.h

uint32_t millis();

I also had to add it to this file
#include "Arduino.h"
The fourth compilation also complained about the same file estimator_22states.h
Commented out getMicros

 // uint64_t tNow = getMicros();
 uint64_t tNow = micros();

Well, the fifth compilation put me in a stupor, because I couldn’t find a solution

.platformio/packages/toolchain-gccarmnoneeabi/bin/../lib/gcc/arm-none-eabi/12.3.1/../../../../arm-none-eabi/lib/thumb/ v7e-m+fp/hard\libc_nano.a(libc_a-openr.o): in function _open_r': openr.c:(.text._open_r+0x10): undefined reference to _open'

I think we're talking about this file:

.platformio\packages\toolchain-gccarmnoneeabi\arm-none-eabi\lib\thumb\v7e-m+fp\hard\libc_nano.a

As far as I understand, the InertialNav code was written almost 10 years ago, during the time of the old C++ standards, and today's gcc arm-none-eabi compilers and toolchain no longer want to process such old implementations.
Perhaps there is already code adapted for platformio, or someone in the future will be puzzled by the same problem and find this issue.
In any case, the discussion will not be useless.
Code below

#include "Arduino.h"

#include "estimator_22states.h"

#include <stdint.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>

bool endOfData = false; //boolean set to true when all files have returned data

// Estimated time delays (msec)
uint32_t msecVelDelay = 230;
uint32_t msecPosDelay = 210;
uint32_t msecHgtDelay = 350;
uint32_t msecRngDelay = 100;
uint32_t msecMagDelay = 30;
uint32_t msecTasDelay = 210;
uint32_t msecOptFlowDelay = 55;

// IMU input data variables
float imuIn;
float tempImu[8];
float IMUtimestamp;
static uint32_t IMUmsec = 0;
static uint64_t IMUusec = 0;

// Magnetometer input data variables
float magIn;
float tempMag[8];
float tempMagPrev[8];
float posNED[3];
float MAGtimestamp = 0;
uint32_t MAGmsec = 0;
uint32_t lastMAGmsec = 0;
bool newDataMag = false;

// AHRS input data variables
float ahrsIn;
float tempAhrs[7];
float tempAhrsPrev[7];
float AHRStimestamp = 0;
uint32_t AHRSmsec = 0;
uint32_t lastAHRStime = 0;
float ahrsEul[3];
float ahrsErrorRP;
float ahrsErrorYaw;
float eulerEst[3]; // Euler angles calculated from filter states
float eulerDif[3]; // difference between Euler angle estimated by EKF and the AHRS solution
float gpsRaw[7];

// ADS input data variables
float adsIn;
float tempAds[10];
float tempAdsPrev[10];
float ADStimestamp = 0;
uint32_t ADSmsec = 0;
uint32_t lastADSmsec = 0;
float Veas;
bool newAdsData = false;
bool newDataGps = false;
bool newOptFlowData = false;

float onboardTimestamp = 0;
uint32_t onboardMsec = 0;
uint32_t lastOnboardMsec = 0;
bool newOnboardData = false;

float onboardPosNED[3];
float onboardVelNED[3];
float onboardLat;
float onboardLon;
float onboardHgt;

uint32_t flowMsec = 0;
uint32_t lastFlowMsec = 0;
bool newFlowData = false;

float flowTimestamp;      // in seconds
float flowRawPixelX;       // in pixels
float flowRawPixelY;       // in pixels
float flowDistance;        // in meters
float flowQuality;   // distance normalized between 0 and 1
float flowSensorId;        // sensor ID
float flowGyroX = 0.0f;
float flowGyroY = 0.0f;
float flowGyroBiasX = 0.0f;
float flowGyroBiasY = 0.0f;

float flowRadX;
float flowRadY;

float flowRawGroundSpeedX;
float flowRawGroundSpeedY;

uint32_t distMsec = 0;
uint32_t lastDistMsec = 0;
bool newDistData = false;

float distTimestamp = 0.0f;
bool distValid = false;
float distGroundDistance;
float distLastValidReading;

// input data timing
uint64_t msecAlignTime;
uint64_t msecStartTime;
uint64_t msecEndTime;

float gpsGndSpd;
float gpsCourse;
float gpsVelD;

bool resetTests;
float dt;

bool timeoutTested;
bool nanTested;
bool hgtTested;

AttPosEKF   *_ekf;

// Data file identifiers
FILE * pImuFile = nullptr;
FILE * pMagFile = nullptr;
FILE * pGpsFile = nullptr;
FILE * pAhrsFile = nullptr;
FILE * pAdsFile = nullptr;
FILE * pStateOutFile = nullptr;
FILE * pGlitchOutFile = nullptr;
FILE * pEulOutFile = nullptr;
FILE * pCovOutFile = nullptr;
FILE * pRefPosVelOutFile = nullptr;
FILE * pVelPosFuseFile = nullptr;
FILE * pMagFuseFile = nullptr;
FILE * pTasFuseFile = nullptr;
FILE * pRngFuseFile = nullptr;
FILE * pOptFlowFuseFile = nullptr;
FILE * pTimeFile = nullptr;
FILE * pGpsRawOUTFile = nullptr;
FILE * pGpsRawINFile = nullptr;
FILE * validationOutFile = nullptr;
FILE * pOnboardPosVelOutFile = nullptr;
FILE * pOnboardFile = nullptr;
FILE * pInFlowFile = nullptr;
FILE * pInDistFile = nullptr;
FILE * pOutFlowFile = nullptr;

FILE * open_with_exit(const char* filename, const char* flags)
{
    FILE *f = fopen(filename, flags);

    if (!f) {
        printf("FAILED TO OPEN FILE: %s\n", filename);
        exit(1);
    }

    return f;
}

int printstates() {
    printf("States:\n");
    unsigned i = 0;
    printf("Quaternion:\n");
    for (; i<4; i++)
    {
        printf(" %e", _ekf->states[i]);
    }
    printf("\n");
    for (; i<4+6; i++)
    {
        printf(" %e", _ekf->states[i]);
    }
    printf("\n");
    for (; i<4+6+6; i++)
    {
        printf(" %e", _ekf->states[i]);
    }
    printf("\n");
    for (; i < sizeof(_ekf->states) / sizeof(_ekf->states[0]); i++)
    {
        printf(" %e", _ekf->states[i]);
    }
    printf("\n");

    return 0;
}

// int main(int argc, char *argv[])
// {
void setup() {
    // Instantiate EKF
    _ekf = new AttPosEKF();

    // open data files
    pImuFile = open_with_exit ("IMU.txt","r");
    pMagFile = open_with_exit ("MAG.txt","r");
    pGpsFile = open_with_exit ("GPS.txt","r");
    pAhrsFile = open_with_exit ("ATT.txt","r");
    pAdsFile = open_with_exit ("NTUN.txt","r");
    pTimeFile = open_with_exit ("timing.txt","r");
    pGlitchOutFile = open_with_exit ("GlitchOffsetOut.txt","w");
    pStateOutFile = open_with_exit ("StateDataOut.txt","w");
    pEulOutFile = open_with_exit ("EulDataOut.txt","w");
    pCovOutFile = open_with_exit ("CovDataOut.txt","w");
    pRefPosVelOutFile = open_with_exit ("RefVelPosDataOut.txt","w");
    pVelPosFuseFile = open_with_exit ("VelPosFuse.txt","w");
    pMagFuseFile = open_with_exit ("MagFuse.txt","w");
    pTasFuseFile = open_with_exit ("TasFuse.txt","w");
    pRngFuseFile = open_with_exit ("RngFuse.txt","w");
    pOptFlowFuseFile = open_with_exit ("OptFlowFuse.txt","w");
    pGpsRawINFile = fopen ("GPSraw.txt","r");
    pGpsRawOUTFile = open_with_exit ("GPSrawOut.txt","w");
    validationOutFile = fopen("ValidationOut.txt", "w");
    pOnboardFile = fopen ("GPOSonboard.txt","r");
    pOnboardPosVelOutFile = open_with_exit ("OnboardVelPosDataOut.txt","w");

    pInFlowFile = fopen ("FLOW.txt","r");
    pInDistFile = fopen ("DIST.txt","r");
    pOutFlowFile = fopen ("FlowRawOut.txt","w");

    printf("Filter start\n");

    memset(gpsRaw, 0, sizeof(gpsRaw));

    // readTimingData();

    printf("First data instances loaded\n");

    dt = 0.0f; // time lapsed since last covariance prediction

    // bool resetTests = false;

    // Test resets
    // if (argc > 1 && (strcmp(argv[1], "reset") == 0)) {
        resetTests = true;
    // }

    timeoutTested = false;
    nanTested = false;
    hgtTested = false;

} // SETUP

    // while (true) {
void loop () {

        // read test data from files for next timestamp
        unsigned nreads = 1;

        // Decide wether to perform any reset tests

        if (resetTests) {

            // Trigger a NaN reset after 25% of the log
            if (!nanTested && (IMUmsec > (msecEndTime - msecStartTime) / 4)) {
                _ekf->states[0] = 0.0f / 0.0f;
                _ekf->states[9] = 0.0f / 0.0f;
                nanTested = true;
                printf("WARNING: TRIGGERING NAN STATE ON PURPOSE!\n");
            }

            // Trigger a timeout at half the log
            if (!timeoutTested && (IMUmsec > (msecEndTime - msecStartTime) / 2)) {
                nreads = 20;
                timeoutTested = true;
                printf("WARNING: TRIGGERING SENSOR TIMEOUT ON PURPOSE!\n");
            }

            // Trigger a NaN altitude at 12.5% of the log
            if (!hgtTested && (IMUmsec > (msecEndTime - msecStartTime) / 8)) {
                _ekf->hgtMea = 0.0f / 0.0f;
                hgtTested = true;
                printf("WARNING: TRIGGERING NaN ALT MEASUREMENT ON PURPOSE!\n");
            }
        }

        // Supporting multiple reads at once to support timeout simulation.
        // The default is however to only read one dataset at a time

        // We need to re-do the dtIMU calculation here so that
        // dtIMU correctly skips if we skip readings.
        uint64_t IMUmsecPrev = IMUmsec;
        for (unsigned i = 0; i < nreads; i++) {
            // readIMUData();
            // readGpsData();
            // readMagData();
            // readAirData();
            // readAhrsData();
            // readOnboardData();
            // readFlowData();
            // readDistData();
        }

        // Apply dtIMU here after 1 or more reads, simulating skipped sensor
        // readings if desired.
        _ekf->dtIMU     = 0.001f*(IMUmsec - IMUmsecPrev);

        if (IMUmsec > msecEndTime || endOfData)
        {
            printf("Reached end at %8.4f seconds (end of logfile reached: %s)", IMUmsec/1000.0f, (endOfData) ? "YES" : "NO");
            // CloseFiles();
            // break;
        }

        if ((IMUmsec >= msecStartTime) && (IMUmsec <= msecEndTime))
        {
            // Initialise states, covariance and other data
            if ((IMUmsec > msecAlignTime) && !_ekf->statesInitialised && (_ekf->GPSstatus == 3))
            {
                if (pGpsRawINFile != nullptr)
                {
                    _ekf->velNED[0] = gpsRaw[4];
                    _ekf->velNED[1] = gpsRaw[5];
                    _ekf->velNED[2] = gpsRaw[6];
                }
                else
                {
                    calcvelNED(_ekf->velNED, gpsCourse, gpsGndSpd, gpsVelD);
                }
                _ekf->InitialiseFilter(_ekf->velNED, _ekf->gpsLat, _ekf->gpsLon, _ekf->gpsHgt, 0.0f);

            } else if ((IMUmsec > msecAlignTime) && !_ekf->statesInitialised) {

                float initVelNED[3];

                initVelNED[0] = 0.0f;
                initVelNED[1] = 0.0f;
                initVelNED[2] = 0.0f;
                _ekf->posNE[0] = posNED[0];
                _ekf->posNE[1] = posNED[1];

                _ekf->InitialiseFilter(initVelNED, 0.0, 0.0, 0.0f, 0.0f);

            } else if (_ekf->statesInitialised) {

                // Run the strapdown INS equations every IMU update
                _ekf->UpdateStrapdownEquationsNED();
                #if 1
                // debug code - could be turned into a filter monitoring/watchdog function
                float tempQuat[4];
                for (uint8_t j=0; j<4; j++) tempQuat[j] = _ekf->states[j];
                _ekf->quat2eul(eulerEst, tempQuat);
                for (uint8_t j=0; j<=2; j++) eulerDif[j] = eulerEst[j] - ahrsEul[j];
                if (eulerDif[2] > M_PI) eulerDif[2] -= 2 * M_PI;
                if (eulerDif[2] < -M_PI) eulerDif[2] += 2 * M_PI;
                #endif
                // store the predicted states for subsequent use by measurement fusion
                _ekf->StoreStates(IMUmsec);
                // Check if on ground - status is used by covariance prediction
                bool onground = (((AttPosEKF::sq(_ekf->velNED[0]) + AttPosEKF::sq(_ekf->velNED[1]) + AttPosEKF::sq(_ekf->velNED[2])) < 4.0f) && (_ekf->VtasMeas < 8.0f));

                _ekf->setOnGround(onground);
                // sum delta angles and time used by covariance prediction
                _ekf->summedDelAng = _ekf->summedDelAng + _ekf->correctedDelAng;
                _ekf->summedDelVel = _ekf->summedDelVel + _ekf->dVelIMU;
                dt += _ekf->dtIMU;
                // perform a covariance prediction if the total delta angle has exceeded the limit
                // or the time limit will be exceeded at the next IMU update
                if ((dt >= (_ekf->covTimeStepMax - _ekf->dtIMU)) || (_ekf->summedDelAng.length() > _ekf->covDelAngMax))
                {
                    _ekf->CovariancePrediction(dt);
                    _ekf->summedDelAng.zero();
                    _ekf->summedDelVel.zero();
                    dt = 0.0f;
                }

                // Set global time stamp used by EKF processes
                _ekf->globalTimeStamp_ms = IMUmsec;
            }

            // Fuse optical flow measurements
            if (newFlowData && _ekf->statesInitialised && _ekf->useOpticalFlow && flowQuality > 0.5 && _ekf->Tnb.z.z > 0.71f)
            {
                // recall states and angular rates stored at time of measurement after adjusting for delays
                _ekf->RecallStates(_ekf->statesAtFlowTime, (IMUmsec - msecOptFlowDelay));
                _ekf->RecallOmega(_ekf->omegaAcrossFlowTime, (IMUmsec - 2*msecOptFlowDelay));

                // Calculate bias errorsfor flow sensor internal gyro
                flowGyroBiasX = 0.999f * flowGyroBiasX + 0.001f * (flowGyroX - _ekf->omegaAcrossFlowTime[0]);
                flowGyroBiasY = 0.999f * flowGyroBiasY + 0.001f * (flowGyroY - _ekf->omegaAcrossFlowTime[1]);

                //use sensor internal rates corrected for bias errors
                _ekf->omegaAcrossFlowTime[0] = flowGyroX - flowGyroBiasX;
                _ekf->omegaAcrossFlowTime[1] = flowGyroY - flowGyroBiasY;

                // calculate rotation matrix
                // Copy required states to local variable names
                float q0 = _ekf->statesAtFlowTime[0];
                float q1 = _ekf->statesAtFlowTime[1];
                float q2 = _ekf->statesAtFlowTime[2];
                float q3 = _ekf->statesAtFlowTime[3];
                float q00 = _ekf->sq(q0);
                float q11 = _ekf->sq(q1);
                float q22 = _ekf->sq(q2);
                float q33 = _ekf->sq(q3);
                float q01 = q0 * q1;
                float q02 = q0 * q2;
                float q03 = q0 * q3;
                float q12 = q1 * q2;
                float q13 = q1 * q3;
                float q23 = q2 * q3;
                _ekf->Tnb_flow.x.x = q00 + q11 - q22 - q33;
                _ekf->Tnb_flow.y.y = q00 - q11 + q22 - q33;
                _ekf->Tnb_flow.z.z = q00 - q11 - q22 + q33;
                _ekf->Tnb_flow.y.x = 2*(q12 - q03);
                _ekf->Tnb_flow.z.x = 2*(q13 + q02);
                _ekf->Tnb_flow.x.y = 2*(q12 + q03);
                _ekf->Tnb_flow.z.y = 2*(q23 - q01);
                _ekf->Tnb_flow.x.z = 2*(q13 - q02);
                _ekf->Tnb_flow.y.z = 2*(q23 + q01);

                // scale from raw pixel flow rate to radians/second
                //float scaleFactor = 0.03f; // best value for quad106.zip data using the 16 mm lens
                //float scaleFactor = 0.06f; // best value for InputFilesPX4_flow.zip data
                //float scaleFactor = 0.882f; // best value for quad123.zip data which outputs flow rates that have already been scaled to rad/sec
                float scaleFactor = 1.000f; // best value for quadOptFlowLogData.zip data which outputs flow rates that have already been scaled to rad/sec
                flowRadX = flowRawPixelX * scaleFactor;
                flowRadY = flowRawPixelY * scaleFactor;

                // calculate motion compensated angular flow rates used for fusion in the main nav filter
                _ekf->flowRadXYcomp[0] = flowRadX/_ekf->flowStates[0] + _ekf->omegaAcrossFlowTime[0];
                _ekf->flowRadXYcomp[1] = flowRadY/_ekf->flowStates[0] + _ekf->omegaAcrossFlowTime[1];

                // these flow rates are not motion compensated and are used for focal length scale factor estimation
                _ekf->flowRadXY[0] = flowRadX;
                _ekf->flowRadXY[1] = flowRadY;

                // perform optical flow fusion
                _ekf->fuseOptFlowData = true;
                _ekf->fuseRngData = false;

                // don't try to estimate focal length scale factor if GPS is not being used or there is no range finder.
                if (_ekf->useGPS && _ekf->useRangeFinder) {
                    _ekf->OpticalFlowEKF();
                }
                _ekf->FuseOptFlow();
                _ekf->fuseOptFlowData = false;

                // estimate speed over ground for cross-check of data (debugging only)
                float tempQuat[4];
                float euler[3];
                for (uint8_t j=0; j<4; j++) tempQuat[j] = _ekf->states[j];
                _ekf->quat2eul(euler, tempQuat);
                float bx = (flowRadX - _ekf->angRate.x) * distLastValidReading;
                float by = (flowRadY - _ekf->angRate.y) * distLastValidReading;
                flowRawGroundSpeedY = cos(euler[2]) * bx + -sin(euler[2]) * by;
                flowRawGroundSpeedX = -(sin(euler[2]) * bx + cos(euler[2]) * by);
            } else {
                _ekf->fuseOptFlowData = false;
            }

            // Fuse Ground distance Measurements
            if (newDistData && _ekf->statesInitialised && _ekf->useRangeFinder)
            {
                if (distValid > 0.0f && _ekf->Tnb.z.z > 0.9f) {
                    distLastValidReading = distGroundDistance;
                    _ekf->rngMea = distGroundDistance;
                    _ekf->fuseRngData = true;
                    _ekf->fuseOptFlowData = false;
                    _ekf->RecallStates(_ekf->statesAtRngTime, (IMUmsec - msecRngDelay));
                    _ekf->OpticalFlowEKF();
                    _ekf->fuseRngData = false;
                }
            }

                // Fuse GPS Measurements
                if (newDataGps)
                {
                    // calculate a position offset which is applied to NE position and velocity wherever it is used throughout code to allow GPS position jumps to be accommodated gradually
                    _ekf->decayGpsOffset();

                    // Convert GPS measurements to Pos NE, hgt and Vel NED
                    if (pGpsRawINFile != nullptr)
                    {
                        _ekf->velNED[0] = gpsRaw[4];
                        _ekf->velNED[1] = gpsRaw[5];
                        _ekf->velNED[2] = gpsRaw[6];
                    }
                    else
                    {
                        calcvelNED(_ekf->velNED, gpsCourse, gpsGndSpd, gpsVelD);
                    }
                    calcposNED(posNED, _ekf->gpsLat, _ekf->gpsLon, _ekf->gpsHgt, _ekf->latRef, _ekf->lonRef, _ekf->hgtRef);

                    if (pOnboardFile != nullptr) {
                        calcposNED(onboardPosNED, onboardLat, onboardLon, onboardHgt, _ekf->latRef, _ekf->lonRef, _ekf->hgtRef);

                    }

                    _ekf->posNE[0] = posNED[0];
                    _ekf->posNE[1] = posNED[1];
                     // set fusion flags
                    _ekf->fuseVelData = true;
                    _ekf->fusePosData = true;
                    // recall states stored at time of measurement after adjusting for delays
                    _ekf->RecallStates(_ekf->statesAtVelTime, (IMUmsec - msecVelDelay));
                    _ekf->RecallStates(_ekf->statesAtPosTime, (IMUmsec - msecPosDelay));
                    // run the fusion step
                    _ekf->FuseVelposNED();
                    printf("FuseVelposNED at time = %e \n", IMUtimestamp);
                }
                else
                {
                    _ekf->fuseVelData = false;
                    _ekf->fusePosData = false;
                }
                calcposNED(posNED, _ekf->gpsLat, _ekf->gpsLon, _ekf->gpsHgt, _ekf->latRef, _ekf->lonRef, _ekf->hgtRef);

                // if (pOnboardFile > 0) {
                if (pOnboardFile != nullptr) {
                    calcposNED(onboardPosNED, onboardLat, onboardLon, onboardHgt, _ekf->latRef, _ekf->lonRef, _ekf->hgtRef);
                }

                _ekf->posNE[0] = posNED[0];
                _ekf->posNE[1] = posNED[1];

                // fuse GPS
                if (_ekf->useGPS && IMUmsec < 1000) {
                    _ekf->fuseVelData = true;
                    _ekf->fusePosData = true;
                    _ekf->fuseHgtData = false;
                    // recall states stored at time of measurement after adjusting for delays
                    _ekf->RecallStates(_ekf->statesAtVelTime, (IMUmsec - msecVelDelay));
                    _ekf->RecallStates(_ekf->statesAtPosTime, (IMUmsec - msecPosDelay));
                    // record the last fix time
                    _ekf->lastFixTime_ms = IMUmsec;
                    // run the fusion step
                    _ekf->FuseVelposNED();
                } else {
                    _ekf->fuseVelData = false;
                    _ekf->fusePosData = false;
                    _ekf->fuseHgtData = false;
                }
            }
            else
            {
                _ekf->fuseVelData = false;
                _ekf->fusePosData = false;
                _ekf->fuseHgtData = false;
            }

            if (newAdsData && _ekf->statesInitialised)
            {
                // Could use a blend of GPS and baro alt data if desired
                _ekf->hgtMea = 1.0f*_ekf->baroHgt + 0.0f*_ekf->gpsHgt - _ekf->hgtRef - _ekf->baroHgtOffset;
                _ekf->fuseVelData = false;
                _ekf->fusePosData = false;
                _ekf->fuseHgtData = true;
                // recall states stored at time of measurement after adjusting for delays
                _ekf->RecallStates(_ekf->statesAtHgtTime, (IMUmsec - msecHgtDelay));
                // run the fusion step
                _ekf->FuseVelposNED();
//                printf("time = %e \n", IMUtimestamp);
            }
            else
            {
                _ekf->fuseVelData = false;
                _ekf->fusePosData = false;
                _ekf->fuseHgtData = false;
            }

            // Fuse Magnetometer Measurements
            if (newDataMag && _ekf->statesInitialised && _ekf->useCompass)
            {
                _ekf->fuseMagData = true;
                _ekf->RecallStates(_ekf->statesAtMagMeasTime, (IMUmsec - msecMagDelay)); // Assume 50 msec avg delay for magnetometer data
                _ekf->magstate.obsIndex = 0;
                _ekf->FuseMagnetometer();
                _ekf->FuseMagnetometer();
                _ekf->FuseMagnetometer();
            }
            else
            {
                _ekf->fuseMagData = false;
            }

            // Fuse Airspeed Measurements
            if (newAdsData && _ekf->statesInitialised && _ekf->VtasMeas > 8.0f && _ekf->useAirspeed)
            {
                _ekf->fuseVtasData = true;
                _ekf->RecallStates(_ekf->statesAtVtasMeasTime, (IMUmsec - msecTasDelay)); // assume 100 msec avg delay for airspeed data
                _ekf->FuseAirspeed();
            }
            else
            {
                _ekf->fuseVtasData = false;
            }

                struct ekf_status_report ekf_report;

                /*
                 *    CHECK IF THE INPUT DATA IS SANE
                 */
                int check = _ekf->CheckAndBound(&ekf_report);

                switch (check) {
                    case 0:
                        /* all ok */
                        break;
                    case 1:
                    {
                        printf("NaN in states, resetting\n");
                        printf("fail states: ");
                        for (unsigned i = 0; i < ekf_report.n_states; i++) {
                            printf("%f ",ekf_report.states[i]);
                        }
                        printf("\n");

                        printf("states after reset: ");
                        for (unsigned i = 0; i < ekf_report.n_states; i++) {
                            printf("%f ",_ekf->states[i]);
                        }
                        printf("\n");
                        break;
                    }
                    case 2:
                    {
                        printf("stale IMU data, resetting\n");
                        break;
                    }
                    case 3:
                    {
                        printf("switching to dynamic state\n");
                        break;
                    }
                    case 4:
                    {
                        printf("excessive gyro offsets\n");
                        break;
                    }
                    case 5:
                    {
                        printf("GPS velocity diversion\n");
                        break;
                    }
                    case 6:
                    {
                        printf("Excessive covariances\n");
                        break;
                    }


                    default:
                    {
                        printf("unknown reset condition\n");
                    }
                }

                if (check) {
                    printf("RESET OCCURED AT %d milliseconds\n", (int)IMUmsec);

                    if (!ekf_report.velHealth || !ekf_report.posHealth || !ekf_report.hgtHealth || ekf_report.gyroOffsetsExcessive) {
                    printf("health: VEL:%s POS:%s HGT:%s OFFS:%s\n",
                        ((ekf_report.velHealth) ? "OK" : "ERR"),
                        ((ekf_report.posHealth) ? "OK" : "ERR"),
                        ((ekf_report.hgtHealth) ? "OK" : "ERR"),
                        ((!ekf_report.gyroOffsetsExcessive) ? "OK" : "ERR"));
                    }

                    if (ekf_report.velTimeout || ekf_report.posTimeout || ekf_report.hgtTimeout || ekf_report.imuTimeout) {
                        printf("timeout: %s%s%s%s\n",
                            ((ekf_report.velTimeout) ? "VEL " : ""),
                            ((ekf_report.posTimeout) ? "POS " : ""),
                            ((ekf_report.hgtTimeout) ? "HGT " : ""),
                            ((ekf_report.imuTimeout) ? "IMU " : ""));
                    }
                }

                // debug output
                //printf("Euler Angle Difference = %3.1f , %3.1f , %3.1f deg\n", rad2deg*eulerDif[0],rad2deg*eulerDif[1],rad2deg*eulerDif[2]);
                if ((IMUmsec >= msecStartTime) && (IMUmsec <= msecEndTime))
                {
                	// WriteFilterOutput();
                }

            // } // WHILE

            // State vector:
            // 0-3: quaternions (q0, q1, q2, q3)
            // 4-6: Velocity - m/sec (North, East, Down)
            // 7-9: Position - m (North, East, Down)
            // 10-12: Delta Angle bias - rad (X,Y,Z)
            // 13: Delta Velocity Z bias -m/s
            // 14-15: Wind Vector  - m/sec (North,East)
            // 16-18: Earth Magnetic Field Vector - milligauss (North, East, Down)
            // 19-21: Body Magnetic Field Vector - milligauss (X,Y,Z)
            // 22: Terrain Vertical Offset - m

            // printf("\n");
            // printf("dtIMU: %8.6f, dt: %8.6f, imuMsec: %u\n", _ekf->dtIMU, dt, IMUmsec);
            // printf("posNED: %8.4f, %8.4f, %8.4f, velNED: %8.4f, %8.4f, %8.4f\n", (double)_ekf->posNED[0], (double)_ekf->posNED[1], (double)_ekf->posNED[2],
            //     (double)_ekf->velNED[0], (double)_ekf->velNED[1], (double)_ekf->velNED[2]);
            // printf("vTAS: %8.4f baro alt: %8.4f\n", _ekf->VtasMeas, _ekf->hgtMea);
            // printf("mag: %8.4f, %8.4f, %8.4f\n", (double)_ekf->magData.x, (double)_ekf->magData.y, (double)_ekf->magData.z);
            // printf("states (quat)        [1-4]: %8.4f, %8.4f, %8.4f, %8.4f\n", (double)_ekf->states[0], (double)_ekf->states[1], (double)_ekf->states[2], (double)_ekf->states[3]);
            // printf("states (vel m/s)     [5-7]: %8.4f, %8.4f, %8.4f\n", (double)_ekf->states[4], (double)_ekf->states[5], (double)_ekf->states[6]);
            // printf("states (pos m)      [8-10]: %8.4f, %8.4f, %8.4f\n", (double)_ekf->states[7], (double)_ekf->states[8], (double)_ekf->states[9]);
            // printf("states (delta ang) [11-13]: %8.4f, %8.4f, %8.4f\n", (double)_ekf->states[10], (double)_ekf->states[11], (double)_ekf->states[12]);
            // printf("states (delta vel) [14]: %8.4ff\n", (double)_ekf->states[13]);
            // printf("states (wind)      [15-16]: %8.4f, %8.4f\n", (double)_ekf->states[14], (double)_ekf->states[15]);
            // printf("states (earth mag) [17-19]: %8.4f, %8.4f, %8.4f\n", (double)_ekf->states[16], (double)_ekf->states[17], (double)_ekf->states[18]);
            // printf("states (body mag)  [20-22]: %8.4f, %8.4f, %8.4f\n", (double)_ekf->states[19], (double)_ekf->states[20], (double)_ekf->states[21]);
            // printf("states (terain offset) [23]: %8.4ff\n", (double)_ekf->states[22]);
            // printf("states: %s %s %s %s %s %s %s %s %s\n",
            //     (_ekf->statesInitialised) ? "INITIALIZED" : "NON_INIT",
            //     (_ekf->onGround) ? "ON_GROUND" : "AIRBORNE",
            //     (_ekf->fuseVelData) ? "FUSE_VEL" : "INH_VEL",
            //     (_ekf->fusePosData) ? "FUSE_POS" : "INH_POS",
            //     (_ekf->fuseHgtData) ? "FUSE_HGT" : "INH_HGT",
            //     (_ekf->fuseMagData) ? "FUSE_MAG" : "INH_MAG",
            //     (_ekf->fuseVtasData) ? "FUSE_VTAS" : "INH_VTAS",
            //     (_ekf->useAirspeed) ? "USE_AIRSPD" : "IGN_AIRSPD",
            //     (_ekf->useCompass) ? "USE_COMPASS" : "IGN_COMPASS");

    delete _ekf;

    printf("\n\nSuccess: Finished processing complete dataset. Text files written.\n");
} // MAIN or LOOP

I did not touch the files, with the exception of millis and micros.
src

P.S.
Assembly by command make also does not work now.

make clean
make

mkdir -p obj
g++ -c -o obj/main_closed_loop_float.o main_closed_loop_float.cpp -I. -D__EXPORT="" -Dnullptr="0" -g -fno-strict-aliasing -fno-strength-reduce -fomit-frame-pointer -funsafe-math-optimizations -Wall -Wextra -Wshadow -Wfloat-equal -Wpointer-arith -Wpacked -Wno-unused-parameter -Werror=format-security -Wmissing-field-initializers -Wfatal-errors -Wno-missing-field-initializers -std=c++11
mkdir -p obj
g++ -c -o obj/estimator_22states.o estimator_22states.cpp -I. -D__EXPORT="" -Dnullptr="0" -g -fno-strict-aliasing -fno-strength-reduce -fomit-frame-pointer -funsafe-math-optimizations -Wall -Wextra -Wshadow -Wfloat-equal -Wpointer-arith -Wpacked -Wno-unused-parameter -Werror=format-security -Wmissing-field-initializers -Wfatal-errors -Wno-missing-field-initializers -std=c++11
In file included from estimator_22states.cpp:40:
estimator_22states.h: In constructor ‘AttPosEKF::AttPosEKF()’:
estimator_22states.h:260:11: warning: ‘AttPosEKF::moCompR_LOS’ will be initialized after [-Wreorder]
260 | float moCompR_LOS; // scaler from sensor gyro rate to uncertainty in LOS rate
| ^~~~~~~~~~~
estimator_22states.h:163:13: warning: ‘int16_t AttPosEKF::gpsGlitchAccel’ [-Wreorder]
163 | int16_t gpsGlitchAccel; // Magnitude of discrepancy between inertial and GPS Horizontal acceleration above which GPS data is ignored : cm/s^2
| ^~~~~~~~~~~~~~
estimator_22states.cpp:57:1: warning: when initialized here [-Wreorder]
57 | AttPosEKF::AttPosEKF() :
| ^~~~~~~~~
In file included from estimator_22states.cpp:40:
estimator_22states.h:166:13: warning: ‘AttPosEKF::gpsPosInnovNSTD’ will be initialized after [-Wreorder]
166 | int8_t gpsPosInnovNSTD; // Number of standard deviations allowed in GPS position innovation consistency check
| ^~~~~~~~~~~~~~~
estimator_22states.h:165:13: warning: ‘int8_t AttPosEKF::gpsVelInnovNSTD’ [-Wreorder]
165 | int8_t gpsVelInnovNSTD; // Number of standard deviations allowed in GPS velocity innovation consistency check
| ^~~~~~~~~~~~~~~
estimator_22states.cpp:57:1: warning: when initialized here [-Wreorder]
57 | AttPosEKF::AttPosEKF() :
| ^~~~~~~~~
estimator_22states.cpp: In member function ‘void AttPosEKF::ZeroVariables()’:
estimator_22states.cpp:3398:11: warning: ‘void* memset(void*, int, size_t)’ clearing an object of non-trivial type ‘struct AttPosEKF::mag_state_struct’; use assignment or value-initialization instead [-Wclass-memaccess]
3398 | memset(&magstate, 0, sizeof(magstate));
| ~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
In file included from estimator_22states.cpp:40:
estimator_22states.h:92:12: note: ‘struct AttPosEKF::mag_state_struct’ declared here
92 | struct mag_state_struct {
| ^~~~~~~~~~~~~~~~
mkdir -p obj
g++ -c -o obj/estimator_utilities.o estimator_utilities.cpp -I. -D__EXPORT="" -Dnullptr="0" -g -fno-strict-aliasing -fno-strength-reduce -fomit-frame-pointer -funsafe-math-optimizations -Wall -Wextra -Wshadow -Wfloat-equal -Wpointer-arith -Wpacked -Wno-unused-parameter -Werror=format-security -Wmissing-field-initializers -Wfatal-errors -Wno-missing-field-initializers -std=c++11
g++ -g -o estimator_closed_loop_test obj/main_closed_loop_float.o obj/estimator_22states.o obj/estimator_utilities.o -I. -D__EXPORT="" -Dnullptr="0" -g -fno-strict-aliasing -fno-strength-reduce -fomit-frame-pointer -funsafe-math-optimizations -Wall -Wextra -Wshadow -Wfloat-equal -Wpointer-arith -Wpacked -Wno-unused-parameter -Werror=format-security -Wmissing-field-initializers -Wfatal-errors -Wno-missing-field-initializers -std=c++11 -lm

Thanks for any hint.

Running the Simulink Model

Hi,

I am not able to run the simulink model. I tried executing RunNavFilterTestHarness21 file in rapid accelerator model, but it throws "Unable to build a standalone executable to simulate the model 'NavFilterTestHarness21' in rapid accelerator mode.". I had the same problem running with RunNavFilterTestHarness24. Could you help me out ?

Thanks.

Regards,
Vinod

Insufficient memory

Hi everyone,

I tried to run the codes following the given instructions on Matlab R2013b. But it said that

>>  RunNavFilterTestHarness24
Making simulation target "NavFilterTestHarness24_sfun", ... 


 
J:\dronematlab\InertialNav\models\slprj\_sfprj\NavFilterTestHarness24\_self\sfun\src>"J:\dronematlab\R2013b\R2013b\sys\lcc64\lcc64\bin\lccmake.exe" -f NavFilterTestHarness24_sfun.lmk  
Error c3_navfiltertestharness24.c: 35132  insufficient memory 
    "J:\dronematlab\R2013b\R2013b\sys\lcc64\lcc64\bin\lcc64.exe" -dll -noregistrylookup  -c -Zp8 -DLCC_WIN64 -DMATLAB_MEX_FILE -DMX_COMPAT_32 -nodeclspec  -I"J:\dronematlab\R2013b\R2013b\sys\lcc64\lcc64\include64"  -I"J:\dronematlab\R2013b\R2013b\extern\include"  -I"J:\dronematlab\R2013b\R2013b\simulink\include"  -I"J:\dronematlab\R2013b\R2013b\rtw\c\src"  -I"J:\dronematlab\R2013b\R2013b\stateflow\c\mex\include"  -I"J:\dronematlab\R2013b\R2013b\stateflow\c\debugger\include"  -I"J:\dronematlab\R2013b\R2013b\toolbox\stateflow\src\sf_runtime\export\include\sf_runtime"  "NavFilterTestHarness24_sfun.c" 
    "J:\dronematlab\R2013b\R2013b\sys\lcc64\lcc64\bin\lcc64.exe" -dll -noregistrylookup  -c -Zp8 -DLCC_WIN64 -DMATLAB_MEX_FILE -DMX_COMPAT_32 -nodeclspec  -I"J:\dronematlab\R2013b\R2013b\sys\lcc64\lcc64\include64"  -I"J:\dronematlab\R2013b\R2013b\extern\include"  -I"J:\dronematlab\R2013b\R2013b\simulink\include"  -I"J:\dronematlab\R2013b\R2013b\rtw\c\src"  -I"J:\dronematlab\R2013b\R2013b\stateflow\c\mex\include"  -I"J:\dronematlab\R2013b\R2013b\stateflow\c\debugger\include"  -I"J:\dronematlab\R2013b\R2013b\toolbox\stateflow\src\sf_runtime\export\include\sf_runtime"  "c1_NavFilterTestHarness24.c" 
    "J:\dronematlab\R2013b\R2013b\sys\lcc64\lcc64\bin\lcc64.exe" -dll -noregistrylookup  -c -Zp8 -DLCC_WIN64 -DMATLAB_MEX_FILE -DMX_COMPAT_32 -nodeclspec  -I"J:\dronematlab\R2013b\R2013b\sys\lcc64\lcc64\include64"  -I"J:\dronematlab\R2013b\R2013b\extern\include"  -I"J:\dronematlab\R2013b\R2013b\simulink\include"  -I"J:\dronematlab\R2013b\R2013b\rtw\c\src"  -I"J:\dronematlab\R2013b\R2013b\stateflow\c\mex\include"  -I"J:\dronematlab\R2013b\R2013b\stateflow\c\debugger\include"  -I"J:\dronematlab\R2013b\R2013b\toolbox\stateflow\src\sf_runtime\export\include\sf_runtime"  "c2_NavFilterTestHarness24.c" 
    "J:\dronematlab\R2013b\R2013b\sys\lcc64\lcc64\bin\lcc64.exe" -dll -noregistrylookup  -c -Zp8 -DLCC_WIN64 -DMATLAB_MEX_FILE -DMX_COMPAT_32 -nodeclspec  -I"J:\dronematlab\R2013b\R2013b\sys\lcc64\lcc64\include64"  -I"J:\dronematlab\R2013b\R2013b\extern\include"  -I"J:\dronematlab\R2013b\R2013b\simulink\include"  -I"J:\dronematlab\R2013b\R2013b\rtw\c\src"  -I"J:\dronematlab\R2013b\R2013b\stateflow\c\mex\include"  -I"J:\dronematlab\R2013b\R2013b\stateflow\c\debugger\include"  -I"J:\dronematlab\R2013b\R2013b\toolbox\stateflow\src\sf_runtime\export\include\sf_runtime"  "c3_NavFilterTestHarness24.c" 
J:\dronematlab\R2013b\R2013b\sys\lcc64\lcc64\bin\lccmake.exe: Error code 1 



Error using RunNavFilterTestHarness24 (line 4)
(SLSF Diagnostic)

And I also tried to run it on R2013a. It said:

>> RunNavFilterTestHarness24
Error using RunNavFilterTestHarness24 (line 4)
**** Error ****

Does anyone know how to solve this? Thanks!

GenerateEquations22states error in 2018b

>> GenerateEquations22states
Error using sym/subexpr (line 27)
Second argument must be a variable or a char array or
string naming a variable.

Error in OptimiseAlgebra (line 17)
    [SymExpOut,SubExpOut]=subexpr(SymExpIn,SubExpIn);

Error in GenerateEquations22states (line 133)
[F,SF]=OptimiseAlgebra(F,'SF');

What representation is expected for gyro measurements?

AttPosEKF::UpdateStrapdownEquationsNED constructs deltaQuat by assuming that correctedDelAng is in an axis-angle representation, where the magnitude of the Vector3f is the angle, and the direction of the vector is the axis of rotation.

I expected that any gyro measurement would instead be represented as Euler angles, that is, three orthogonal measurements of rotation rate. That's not algebraically equivalent to the axis-angle representation, as far as I can tell, except when any two of the components are 0.

What's the justification for this approximation? Is it expected that UAVs will usually only rotate around one body-frame axis at a time? Or perhaps you determined experimentally that the difference is small in practice?

The computation of summedDelAng has a related issue: summing rotation rates component-wise is only valid as long as all rotation is around the same gyro axis. Again, what's the intuition behind this?

calc_truth_deriv

What's the meaning of line 46 in the calc_truth_deriv.m?
acc_deriv = 10*w*sin(wt) - 10*Acc - 10*Vel;
Is this a formula like F=macc+cvel+k*pos?But why it is 10?@priseborough
Thanks.

how to understand the code

Hi, priseborough! I am really grateful that you released the code about inertial navigation. But I can't understand the details. Could you please add some doc about the theory! Any help would be really appreciated!

NavFilterTestHarness24': Undefined function or variable 'startTime'

Hi
I am getting the error below when trying to run NavFilterTestHarness24.mdl in Matlab 2015a. I have followed the steps in the readme. Any clue why its happening ?

Thank you

Error evaluating expression 'startTime' for 'StartTime' specified in the Configuration Parameters dialog for block diagram 'NavFilterTestHarness24': Undefined function or variable 'startTime'

NavFilterTestHarness24:

Hi Dear
I have followed the steps in the readme , but get the error below when running NavFilterTestHarness24.mdl in Matlab 2012a. How can I slove it?

Thank you

Unable to build a standalone executable to simulate the model 'NavFilterTestHarness24' in rapid accelerator mode. This could be because the model
has blocks that either do not support code generation or generate code only for specific hardware targets. Models with such blocks cannot be
simulated in rapid accelerator mode.
The rtwbuild error (message and identifier shown below) may identify these blocks.
message: Error(s) encountered while building model "NavFilterTestHarness24"
identifier: 'RTW:buildProcess:fatalBuildError'

Caused by:
Error using RunNavFilterTestHarness24 (line 4)
Error(s) encountered while building model "NavFilterTestHarness24"

attitude plot data source

The python plot show two euler angle sources : onboard and estimate. The estimate data comes from EKF, onboard data comes from ATT.txt, I thought the euler angle onboard data should come from the primary calculation of _ekf->UpdateStrapdownEquationsNED() function, but this function WriteFilterOutput() shows onboard data is from ATT.txt, so I am confused, what sensor data could this ATT.txt represents?

EFK ? TestHarness24

Hello
I would like to use the EKF for offline wind calculation. I already use the pixhawk on line with good result compare to meteorology system. In the present case i want to use it but the calcul is not good. I try to use NavfilterTestHarnes24. I use LoadMissionPlannerPlaneFlashData model to create my data vector
but what is :
msecVelDelay = 230;
msecPosDelay = 210;
msecHgtDelay = 350;
msecMagDelay = 30;
msecTasDelay = 210;

the vector time is in milliseconde or seconde

at the end i just want wind calculation form EKF2

Thanks for your help
Patrice

Using This Code for AHRS or VRU applications

Hi
Thanks for your impressive code.
I want to use your code as a VRU estimation code. I mean, just estimating roll and pith angles.
So in this application, the only observations are delta_angles and aelta_velocities. but in your code both of them are variable parameters and not observations. So i don't know how to run correction step in kalman filter. because there is no observation.
Can anyone tell me that, how to use this algorithm in such situations?

Incremental initialization

It would be beneficial if the C code would support a partial initialization without GPS lock to provide a somewhat accurate attitude and altitude estimate in absence of GPS. This will of course make it weak w.r.t. centrifugal forces and other accelerations during that phase, but will give a better pilot experience when handling the craft and would allow indoor tests of the attitude control response.

About understanding ekf2 code of px4

Hello,
I am a graduate student. I have just started working on ekf2 (https://github.com/PX4/Firmware/blob/master/src/modules/ekf2/ekf2_main.cpp) code of px4 firmware. Kindly can you help me to understand it.
I am currently referring to this site http://ardupilot.org/dev/docs/ekf2-estimation-system.html but not sure that the same code is written in cpp as ekf2 for px4 firmware.
Please reply when you get time if possible provide me literature which explains about ekf2 code.
Thanks in advance.

Mathmatical Equations of InertialNav filter

First of all,
I would like to appreciate your matlab and c++ code of InertialNav.
It helps me a lot.

Is there any Mathmatical Equations of GenerateEquations24states.m file and simulink file?
(Any reference journal, conference paper, or any document that describes the Mathmatical Equations of navigation filter would be the best)

Best Regards

Align time

Hi priseborough
Iam trying to insert my data into your code but it come out an with this error

Error in LoadNavFilterTestData (line 77)
alignTime = min(IMUtime(IMUframe>GPSframe(find(GndSpd >8, 1 ))))- 10 ;

May i know on how to align time gps and imu?Or do you have any idea on this issue?

Thanks for cooperation and time

Idea of Adding a Second Optical Flow Sensor to Recover Roll and Pitch Rates

I am thinking of having two optical flow sensors on a multicopter, with one canted away from nadir into X and the other into Y. Because the coefficients of the two sensors from translation and rotation rates are different, it should be possible to recover roll and pitch rates as well as translation from optical sensors and improve hovering performance. My question is, is this a worthwhile project? I would like to learn more about multicopter navigation and the ArduPilot codebase, and so I’m looking for a project to take on.

Arducopter uses EKF or CF

Dear Concerned,

Many thanks for the complete and complete detailed code.

I have some doubts. The paper "Lim, H.; Park, J.; Lee, D.; Kim, H.J. Build your own quadrotor: Open-source projects on unmanned aerial vehicles. IEEE Robotics & Automation Magazine 2012, 19, 33–45; DOI: 10.1109/MRA.2012.2205629" says that Arducopter is based on Non-linear Complementary Filter, however, your page says it's based on EKF. Could you please clarify/ give more information if from any new version the EKF is been implemented?
Thanks in advance.

Stuck in [PP,SPP]=OptimiseAlgebra(PP,'SPP');

Matlab R2021a, run m files "GenerateNavFilterEquations.m", stuck in [PP,SPP]=OptimiseAlgebra(PP,'SPP'); I have spend about 30hours run this code. The index loop to 7800, and still not finish. I have updated the function "OptimiseAlgebra". Any one known where is wrong?
2023-02-09 22 42 32

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.