【51单片机快速入门指南】4.4.3:Madgwick AHRS 九轴姿态融合获取四元数、欧拉角

最新更新时间:2022-07-15来源: csdn关键字:51单片机  AHRS  四元数  欧拉角 手机看文章 扫描二维码
随时随地手机看文章

STC15F2K60S2 22.1184MHz

Keil uVision V5.29.0.0

PK51 Prof.Developers Kit Version:9.60.0.0

上位机:Vofa+ 1.3.10


移植自AHRS —— LOXO,算法作者:SOH Madgwick


传感器的方向

在这里插入图片描述

源码

       所用MCU为STC15F2K60S2 使用内部RC时钟,22.1184MHz


       stdint.h见【51单片机快速入门指南】1:基础知识和工程创建

       软件I2C程序见【51单片机快速入门指南】4: 软件 I2C

       串口部分见【51单片机快速入门指南】3.3:USART 串口通信

       MPU6050驱动程序见【51单片机快速入门指南】4.3: I2C读取MPU6050陀螺仪的原始数据

       HMC5883L/QMC5883L驱动程序见【51单片机快速入门指南】4.4:I2C 读取HMC5883L / QMC5883L 磁力计

       磁力计的椭球拟合校准见【51单片机快速入门指南】4.4.1:python串口接收磁力计数据并进行最小二乘法椭球拟合


       beta要按需调整,我这里取1.0


Madgwick_9.c

//=====================================================================================================

//

// Implementation of Madgwick's IMU and AHRS algorithms.

// See: http://www.x-io.co.uk/node/8#open_source_ahrs_and_imu_algorithms

//

// Date Author          Notes

// 29/09/2011 SOH Madgwick    Initial release

// 02/10/2011 SOH Madgwick Optimised for reduced CPU load

// 19/02/2012 SOH Madgwick Magnetometer measurement is normalised

//

//=====================================================================================================


//---------------------------------------------------------------------------------------------------

// Header files

#include

#include "MPU6050.h"


//---------------------------------------------------------------------------------------------------

// Definitions


#define beta 1.0f // 2 * proportional gain (Kp)


//---------------------------------------------------------------------------------------------------

// Variable definitions


float q0 = 1.0f, q1 = 0.0f, q2 = 0.0f, q3 = 0.0f; // quaternion of sensor frame relative to auxiliary frame

float Pitch = 0.0f, Roll = 0.0f, Yaw = 0.0f;


//====================================================================================================

// Functions


float sampleFreq = 1;

float GYRO_K = 1;


void MPU6050_Madgwick_Init(float loop_ms)

{

sampleFreq = 1000. / loop_ms; //sample frequency in Hz

switch((MPU_Read_Byte(MPU_GYRO_CFG_REG) >> 3) & 3)

{

case 0:

GYRO_K = 1./131/57.3;

break;

case 1:

GYRO_K = 1./65.5/57.3;

break;

case 2:

GYRO_K = 1./32.8/57.3;

break;

case 3:

GYRO_K = 1./16.4/57.3;

break;

}

}


//---------------------------------------------------------------------------------------------------

// Fast inverse square-root

// See: http://en.wikipedia.org/wiki/Fast_inverse_square_root


float invSqrt(float x) 

{

float halfx = 0.5f * x;

float y = x;

long i = *(long*)&y;

i = 0x5f3759df - (i>>1);

y = *(float*)&i;

y = y * (1.5f - (halfx * y * y));

return y;

}


//---------------------------------------------------------------------------------------------------

// AHRS algorithm update

//---------------------------------------------------------------------------------------------------

// IMU algorithm update


void MadgwickAHRSupdate_6(float gx, float gy, float gz, float ax, float ay, float az) 

