Estimating Movement Smoothness From Inertial Measurement Units

Inertial measurement units (IMUs) are increasingly used to estimate movement quality and quantity to the infer the nature of motor behavior. The current literature contains several attempts to estimate movement smoothness using data from IMUs, many of which assume that the translational and rotational kinematics measured by IMUs can be directly used with the smoothness measures spectral arc length (SPARC) and log dimensionless jerk (LDLJ-V). However, there has been no investigation of the validity of these approaches. In this paper, we systematically evaluate the use of these measures on the kinematics measured by IMUs. We show that: (a) SPARC and LDLJ-V are valid measures of smoothness only when used with velocity; (b) SPARC and LDLJ-V applied on translational velocity reconstructed from IMU is highly error prone due to drift caused by integration of reconstruction errors; (c) SPARC can be applied directly on rotational velocities measured by a gyroscope, but LDLJ-V can be error prone. For discrete translational movements, we propose a modified version of the LDLJ-V measure, which can be applied to acceleration data (LDLJ-A). We evaluate the performance of these measures using simulated and experimental data. We demonstrate that the accuracy of LDLJ-A depends on the time profile of IMU orientation reconstruction error. Finally, we provide recommendations for how to appropriately apply these measures in practice under different scenarios, and highlight various factors to be aware of when performing smoothness analysis using IMU data.

In this supplementary document, all equations and figures are referenced with Roman numerals. Any reference using the Arabic numerals refers to contect in the main document.
Appendix A Kinematics of a segment moving relative to another Translational motion -The linear velocity of S relative to B (Fig. 2 in the main text) is given by: Rotational motion -The angular velocity of S with respect to B is given by: Combining Eqs. VI and 4, we can obtain B ω BS from the IMU's gyroscope signal S w Appendix B

Analysis of the properties of LDLJ-A
In this paper, we defined a modified version of LDLJ-V derived from movement acceleration (LDLJ-A or λ a L , Eq. 16). λ a L and λ v L (Eq. 2) differ only in their scaling factor (t 2 −t 1 ) , respectively, since Thus, This relationship is non-linear due to the 2 ln a peak v peak factor. Consistency of LDLJ-A. The consistency of the LDLJ-V and LDLJ-A measures were evaluated using simulated 1D movements, which were generated as a sum of a varying number of submovements with varying inter-submovement intervals. The movement velocity and acceleration profiles of these 1D movements were of the form, where v (·) and a (·) are the velocity and acceleration of the movement, respectively, ψ (t) = e −25(t−0.5) 2 is a smooth submovement, δT i ∈ 0.001i i ∈ N and 0 ≤ i ≤ 2000 are the intersubmovement intervals, and N s ∈ {2, 4, 8} is the number of submovements. The smoothness values of this set of movements was evaluated using the LDLJ-V and LDLJ-A applied to the velocity and acceleration data respectively.  Figure II: Scatter plot of λ v L and λ a L for a set of movements with varying number of submovements. The scatter plot shows that there is a fair correlation between the smoothness values reported by the two movements on the same movement data, but there is no one-to-one relationship.
λ v L and λ a L have roughly similar responses to changes in δT i , and show approximately monotonic responses to the change in submovement characteristics (Fig. I). However, in general, for the same movement, smoothness values λ v L and λ a L will not be the same. LDLJ-A and LDLJ-V cannot be interchanged. Although λ v L and λ a L do not provide the same smoothness values, they could still be used interchangeably if they would report the same relative smoothness between any two movements. Mathematically, if there is be an isomorphism between the measures, that is, given movements M 1 and M 2 with velocity profiles v 1 and v 2 , respectively, then To check if such an isomorphism exists between λ v L and λ a L , we simulated a set of discrete, minimum jerk, 3D movements through via-points. The generation of these movements was posed as a constrained optimization problem and was solved using the CVXPY solver in Python, similar to the approach used in [2]: where x is a stacked vector of positions corresponding to the movement trajectory, D is the matrix performing the triple derivative to obtain movement jerk from x, and A and b represent the spatio-temporal constraints corresponding to initial, final, and via-points position, and the zero initial, and final velocity and acceleration constraints. The solution x m obtained using this procedure would possess the minimum squared jerk Dx 2 2 for the given spatial and temporal constraints.
We simulated movements with initial position x i = 0 0 0 , final position x f = 0 1 0 , and duration T = 1s. A set of N via-point spatio-temporal constraints {x v n , t n } N n=1 were randomly generated, where the n th via-point x v n is traversed at time t n . A set X m of 251 minimum jerk trajectories were generated, corresponding to: (a) 25 different spatio-temporal constraints for each of 10 different values of N ; and (b) one movement without any via-points. The sampling period was δt s = 0.001s, resulting in 1000 sample points per movement.
The smoothness of the movements in X m estimated using λ v L and λ a L were fairly well correlated (Fig. II), but not isomorphic. This implies that the two measures have slightly different structures, and one cannot be exactly obtained from the other.

Appendix C SPARC on translational velocity and acceleration data
The SPARC measure was also applied to the simulated data (from Appendix B) to evaluate the measure's consistency when estimated from movement velocity or acceleration (Fig. III). Similar to [1], there was a reasonably consistent response when using velocity data ( Fig. III(a)). On the other hand, when movement acceleration is used (Fig. III(b)), it does not result in a consistent response (i.e., a monotonic response to change in inter-submovement interval). This difference in behavior is due to the definition of SPARC: it was inspired by the idea of movement planning and execution by combining submovements. SPARC measures the level of dispersion of a signal in time (e.g., relative periods of movement and arrest). This measure of temporal dispersion acts as a reasonable measure of smoothness only when it is applied on the velocity, because velocity dispersion is strongly related to movement intermittency. Temporal dispersion in all other kinematic variables do not provide a unique measure of intermittency.

