没有一定的数学和控制理论基础 对卡尔曼滤波的理解还是有一定困难的;
X(k)=AX(k-1)+BU(k)+W(k) 状态反馈矩阵
Z(k)=HX(k)+V(k) 观测阵 现控知识
w v 为过程噪声和测量噪声
卡尔曼滤波主要五个核心公式
(1) X(k|k-1)=AX(k-1|k-1)+BU(k)
(2)P(k|k-1)=AP(k-1|k-1)A’+Q
(3)X(k|k)=X(k|k-1)+Kg(k)(Z(k)-HX(k|k-1))
(4)Kg(k)=P(k|k-1)H’/(HP(k|k-1)H’+R)
(5)P(k|k)=(I-Kg(k)H)P(k|k-1)
其中 X(k|k-1)为上一期状态预测结果 X(k-1|k-1)上一次状态最优结果 U(k)为现在状态的控制量
P(k|k-1)是X(k|k-1)对应的协方差 Q R为过 观噪声协方差
float kalmanFilter(float inData)
{
static float lastData=0;
static float p, q, r, kGain=0;//p不wei0 r q 自己估
p = p+q;
kGain = p/(p+r); //卡尔曼增益
inData = lastData+(kGain*(inData-lastData)); //最优解
p = (1-kGain)*p; //协方差
lastData = inData;
return inData;
}
N=500;
w(1)=0;%过程噪声
w=randn(1,N); %w为1行N列的矩阵,且服从N(1,N)的正态分布
x(1)=0;
a=1; %房间温度我们认为不变 自己预测
for k=2:N;
x(k)=a*x(k-1)+w(k-1);%系统趋势值 无控制量 白噪声为wk
end
V=randn(1,N);
q1=std(V);
Rvv=q1.^2%系统测量协方差 r
q2=std(x);
Rxx=q2.^2;
q3=std(w);
Rww=q3.^2 %系统过程协方差 q
c=0.2;
Y=c*x+V; %测量值 实际上的测量值会绕着真值波动
p(1)=10;
s(1)=1;
for t=2:N;
p1(t)=a.^2*p(t-1)+Rww;%用上次协方差计算
b(t)=c*p1(t)/(c.^2*p1(t)+Rvv);% 卡尔曼增益
s(t)=a*s(t-1)+b(t)*(Y(t)-a*c*s(t-1));%此次最优化估算结果
p(t)=p1(t)-c*b(t)*p1(t);%更新本次协方差
end
t=1:N;
plot(t,s,'r',t,Y,'g',t,x,'k');
legend('Kalman滤波结果','观测值','系统值');
r参数调整滤波后的曲线与实测曲线的相近程度,r越小越接近。
q参数调滤波后的曲线平滑程度,q越小越平滑。
以下三类
r=RVV;
r=RVV*10;
r=RVV/10
github 上的用卡尔曼姿态解算
#include "Wire.h"
#include "I2Cdev.h"
#include "MPU6050.h"
MPU6050 accelgyro;
unsigned long now, lastTime = 0;
float dt; //微分时间
int16_t ax, ay, az, gx, gy, gz; //加速度计陀螺仪原始数据
float aax=0, aay=0,aaz=0, agx=0, agy=0, agz=0; //角度变量
long axo = 0, ayo = 0, azo = 0; //加速度计偏移量
long gxo = 0, gyo = 0, gzo = 0; //陀螺仪偏移量
float pi = 3.1415926;
float AcceRatio = 16384.0; //加速度计比例系数
float GyroRatio = 131.0; //陀螺仪比例系数
uint8_t n_sample = 8; //加速度计滤波算法采样个数
float aaxs[8] = {0}, aays[8] = {0}, aazs[8] = {0}; //x,y轴采样队列
long aax_sum, aay_sum,aaz_sum; //x,y轴采样和
float a_x[10]={0}, a_y[10]={0},a_z[10]={0} ,g_x[10]={0} ,g_y[10]={0},g_z[10]={0}; //加速度计协方差计算队列
float Px=1, Rx, Kx, Sx, Vx, Qx; //x轴卡尔曼变量
float Py=1, Ry, Ky, Sy, Vy, Qy; //y轴卡尔曼变量
float Pz=1, Rz, Kz, Sz, Vz, Qz; //z轴卡尔曼变量
void setup()
{
Wire.begin();
Serial.begin(115200);
accelgyro.initialize(); //初始化
unsigned short times = 200; //采样次数
for(int i=0;i<times;i++)
{
accelgyro.getMotion6(&ax, &ay, &az, &gx, &gy, &gz); //读取六轴原始数值
axo += ax; ayo += ay; azo += az; //采样和
gxo += gx; gyo += gy; gzo += gz;
}
axo /= times; ayo /= times; azo /= times; //计算加速度计偏移
gxo /= times; gyo /= times; gzo /= times; //计算陀螺仪偏移
}
void loop()
{
unsigned long now = millis(); //当前时间(ms)
dt = (now - lastTime) / 1000.0; //微分时间(s)
lastTime = now; //上一次采样时间(ms)
accelgyro.getMotion6(&ax, &ay, &az, &gx, &gy, &gz); //读取六轴原始数值
float accx = ax / AcceRatio; //x轴加速度
float accy = ay / AcceRatio; //y轴加速度
float accz = az / AcceRatio; //z轴加速度
aax = atan(accy / accz) * (-180) / pi; //y轴对于z轴的夹角
aay = atan(accx / accz) * 180 / pi; //x轴对于z轴的夹角
aaz = atan(accz / accy) * 180 / pi; //z轴对于y轴的夹角
aax_sum = 0; // 对于加速度计原始数据的滑动加权滤波算法
aay_sum = 0;
aaz_sum = 0;
for(int i=1;i<n_sample;i++)
{
aaxs[i-1] = aaxs[i];
aax_sum += aaxs[i] * i;
aays[i-1] = aays[i];
aay_sum += aays[i] * i;
aazs[i-1] = aazs[i];
aaz_sum += aazs[i] * i;
}
aaxs[n_sample-1] = aax;
aax_sum += aax * n_sample;
aax = (aax_sum / (11*n_sample/2.0)) * 9 / 7.0; //角度调幅至0-90°
aays[n_sample-1] = aay; //此处应用实验法取得合适的系数
aay_sum += aay * n_sample; //本例系数为9/7
aay = (aay_sum / (11*n_sample/2.0)) * 9 / 7.0;
aazs[n_sample-1] = aaz;
aaz_sum += aaz * n_sample;
aaz = (aaz_sum / (11*n_sample/2.0)) * 9 / 7.0;
float gyrox = - (gx-gxo) / GyroRatio * dt; //x轴角速度
float gyroy = - (gy-gyo) / GyroRatio * dt; //y轴角速度
float gyroz = - (gz-gzo) / GyroRatio * dt; //z轴角速度
agx += gyrox; //x轴角速度积分
agy += gyroy; //x轴角速度积分
agz += gyroz;
/* kalman start */
Sx = 0; Rx = 0;
Sy = 0; Ry = 0;
Sz = 0; Rz = 0;
for(int i=1;i<10;i++)
{ //测量值平均值运算
a_x[i-1] = a_x[i]; //即加速度平均值
Sx += a_x[i];
a_y[i-1] = a_y[i];
Sy += a_y[i];
a_z[i-1] = a_z[i];
Sz += a_z[i];
}
a_x[9] = aax;
Sx += aax;
Sx /= 10; //x轴加速度平均值
a_y[9] = aay;
Sy += aay;
Sy /= 10; //y轴加速度平均值
a_z[9] = aaz;
Sz += aaz;
Sz /= 10;
for(int i=0;i<10;i++)
{
Rx += sq(a_x[i] - Sx);
Ry += sq(a_y[i] - Sy);
Rz += sq(a_z[i] - Sz);
}
Rx = Rx / 9; //得到方差
Ry = Ry / 9;
Rz = Rz / 9;
Px = Px + 0.0025; // 0.0025在下面有说明...
Kx = Px / (Px + Rx); //计算卡尔曼增益
agx = agx + Kx * (aax - agx); //陀螺仪角度与加速度计速度叠加
Px = (1 - Kx) * Px; //更新p值
Py = Py + 0.0025;
Ky = Py / (Py + Ry);
agy = agy + Ky * (aay - agy);
Py = (1 - Ky) * Py;
Pz = Pz + 0.0025;
Kz = Pz / (Pz + Rz);
agz = agz + Kz * (aaz - agz);
Pz = (1 - Kz) * Pz;
/* kalman end */
Serial.print(agx);Serial.print(",");
Serial.print(agy);Serial.print(",");
Serial.print(agz);Serial.println();