What units does Matlab's function angvel produce?

Question

Matlab introduced a function angvel since 2020. According to the documentation, the function is supposed to compute angular velocity tensor from a sequence of quaternions. The details seem puzzling. No units are mentioned in the documentation; explicit passing of timestep dt is required; check-runs on elementary rotations produce mixed results.

Ex 1. Specifying 0, 90, and 180 degree rotations about z-axis

R1 = [1 0 0; 0 1 0; 0 0 1]; 
R2 = [0 -1 0; 1 0 0; 0 0 1];
R3 = [-1 0 0; 0 -1 0; 0 0 1];

with the constant time step of 1 second, I would expect numeric angular velocities (deg/s) 0, 90, and 90, or (rad/s) 0, 1.58, 1.58, or (Pi-s/s) 0, 0.5 0.5, or (Rev/s) 0, 0.25 0.25. Yet, the following expression

angvel(quaternion(rotm2quat(cat(3,R1,R2,R3))),1,"point")

produces

         0         0         0
         0         0    1.4142
         0         0    1.4142

Why is the computed velocity sqrt(2)?..

But even further, let us specify 0, 180, and 360 degree rotations instead. I would expect the velocity to be double of the previous.

R1 = [1 0 0; 0 1 0; 0 0 1];
R2 = [-1 0 0; 0 -1 0; 0 0 1];
R3 =  [1 0 0; 0 1 0; 0 0 1];

But the output is instead

     0     0     0
     0     0     2
     0     0    -2

I understand that there might be a jump caused by the sign-convention in case of 180-degree rotation. But why is there such a weird scaling from sqrt(2) to 2?

P.S.1 changing "frame" to "point" convention does not affect the result. P.S.2 specifying angles smaller than 90 results in "correctly" scaled velocities, e.g. 45-degree rotation produces velocity 0.707

I would like to figure out how to numerically estimate angular speed (not velocities, but presumably the square root of sum of their squares) of a rotation specified as a Nx3x3 sequence of rotation matrices.

Answer 1

Although the method is not mentioned in the doc, all the testing I have done indicates that the angvel() function is a small angle approximation calculator for the angular velocity, not an actual angular velocity calculator that uses pure rotation principle calculations. Alas, the MATLAB doc should prominently mention this very important point, but it doesn't! You are comparing apples to oranges here. Your assumptions about what the actual angular velocity rates are is entirely correct, and if you want actual angular velocities you will need to write your own code for this ... but that is not what the angvel() function is calculating. The angvel() function essentially assumes a numerical Euler step between two quaternions and backs out the approximate numerical rate from that. This is not at all appropriate for large angles, but I will demonstrate what is going on in angvel():

Start with the derivative of an arbitrary quaternion q using appropriate convention:

qdot = (1/2) * q * w

Solving for the angular velocity w gives

w = 2 * q^-1 * qdot

Then suppose you had two different quaternions q1 and q2 and wanted to know the angular velocity that would result in q1 being rotated into q2 in some delta time dt. A crude numerical approximation (e.g., Euler step) ONLY APPROPRIATE FOR SMALL ANGLES to qdot could be:

qdot = deltaq / dt = (q2 - q1) / dt

So you would get the following approximate expression for the angular velocity w, letting q = q1 for this calculation:

w = 2 * q1^-1 * (q2 - q1) / dt

Now just plug in your quaternions and see the results:

>> R1 = [1 0 0; 0 1 0; 0 0 1];
>> R2 = [0 -1 0; 1 0 0; 0 0 1];
>> R3 = [-1 0 0; 0 -1 0; 0 0 1];
>> q = rotm2quat(cat(3,R1,R2,R3));
>> Q = quaternion(q) % DCM converted to quaternion type
Q = 
  3×1 quaternion array
           1 +       0i +       0j +       0k
     0.70711 +       0i +       0j + 0.70711k
           0 +       0i +       0j +       1k
>> 2*acos(q(:,1)) % The actual cumulative angles match expectations
ans =
         0
    1.5708
    3.1416
