#include <stdio.h>


float Q_angle=0.001;//过程噪声的协方差
float Q_gyro=0.003;//0.003 过程噪声的协方差 过程噪声的协方差为一个一行两列矩阵
float R_angle=0.5;// 测量噪声的协方差 既测量偏差
float dt=0.005;//注意：dt的取值为kalman滤波器采样时间;
float P[2][2] = {{ 1, 0 },{ 0, 1 }};
float Angle,Angle_dot;
float Pdot[4] ={0,0,0,0};
const char C_0 = 1;
float Q_bias, Angle_err, PCt_0, PCt_1, E, K_0, K_1, t_0, t_1;

/**************************************************************************
函数功能：简易卡尔曼滤波
入口参数：加速度、角速度
返回  值：无
**************************************************************************/
void Kalman_Filter(float Angle_m,float Gyro_m) //gyro_m:gyro_measure
{
    Angle+=(Gyro_m-Q_bias) * dt;//先验估计误差协方差的微分
    Pdot[0]=Q_angle - P[0][1] - P[1][0];// Pk-先验估计误差协方差的微分
    Pdot[1]=- P[1][1];
    Pdot[2]=- P[1][1];
    Pdot[3]=Q_gyro;
    
    P[0][0] += Pdot[0] * dt;// Pk- 先验估计误差协方差微分的积分
    P[0][1] += Pdot[1] * dt;// = 先验估计误差协方差
    P[1][0] += Pdot[2] * dt;
    P[1][1] += Pdot[3] * dt;
    
    Angle_err = Angle_m - Angle;//zk-先验估计
    
    PCt_0 = C_0 * P[0][0];
    PCt_1 = C_0 * P[1][0];
    
    E = R_angle + C_0 * PCt_0;
    
    K_0 = PCt_0 / E;//Kk
    K_1 = PCt_1 / E;
    
    t_0 = PCt_0;
    t_1 = C_0 * P[0][1];
    
    P[0][0] -= K_0 * t_0;//后验估计误差协方差
    P[0][1] -= K_0 * t_1;
    P[1][0] -= K_1 * t_0;
    P[1][1] -= K_1 * t_1;
    
    Angle += K_0 * Angle_err;//后验估计
    Q_bias += K_1 * Angle_err;//后验估计
    Angle_dot = Gyro_m-Q_bias;//输出值(后验估计)的微分 = 角速度
}