{

float recipNorm;

float s0, s1, s2, s3;

float qDot1, qDot2, qDot3, qDot4;

float _2q0, _2q1, _2q2, _2q3, _4q0, _4q1, _4q2 ,_8q1, _8q2, q0q0, q1q1, q2q2, q3q3;


//将陀螺仪AD值转换为 弧度/s

gx = gx * GYRO_K;

gy = gy * GYRO_K;

gz = gz * GYRO_K;


// Rate of change of quaternion from gyroscope

qDot1 = 0.5f * (-q1 * gx - q2 * gy - q3 * gz);

qDot2 = 0.5f * (q0 * gx + q2 * gz - q3 * gy);

qDot3 = 0.5f * (q0 * gy - q1 * gz + q3 * gx);

qDot4 = 0.5f * (q0 * gz + q1 * gy - q2 * gx);


// Compute feedback only if accelerometer measurement valid (avoids NaN in accelerometer normalisation)

if(!((ax == 0.0f) && (ay == 0.0f) && (az == 0.0f))) {


// Normalise accelerometer measurement

recipNorm = invSqrt(ax * ax + ay * ay + az * az);

ax *= recipNorm;

ay *= recipNorm;

az *= recipNorm;   


// Auxiliary variables to avoid repeated arithmetic

_2q0 = 2.0f * q0;

_2q1 = 2.0f * q1;

_2q2 = 2.0f * q2;

_2q3 = 2.0f * q3;

_4q0 = 4.0f * q0;

_4q1 = 4.0f * q1;

_4q2 = 4.0f * q2;

_8q1 = 8.0f * q1;

_8q2 = 8.0f * q2;

q0q0 = q0 * q0;

q1q1 = q1 * q1;

q2q2 = q2 * q2;

q3q3 = q3 * q3;


// Gradient decent algorithm corrective step

s0 = _4q0 * q2q2 + _2q2 * ax + _4q0 * q1q1 - _2q1 * ay;

s1 = _4q1 * q3q3 - _2q3 * ax + 4.0f * q0q0 * q1 - _2q0 * ay - _4q1 + _8q1 * q1q1 + _8q1 * q2q2 + _4q1 * az;

s2 = 4.0f * q0q0 * q2 + _2q0 * ax + _4q2 * q3q3 - _2q3 * ay - _4q2 + _8q2 * q1q1 + _8q2 * q2q2 + _4q2 * az;

s3 = 4.0f * q1q1 * q3 - _2q1 * ax + 4.0f * q2q2 * q3 - _2q2 * ay;

recipNorm = invSqrt(s0 * s0 + s1 * s1 + s2 * s2 + s3 * s3); // normalise step magnitude

s0 *= recipNorm;

s1 *= recipNorm;

s2 *= recipNorm;

s3 *= recipNorm;


// Apply feedback step

qDot1 -= beta * s0;

qDot2 -= beta * s1;

qDot3 -= beta * s2;

qDot4 -= beta * s3;

}


// Integrate rate of change of quaternion to yield quaternion

q0 += qDot1 * (1.0f / sampleFreq);

q1 += qDot2 * (1.0f / sampleFreq);

q2 += qDot3 * (1.0f / sampleFreq);

q3 += qDot4 * (1.0f / sampleFreq);


// Normalise quaternion

recipNorm = invSqrt(q0 * q0 + q1 * q1 + q2 * q2 + q3 * q3);

q0 *= recipNorm;

q1 *= recipNorm;

q2 *= recipNorm;

q3 *= recipNorm;


Pitch = asin(-2.0f * (q1*q3 - q0*q2))* 57.3f;

Roll = atan2(q0*q1 + q2*q3, 0.5f - q1*q1 - q2*q2) * 57.3f;

Yaw = atan2(q1*q2 + q0*q3, 0.5f - q2*q2 - q3*q3)* 57.3f;

}


void MadgwickAHRSupdate_9(float gx, float gy, float gz, float ax, float ay, float az, float mx, float my, float mz) 

