5iMX宗旨:分享遥控模型兴趣爱好

5iMX.com 我爱模型 玩家论坛 ——专业遥控模型和无人机玩家论坛(玩模型就上我爱模型,创始于2003年)
楼主: dongfang
打印 上一主题 下一主题

开源!我的捷联系统的四元数法姿态算法和程序

[复制链接]
61
发表于 2012-1-28 02:36 | 只看该作者
最近有时间也玩一下姿态解算。楼主看下我的理解正确么
Quaternion nowQ;                    //获得原始姿态的四元数  [a]
void loop() {                       //进入循环
  int intx, inty,intz;
  float  pitch,roll,yaw;

  gyro.ReadGyroOutCalibrated_Radian(&pitch, &roll, &yaw);   //读取陀螺输出。并且校准(弧度)

  EulerAngle dt;                                     //获得dt。就是采样刷新time
  dt.fRoll=roll;                                     //计算roll的变化量 (gyrox*dt)
  dt.fPitch=pitch;                                     //计算pitch的变化量(gyroy*dt)
  dt.fYaw=-yaw;                                        //计算yaw的变化量(有反向,和安装方向有关)(gyroz*dt)

  Quaternion dQ;                                        //计算当前姿态的四元数
  dQ.FromEulerAngle(dt);                                //把dt的欧拉转换为四元数   
  nowQ=nowQ.Multiply(dQ);                                //做四元数的乘法运算,得到新的四元数(nowQ=nowQ*dQ(是dt转换得到))[a=a*b]


  count++;                                              //计数器++
  if (count>1000){                                       //如果count大于1000则将nowQ四元转换为欧拉角
    EulerAngle nowG=nowQ.ToEulerAngle();

    Serial.print(nowG.fRoll/3.1415926535*180,11);//横滚
    Serial.print(",");                                            //将弧度转换为度,串口输出
    Serial.print(nowG.fPitch/3.1415926535*180,11);//俯仰
    Serial.print(",");
    Serial.print(nowG.fYaw/3.1415926535*180,11);//偏航
    Serial.print(",");
    Serial.print(nowQ.getS(),11);//偏航
    Serial.println();
    count=0;
  }
}
有点迷惑的是这里:
Quaternion dQ;                                        //计算当前姿态的四元数
  dQ.FromEulerAngle(dt);                                //把dt的欧拉转换为四元数   
既然后面做四元乘法用的是a=a*b。那Quaternion dQ;       //计算当前姿态的四元数
有什么用呢。
后面的计算就是将上一次的nowq作为原始姿态和更新的dt做乘法。这样就用三轴陀螺计算出姿态,如果要做acc的融合的话就另说了。
找到了一个算法:
//=====================================================================================================
// IMU.c
// S.O.H. Madgwick
// 25th September 2010
//=====================================================================================================
// Description:
//
// Quaternion implementation of the 'DCM filter' [Mayhony et al].
//
// User must define 'halfT' as the (sample period / 2), and the filter gains 'Kp' and 'Ki'.
//
// Global variables 'q0', 'q1', 'q2', 'q3' are the quaternion elements representing the estimated
// orientation.  See my report for an overview of the use of quaternions in this application.
//
// User must call 'IMUupdate()' every sample period and parse calibrated gyroscope ('gx', 'gy', 'gz')
// and accelerometer ('ax', 'ay', 'ay') data.  Gyroscope units are radians/second, accelerometer
// units are irrelevant as the vector is normalised.
//
//=====================================================================================================

//----------------------------------------------------------------------------------------------------
// Header files

#include "IMU.h"
#include <math.h>

//----------------------------------------------------------------------------------------------------
// Definitions

#define Kp 2.0f                        // proportional gain governs rate of convergence to accelerometer/magnetometer
#define Ki 0.005f                // integral gain governs rate of convergence of gyroscope biases
#define halfT 0.5f                // half the sample period

//---------------------------------------------------------------------------------------------------
// Variable definitions

float q0 = 1, q1 = 0, q2 = 0, q3 = 0;        // quaternion elements representing the estimated orientation
float exInt = 0, eyInt = 0, ezInt = 0;        // scaled integral error

//====================================================================================================
// Function
//====================================================================================================