>> angvel(Q,1,"point") % Small angle approximation (inappropriate for large angles!!!)
ans = 3x3
         0         0         0
         0         0    1.4142
         0         0    1.4142
>> dt = 1;
>> 2 * conj(Q(1)) * (Q(2) - Q(1)) / dt % Manual small angle approximation calculation
ans = 
  quaternion
    -0.58579 +       0i +       0j +  1.4142k
>> 2 * conj(Q(2)) * (Q(3) - Q(2)) / dt % Manual small angle approximation calculation
ans = 
  quaternion
    -0.58579 +       0i +       0j +  1.4142k

Then just pick off the vector part for the answer (the last three elements since MATLAB is scalar-first) and lo and behold it matches angvel(). Of course, this is a very lousy result because we used quaternions representing large angular movement in a small angle approximation formula. Bottom line, don't use angvel() for these calculations.

If you want angular velocity, just use the usual form to get the relative rotational quaternion between two quaternions and work from that (with appropriate divide by 0 checks not shown). E.g.:

>> qrot = compact(conj(Q(1)) * Q(2));
>> e = qrot(2:4); % Vector part
>> norme = norm(e);
>> ang = 2 * atan2(norme,qrot(1));
>> w = (ang / dt) * (e / norme)
w = 1x3
     0     0   1.5708

If you just want angular speed then ang/dt will suffice. E.g., here you could use acos():

>> 2 * acos(min(max(qrot(1),-1),1)) / dt
  ans = 1.5708

I am not a huge fan of using acos() without argument protection against it being slightly larger than 1 in magnitude due to numerical effects, hence the extra code. Typically I just use the atan2() formula which avoids this concern.

See this MATLAB Answers link for a related post: https://www.mathworks.com/matlabcentral/answers/564368-whats-difference-beetween-angvel-and-rotvec?s_tid=srchtitle

Another related link, which points out that angvel( ) uses a small angle approximation that is different from the method used above: https://www.mathworks.com/matlabcentral/answers/2153395-what-units-does-matlab-s-function-angvel-produce?s_tid=srchtitle

Answer 2

The angvel function calculates angular velocity in radians per second (rad/s) however, the way angular velocity is derived from quaternions or rotation matrices can lead to values that are not intuitive, especially with how rotation transformations work.

The reason you're seeing values like sqrt(2) and 2 instead of the expected values can be due to reasons such as

• Angular velocity is derived from the change in orientation over time; when you pass a sequence of quaternions to angvel, it computes the angular velocity at each step in time.
• Quaternions represent rotations in non-linear ways. Computing angvel can include issues with scaling when rotations approach certain limits. This is why quaternions interpolate rotations and handle orientation.
• As already mentioned, scaling anomalies when you're observing values like sqrt(2) is related to why angular velocity vectors are not purely aligned with the axis of rotation. The calculation yields a vector where the magnitude represents the angular speed, but the vector’s components can result in seemingly unexpected numbers.

This is why a 90 degree rotation might have an angular velocity magnitude that is sqrt(2), not exactly 1 rad/s as one should expect.

• A Rotation of 180 degrees is a special case because it can represent the same rotation using different signs (+/-) and can lead to flips in the sign of angvel.

In your example, the angular velocities you might expect are based on simple geometric reasoning. However, the way quaternion computes results in angular velocities includes unintuitive factors, since how quaternions encode rotational changes which MATLAB computes the derivative over time. The negative sign in your 0, 180 and 360 degree rotations represents a flip in orientation.

If you're interested in calculating the angular speed rather than the full angular velocity vector you can use:

% Compute angular speed from angular velocity vectors
ang_speed = sqrt(sum(omega.^2, 2));  % norm of each angular velocity vector

If omega is the angular velocity matrix returned by angvel, it will give you the scalar angular speed (rad/s).

Useful Learning Materials:

Rotations, Orientation, and Quaternions - MATLAB

rotm2quat - MATLAB

Quaternions and Rotation Sequences - Jack B. Kuipers

Math Stack Exchange - STACK

Quaternions and 3d rotation, explained interactively - 3Blue1Brown