Numerical differentiation on circular data
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.