File content as of revision 3:5358a691a100:

#include "mbed.h"
#include "MPU6050.h"
#include "HMC5883L.h"
#include "LPS25H.h"
#include "GMS6_CR6.h"
#include "ErrorLogger.h"
#include "Vector.h"
#include "Matrix.h"
#include "Vector_Matrix_operator.h"
#include "myConstants.h"
#include "Log.h"

/********** private define    **********/
#define TRUE    1
#define FALSE   0
/********** private macro     **********/

/********** private typedef   **********/

/********** private variables **********/
I2C         i2c(PB_9, PB_8);      // I2Cポート
MPU6050     mpu(&i2c);                  // 加速度・角速度センサ
HMC5883L    hmc(&i2c);                  // 地磁気センサ
LPS25H      lps(&i2c);                  // 気圧センサ
Serial      gps(PA_11, PA_12);                // GPS通信用シリアルポート
Serial      pc(SERIAL_TX, SERIAL_RX);   // PC通信用シリアルポート
GMS6_CR6    gms(&gps, &pc);             // GPS
Ticker      INT_timer;                      // 割り込みタイマー
//Log         logger(PA_9, PA_10, PB_5, PB_4, PB_3, PB_10, "sd");    // ロガー（microSD、XBee）

const float dt = 0.1f;           // 割り込み周期（s）

int lps_cnt = 0;                    // 気圧センサ読み取りカウント
uint8_t INT_flag = TRUE;            // 割り込み可否フラグ
Vector acc(3);                       // 加速度(m/s^2)
Vector gyro(3);                      // 角速度(deg/s)
Vector geomag(3);                    // 地磁気(?)
float press;                        // 気圧(hPa)

Vector g(3);                        // 重力ベクトル
//Vector n(3);                        // 地磁気ベクトル

/* ----- Kalman Filter ----- */
Vector pri_x(6);
Matrix pri_P(6, 6);
Vector post_x(6);
Matrix post_P(6, 6);
Matrix F(6, 6), H(3, 6);
Matrix R(6, 6), Q(3, 3);
Matrix I(6, 6);
Matrix K(6, 3);
Matrix S(3, 3), inv(3, 3);
/* ----- ------------- ----- */

Timer timer;

char data[1024] = {};

/********** private functions **********/
void KalmanInit();
void KalmanUpdate();
void INT_func();                    // 割り込み用関数
void toString(Matrix& m);
void toString(Vector& v);

/********** main function     **********/

int main() {
    
    i2c.frequency(400000);          // I2Cの通信速度を400kHzに設定
    
    if(!mpu.init()) AbortWithMsg("mpu6050 Initialize Error !!");        // mpu6050初期化
    if(!hmc.init()) AbortWithMsg("hmc5883l Initialize Error !!");       // hmc5883l初期化
    
    //char* title = "log data\r\n";                                       // ログのタイトル行
    //if(!logger.initialize_sdlog(title)) return 0;                       // ログファイル初期設定
    
    KalmanInit();
    
    INT_timer.attach(&INT_func, dt);  // 割り込み有効化（Freq = 0.01fなので、10msおきの割り込み）
    
    //重力ベクトルの初期化
    g.SetComp(1, 0.0f);
    g.SetComp(2, 0.0f);
    g.SetComp(3, 1.0f);
    
    /* ---------- ↓↓↓ ここからメインループ ↓↓↓ ---------- */
    while(1) {
        
        // 0.1秒おきにセンサーの出力をpcへ出力
        wait(0.1f);
        
        INT_flag = FALSE;           // 割り込みによる変数書き換えを阻止
        
        // センサ類の全出力値をPCに送信
        /*
        pc.printf("%.3f,", acc.GetComp(1));
        pc.printf("%.3f,", acc.GetComp(2));
        pc.printf("%.3f,", acc.GetComp(3));
        pc.printf("%.3f,", gyro.GetComp(1));
        pc.printf("%.3f,", gyro.GetComp(2));
        pc.printf("%.3f,", gyro.GetComp(3));
        pc.printf("%.3f,", geomag.GetComp(1));
        pc.printf("%.3f,", geomag.GetComp(2));
        pc.printf("%.3f\r\n", geomag.GetComp(3));
        */
        //pc.printf("%.3f\r\n", press);
        
        
        pc.printf("%.3f,", geomag.GetComp(1));
        pc.printf("%.3f,", geomag.GetComp(2));
        pc.printf("%.3f,", geomag.GetComp(3));
        pc.printf("%.3f,", post_x.GetComp(1));
        pc.printf("%.3f,", post_x.GetComp(2));
        pc.printf("%.3f\r\n", post_x.GetComp(3));
        
        /*
        sprintf(data, "%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f\r\n", 
                acc.GetComp(1), acc.GetComp(2), acc.GetComp(3), 
                gyro.GetComp(1), gyro.GetComp(2), gyro.GetComp(3), 
                geomag.GetComp(1), geomag.GetComp(2), geomag.GetComp(3), press);
        logger.puts(data);
        */
        INT_flag = TRUE;            // 割り込み許可
        
    }
    
    /* ---------- ↑↑↑ ここまでメインループ ↑↑↑ ---------- */
}

