IMU Attitude Estimation from Angular Velocity rotates by roughly half the expected value. Why?

Question

I'm working on getting the orientation of an IMU in space, following the instructions from https://ahrs.readthedocs.io/en/latest/filters/madgwick.html. I can already determine the orientation on the basis of the magnetometer (this suffers badly from very fast movements though, so I want to go the gyro route).

The rotation basically does what it is supposed to do (I use the rotation matrix from the quaternion to render a colored box in real time), but at a much lower rate. In essence, the rotation that results from my code is approximately (but not quite) 1/2 as fast as the one in the real world.

Simply multiplying the angle with 2

• does not quite solve the problem (there still results a certain lag)
• does not help me in understanding the problem

Is what I am doing conceptually sound? Is there a major issue I just don't see? As I need this stuff for my master's thesis (in exercise science), any help is greatly appreciated.

UPDATE: I have implemented the advice by @JamesTursa and switched to Quaternion.CreateFromAxisAngle() instead of Quaternion.CreateFromYawPitchRoll(). Unfortunately, this did not solve the issue. I have also implemented a unit test for the Quaternion class to check whether or not it behaves as expected. It does, so I doubt there's an error with that.

After reading through Efficient quaternion angular velocity, I have also tried using the half angle, but this does exactly what I supposed and simply slows the rotation further.

So, I am still stuck with a rotation that basically does approximately (but not quite) half of what it is supposed to. Again, any help is eppreciated.

/// <summary>
        /// Numerically integrate angular velocities in x,y,z to get orientation in x,y,z
        /// Does not work. Rotation happens at around (but not quite) half the rate as in the real world.
        /// </summary>
        /// <param name="dt">time in ms since last invocation</param>
        private void UpdateAngularVelocities(double dt)
        {
            // rotation threshold (°/sec) to minimize noise
            // set to 0 for the time being
            double threshold = 0.0; 

            // iterative quaternion multiplication based on madgwick filter
            // https://ahrs.readthedocs.io/en/latest/filters/madgwick.html
            // keep an AngularOrientation Vector for illustration purposes

            // apply rotation velocity threshold
            // ie, discard any rotation that is slower than the defined threshold
            var wtx = (Math.Abs(w[0]) > threshold) ? w[0] : 0.0;
            var wty = (Math.Abs(w[1]) > threshold) ? w[1] : 0.0;
            var wtz = (Math.Abs(w[2]) > threshold) ? w[2] : 0.0;

            // rotation velocities, stored as a Vector for easier use
            // theses will be used to onstruct the Quaternion
            Vector3 deltaTRotation = new Vector3((float)ConvertToRadians(wtx), (float)ConvertToRadians(wty), (float)ConvertToRadians(wtz));
            Vector3 deltaAxis = deltaTRotation;

            // scale rotation with time
            // dt is supplied in ms, we need s to get from °/s to °
            // use half angle because of
            // https://stackoverflow.com/questions/24197182/efficient-quaternion-angular-velocity/24201879#24201879
            deltaTRotation.X *= (float)(dt / 1000.0) * 0.5f;
            deltaTRotation.Y *= (float)(dt / 1000.0) * 0.5f;
            deltaTRotation.Z *= (float)(dt / 1000.0) * 0.5f;


            // build quaternion from scaled rotations
            // use Quaternion.CreateFromAxisAngle() to get simultaneous angles
            // https://stackoverflow.com/questions/72578478/imu-orientation-from-gyro
            // Quaternion qDelta = Quaternion.CreateFromYawPitchRoll(deltaTRotation.Y, deltaTRotation.X, deltaTRotation.Z);
            Quaternion qDelta = Quaternion.CreateFromAxisAngle(Vector3.Normalize(deltaAxis), deltaTRotation.Length());

            // deal with NaN, in case no Rotation takes place and hence, ||qCurrent|| = 0 which makes Normalize throw NaN
            if (qDelta.Length() > 0.0)
            {
                // combine new rotation with previous rotation
                qRotation = qDelta * qRotation;
                //AdditiveOrientation += 0.5f * Vector4.Normalize(new Vector4(deltaTRotation.X, deltaTRotation.Y, deltaTRotation.Z, deltaTRotation.W)) * qDelta.Length();
            }

            // https://ahrs.readthedocs.io/en/latest/filters/madgwick.html suggests an addition of the old state and the new state.
            // this results in the quaternion not being a unit quaternion / versor anymore. Ditch the addition and stick with simple multiplication
            //qRotation = qRotation + Quaternion.Multiply(qCurrent * qRotation, (float)(0.5));// * dt / 1000.0)); 
            //qRotation = Quaternion.Normalize(qRotation);


            // rotate magnetic orientation for illustration purposes
            // always start rotation from (arbitrary) [1,0,0] and apply most recent quaternion
            // does not work yet
            // rotates roughly half the angle it should
            AngularOrientation = Vector4.Transform(new Vector4(1.0f,0.0f,0.0f,1.0f), qRotation);

            // save angular velocities for later use
            OmegaXData.Add(deltaTRotation.X);
            OmegaYData.Add(deltaTRotation.Y);
            OmegaZData.Add(deltaTRotation.Z);
        }