*2.3. Algorithm*

The algorithm flowchart and sample data are presented in Figure 2. The input data consisted of the three acceleration components measured by a single sensor (Figure 2A). The analysis was applied on windowed finite time series of the acceleration module (Figure 2B) and of single acceleration components.

**Figure 2.** Method data flow: (**A**) raw recording from triaxial accelerometer; (**B**) window of acceleration module; (**C**) autocorrelation function computed on standardized windowed data (the cross marks the autocorrelation peak which corresponds to the window fundamental period); (**D**) the regularity time profile (the dotted line represents the window average value considered as a regularity index (RI)); (**E**) the fundamental period time profile (the dotted line represents the window average value considered as the period index (PI)).

The duration (*n* seconds, or alternatively *N* samples) of the signal window was fixed, on each trial, to three times the fundamental period preliminarily quantified on the entire one-minute recording (this duration ranged from about 2 s for the fastest running to about 4 s for the slower walking). Such a choice was sufficient for a robust computation of the autocorrelation function and avoided regularity variations across a higher number of strides may have affected the index estimation [31].

The windowed data were first standardized (removed average and then divided by standard deviation) and then an autocorrelation analysis was applied, which is equivalent to an autocovariance analysis since a unitary covariance is obtained when correlating the signal with itself at zero lag. Thus, a unitary value of the autocorrelation function reflects a perfect regularity, while a null value marks the absence of any cyclic component, i.e., the absence of any regularity.

The autocorrelation sequence for a given data window (Figure 2C), which is a function of the sample lag *m*, was implemented as follows:

$$\mathcal{R}\_{\text{unbiased}}[m] = \frac{1}{N - m} \sum\_{n=0}^{N-m-1} y[n+m]y[n]. \tag{3}$$

where *R* is the autocorrelation function in the unbiased form, *y* is the analyzed variable, *N* is the number of data samples in the considered time window, and *m* is the number of samples quantifying the lag time of the autocorrelation function.

Notably, the previous equation includes a coefficient corrected for *m*, to compensate for the bias resulting from the *m* number of samples associated with the considered time lag (in fact, the larger the lag, the fewer samples support the computation) [29]. Negative lag values were not considered since they did not provide additional information, as autocorrelation is an even function.

Window periodicity, i.e., the fundamental period duration of the window data, corresponds to the time of the autocorrelation peak (Figure 2C). Window regularity corresponds to the peak value; the closer the input data are to a perfect periodic behavior, the closer the window regularity is to a unitary value.

The previous analysis was then repeated for the subsequent data window which was time-shifted from the previous window by 0.1 s. This implies that two subsequent data windows, for locomotion recordings, may overlap for more than 95% of their data. In order to speed up and increase algorithm robustness, when considering the next data window, the previous identification of the autocorrelation peak was considered as the starting seed in the search of the actual peak of the autocorrelation function, thus substantially avoiding the identification of peaks related to higher harmonics of the fundamental period.

The regularity index (RI) and period index (PI) were computed as the average values of the sequences of window regularity (Figure 2D) and periodicity (Figure 2E), respectively.

All analyses were implemented in Matlab (version R2017b, The MathWorks Inc, USA).