void KalmanInit() {
    
    // 誤差共分散行列の値を決める（対角成分のみ）
    float alpha_R = 60.0f;
    float alpha_Q = 100.0f;
    R *= alpha_R;
    Q *= alpha_Q;
    
    // 状態方程式のヤコビアンの初期値を代入（時間変化あり）
    float f[36] = {
        1.0f, gyro.GetComp(3)*dt, -gyro.GetComp(2)*dt, 0.0f, 0.0f, 0.0f, 
        -gyro.GetComp(3)*dt, 1.0f, gyro.GetComp(1)*dt, 0.0f, 0.0f, 0.0f, 
        gyro.GetComp(2)*dt, -gyro.GetComp(1)*dt, 1.0f, 0.0f, 0.0f, 0.0f, 
        0.0f, 0.0f, 0.0f, 1.0f, 0.0f, 0.0f, 
        0.0f, 0.0f, 0.0f, 0.0f, 1.0f, 0.0f,  
        0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 1.0f
    };
    
    F.SetComps(f);
    
    // 観測方程式のヤコビアンの値を設定（時間変化無し）
    float h[18] = {
        1.0f, 0.0f, 0.0f, 1.0f, 0.0f, 0.0f,  
        0.0f, 1.0f, 0.0f, 0.0f, 1.0f, 0.0f,  
        0.0f, 0.0f, 1.0f, 0.0f, 0.0f, 1.0f
    };
    
    H.SetComps(h);
}

void KalmanUpdate() {
    // ヤコビアンの更新
    float f[36] = {
        1.0f, gyro.GetComp(3)*dt, -gyro.GetComp(2)*dt, 0.0f, 0.0f, 0.0f, 
        -gyro.GetComp(3)*dt, 1.0f, gyro.GetComp(1)*dt, 0.0f, 0.0f, 0.0f, 
        gyro.GetComp(2)*dt, -gyro.GetComp(1)*dt, 1.0f, 0.0f, 0.0f, 0.0f, 
        0.0f, 0.0f, 0.0f, 1.0f, 0.0f, 0.0f, 
        0.0f, 0.0f, 0.0f, 0.0f, 1.0f, 0.0f,  
        0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 1.0f
    };
    
    F.SetComps(f);
    
    // 事前推定値の更新
    pri_x = F * post_x;
    // 事前誤差分散行列の更新
    pri_P = F * post_P * F.Transpose() + R;
    
    // カルマンゲインの計算
    S = Q + H * pri_P * H.Transpose();
    float det;
    if((det = S.Inverse(inv)) >= 0.0f) {
        pc.printf("E:%.3f\r\n", det);
        return;     // 万が一、逆行列が見つからなかった場合は前回の推定値を保持して終了
    }
    K = pri_P * H.Transpose() * inv;    
    
    // 事後推定値の更新
    post_x = pri_x + K * (geomag - H * pri_x);
    // 事後誤差分散行列の更新
    post_P = (I - K * H) * pri_P;
}

void INT_func() {
    if(INT_flag != FALSE) {
        
        timer.reset();
        timer.start();
        // センサーの値を更新
        mpu.read();
        hmc.read();
        
        for(int i=0; i<3; i++) {
            acc.SetComp(i+1, (float)mpu.data.value.acc[i] * ACC_LSB_TO_G);
            gyro.SetComp(i+1, (float)mpu.data.value.gyro[i] * GYRO_LSB_TO_DEG * DEG_TO_RAD);
            geomag.SetComp(i+1, (float)hmc.data.value[i] * MAG_LSB_TO_GAUSS);
        }
        
        Vector delta_g = Cross(gyro, g);                            // Δg = ω × g
        g = 0.9f * (g - delta_g * dt) + 0.1f * acc.Normalize();   // 相補フィルタ
        g = g.Normalize();
        
        
        KalmanUpdate();
        
        
        // LPS25Hによる気圧の取得は10Hz
        lps_cnt = (lps_cnt+1)%10;
        if(lps_cnt == 0) {
            press = (float)lps.pressure() * PRES_LSB_TO_HPA;
        }
        
        timer.stop();
        //pc.printf("%d(us)\r\n", timer.read_us());
    }
}

void toString(Matrix& m) {
    
    for(int i=0; i<m.GetRow(); i++) {
        for(int j=0; j<m.GetCol(); j++) {
            pc.printf("%.6f\t", m.GetComp(i+1, j+1));
        }
        pc.printf("\r\n");
    }
    
}

void toString(Vector& v) {
    
    for(int i=0; i<v.GetDim(); i++) {
        pc.printf("%.6f\t", v.GetComp(i+1));
    }
    pc.printf("\r\n");
    
}