{

float recipNorm;

float s0, s1, s2, s3;

float qDot1, qDot2, qDot3, qDot4;

float hx, hy;

float _2q0mx, _2q0my, _2q0mz, _2q1mx, _2bx, _2bz, _4bx, _4bz, _2q0, _2q1, _2q2, _2q3, _2q0q2, _2q2q3, q0q0, q0q1, q0q2, q0q3, q1q1, q1q2, q1q3, q2q2, q2q3, q3q3;


// Use IMU algorithm if magnetometer measurement invalid (avoids NaN in magnetometer normalisation)

if((mx == 0.0f) && (my == 0.0f) && (mz == 0.0f)) 

{

MadgwickAHRSupdate_6(gx, gy, gz, ax, ay, az);

return;

}


//将陀螺仪AD值转换为 弧度/s

gx = gx * GYRO_K;

gy = gy * GYRO_K;

gz = gz * GYRO_K;


// Rate of change of quaternion from gyroscope

qDot1 = 0.5f * (-q1 * gx - q2 * gy - q3 * gz);

qDot2 = 0.5f * (q0 * gx + q2 * gz - q3 * gy);

qDot3 = 0.5f * (q0 * gy - q1 * gz + q3 * gx);

qDot4 = 0.5f * (q0 * gz + q1 * gy - q2 * gx);


// Compute feedback only if accelerometer measurement valid (avoids NaN in accelerometer normalisation)

if(!((ax == 0.0f) && (ay == 0.0f) && (az == 0.0f))) 

{


// Normalise accelerometer measurement

recipNorm = invSqrt(ax * ax + ay * ay + az * az);

ax *= recipNorm;

ay *= recipNorm;

az *= recipNorm;   


// Normalise magnetometer measurement

recipNorm = invSqrt(mx * mx + my * my + mz * mz);

mx *= recipNorm;

my *= recipNorm;

mz *= recipNorm;


// Auxiliary variables to avoid repeated arithmetic

_2q0mx = 2.0f * q0 * mx;

_2q0my = 2.0f * q0 * my;

_2q0mz = 2.0f * q0 * mz;

_2q1mx = 2.0f * q1 * mx;

_2q0 = 2.0f * q0;

_2q1 = 2.0f * q1;

_2q2 = 2.0f * q2;

_2q3 = 2.0f * q3;

_2q0q2 = 2.0f * q0 * q2;

_2q2q3 = 2.0f * q2 * q3;

q0q0 = q0 * q0;

q0q1 = q0 * q1;

q0q2 = q0 * q2;

q0q3 = q0 * q3;

q1q1 = q1 * q1;

q1q2 = q1 * q2;

q1q3 = q1 * q3;

q2q2 = q2 * q2;

q2q3 = q2 * q3;

q3q3 = q3 * q3;


// Reference direction of Earth's magnetic field

hx = mx * q0q0 - _2q0my * q3 + _2q0mz * q2 + mx * q1q1 + _2q1 * my * q2 + _2q1 * mz * q3 - mx * q2q2 - mx * q3q3;

hy = _2q0mx * q3 + my * q0q0 - _2q0mz * q1 + _2q1mx * q2 - my * q1q1 + my * q2q2 + _2q2 * mz * q3 - my * q3q3;

_2bx = sqrt(hx * hx + hy * hy);

_2bz = -_2q0mx * q2 + _2q0my * q1 + mz * q0q0 + _2q1mx * q3 - mz * q1q1 + _2q2 * my * q3 - mz * q2q2 + mz * q3q3;

_4bx = 2.0f * _2bx;

_4bz = 2.0f * _2bz;


// Gradient decent algorithm corrective step

s0 = -_2q2 * (2.0f * q1q3 - _2q0q2 - ax) + _2q1 * (2.0f * q0q1 + _2q2q3 - ay) - _2bz * q2 * (_2bx * (0.5f - q2q2 - q3q3) + _2bz * (q1q3 - q0q2) - mx) + (-_2bx * q3 + _2bz * q1) * (_2bx * (q1q2 - q0q3) + _2bz * (q0q1 + q2q3) - my) + _2bx * q2 * (_2bx * (q0q2 + q1q3) + _2bz * (0.5f - q1q1 - q2q2) - mz);

s1 = _2q3 * (2.0f * q1q3 - _2q0q2 - ax) + _2q0 * (2.0f * q0q1 + _2q2q3 - ay) - 4.0f * q1 * (1 - 2.0f * q1q1 - 2.0f * q2q2 - az) + _2bz * q3 * (_2bx * (0.5f - q2q2 - q3q3) + _2bz * (q1q3 - q0q2) - mx) + (_2bx * q2 + _2bz * q0) * (_2bx * (q1q2 - q0q3) + _2bz * (q0q1 + q2q3) - my) + (_2bx * q3 - _4bz * q1) * (_2bx * (q0q2 + q1q3) + _2bz * (0.5f - q1q1 - q2q2) - mz);

s2 = -_2q0 * (2.0f * q1q3 - _2q0q2 - ax) + _2q3 * (2.0f * q0q1 + _2q2q3 - ay) - 4.0f * q2 * (1 - 2.0f * q1q1 - 2.0f * q2q2 - az) + (-_4bx * q2 - _2bz * q0) * (_2bx * (0.5f - q2q2 - q3q3) + _2bz * (q1q3 - q0q2) - mx) + (_2bx * q1 + _2bz * q3) * (_2bx * (q1q2 - q0q3) + _2bz * (q0q1 + q2q3) - my) + (_2bx * q0 - _4bz * q2) * (_2bx * (q0q2 + q1q3) + _2bz * (0.5f - q1q1 - q2q2) - mz);

s3 = _2q1 * (2.0f * q1q3 - _2q0q2 - ax) + _2q2 * (2.0f * q0q1 + _2q2q3 - ay) + (-_4bx * q3 + _2bz * q1) * (_2bx * (0.5f - q2q2 - q3q3) + _2bz * (q1q3 - q0q2) - mx) + (-_2bx * q0 + _2bz * q2) * (_2bx * (q1q2 - q0q3) + _2bz * (q0q1 + q2q3) - my) + _2bx * q1 * (_2bx * (q0q2 + q1q3) + _2bz * (0.5f - q1q1 - q2q2) - mz);

recipNorm = invSqrt(s0 * s0 + s1 * s1 + s2 * s2 + s3 * s3); // normalise step magnitude

s0 *= recipNorm;

s1 *= recipNorm;

s2 *= recipNorm;

s3 *= recipNorm;


// Apply feedback step

qDot1 -= beta * s0;

qDot2 -= beta * s1;

qDot3 -= beta * s2;

qDot4 -= beta * s3;

}


// Integrate rate of change of quaternion to yield quaternion

q0 += qDot1 * (1.0f / sampleFreq);

q1 += qDot2 * (1.0f / sampleFreq);

q2 += qDot3 * (1.0f / sampleFreq);

q3 += qDot4 * (1.0f / sampleFreq);


// Normalise quaternion

recipNorm = invSqrt(q0 * q0 + q1 * q1 + q2 * q2 + q3 * q3);

q0 *= recipNorm;

q1 *= recipNorm;

q2 *= recipNorm;

q3 *= recipNorm;


Pitch = asin(-2.0f * (q1*q3 - q0*q2))* 57.3f;

Roll = atan2(q0*q1 + q2*q3, 0.5f - q1*q1 - q2*q2) * 57.3f;

Yaw = atan2(q1*q2 + q0*q3, 0.5f - q2*q2 - q3*q3)* 57.3f;

}


