1
$\begingroup$

I have multiple time series of position data (x, y, z) and orientation data expressed as Euler/Tait-Bryan angles (yaw, pitch, roll) obtained from a head tracking device. I need to obtain velocity information from this data to understand the average and maximum head speed of the users.

Since a first-order differencing amplifies high-frequency noise, as discussed in this answer, I consider applying smoothing and particularly, a Savitzky-Golay filter, to get less noisy estimates of the velocity. For the position components, e.g. $x$, I can compute the smoothed derivative easily using the Scipy's implementation:

vx = savgol_filter(x, window_length=29, polyorder=4, deriv=1, delta=dt)

where dt=0.005s is the time interval between my samples.

However, I'm not certain how I should compute the angular velocities from the Euler angles because the data contains discontinuities at the boundaries $\pm{180}$ degrees. Consider for example this portion from the yaw component:

import numpy as np
alpha = np.array([-178.06, -178.48, -178.91, -179.37, -179.83, 179.72, 179.3, 178.88, 178.46, 178.04])

As a workaround, I first computed the first-difference and set the erroneous values to 0 that are caused by the discontinuities (I didn't bother to compute the correct values for these because the effect on average/maximum velocity is negligible given that I have thousands of samples). Then I applied the SG filter to the corrected first-difference and scaled by my sampling time. See below:

alpha_diff = np.diff(alpha, 1)
alpha_diff[np.abs(alpha_diff) > 180] = 0
v_alpha = savgol_filter(alpha_diff, window_length=29, polyorder=4)
v_alpha *= 1/dt

Is this a valid approach to smooth circular/spherical data like Euler angles that continuities discontinuities at the boundaries? I'm just not sure if this workaround conceptually makes sense.

Improve this question
$\endgroup$

1 Answer 1

Reset to default
2
$\begingroup$

I'm not certain how I should compute the angular velocities from the Euler angles because the data contains discontinuities at the boundaries ±180 degrees.

Welcome to the wonderful world of trying to comb the hairy ball flat. Do a web search on "Hairy ball theorem", then again on "hairy ball theorem" and "Euler angles". You'll wish you'd never heard of 3D rotations.

Basically -- you can't get there from here, at least not in any direct manner.

My knee-jerk reaction to this would be to convert rotations to quaternions, then subtract**, then derotate. I.e. $\hat\omega \simeq 2\frac{q_n - q_{n-1}}{T_s}$ where $q_n$ is the $n^{th}$ quaternion, $T_s$ is your sampling interval, and $\hat\omega$ is your estimated rate in the head coordinate system. To get this rotation into your "home" coordinate system you'd derotate it by $q_n$: $\hat\omega_{h} = q_n^* \hat\omega = 2\frac{1- q_n^* q_{n-1}}{T_s}$. Check my math here! You should get a rotational speed, in your desired coordinate system, with $\mathcal{Re}(\hat\omega_h) \simeq 0$ -- but the signs may be wrong, or some other stupidity, because I haven't checked the math myself before typing this.

You'll need a de-jumping algorithm (I don't know the common name) because quaternions avoid the hairy ball singularity in part by representing 3D rotations as half-rotations*. This means that going from a known 3D rotation with three elements has, of necessity, a $\pm360^\circ$ ambiguity in the final result. This ambiguity doesn't matter at all in the 3D world, but it shows as a $180^\circ$ ambiguity in the quaternion itself. You need to avoid this by calculating your quaternion from your Euler angles, then multiplying the result by $-1$ as necessary to avoid big jumps.

There may be some method to get the same thing that involves less calculation. It would have to involve calculating the least-magnitude rotation matrix necessary to make one set of Euler angles rotate into a second. I'd look for it if I were you, but I don't know of it off hand.

* If you have a rotation quaternion $q$ and a vector $\vec v$, you rotate with $\vec v_{rot} = q \vec v q^*$ where $q^*$ is conjugation. Each multiplication rotates by half the angle described by $q$, so $q$ can rotate a vector by up to $720^\circ$ -- with the obvious $360^\circ$ ambiguity.

** This used to say "simply subtract", then someone asked a question, then I realized things were more complicated.

Improve this answer
$\endgroup$
11
  • 1
    $\begingroup$ de-jumping: when it's used on phase signals, we often speak of "phase un-wrapping", so "un-wrapping" does sound like a candidate. $\endgroup$
    – Marcus Müller
    May 20, 2020 at 18:11
  • $\begingroup$ The popularity of quaternions, or rather the Clifford algebra of the double cover group of rotations SU(2), for simple rotations in computer graphics eludes me to this day. Clifford algebras, and geometric algebras in general, have their use, but you can have what you are describing much simpler. Rotations in 3-dimensional real space are represented by the rigid rotation group SO(3), which has a 3-dimensional tangent space. This allows you to pick three orthogonal (anti-Hermitian) generators $X,Y,Z$ and write any rotation as $\exp(x*X+y*Y+z*Z)$ with the matrix exponential. $\endgroup$
    – Jazzmaniac
    May 20, 2020 at 18:53
  • $\begingroup$ @Jazzmaniac I'm not familiar with that math -- but your example suggests that the actual mechanics may be more computationally intensive. One of the nice things about quaternions is that while the math is wacky, the arithmetic is straightforward to code, and doesn't use more than multiplies and additions. $\endgroup$
    – TimWescott
    May 20, 2020 at 18:56
  • $\begingroup$ Consequently, you can decompose any rotation matrix into a generator by taking the matrix logarithm. Decomposing the result into the generator basis then gives you three real group coordinates $x,y,z$. The inner product required for the decomposition is the usual inner product over matrices (the trace over the adjoint product). These coordinates can be treated like a normal signal of independent components and therefore filtered and afterwards used to create the rotation matrix again using the matrix exponential. $\endgroup$
    – Jazzmaniac
    May 20, 2020 at 18:57
  • $\begingroup$ @TimWescott the computations required for the Lie-group calculations are just the same as for the Quaternions. Instead of the quaternion algebra you have matrix multiplications. And the matrix exponential and logarithm replace the trigonometry that you need for the conversion from Euler angles to quaternions and back. The mathematical framework of the Lie group is a lot more straight forward for this kind of application, however. $\endgroup$
    – Jazzmaniac
    May 20, 2020 at 18:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.