void IMUupdate(float gx, float gy, float gz, float ax, float ay, float az) {
        float norm;
        float vx, vy, vz;
        float ex, ey, ez;         
       
        // normalise the measurements
        norm = sqrt(ax*ax + ay*ay + az*az);      
        ax = ax / norm;
        ay = ay / norm;
        az = az / norm;      
       
        // estimated direction of gravity
        vx = 2*(q1*q3 - q0*q2);
        vy = 2*(q0*q1 + q2*q3);
        vz = q0*q0 - q1*q1 - q2*q2 + q3*q3;
       
        // error is sum of cross product between reference direction of field and direction measured by sensor
        ex = (ay*vz - az*vy);
        ey = (az*vx - ax*vz);
        ez = (ax*vy - ay*vx);
       
        // integral error scaled integral gain
        exInt = exInt + ex*Ki;
        eyInt = eyInt + ey*Ki;
        ezInt = ezInt + ez*Ki;
       
        // adjusted gyroscope measurements
        gx = gx + Kp*ex + exInt;
        gy = gy + Kp*ey + eyInt;
        gz = gz + Kp*ez + ezInt;
       
        // integrate quaternion rate and normalise
        q0 = q0 + (-q1*gx - q2*gy - q3*gz)*halfT;
        q1 = q1 + (q0*gx + q2*gz - q3*gy)*halfT;
        q2 = q2 + (q0*gy - q1*gz + q3*gx)*halfT;
        q3 = q3 + (q0*gz + q1*gy - q2*gx)*halfT;  
       
        // normalise quaternion
        norm = sqrt(q0*q0 + q1*q1 + q2*q2 + q3*q3);
        q0 = q0 / norm;
        q1 = q1 / norm;
        q2 = q2 / norm;
        q3 = q3 / norm;
}

//====================================================================================================
// END OF CODE
//====================================================================================================
这个四元貌似有acc的融合功能???请指教,谢谢

[ 本帖最后由 峰回路转 于 2012-1-28 02:38 编辑 ]

欢迎继续阅读楼主其他信息

62
发表于 2012-1-28 07:38 | 只看该作者
楼主没有考虑线加速度对角速度的影响吧?
63
 楼主| 发表于 2012-1-28 10:51 | 只看该作者
