C++演示IMU滤波器:互补、卡尔曼与Mahony/Madgwick对比
2026/9/15 7:53:54 网站建设 项目流程

简介:面向嵌入式、机器人、无人机及自动驾驶学习者的IMU滤波算法演示项目,用C++实现互补滤波器、卡尔曼滤波器以及Mahony/Madgwick三类常用姿态解算算法,可对加速度计与陀螺仪数据进行实时融合,直接服务于姿态估计与运动控制应用开发。代码框架开放,内部整合了Qt和QCustomPlot绘图组件,能够将原始传感器数据与三种滤波结果同步显示在界面上,便于对比不同算法在动态响应、抗噪性能和稳定性上的差异,对理解滤波器调参规律和选型很有帮助。压缩包共32个文件,主要由dll动态库、h头文件、cpp源文件及pro/ui工程文件构成,附带的README文档说明了使用方法和参考来源,整体仅545KB,轻量易读。资源还整理了大量外部教程与算法详解,覆盖从基础IMU指南到DCM、数据融合等进阶主题,可帮助系统梳理IMU数据处理知识。已有2007人学习下载,适合嵌入式爱好者、算法工程师和学生作为学习与二次开发的参考。

1. 为什么需要一个演示 IMU 滤波器的 C++ 项目

第一次把六轴 IMU 数据接进上位机时,大多数人都会卡在同一步:陀螺仪积分出来的角度三分钟后漂出好几度,加速度计在跑步机上抖得像噪声发生器,随手抄来的融合公式又看不懂四元数。这个标题所描述的“用 C++ 演示 IMU 滤波器”项目,本质上是把姿态估计里最常用的三条路线——互补滤波、卡尔曼滤波、Mahony 和 Madgwick 滤波——放进同一份工程里,让你先用仿真数据把算法跑明白,再决定哪条路线适合自己手头的传感器和算力。目标读者不是算法研究员,而是要做运动检测、平衡控制、AR 姿态标定的一线开发:你不需要推公式推到呕心沥血,但至少要能区分“这个滤波器的滞后是哪来的”和“参数往哪个方向调”。

2. 滤波前先理解 IMU 数据与三种算法各自的融合思路

把代码之前,先花三节时间把传感器模型讲清楚。因为这一项目里所谓“演示”,最容易被误解成“拿现成公式调用一遍”;真正有价值的其实是三个滤波器对同一个观测噪声的不同假设。

2.1 IMU 的原始误差从哪来:陀螺仪积分漂移与加速度计非重力干扰

陀螺仪输出的是角速度,姿态必须通过对角速度做数值积分才能得到。积分会把零偏积累起来,一个稳定为 0.01 rad/s 的常值零偏,在一分钟后就会造成约 0.6 弧度的姿态偏差,换算成角度大约 34 度。这就是“飘”的根源。加速度计则相反,它直接测量比力,静止时重力方向可以从中解算出俯仰和横滚角,但它对线性加速度没有区分能力,运动越剧烈,水平加速度混入越严重。

两者误差特性互补,所以需要融合。一种融合思路在频域看最直观:陀螺仪积分后的姿态在低频段可信度高,但高频漂移大;加速度计解算的姿态在高频噪声大,长期却没有漂移。取长补短的直接手段可以是一个高通滤波和一个低通滤波的组合,这就是互补滤波的思想雏形。在 C++ 里演示时,我不会拿真实传感器直接跑,而是先建立一个带噪声和零偏的合成信号源,让三个滤波器面对完全相同的输入;这样对比出来的“谁更准”才有意义。

