What units does Matlab's function angvel produce?

6 次查看（过去 30 天）
Adan
Adan 2024-9-17
编辑： James Tursa 2024-9-19
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.
  4 个评论
显示 2更早的评论
James Tursa
James Tursa 2024-9-18
编辑：James Tursa 2024-9-18
@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 2024-9-19
编辑：James Tursa 2024-9-19
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.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

Paul
Paul 2024-9-18
Ran in:
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.
ScottB
ScottB 2024-9-17
编辑：ScottB 2024-9-17
  1 个评论
Paul
Paul 2024-9-17
编辑：Paul 2024-9-18
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Axes Transformations 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by