Numerical differentiation on circular data 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) &gt; 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.
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.