struct ImuSample { float gx, gy, gz; // 角速度,单位 rad/s float ax, ay, az; // 比力,单位 m/s^2 }; // 合成噪声模型:固定零偏 + 高斯白噪声 ImuSample makeSyntheticSample(float t, float bias, float noiseStd) { ImuSample s; s.gx = 0.3f * sinf(2.0f * 3.14159f * 0.5f * t) + bias + noise(); s.gy = -0.2f * cosf(2.0f * 3.14159f * 0.3f * t) + bias + noise(); s.gz = 0.0f; // 加速度计由理论姿态投影到载体坐标,再加白噪声 // pitch = 0.5sin(2pi*0.5t),roll = 0.3cos(2pi*0.3t) return s; }

这段代码里noise()可以是任何 Box-Muller 变换实现。注意我刻意把零偏设成常值,因为它能直观展示卡尔曼滤波和 Mahony 积分项对零偏的补偿能力;高斯噪声则用来模拟加速度计的高频抖动。

2.2 互补滤波器:为什么 alpha 公式比“陀螺仪占 98%”更值得记

互补滤波的实现者常把输出写成angle = alpha * (angle + gyro*dt) + (1-alpha) * accAngle,然后告诉你 alpha 取 0.98。这个公式没错,但它遗漏了一个关键量:时间常数。

对一阶互补滤波而言,alpha 与采样周期 dt 共同决定融合截止频率,关系如式tau = (alpha * dt) / (1 - alpha)。假设控制周期 10ms,alpha=0.98,tau=0.49 秒;这意味着来自陀螺仪的积分会在 0.49 秒的时间尺度上被加速度计“拉回”。把 alpha 提到 0.995,tau 就变成 1.99 秒,短期动态保留得更好,但动态加速度干扰也更容易混入。所以正确调参逻辑不是“越大越好”,而是“先定 tau,再反算 alpha”。

这个结论在演示项目中最适合用一张参数表体现:

目标场景建议 tau10ms 周期对应 alpha后果
静态倾角测量0.1~0.3s0.90~0.97响应快,静止时抖动明显
无人机姿态控制0.5~1.0s0.98~0.99动态下更平滑,滚转略滞后
低帧率处理(50Hz)1.0~2.0s0.98~0.99积分噪声被拉回力度减小

看这张表时会发现一个规律:alpha 对采样率极其敏感,如果直接把 100Hz 调好的 0.98 拿到 20Hz 系统里,等效时间常数会变化,所以项目代码里应该用alpha = tau / (tau + dt)反向计算,而不是把 alpha 写死成常量。做这个项目时,我会在配置文件中只暴露 tau,所有滤波器内部统一换算,这样三个算法在对比时至少处在同一时间尺度上。

2.3 卡尔曼滤波在姿态问题中的状态空间假设:它到底在估计什么

卡尔曼滤波在姿态估计中经常被简化成一维角度滤波器,这是演示项目最常用的入口,也是最容易误解的地方。完整的三维姿态卡尔曼需要 6 到 7 维状态量,同时估计角度和陀螺仪零偏,推导与调参会迅速复杂化,对演示来说并不友好。

一维模型的状态向量是x = [angle, bias]^T,过程方程写为:

angle_new = angle_old + (gyro - bias) * dt bias_new = bias_old

其中角速度被视为系统输入,加速度计解算角度被视为观测。这样做的前提假设是:角速度短期内的积分是可用的,而 bias 是缓变常数。卡尔曼滤波的贡献在于,它能根据噪声统计特性动态调整“相信陀螺仪还是相信加速度计”,这比固定 alpha 的互补滤波更灵活,但代价是必须给出合理的过程噪声 Q 和观测噪声 R。

2.4 Mahony&Madgwick 滤波器的关键差异:叉积修正与梯度下降

Mahony 和 Madgwick 虽然经常被并列提起,但它们对“加速度计如何修正陀螺仪”的建模完全不同。Mahony 把问题看成 PI 控制器:用加速度计实际观测值与姿态旋转矩阵推算出的理论重力方向做叉积,叉积大小近似等于角度误差,将误差通过kpki反馈到陀螺仪。Madgwick 则把它看成优化问题:构造一个关于四元数的目标函数,计算其梯度方向,然后沿着梯度反方向走一小步来纠正四元数,这一步由beta控制。

换句话说,Mahony 的反馈增益只有两个数,物理意义明确;Madgwick 的beta代表对测量误差的信任程度,值越大纠正越快,但会把加速度计噪声直接引入姿态。两者在代码结构上非常相似,都依赖同一个四元数微分方程:

q_dot = 0.5 * q ⊗ omega

差别仅在于把“误差修正项”以不同方式叠加进q_dot。演示项目里如果把这两种都实现,就能直观看到:相同参数、相同噪声下,Mahony 更容易调平,Madgwick 在动态激励更剧烈时表现更稳。这一点放到后面调参部分再展开。

3. 在 C++ 里落地互补、卡尔曼与 Mahony&Madgwick 的最小实现

这一章直接面对标题里“C++_代码_相关文件”这个检索意图。下面给出的实现刻意保持最小依赖,不引入 Eigen,不加线程,目的是让你能照抄到自己的工程里改。我会先定义公共数据结构,然后逐个说明每个滤波器在回调里应如何更新。

3.1 定义一个公共四元数结构和采样周期约定

所有滤波器输出统一用四元数,因为欧拉角存在万向锁和多值性问题;调试时如果需要看角度,再在最后一层把四元数转成欧拉角。四元数约定为w, x, y, z,单位化,旋转顺序只影响矩阵展开细节,在一个项目内保持一致即可。

struct Quat { float w, x, y, z; void normalize() { float n = sqrtf(w*w + x*x + y*y + z*z); if (n > 1e-8f) { w /= n; x /= n; y /= n; z /= n; } } };

所有滤波器都遵守同一个接口:传入陀螺仪角速度、加速度计比力、dt,更新内部姿态。不把 dt 设为常量,是为了让回放离线日志时能处理变采样率,也便于测试不同频率对同一个 alpha/beta 的影响。

3.2 互补滤波器:面向欧拉角的最简实现与 tau 换算

下面的类保留了最常见的欧拉角形式,因为互补滤波通常被人用于直接输出俯仰和滚转角。实际使用中我建议把它输出的欧拉角转成四元数,避免后续控制链路出现角度跳变。

class ComplementaryFilter { public: void init(float tau, float dt) { float ts = dt; // 采样周期 alpha_ = tau / (tau + ts); // 统一按时间常数配置 } void update(float gx, float gy, float accPitch, float accRoll) { pitch_ += gx * dt_; roll_ += gy * dt_; pitch_ = alpha_ * pitch_ + (1.0f - alpha_) * accPitch; roll_ = alpha_ * roll_ + (1.0f - alpha_) * accRoll; } private: float alpha_ = 0.98f; float dt_ = 0.01f; float pitch_ = 0.0f, roll_ = 0.0f; };

这段代码里,alpha_ * pitch前一半在“信任陀螺仪的积分”,(1-alpha) * accPitch在“信任加速度计的绝对角度”。注意陀螺仪积分前没有手动减零偏,所以这个版本的互补滤波无法消除零偏漂移;用于演示时在仿真里加零偏,自然能看到它输出在长时间尺度上偏离真值。把alpha_换成alpha_ = tau / (tau + dt_)后,即使采样率变化,频率特性也能保持一致。

3.3 一维卡尔曼滤波:用两个状态量同时估计角度与零偏

这里给出一个可以在演示项目中跑起来的一维卡尔曼实现。状态量为角度和陀螺仪零偏,输入是陀螺仪 raw 值,观测是加速度计解算角。由于零偏被显式建模,它的长期漂移抑制能力会明显好于上面互补滤波版本。

class KalmanAngle { public: KalmanAngle() { P_ = { {1, 0}, {0, 1} }; } void setAngleNoise(float angleNoise, float biasNoise) { Q_ = { {angleNoise, 0}, {0, biasNoise} }; } void setMeasureNoise(float r) { R_ = r; } float update(float gyroRate, float accAngle, float dt) { // 预测:使用上一时刻的零偏估计修正角速度 angle_ += (gyroRate - bias_) * dt; P_[0][0] += dt * (dt * Q_[1][1] + Q_[0][0]); P_[0][1] -= dt * Q_[1][1]; P_[1][0] -= dt * Q_[1][1]; P_[1][1] += dt * Q_[1][1]; // 更新:卡尔曼增益由预测协方差与观测噪声决定 float S = P_[0][0] + R_; float K0 = P_[0][0] / S; float K1 = P_[1][0] / S; float y = accAngle - angle_; // 新息 angle_ += K0 * y; bias_ += K1 * y; float P00 = P_[0][0], P01 = P_[0][1], P10 = P_[1][0], P11 = P_[1][1]; P_[0][0] = (1 - K0) * P00; P_[0][1] = (1 - K0) * P01; P_[1][0] = -K1 * P00 + P10; P_[1][1] = -K1 * P01 + P11; return angle_; } private: float angle_ = 0.0f, bias_ = 0.0f; float R_ = 0.1f; float Q_[2][2]; float P_[2][2]; };

这段代码的关键点在于:P_[0][1]P_[1][0]不再是对称零矩阵,因为角速度与零偏通过积分耦合在一起。Q_[1][1]对应“零偏随机游走强度”,它决定滤波器能不能跟随缓慢变化的零偏:设为 1e-6 时滤波器认为零偏基本不变,设到 1e-4 时则对零偏变化更敏感,但也更容易被噪声带偏。演示时把这个参数漏掉,是“卡尔曼完全不服”的最常见原因。

3.4 Mahony 与 Madgwick 的四元数核心代码

Mahony 的更新函数短小精悍,核心就是叉积修正加上 PI 控制。为了便于阅读,我把积分误差作为静态量保留在函数外。

class MahonyFilter { public: void update(Quat& q, float gx, float gy, float gz, float ax, float ay, float az, float dt) { float norm = sqrtf(ax*ax + ay*ay + az*az); if (norm < 1e-6f) return; ax /= norm; ay /= norm; az /= norm; // 当前四元数下的理论重力方向(载体坐标) float q0=q.w, q1=q.x, q2=q.y, q3=q.z; float vx = 2.0f*(q1*q3 - q0*q2); float vy = 2.0f*(q0*q1 + q2*q3); float vz = q0*q0 - q1*q1 - q2*q2 + q3*q3; // 加速度计观测量与理论值之间的叉积误差 float ex = ay*vz - az*vy; float ey = az*vx - ax*vz; float ez = ax*vy - ay*vx; exInt_ += ex * ki_ * dt; eyInt_ += ey * ki_ * dt; ezInt_ += ez * ki_ * dt; gx += kp_*ex + exInt_; gy += kp_*ey + eyInt_; gz += kp_*ez + ezInt_; // 四元数一阶积分 q.w += 0.5f*dt * (-q1*gx - q2*gy - q3*gz); q.x += 0.5f*dt * ( q0*gx - q3*gy + q2*gz); q.y += 0.5f*dt * ( q3*gx + q0*gy - q1*gz); q.z += 0.5f*dt * (-q2*gx + q1*gy + q0*gz); q.normalize(); } private: float kp_ = 0.5f, ki_ = 0.0f; float exInt_ = 0.0f, eyInt_ = 0.0f, ezInt_ = 0.0f; };

代码中kp_是比例增益,ki_是积分增益。只设比例项也能跑,零偏导致的稳态误差会残留在姿态中;把ki_设到kp_的十分之一量级,长时间漂移通常能压到较低水平,但过大的ki_容易引起低频振荡,表现是姿态在静止时缓慢上下摆动。

Madgwick 的加速度计部分通常写成梯度下降步长的形式。把理论目标函数f_g与雅可比矩阵J_g合并计算,就得到下面的梯度步s

class MadgwickFilter { public: void update(Quat& q, float gx, float gy, float gz, float ax, float ay, float az, float beta, float dt) { float norm = sqrtf(ax*ax + ay*ay + az*az); if (norm < 1e-6f) return; ax /= norm; ay /= norm; az /= norm; float q0=q.w, q1=q.x, q2=q.y, q3=q.z; float f1 = 2.0f*(q1*q3 - q0*q2) - ax; float f2 = 2.0f*(q0*q1 + q2*q3) - ay; float f3 = 1.0f - 2.0f*(q1*q1 + q2*q2) - az; float s0 = -2.0f*q2*f1 + 2.0f*q3*f2 - 2.0f*q0*f3; float s1 = 2.0f*q3*f1 + 2.0f*q0*f2 - 4.0f*q1*f3; float s2 = -2.0f*q0*f1 + 2.0f*q3*f2 - 4.0f*q2*f3; float s3 = 2.0f*q1*f1 + 2.0f*q2*f2; float snorm = sqrtf(s0*s0 + s1*s1 + s2*s2 + s3*s3); s0 /= snorm; s1 /= snorm; s2 /= snorm; s3 /= snorm; // 梯度下降修正项与陀螺仪积分项合成 q.w += 0.5f*dt * (-q1*gx - q2*gy - q3*gz) - beta*dt*s0; q.x += 0.5f*dt * ( q0*gx + q2*gz - q3*gy) - beta*dt*s1; q.y += 0.5f*dt * ( q0*gy - q1*gz + q3*gx) - beta*dt*s2; q.z += 0.5f*dt * ( q0*gz + q1*gy - q2*gx) - beta*dt*s3; q.normalize(); } };

这里s的含义是“让目标函数下降最快的单位方向”,beta决定沿该方向前进的速度。beta = 0时 Madgwick 退化为纯陀螺仪积分;beta太大则等于把加速度计噪声原样卷积进四元数。与原版论文不同,我故意没写磁力计分支,因为在大多数视频演示项目里,磁力计数据反而会引入更多环境干扰。

4. 用合成数据测试:噪声参数、初始值与误差评估怎么做

第 3 章代码能跑起来之后,真正决定这个演示项目能不能给人信服结果的是测试方法。循环播放理想正弦信号没有意义:必须让加速度计加入噪声、陀螺仪加入零偏,再对比三种估计和真值之间的误差。

4.1 合成测试信号怎么生成,才能看出三种滤波器的差异

测试信号我建议包含三段:静止 10 秒,正弦摆动 20 秒,阶跃突变 5 秒。静止段暴露零偏漂移;正弦段暴露动态滞后,正弦频率应设在 0.2~1 Hz;阶跃段暴露超调。加速度计噪声标准差设为 0.15 m/s²,陀螺仪零偏设为 0.02 rad/s,这个量级接近常见低成本 IMU 的表现。

// 按时间生成理论欧拉角,再投影成加速度计真值 void generateTruthAndRaw(float t, float& pitch, float& roll, float& ax, float& ay, float& az) { if (t < 10.0f) { pitch = 0.0f; roll = 0.0f; } else if (t < 30.0f) { pitch = 0.35f * sinf(2 * M_PI * 0.5f * (t - 10)); roll = 0.25f * sinf(2 * M_PI * 0.3f * (t - 10)); } else { pitch = 0.4f; roll = -0.3f; } // 世界系重力 [0,0,g] 经旋转矩阵转到载体系 // 这里演示用简化的欧拉角正变换 ax = -sinf(pitch) * 9.81f + noiseStd * gaussian(); ay = cosf(pitch) * sinf(roll) * 9.81f + noiseStd * gaussian(); az = cosf(pitch) * cosf(roll) * 9.81f + noiseStd * gaussian(); }

注意上面把noiseStd * gaussian()放在ax/ay/az上,但没在陀螺仪上模拟运动过程,而是直接计算真值角速度再叠加零偏,这样后续评估时才有一条绝对可靠的真值曲线。实际项目里陀螺仪信号也应由姿态微分得到,否则“零偏补偿效果”无从体现。

4.2 三个滤波器的关键参数怎么调:一组能用的初始值

没有磁力计参与时,各滤波器的调参入口完全不同。互补滤波只有一个tau;卡尔曼至少四个参数Q_angle/Q_bias/R_;Mahony 两个增益kp/ki;Madgwick 只有一个beta。下面这组数值适合 100Hz、中等噪声,可以作为调试起点,但不应该直接搬进成品。

滤波器参数初始建议值调参方向
互补滤波tau0.5~1.0s静止噪声大则增大 tau,动态滞后大则减小 tau
卡尔曼Q_angle1e-4如果跟踪慢,减小;如果噪声大,增大
卡尔曼Q_bias1e-5零偏漂移补偿不足则增大
卡尔曼R0.1加速度计噪声越大,R 应越大
Mahonykp0.5动态响应慢则增大,噪声放大则减小
Mahonyki0.05长时间漂移大则增大,但不建议超过 kp/5
Madgwickbeta0.1静止抖动大则减小,阶跃恢复慢则增大

调参口诀只有两句:优先满足“静止不飘”,再压“动态滞后”。如果静止时噪声已经很明显,无论哪个滤波器,第一反应都应该是增加对加速度计的信任度权重,而不是去动陀螺仪参数。

4.3 从原始数据到误差评估:输出 CSV 并计算 RMS、最大偏差和滞后

三个滤波器输出后,把四元数转成欧拉角,与真值做差。演示项目里最容易量化的三个指标是:静态 RMS、动态 RMS、阶跃后 90% 恢复时间。

