Calculate a rotational speed for an object spinning in 3 axes

Question

I've got an object spinning in 3 axes and I'm tracking it with a motion capture system. For each timepoint, depending on how I export the data, I either get 4 columns of data for the quaternion rotation or 3 columns for the rotation around the global XYZ axes.

For the XYZ data, I can calculate the rotational velocity in each axis but I'm not sure how to combine this data into a speed. I've tried using the sqrt(dx^2 + dy^2 + dz^2) formula you would use for combining linear velocities into a speed but when I picture the problem this doesn't seem like the right approach.

For the quaternion data I used the pyquaternion library in python to calculate the angle between consecutive data points then calculated the rate of change in the angle. For the angle calculation I used the approach given here; I've copied the approach below in case the link breaks. This seems like a better approach than my attempt with the XYZ data but I'd never heard of quaternions until today so I'm doubting my approach.

Assuming these represent attitude rotations from one coordinate frame to another, if you are simply asking what is the minimum rotation to take you from one quaternion to the other, you simply multiply one quaternion by the conjugate of the other and then pick off the rotation angle of the resulting quaternion.
But we really need to know what these quaternions represent, and what angle you are trying to recover, before we know what you want.
E.g., suppose x and y represent ECI->BODY rotation quaternions, and you want to know the minimum rotation angle that would take you from the x BODY position to the y BODY position. Then you could do this:
>> x = [ 0.968, 0.008, -0.008, 0.252]; x = x/norm(x); % ECI->BODY1
>> y = [ 0.382, 0.605,  0.413, 0.563]; y = y/norm(y); % ECI->BODY2
>> z = quatmultiply(quatconj(x),y) % BODY1->BODY2
>> z = 0.5132    0.6911    0.2549    0.4405
>> a = 2*acosd(z(4)) % min angle rotation from BODY1 to BODY2
>> a = 127.7227
But, again, these calculations are dependent on how I have the quaternions defined. Your specific case may be different. The scalar part is the 4th element.

Answer 1

I believe you are right that the linear velocity formula is incorrect - AIUI, the correct method would be to combine the three axis-aligned rotations over your timestep, calculate the net rotation (which may have an unaligned axis) and then retrieve the rotation angle of this single rotation. However, this sounds exactly the same concept as what you have with the quarternions, and I expect that the matrix arithmetic involved in doing the same with three scalars would be less convienient than what you already have.

Answer 2

You have an object with orientation quaternion $\hat{q} = \pmatrix{ \vec{q}_v \\ q_s } $ where $\vec{q}_v$ is the vector part of the quaternion and $q_s$ the scalar part.

If this object is rotating with $\vec{\omega} = \pmatrix{\omega_x \\ \omega_y \\ \omega_z}$ then the orientation is going to have the following derivative.

$$ \tfrac{\rm d}{{\rm d}t} \hat{q} = \tfrac{1}{2} \pmatrix{\vec{\omega} \\ 0} \otimes \pmatrix{ \vec{q}_v \\ q_s} $$ where $\otimes$ is the quaternion product. The above is expanded into the following

$$ \pmatrix{ \dot{\vec{q}}_v \\ \dot{q}_s } = \tfrac{1}{2} \pmatrix{q_s \vec{\omega} + \vec{\omega}\times \vec{q}_v \\ -\vec{q}_v \cdot \vec{\omega}} $$

Now to compute the reverse, when an orientation $\hat{q}$ and its rate $\dot{\hat{q}}$ we are given and we calculate the rotational velocity

$$ \vec{\omega} = 2 \left( q_s \dot{\vec{q}}_v - \vec{q}_v \dot{q}_s + \vec{q}_v \times \dot{\vec{q}}_v \right) $$

and rotational speed

$$ \omega = \| \vec{\omega} \| = \sqrt{ \omega_x ^2 + \omega_y^2 + \omega_z^2} $$

The above is the instantaneous rotational speed, given a quaternion rate.

But in your case, you have two quaternions and I assume a time step $h$ and you might want the average rotational speed between the two orientations.

The interpolation between $\hat{q}_1$ and $\hat{q}_2$ is done mathematically with the slerp function.

$$ \hat{q}(t) = \left( \frac{ \sin\left( (1-\tfrac{t}{h}) \theta \right)}{ \sin \left( \theta \right) } \right) \hat{q}_1 + \left( \frac{ \sin\left( \tfrac{t}{h} \theta \right)}{ \sin \left( \theta \right) } \right) \hat{q}_2$$

where the relative rotation angle is $\theta$ is the angle encoded in the relative rotation $\hat{q}_h = \hat{q}_1^{-1} \otimes \hat{q}_2$.

The angle is calculated from $\hat{q}_h = \pmatrix{ \vec{q}_h \\ q_h }$ as $$ \theta = 2 \tan^{-1} \left( \frac{ \| \vec{q}_h \| }{q_h} \right)$$ $$ \hat{\rm axis} = \frac{ \vec{q}_h }{ \sin( \theta/2) }$$

If that is the case, then the average speed is just $\boxed{\dot{\theta} = \frac{\theta}{h}}$.