What units does Matlab's function angvel produce?

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")

produce

         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.

[James Tursa]

Cross posted here (with the small angle approximation answer):

https://stackoverflow.com/questions/78984948/what-units-does-matlabs-function-angvel-produce

Also this was previously discussed in this 2020 thread (also with the small angle approximation answer):

https://www.mathworks.com/matlabcentral/answers/564368-whats-difference-beetween-angvel-and-rotvec?s_tid=srchtitle

[Paul]

@James Tursa

After reading both of those links, I just want to point out that angvel does not do an arithmetic subtraction of quaternion elements.

It just computes the quaternion from one frame to the next as (assuming frame rotations)

qdiff = q2*conj(q1)

and then computes the "angular velocity" from qdiff as shown below.

Seems like it would have been very simple to compute angvel correctly under what appears to be the underlying assumption that the angular velocity is a constant vector over the interval.

The fact that angvel doesn't even explain what it means by "angular velocity" is not good (to state it nicely).

[James Tursa]

@Paul So a different small angle approximation, I guess (possibly mathematically equivalent?). I don't have a MATLAB version with angvel( ) so couldn't look at the implementation, but based on my testing it seemed pretty obvious something like that was going on. Nevertheless, I don't know when I would ever want to use such a function when exact rotational methods are available and not difficult to implement.

[Paul]

Fair enough on the small angle approximation, but I don't think exactly mathematically equivalent. I would never want to use such a function and am seriously wondering why the developers did what they did and in what context they thought it would be more useful than an exact method, which I don't think is that much harder to implement than what they did anyway.

Answer 1

[Paul]

angvel is fundamentally flawed.

Here's the example (truncated) from the doc

eulerAngles = [(0:10:40).',zeros(numel(0:10:40),2)];
q = quaternion(eulerAngles,'eulerd','ZYX','frame');
dt = 1;
format long e
av = angvel(q,dt,'frame') % units in rad/s

av = 5×3

0 0 0 0 0 1.743114854953163e-01 0 0 1.743114854953164e-01 0 0 1.743114854953163e-01 0 0 1.743114854953164e-01

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

That third component should be

10*pi/180/dt %?

ans = 
     1.745329251994329e-01

Inspecting the code for angvel shows that it intends for the direction of AV(k,:) to be along the axis of rotation for the rotation from qi(k-1,:) to qi(k,:) (resolved in the k-1 frame) and the magnitude of AV(k,:) is equal to the angle of that rotation divided by dt.

However the code is incorrect. In terms of the angle (phi) and axis of rotation (u, a unit vector), the quaternion (expressed as a vector) is:

q = [cos(phi/2), sin(phi/2)*u]

What angvel does is effectively (to within a sign): av = 2/dt*q(2:4).

In the example above, that would be

2/dt*sind(10/2)*[0 0 1]

ans = 1×3

0 0 1.743114854953163e-01

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

which is exactly what was produced by angvel (to rounding in the last digit).

So the code is relying on the small angle approximation:

2*sin(phi/2) almost= phi.

Don't know why it needs any approximation, which of course gets worse as the angle gets bigger as in the code in the question

eulerAngles = [(0:90:90).',zeros(numel(0:90:90),2)];
q = quaternion(eulerAngles,'eulerd','ZYX','frame');
dt = 1;
av = angvel(q,dt,'frame') % units in rad/s

av = 2×3

0 0 0 0 0 1.414213562373095e+00

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

We see the function returned

2/dt*sind(90/2)

ans = 
     1.414213562373095e+00

instead of

90*pi/180/dt  % pi/2 rad/sec

ans = 
     1.570796326794897e+00

@Adan

I suggest you file a bug report. If you do, please comment back here with a summary of the response from Tech Support.

Answer 2

[ScottB]

Documentation asserts that the results are in radians per second.

https://www.mathworks.com/help/robotics/ref/quaternion.angvel.html?searchHighlight=angvel&s_tid=srchtitle_support_results_1_angvel

[Paul]

The units of the output are only shown in the Examples section. It wouild be better if the units were stated in the Output Arguments section in the description of AV.