  • 静态 RMS 评估滤波后的残余噪声水平,单位度。
  • 动态 RMS 评估跟踪误差,包含滞后导致的相位落后。
  • 90% 恢复时间评估滤波器对突变的响应速度,单位秒。

这部分不需要实时可视化,先把时间戳、真值、三个滤波器的姿态写入 CSV,然后直接用 Python 的 pandas 统计。我给项目配套脚本,让结果可重复输出:

fprintf(log, "%f,%f,%f,%f,%f\n", t, pitchTruth, compPitch, kalmanPitch, mahonyPitch);
import pandas as pd df = pd.read_csv("result.csv", header=0, names=["t","truth","comp","kalman","mahony"]) err_comp = (df["comp"] - df["truth"]).abs() err_kalman = (df["kalman"] - df["truth"]).abs() err_mahony = (df["mahony"] - df["truth"]).abs() print("互补 RMS:", err_comp.std()) print("卡尔曼 RMS:", err_kalman.std()) print("Mahony RMS:", err_mahony.std())

这里不直接打印均值,是因为静态段与动态段混在一起时均值会把正负误差抵消;标准差能同时反映噪声和平滑程度,静态段与动态段切分开后分别计算更有意义。把这两段分开看,你会得到两条结论:静止段卡尔曼最容易把 RMS 压到 0.1 度以下;动态段 Madgwick 的峰值误差通常比互补滤波小一个量级,但代价是启动阶段有一段明显的收敛过程。

4.4 可视化验证时重点看什么:两条曲线和一段误差曲线

四个滤波器同时输出时,图标不要画得太满,否则看不出差异。我建议一张图只画“真值 + 一个滤波器”,另一张图叠加所有滤波器的误差绝对值。关注三个位置:启动前 2 秒、正弦换向最高点、阶跃发生瞬间。启动段的差异暴露滤波器初值收敛能力;正弦段的峰值误差暴露滞后;阶跃段的超调量暴露阻尼特性。数据采样频率越高,这三个位置的曲线越可靠,建议至少保存 50 Hz,评估时用完整数据长度。

5. 把 Mahony&Madgwick 项目接入真实系统前的 4 个检查项与回放技巧

演示代码和真实嵌入式代码之间隔着一段距离,这里把最容易踩的四个坑列出来。检查项一:轴向一致性。芯片坐标系与算法假设坐标系不一致是头号问题,陀螺仪正方向、加速度计正方向、四元数旋转方向三者必须逐一核对,否则滤波收敛后会停在错误姿态上,且表现为静态误差而不是噪声。检查项二:采样间隔不能依赖sleep,要用时间戳差分;把固定 dt 假设延迟到真实数据上,会引入无法追踪的随机滞后。

检查项三:磁力计要不要加,要依据场景判断。在空中、在室内金属环境做姿态估计,磁力计信号质量往往比陀螺仪零偏更糟;没有磁力计的 Mahony 与 Madgwick 形式已经够用。若后续必须加磁力计,应先把偏置标定完成,再改进f_mJ_m分支,不要只往现有四元数上叠加磁力计观测。检查项四:参数不能只标定一次。alpha、beta、kp 都应该做成可在运行时调整的变量,用日志回放的方式离线调好,再固化到嵌入式端。

最后一个技巧是回放验证。演示项目最常见的价值反而不是“现场跑”,而是把之前录好的传感器日志喂进滤波器,反复比较改动前后差异。在 C++ 里实现回放只需让滤波器更新函数从std::istream读取,不与时间强绑定:

std::ifstream log("imu_log.csv"); float t, gx, gy, gz, ax, ay, az; while (log >> t >> gx >> gy >> gz >> ax >> ay >> az) { float dt = t - lastTime; filterMahony.update(q, gx, gy, gz, ax, ay, az, dt); // 也可以在这里不打印数据,直接累计误差统计量 }

回放时只需把同一份日志分别跑三次互补、卡尔曼、Mahony 与 Madgwick,固定输入、固定初值,比较结果才有说服力。实际项目中我一般会再把“上次现场记录的真实传感器噪声变成合成信号源”这一套流程也做进 CMake target,这样每次调完参数都可以无传感器回归验证。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询