What units does Matlab's function angvel produce?

16 views (last 30 days)
Adan
Adan on 17 Sep 2024
Edited: Adan on 5 Dec 2024 at 8:07
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 games like Skullgirls mod.
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.
  4 Comments
James Tursa
James Tursa on 18 Sep 2024
Edited: James Tursa on 18 Sep 2024
@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
Paul on 19 Sep 2024
Edited: James Tursa on 19 Sep 2024
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.

Sign in to comment.

Answers (2)

Paul
Paul on 18 Sep 2024
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
I suggest you file a bug report. If you do, please comment back here with a summary of the response from Tech Support.

ScottB
ScottB on 17 Sep 2024
Edited: ScottB on 17 Sep 2024
  1 Comment
Paul
Paul on 17 Sep 2024
Edited: Paul on 18 Sep 2024
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.

Sign in to comment.

Tags

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!