原帖由 峰回路转 于 2012-1-28 02:36 发表
最近有时间也玩一下姿态解算。楼主看下我的理解正确么
Quaternion nowQ;                    //获得原始姿态的四元数  [a]
void loop() {                       //进入循环
  int intx, inty,intz;
  float  pi ...


你好,Quaternion dQ,也就是a=a*b中的那个b,是自上一个姿态计算后,发生旋转(增量delta)的角度,可以从角速度读数乘以增量时间得到dt(这里的dt不是时间,是增量角度),再转换成四元数得到。

你的那段带ACC融合的代码看上去不错,但肯定不是标准的四元数算法,因为没看到三角函数,有机会我试一试。
64
 楼主| 发表于 2012-1-28 10:52 | 只看该作者
原帖由 oppop 于 2012-1-28 07:38 发表
楼主没有考虑线加速度对角速度的影响吧?

不是很明白线加速度对角速度的影响指的是什么,我的代码没有集成acc是事实。
65
发表于 2012-1-28 11:51 | 只看该作者
:em15:
66
 楼主| 发表于 2012-1-28 12:37 | 只看该作者
找到了58楼代码的出处了,代码还在更新中,现在已经支持9dof了,

作者的主页在:
http://code.google.com/p/imumargalgorithm30042010sohm/
http://www.x-io.co.uk/node/8#open_source_imu_and_ahrs_algorithms

c代码在此:http://www.x-io.co.uk/res/sw/madgwick_algorithm_c.zip

姿态解算这么高深的问题,有看得懂的代码,感觉开源真是太了不起了。
67
发表于 2012-1-28 19:23 | 只看该作者
原帖由 dongfang 于 2012-1-28 10:51 发表


你好,Quaternion dQ,也就是a=a*b中的那个b,是自上一个姿态计算后,发生旋转(增量delta)的角度,可以从角速度读数乘以增量时间得到dt(这里的dt不是时间,是增量角度),再转换成四元数得到。

你的那段带AC ...

呵呵。那我就明白了。关于四元的公式用这个就可以了吧:void Quaternion::FromEulerAngle(const EulerAngle &ea)
{
        float fCosHRoll = cos(ea.fRoll * .5f);
        float fSinHRoll = sin(ea.fRoll * .5f);
        float fCosHPitch = cos(ea.fPitch * .5f);
        float fSinHPitch = sin(ea.fPitch * .5f);
        float fCosHYaw = cos(ea.fYaw * .5f);
        float fSinHYaw = sin(ea.fYaw * .5f);

        w = fCosHRoll * fCosHPitch * fCosHYaw + fSinHRoll * fSinHPitch * fSinHYaw;
        x = fCosHRoll * fSinHPitch * fCosHYaw + fSinHRoll * fCosHPitch * fSinHYaw;
        y = fCosHRoll * fCosHPitch * fSinHYaw - fSinHRoll * fSinHPitch * fCosHYaw;
        z = fSinHRoll * fCosHPitch * fCosHYaw - fCosHRoll * fSinHPitch * fSinHYaw;
}
这个函数输入gyrox y z输出wxyz也就是q0 q1 q2 q3 对四元的概念还是有些模糊,看来得好好看看角度,弧度,欧拉,四元之间的关系。
68
发表于 2012-1-28 19:25 | 只看该作者
原帖由 dongfang 于 2012-1-28 12:37 发表
找到了58楼代码的出处了,代码还在更新中,现在已经支持9dof了,

作者的主页在:
http://code.google.com/p/imumargalgorithm30042010sohm/
http://www.x-io.co.uk/node/8#open_source_imu_and_ahrs_algorithms ...

这个代码楼主有没有测试过,他融合了三轴加速度和三轴磁场,期待详解。。。。
69
 楼主| 发表于 2012-1-28 22:34 | 只看该作者
奇怪,我刚看到了h.e.u.y.c.k的回复,刚想说两句,突然不见了,是被我自己删了么?
70
发表于 2012-1-28 23:11 | 只看该作者
顶起 楼主是偶像
71
发表于 2012-1-29 09:13 | 只看该作者
做一个46656000行的查表法如何?:em15:
72
发表于 2012-1-29 09:22 | 只看该作者

标题

原帖由 dongfang 于 2012-1-28 22:34 发表
奇怪,我刚看到了h.e.u.y.c.k的回复,刚想说两句,突然不见了,是被我自己删了么?

不想发帖了,免得有人jjww。最近讲啥都有人跟。
你看过回帖了解就好了
73
发表于 2012-1-29 18:45 | 只看该作者
原帖由 dongfang 于 2012-1-28 10:52 发表

不是很明白线加速度对角速度的影响指的是什么,我的代码没有集成acc是事实。


由于MEMS陀螺的固有特性,每个坐标上的角速度输出值会受到线加速度的影响,越便宜的陀螺仪这种影响越明显,即,Linear Acceleration Effect on Bias。
换句话说,低端的陀螺仪对陀螺本身的振动非常敏感。高端的陀螺仪贵的原因一方面是固有噪声低和零漂稳定性低,另一方面则是对加速度的敏感度低。因为高端的陀螺仪每个轴向上都会做两对相同的结构来构成差分的形式,来消除振动带来的影响。这样做的话,成本自然就高了。
一般的、起步比较晚的MEMS陀螺仪生产厂商,一般是不会提Linear Acceleration Effect on Bias这个指标的,因为他们这个指标都很差。但是在实际运用中,振动是不可避免的,振动带来的误差往往占很大比例。这就是为什么有些飞控的内部的IMU要加减振垫子的原因。
74
 楼主| 发表于 2012-1-29 20:58 | 只看该作者
原帖由 oppop 于 2012-1-29 18:45 发表


由于MEMS陀螺的固有特性,每个坐标上的角速度输出值会受到线加速度的影响,越便宜的陀螺仪这种影响越明显,即,Linear Acceleration Effect on Bias。
换句话说,低端的陀螺仪对陀螺本身的振动非常敏感。高端的 ...

首次听说,不知道实际影响有多大,好像感觉不出来么。
75
 楼主| 发表于 2012-1-29 21:34 | 只看该作者
原帖由 (涉嫌广告已屏蔽) 于 2012-1-29 09:22 发表

不想发帖了,免得有人jjww。最近讲啥都有人跟。
你看过回帖了解就好了

唉,人怕出名猪怕胖,就讨论些技术算了,别扯那些是非了。算法上荒废很久了,已经遇到了AVR的运算极限,STM32也没搞定(本身就是新学,IIC没搞定就已经损失了一个测试板)。

最近各种分心的事很多,加上在研究四轴,感觉精力不济。
你要是有时间的话,我建议你的飞控搞一个飞行测试模式,就是依次测试:
测试SETUP: 测试在100米左右安全微风空域进行
TEST 1: 测试油门到达指定速度SPEED的理想PID参数,参数能和SPEED成函数关系最好。我在xplan的测试中已经发现I*I积分是和SPEED直接关系的。
TEST 2: 测试各SPEED下保持平衡飞行的ROLL的理想PID参数。
TEST 3: 测试各SPEED下保持高度、平衡飞行的PITCH的理想PID参数。
TEST 4: 测试各SPEED下保持高度、平衡飞行,仅使用方向舵转弯,各个转弯力度飞行的方向舵的理想PID参数。
TEST 5: 测试各SPEED下保持高度、平衡飞行,组合使用方向舵和副翼转弯,各个转弯力度飞行的方向舵、副翼、升降舵的理想PID参数。
测试错误控制:一旦飞机失控,高度降低20米以上,立刻使用FAILSAFE参数,救回。

我在xplan上完成过一部分测试,也发现了一些有趣的规律,只是程序写的很烂,每次测试都是直接修改了前一次的测试代码,所以没留下来什么能共享的。
76
发表于 2012-1-29 22:56 | 只看该作者
呵呵。各位都是高手,最近假期较长,想好好研究下,平台用的是fy20的平衡仪,上面有stm32f103brt+adxl335+ly530+lpr530a。需要的东西都全了。。把原来的mcu拆了保存好。万一以后想玩平衡仪了还可以装回去,装上一片新的。加装jtag(swd)接口。这个本身串口是引出的。加装hp03气压传感器。stm32没有eeprom,装了一片24c02.试过吧kk_c移植到这个平台上,测试可以工作。但是加速度计等都没用上,这次乘着长假玩玩。。。

[ 本帖最后由 峰回路转 于 2012-1-29 23:03 编辑 ]

1.jpg (45.8 KB, 下载次数: 40)

1.jpg

2.jpg (63.21 KB, 下载次数: 38)

2.jpg
77
发表于 2012-1-29 23:47 | 只看该作者
不懂:em22: 帮顶
78
发表于 2012-1-30 02:35 | 只看该作者
关于欧拉角,gyro输出值积分得到的就是欧拉角么?另外计算时输入的gyro的单位是mv/度/s ? acc的是mv/g
也就是说用gyro积分计算角度偏移量 gyrox_dt=gyrox*dt
用acc计算当前欧拉角 Roll= atan2(ax,az);
                Pitch= atan2(-ay,az);

[ 本帖最后由 峰回路转 于 2012-1-30 03:45 编辑 ]
79
 楼主| 发表于 2012-1-30 22:41 | 只看该作者
原帖由 峰回路转 于 2012-1-30 02:35 发表
关于欧拉角,gyro输出值积分得到的就是欧拉角么?另外计算时输入的gyro的单位是mv/度/s ? acc的是mv/g
也就是说用gyro积分计算角度偏移量 gyrox_dt=gyrox*dt
用acc计算当前欧拉角 Roll= atan2(ax,az);
                Pitch=  ...

gyro输出值*间隔时间=欧拉角。gyro的单位要查芯片手册,数字芯片不用mv。
80
发表于 2012-1-30 23:04 | 只看该作者
原帖由 dongfang 于 2012-1-29 21:34 发表

唉,人怕出名猪怕胖,就讨论些技术算了,别扯那些是非了。算法上荒废很久了,已经遇到了AVR的运算极限,STM32也没搞定(本身就是新学,IIC没搞定就已经损失了一个测试板)。

最近各种分心的事很多,加上在研究四 ...

最近研究XPLANE、APM SIM还有MAVLINK,得把空速管什么的数据都弄出来才好仿真。
xplane可是个好东西啊,可惜只有一架PT60。。。找不到其他模型了,不知道你用什么模型模拟的?
APM功能很强大,正在研究,不过是C++写的,看得我一头雾水。最近又有新版APM2出来,好像只卖200刀,也可以研究一下。
MAVLINK也是个好东西,支持MAVLINK的GCS一大堆,而且有免费库可以用,还是C代码的。
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

关闭

【站内推荐】上一条 /1 下一条

快速回复 返回顶部 返回列表