The documentation is a bit incorrect or not detailed enough. Here’s what is going on:
The gyroscope signal is integrated to create a “delta angle,” i.e. a change in angle. The delta angle bias is tracked as a state and subtracted from the delta angle. The delta angle bias is a result of the gyroscope bias.
This delta angle minus bias quantity is treated as a rotation vector and converted to a delta quaternion. That is the key piece that is missing from these equations. The current orientation estimate (states 1:4, representing a quaternion) is multiplied by this delta quaternion for a new orientation quaternion estimate.
Also not shown in the equations is that this new orientation quaternion estimate is then normalized and the real part is forced to be positive.
It’s hard to capture all this in simple equations, but I’ll make sure the documentation is updated to reflect these details.