//====================================================================================================

// END OF CODE

//====================================================================================================

[1] [2]
关键字:51单片机  AHRS  四元数  欧拉角 编辑:什么鱼 引用地址:【51单片机快速入门指南】4.4.3:Madgwick AHRS 九轴姿态融合获取四元数、欧拉角

上一篇:【51单片机快速入门指南】4.5:I2C 与 TCA6416实现双向 IO 扩展
下一篇:【51单片机快速入门指南】4.4.2:Mahony AHRS 九轴姿态融合获取四元数、欧拉角

推荐阅读

用AT89C51单片机显示倒计时程序
;可设定时间的倒计时定时器,可选择5/15/20/30/35/45/50分钟倒计时;倒计时时间由四位拨码开关的2/3/4位来控制,;第2位表示5分钟,第3位表示15分钟,第4位表示30分钟,;通过不同的组合可以产生5/15/20/30/35/45/50分钟倒计时;P1.0口的外接的发光二极管为状态LED,定时未开始时LED常亮,定时过程中LED闪烁;K1为开始按钮,K2为停止按钮适用STM8S/STM8L/STM8A N76E003 脱机编程器/烧录器/下载器/SP_00【包邮】m.tb.cn/h.UlXVKiOa_bit equ 20h ;数码管个位数存放内存位置b_bit equ 21h ;数码管十位数存放内存位置temp eq
发表于 2023-01-13
51单片机数码管静态显示和动态显示原理及实验 夜猫子
数码管多位数码管,即是两个或两个以上单个数码管并列集中在一起形成一体的数码管。当多位一体时,它们内部的公共端是独立的,而负责显示什么数字的段线全部是连接在一起的,独立的公共端可以控制多位一体中的哪一位数码管点亮,而连接在一起的段线可以控制这个能点亮数码管亮什么数字,通常我们把公共端叫做“位选线”,连接在一起的段线叫做“段选线”有了这两个线后,通过单片机及外部驱动电路就可以控制任意的数码管显示。一般一位数码管有10个引脚,二位数码管也是10个引脚,四位数码管是12个引脚。为了更方便区分段选和位选,请看下原理图:如图为两个4位一体的数码管,可以看到与8个com相连的是两个数码管的位选,位选与引脚相连,所以位选控制那个灯亮。段选可以看到a
发表于 2023-01-12
<font color='red'>51单片机</font>数码管静态显示和动态显示原理及实验 夜猫子
51单片机独立按键和矩阵按键实现
独立按键实验按键是一种电子开关,使用时轻轻按开关按钮就可使开关接通,当松开手时,开关断开。我们开发板上使用的按键及内部简易图如下图所示管脚与管脚之间(注意是距离)距离长的是导通状态,短的是接通状态。通常的按键所用开关为机械弹性开关,当机械触点断开、闭合时,电压信号如下图所示:如图所示,按键闭合式不会立刻稳定的接通,断开时也不会一下子断开,会伴随一些抖动。抖动的时间长短有按键特性决定,一般为5Ms到10ms.按键抖动会引起按键被误读多次。为了确保 CPU 对按键的一次闭合仅作一次处理,必须进行消抖。消抖消抖可分为硬件消抖和软件消抖。为了使电路更加简单,通常采用软件消抖。一般来说一个简单的按键消抖就是先读取按键的状态, 如果得到按键按下
发表于 2023-01-12
51单片机8*8点阵原理及实现
LED点阵(8*8)LED 点阵是由发光二极管排列组成的显示器件,在我们日常生活的电器中随处可见,被广泛应用于汽车报站器,广告屏等。通常应用较多的是 8* 8 点阵,然后使用多个 8 * 8 点阵可组成不同分辨率的 LED点阵显示屏,比如 16* 16 点阵可以使用 4 个 8* 8 点阵构成。因此理解了 8* 8LED点阵的工作原理,其他分辨率的 LED 点阵显示屏都是一样的。这里以 8* 8LED 点阵来做介绍。发光原理8* 8 点阵共由 64 个发光二极管组成,且每个发光二极管是放置在行线和列线的交叉点上,当对应的某一行置 1 电平(行所接的是二极管的阳极,所以为高电平),某一列置 0 电平(列所接的是二极管的阴极极,所以为低
发表于 2023-01-12
51单片机中断基本概念
问题引入在了解基本概念之前,先看三个问题:1.你想使用的中断是哪个?2.你所希望的触发条件是什么?3.你希望在中断之后做什么?可以边看边思考,文章最后给出答案中断概念为什么引入中断?中断是为使单片机具有对外部或内部随机发生的事件实时处理而设置的,中断功能的存在,很大程度上提高了单片机处理外部或内部事件的能力。中断系统特点:①分时操作。CPU 可以分时为多个 I/O 设备服务,提高了计算机的利用率;②实时响应。CPU 能够及时处理应用系统的随机事件,系统的实时性大大增强;③可靠性高。CPU 具有处理设备故障及掉电等突发性事件能力,从而使系统可靠性提高中断过程对于单片机来讲,中断是指CPU在处理某一时间A时,发生了另一事件B请求CPU立
发表于 2023-01-12
51单片机外部中断点亮LED
外部中断软件设计原理中断发生的三个条件①中断源有中断请求;②此中断源的中断允许位为 1;③CPU 开中断(即 EA=1)。比如我们配置外部中断 0,对应的配置程序如下:EA=1;//打开总中断开关EX0=1;//开外部中断 0IT0=0/1;//设置外部中断的触发方式(下降沿触发)如果要配置的是外部中断 1,只需将 EX0 改为 EX1,IT0 改为 IT1在编写程序时通常我们会将外部中断的配置放到一个自定义函数内便于管理维护。如下伪代码所示:void Int0Init(){//设置 INT0IT0=1;//边沿触发方式(下降沿)EX0=1;//打开 INT0 的中断允许。EA=1;//打开总中断}/*当触发中断后即会进入中断服务函
发表于 2023-01-12
小广播
设计资源 培训 开发板 精华推荐

何立民专栏 单片机及嵌入式宝典

北京航空航天大学教授,20余年来致力于单片机与嵌入式系统推广工作。

换一换 更多 相关热搜器件
电子工程世界版权所有 京B2-20211791 京ICP备10001474号-1 电信业务审批[2006]字第258号函 京公网安备 11010802033920号 Copyright © 2005-2023 EEWORLD.com.cn, Inc. All rights reserved