What units does Matlab's function angvel produce?
16 views (last 30 days)
Show older comments
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
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
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.
Answers (2)
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
That third component should be
10*pi/180/dt %?
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]
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
We see the function returned
2/dt*sind(90/2)
instead of
90*pi/180/dt % pi/2 rad/sec
I suggest you file a bug report. If you do, please comment back here with a summary of the response from Tech Support.
0 Comments
See Also
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!