Appendix D
Calculating magnitude of linear jerk from IMU signals without orientation estimation The derivative of the accelerometer signal S a is given by (from Eqs. 3 and 4): where the matrix [ S ω SO ] × is skew-symmetric. Because the operator 2 is rotation invariant, the magnitude of the linear jerk is given by It is important to note that, in practice, the term Sȧ is computed from numerical differentiation. Subtracting the correction term S a × S w, which is derived from the raw IMU signals, is affected by misalignment in time and different sampling of peaks. In simulated data, we found that, in some cases, this can lead to erroneous estimation of the translational jerk. Further investigation is needed before general recommendations can be provided for this approach.

Calculating magnitude of angular jerk from IMU signals without rotation estimation
Similar to translational jerk, the absolute angular jerk can be computed directly from IMU data. The second derivative of the gyroscope signal S w is given by: Thus, the magnitude of the absolute angular jerk can be calculated as: In practice, numerical differentiation would be necessary to obtain Sẅ and Sẇ , which could result in to erroneous estimates. Indeed, the formulation in Eq. XII was found to be sensitive to practical implementation issues (e.g. choice of numerical differentiation methods, sampling frequency, etc.). Therefore, further investigation is needed before general recommendations can be made for this formulation to be used in movement smoothness analyses.

Errors in orientation reconstruction δR result in errors in estimated movement acceleration
Ox OS = δR Oẍ OS + (δR − I) O g. These errors are reflected in the jerk integralÎ j and the peak accelerationâ peak , thus affecting LDLJ-A. In summary: • When δR is constant over time, then there is zero error in the smoothness estimate (Eqs. XIV and XVII).
• When δR = 0, movements with larger acceleration/jerk are likely to have smaller errors in their smoothness estimates (Eqs. XV-XVI and Eqs. XIX-XX).
• Small values ofδ R result in small errors in the jerk integral term.
• When ∆ δR (Eq. 20) is small, errors in peak acceleration are likely to be small.
In the rest of this section, we provide the mathematical details supporting these statements. Error in the jerk integral term. Consider a movement M with true acceleration Oẍ OS and reconstructed acceleration Ox OS . From Eq. 14, the derivative of the reconstructed acceleration is given by O. ..
Thus, the difference between the true and the reconstructed squared jerk magnitudes consists of: : the rate of change of reconstructed acceleration due to time varying rotational error δR.
Inner product between the rotational error δR O ...
x OS and the time-varying rotational errorδ R ( Oẍ OS + O g). Let δJ be the relative error between the mean squared magnitudes of the true and reconstructed jerk magnitudes. Using Eq. XIII: Thus, when the reconstruction error is fixed and does not vary over time, i.e.δ R = 0, then δJ = 0. When the orientation reconstruction error tends towards zero, thenδ R → 0 (assuming no discontinuities in δR), and it can be shown that δJ tends to 0, i.e. lim δR→I δJ = 0. However, when there is time-varying reconstruction error δR = I andδ R = 0, then δJ will be non-zero and its magnitude will be determined by acceleration and jerk magnitudes of the movement: • When Oẍ OS and O ...
x OS are small, δJ is dominated by gravity: δJ can be quite large depending on O ...
x OS 2 2 . Thus, incorrectly reconstructed, slow movements can have large errors in δJ and thus in the LDLJ-A smoothness estimates.
• Movements with larger acceleration/jerk reduce the effect of gravity. Let δJ in Eq. XIV be the error for a movement M measured using an IMU. If the overall amplitude of movement M was larger (no change in duration), such that the acceleration and jerk of this larger movement are A Oẍ OS and A O ...
x OS , respectively, where A > 1 is the scale of the movement. The relative error in the jerk magnitude δJ A will then be, • Shorter movements will have lower error in the jerk integral term. For two movements with the same average relative error in the jerk magnitude, the shorter movement will have an overall lower error because of shorter period of integration.
Error in the peak acceleration term. The value of a peak is computed from the magnitude of the mean subtracted reconstructed acceleration (Eq. 16). The ratio between the magnitude of the reconstructed and the true acceleration is: (XVII) Unlike δJ, δa does not depend onδ R and only on δR. δa will be closed to one as long as the dynamic range of the reconstruction error signal ∆R is close to 0, which implies that δ δR at different times points are similar to each other, i.e.
When ∆ δR ≈= 0 =⇒ δR ≈ R, where R is the fixed "average" orientation reconstruction error, then The contribution of the gravity terms to δa are reduced by a factor of A than that of Eq. XVII. A similar effect can be shown for faster movements.
Demonstration of the effects of reconstruction on LDLJ-A. We demonstrate the effect discussed in this section through a simple simulation. We simulated a minimum jerk movement of 30cm amplitude and 0.5s duration:  where • A β controls the dynamic range of the rotational error ∆ δR .
• A β and f β determined the rate of change of rotational errorδ R (t).
The magnitude of the reconstructed acceleration is affected only by A β and not f β . On the other hand, the size of the errors in the jerk magnitude is related to the rate of change of orientation reconstruction error (Fig. IV(c)). Note thatβ (t) has the same peak value when A β = 5°, f β = 16Hz and A β = 10°, f β = 8Hz. Finally, the error in the jerk magnitude increases with the rate of the change of orientation reconstruction error (Fig. IV(d)).