# What units does Matlab's function angvel produce?

7 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.

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
on 19 Sep 2024

Edited: James Tursa
on 19 Sep 2024

### 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

### Categories

### Products

### Community Treasure Hunt

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

Start Hunting!