首先是上一次更新到这一次更新的旋转矢量。这里存在圆锥误差，可以理解为有限大角增量不是矢量导致的转动不可交换误差。
使用N子样误差补偿，即在第m-1次更新到第m次更新种，陀螺仪等时间读取了N个值。这里一般不使用角速度，而是角增量，因为光纤陀螺等直接输出的是角增量，而MEMS这些均值平滑后输出角增量和角速度差别不大。
设更新周期是T，t从上一次更新开始计时，则理论上是：
$\Delta\vec{\theta}_{mi} = \int^{\frac{i}{N}T}_{\frac{i-1}{N}T}\vec{\omega}(t)dt$（A1）
实际上不可能做到连续的积分，最简单的长方形积或者梯形积分即可，我这个代码就是前者。（当然下一代还是计划搞成梯形积分）这里带来的不可交换误差会在后面求旋转矢量时补偿。
即：
$\Delta\vec{\theta}_{mi} = \frac{1}{N}T\vec{\omega}_{i}$ （A2）
$\Delta\vec{\theta}_{mi} = \frac{1}{2N}T(\vec{\omega}_{i-1}+\vec{\omega}_{i})$ （A3）
第m-1次更新到第m次更新过程中，第i次（i<N）角增量记为$\Delta\vec{\theta}_{mi}$，旋转矢量记作$\vec{\phi}_m$，这里给出在锥运动条件下优化的N子样算法。（还有在多项式角速度假设下优化的。一般认为前者模型更接近现实且更准确，但是也没法一般性证明，我可以给个例子，还是这篇重刚体定点转动的两种可解情况的求解——欧拉-潘索情形与拉格朗日-泊松情形，结果是自由定点转动下角速度是雅可比椭圆函数形式，而三角函数作为其特殊情况，显然拟合得更好些。而且探空火箭一般有两个相等的主轴惯量，文中证明了无外力矩的情况下精确满足锥运动条件。不过探空火箭绕质心的定点转动无法忽略空气阻力，这些嘛……也没办法严谨处理……只能用不完整的经验归纳）
$\vec{\phi}_m = \sum^\limits{N}_{i=1}\Delta\vec{\theta}_{mi}+(\sum^\limits{N-1}_{i=1}K_{N-i}\Delta\vec{\theta}_{mi})\times\vec{\theta}_{mN}$（A4）
其中：
一子样：无修正项
二子样：$K_1 = \frac{2}{3}$
三子样：$K_1 = \frac{27}{20},K_2 = \frac{9}{20}$
四字样：$K_1 = \frac{214}{105},K_2 = \frac{92}{105},K_3 = \frac{54}{105}$
显式解析式过于复杂，一般更高阶求解都是靠解方程得数值解的，具体见《捷联惯导算法与组合导航原理》。另外也有更优化的联系历史数据的算法，但是我觉得低成本MEMS根本没有必要上。
很容易看出，（A1）式等号右侧第一项即所有角增量求和，第二项即为圆锥误差修正。可以看出，角增量趋于零的时候，总的角增量就趋于旋转矢量。这个算法只能保证在特定的模型下，N趋于无穷时结果趋于真实的旋转矢量。现实确实最接近锥运动模型，但是毕竟不是，而且还有噪声干扰。综合考虑精度和计算速度，实际应用上一般取N=2、3、4。4子样就足够高了，对于低成本MEMS，实际上2子样的算法误差就足够小了（一般保证算法误差小于仪器误差的5%）。
注：角度使用弧度制。
这个公式实际上有不同的表述形式，由于一般读取是等时的，所以都可以化简成式（A4）。
接下来我直接给出1到4子样的公式。
一子样：$\vec{\phi}_m = \Delta\vec{\theta}_{m1}$（A5）
二子样：$\vec{\phi}_m = (\Delta\vec{\theta}_{m1}+\Delta\vec{\theta}_{m2})+\frac{2}{3}\Delta\vec{\theta}_{m1}\times\Delta\vec{\theta}_{m2}$ （A6）
三子样：$\vec{\phi}_m = (\Delta\vec{\theta}_{m1}+\Delta\vec{\theta}_{m2}+\Delta\vec{\theta}_{m3})+(\frac{9}{20}\Delta\vec{\theta}_{m1}+\frac{27}{20}\Delta\vec{\theta}_{m2})\times\Delta\vec{\theta}_{m3}$ （A7）
四子样：$\vec{\phi}_m = (\Delta\vec{\theta}_{m1}+\Delta\vec{\theta}_{m2}+\Delta\vec{\theta}_{m3}+\Delta\vec{\theta}_{m4})+(\frac{54}{105}\Delta\vec{\theta}_{m1}+\frac{92}{105}\Delta\vec{\theta}_{m2}+\frac{214}{105}\Delta\vec{\theta}_{m3})\times\Delta\vec{\theta}_{m4}$ （A8）
好了，那A部分完了，开发时可选（A2）到（A4），建议用二子样，直接翻译成代码就行。
推导和误差分析见《捷联惯导算法与组合导航原理》