Here, I will assume that the submovements are minimum jerk submovements, that is, they minimize mean squared jerk (other options are available). These submovements are described as 5th order polynomials with 3 free parameters (which is the minimum possible) - amplitude, duration, and starting time. The equation for the position (relative to the starting point) and the velocity are given by:

where

*D*is the duration,*A*is the amplitude and*t0*the starting time, and*t0*< t < t0 + D. The Matlab function MJxy.m computes this velocity profile.
We can have more than one submovement being executed at a time, in which case we just add the change of positions of both submovements. Our goal is to take a trajectory and decompose it one or more submovements. To do this, we need to find the 3 parameters (or 4 if we are working in 2D, where we have an extra amplitude parameter) that describe each submovement (and decide how many submovements we need for a particular movement). To decide how good a reconstruction is, we need to define a cost, which is based on the difference between the observed trajectory and the reconstructed trajectory. Here will we used the squared difference (MSE). Minimising the MSE is equivalent to minimising the root mean squared error (RMSE), we use the MSE for computational reasons. Once we have the cost, we need to vary the parameters to find the minimum cost (i.e., best fit).

Rather than just minimising the difference between the observed and predicted velocity profiles, we also include a term for the tangential velocity, which we calculate for each submovement and then sum. This is to prevent the algorithm finding two submovements which start at approximately the same time, but have approximately equal and opposite amplitudes. While the resultant velocity will be approximately zero (and so won't affect the MSE), from a theoretical point of view, these seem highly unlikely. So our cost function E is:

To minimise this function in Matlab, we can use fmincon (part of the Optimization toolbox). We use this rather than fminsearch so that we can set constraints on some of the variables, in particular:

- 0 <= t0 < tf - 0.167
- 0.167 < D < 1

These enforce that submovements need to have a duration of at least 167 ms. We also constrain the amplitudes to reasonable values. After running the optimization and getting back the parameters and the error, we need to decide how many submovements we want. I usually just use a threshold (0.02), and select the minimum number of submovements that give an error of less than 0.02. Before running the optimization, we need to select a starting "guess". Here, I randomly select a starting point within the legal ranges for each parameter. However, we may end up in a local minimum. To minimize the chance of this occurring, we start from 10 different starting points, as suggested in Rohrer & Hogan (2006).

The function decompose.m performs the decomposition (calls fmincon, etc). I have provided an example (example.m) to demonstrate how to use the code (which needs a sample data file: sampledata.mat). The output of the example program is shown below:

In this case, two submovements were needed to get the reconstruction error less than 0.02. In this case, this consisted of a small submovement to the right (in green, starting at 378 ms with an amplitude of [0.088,0.125] and duration of 345 ms), followed by another submovement at 551 ms (in red, with an amplitude of [-0.432,0.346], and a duration of 689ms). Note that these numbers which describe the submovements are taken from the best fit parameters.

All files needed to run the demo can be downloaded here in a zip file.