Ambulatory Assessment of Instantaneous Velocity during Walking Using Inertial Sensor Measurements

A novel approach for estimating the instantaneous velocity of the pelvis during walking was developed based on Inertial Measurement Units (IMUs). The instantaneous velocity was modeled by the sum of a cyclical component, decomposed in the Medio-Lateral (ML), VerTical (VT) and Antero-Posterior (AP) directions, and the Average Progression Velocity (APV) over each gait cycle. The proposed method required the availability of two IMUs, attached to the pelvis and one shank. Gait cycles were identified from the shank angular velocity; for each cycle, the Fourier series coefficients of the pelvis and shank acceleration signals were computed. The cyclical component was estimated by Fourier-based time-integration of the pelvis acceleration. A Bayesian Linear Regression (BLR) with Automatic Relevance Determination (ARD) predicted the APV from the stride time, the stance duration, and the Fourier series coefficients of the shank acceleration. Healthy subjects performed tasks of Treadmill Walking (TW) and Overground Walking (OW), and an optical motion capture system (OMCS) was used as reference for algorithm performance assessment. The widths of the limits of agreements (±1.96 standard deviation) were computed between the proposed method and the reference OMCS, yielding, for the cyclical component in the different directions: ML: ±0.07 m/s (±0.10 m/s); VT: ±0.03 m/s (±0.05 m/s); AP: ±0.06 m/s (±0.10 m/s), in TW (OW) conditions. The ARD-BLR achieved an APV root mean square error of 0.06 m/s (0.07 m/s) in the same conditions.


Introduction
An important goal of locomotion is the displacement of the Body Center of Mass (BCOM), whose location in humans is somewhere within the pelvis. During walking, the pelvis moves in the three dimensional (3D) space, showing, relative to the mean forward velocity, speed fluctuations in the direction of progression [1].
One approach to investigate the behavior of the BCOM during gait consists of attaching a marker to the sacrum bone and tracking its motion using an Optical Motion Capture System (OMCS)-the sacral marker method [2]. In spite of being less accurate, this method is appealing for the ease of implementation, as compared with other approaches, such as the segmental analysis and the force plate methods [3]. However, a feature common to all these methods is that their application is restricted to gait labs [4]. Additionally, unless a treadmill is used to reproduce walkway conditions, the capture of a large sample of consecutive gait cycles is precluded [4]. Unobtrusive approaches to estimating the BCOM motion are thus needed, e.g., in the ambulatory assessment of external work during gait [5,6].
In this paper, we report and evaluate a novel inertial sensor-based algorithm that estimates the instantaneous velocity of an Inertial Measurement Unit (IMU) attached to the pelvis (pelvis IMU) during gait. For cyclical motions such as bipedal locomotion, the instantaneous velocity can be modeled by the sum of the 3D cyclical component of velocity and the velocity of forward motion averaged detecting the underlying periodicity of gait data and performs data-driven windowing without loss of time resolution [17].

Datasets
In this paper, we revisited two datasets available from our previous research. The Treadmill Walking (TW) dataset collects IMU and OMCS data from healthy subjects that performed trials of treadmill walking at preset speeds [10] ( Table 1). The Overground Walking (OW) dataset collects IMU and OMCS data from healthy subjects that walked along a "figure-of-eight" pathway at their preferred speed [24] (Table 2). rest period of 5 s with the participants standing still in their upright posture before the start of each trial * Although available, data from these sensors were not used in this paper. Table 2. Setup and experimental protocol for the OW dataset [24].
In both datasets, the experimental setup required two IMUs that were attached to the lumbar spine (L5 level) and to the shank (above the right malleolus). Clusters of four retro-reflective markers were mounted on a small plastic support rigidly attached to each IMU. In the case of the TW dataset, three additional retro-reflective markers were available in correspondence of the heel, the first metatarsal head and the fifth distal phalanx of the instrumented foot. Standard precautions were taken for aligning the IMU axes approximately along the anatomical directions and securing the IMUs in place firmly [10,24].

IMU and OMCS Data Pre-Processing
The procedures of IMU sensor calibration, frame registration, synchronization and conditioning of IMU and OMCS data streams were the same in the construction of both datasets [24]. The Cartesian coordinate systems fixed with the OMCS and an IMU were denoted, respectively, as the Global Earth-fixed Frame (with one axis aligned with Gravity) (GGF) and the Unit Local Frame (ULF). The orientation of the IMU axes were X: Antero-Posterior (AP) and positive forward; Y: Medio-Lateral (ML) and positive to the right; and Z: VerTical (VT) aligned with the direction of gravity and positive downwards, Figure 1.  Virtual markers were created in correspondence of the center of the tri-axial accelerometer chips, where the origins of the ULFs were located, by using the positions of the four retro-reflective markers associated to each IMU (Tables 1 and 2).

Mathematical Processing
During steady-state level locomotion it is reasonable to assume that limb, or trunk displacement data can be modeled as quasi-periodic functions of time. Body segment displacement data have been obtained that can be accurately described through Fourier series that contain up to M = 6 significant harmonics [4,7]: where the DC component accounts for the effect of forward motion and the period is the stride time. In this paper, the stride time is identified by the time elapsed between successive contacts of the same foot with the ground. The Fourier series can also be represented in the so-called phase-angle form: Suppose that ( ) describes the instantaneous linear velocity of a point on the human body in the progression direction. The DC component is identified with the APV of this point, with the cyclical component of velocity being described by the partial sum series in Equations (1) and (2). Figure 2 reports the block diagram of the method proposed in this paper to estimate the instantaneous velocity of the pelvis. All computations were performed using functions written in Virtual markers were created in correspondence of the center of the tri-axial accelerometer chips, where the origins of the ULFs were located, by using the positions of the four retro-reflective markers associated to each IMU (Tables 1 and 2).

Mathematical Processing
During steady-state level locomotion it is reasonable to assume that limb, or trunk displacement data can be modeled as quasi-periodic functions of time. Body segment displacement data have been obtained that can be accurately described through Fourier series that contain up to M = 6 significant harmonics [4,7]: where the DC component c 0 accounts for the effect of forward motion and the period T is the stride time. In this paper, the stride time is identified by the time elapsed between successive contacts of the same foot with the ground. The Fourier series can also be represented in the so-called phase-angle form: Suppose that f (t) describes the instantaneous linear velocity of a point on the human body in the progression direction. The DC component c 0 is identified with the APV of this point, with the cyclical component of velocity being described by the partial sum series in Equations (1) and (2). Figure 2 reports the block diagram of the method proposed in this paper to estimate the instantaneous velocity of the pelvis. All computations were performed using functions written in MATLAB (The MathWorks, Natick, MA, USA).

Gait Phases Segmentation
The ML component of the angular velocity measured by the shank gyroscope was used to perform the gait phases segmentation. A four-state left-right Hidden Markov Model (HMM) was developed with two-dimensional emission vectors that included the measured angular velocity and its first time difference [25]. The emissions were modeled as Gaussian mixtures with three modes. The gait phases that were paired to model states were defined by consecutive occurrences of Foot Strike (FS), Flat Foot (FF), Heel Off (HO) and Toe Off (TO).
IMU and OMCS data from the TW dataset were used for HMM training. The data from the foot markers were used to generate the reference data for the FS, FF, HO and TO gait events in the k-th gait cycle, namely ( ), ( ), ( ), ( ). Model parameters were trained in a supervised way, and testing was done using a Leave-One-Out-Subject (LOSO) validation study: one subject was tested with a model learnt using data from all other subjects; the training-testing procedure was then repeated, with all subjects being used once for testing.

Target Extraction
The nominal speed of the treadmill machine was the reference chosen for the APV (TW-APVREF) in the TW dataset. In the case of the OW dataset a different procedure was needed. A task-specific HMM was not trained, since, in contrast with the TW dataset, OMCS foot-marker data were not available (Table 2); rather, we used the HMM structure learnt using the TW dataset. The X and Y coordinates of the virtual marker at the pelvis (shank) were submitted to a first central difference method for numerical differentiation, after being low-pass filtered using a forward-backward second-order Butterworth filter (cut-off frequency: 3 Hz). The norm of the instantaneous velocity in the horizontal plane was then computed and averaged for each detected gait cycle, yielding the pelvis (shank) OW-APVREF.

Strap-Down Rotation
The Extended Kalman Filter (EKF) discussed in [24] was used to estimate the quaternion from the pelvis ULF to the GGF. The estimated quaternion was used to perform the strap-down rotation of the accelerometer output, so as to obtain the linear acceleration of the pelvis IMU. In the EKF design the magnetic sensor measurements were dismissed, although they were potentially helpful to stabilize the yaw estimate; this was done to meet the experimental setups that are typically used in studying the kinematics of pelvic motion [10].

Fourier Analysis for Estimating the Cyclical Component
Data from the pelvis IMU were submitted to Fourier analysis at each gait cycle (stride) [10]. In

Gait Phases Segmentation
The ML component of the angular velocity measured by the shank gyroscope was used to perform the gait phases segmentation. A four-state left-right Hidden Markov Model (HMM) was developed with two-dimensional emission vectors that included the measured angular velocity and its first time difference [25]. The emissions were modeled as Gaussian mixtures with three modes. The gait phases that were paired to model states were defined by consecutive occurrences of Foot Strike (FS), Flat Foot (FF), Heel Off (HO) and Toe Off (TO).
IMU and OMCS data from the TW dataset were used for HMM training. The data from the foot markers were used to generate the reference data for the FS, FF, HO and TO gait events in the k-th gait cycle, namely t FS (k), t FF (k), t HO (k), t TO (k). Model parameters were trained in a supervised way, and testing was done using a Leave-One-Out-Subject (LOSO) validation study: one subject was tested with a model learnt using data from all other subjects; the training-testing procedure was then repeated, with all subjects being used once for testing.

Target Extraction
The nominal speed of the treadmill machine was the reference chosen for the APV (TW-APV REF ) in the TW dataset. In the case of the OW dataset a different procedure was needed. A task-specific HMM was not trained, since, in contrast with the TW dataset, OMCS foot-marker data were not available (Table 2); rather, we used the HMM structure learnt using the TW dataset. The X and Y coordinates of the virtual marker at the pelvis (shank) were submitted to a first central difference method for numerical differentiation, after being low-pass filtered using a forward-backward second-order Butterworth filter (cut-off frequency: 3 Hz). The norm of the instantaneous velocity in the horizontal plane was then computed and averaged for each detected gait cycle, yielding the pelvis (shank) OW-APV REF .

Strap-Down Rotation
The Extended Kalman Filter (EKF) discussed in [24] was used to estimate the quaternion from the pelvis ULF to the GGF. The estimated quaternion was used to perform the strap-down rotation of the accelerometer output, so as to obtain the linear acceleration of the pelvis IMU. In the EKF design the magnetic sensor measurements were dismissed, although they were potentially helpful to stabilize the Sensors 2016, 16, 2206 6 of 17 yaw estimate; this was done to meet the experimental setups that are typically used in studying the kinematics of pelvic motion [10].

Fourier Analysis for Estimating the Cyclical Component
Data from the pelvis IMU were submitted to Fourier analysis at each gait cycle (stride) [10]. In short, each linear acceleration component (namely: ML, VT, AP) of the k-th stride was analyzed and the Fourier series coefficients were computed up to M = 6: where T k is the duration of the k-th stride. The analytical integration of Equation (3) where a ML0 (k), a VT0 (k) and a AP0 (k) were set to zero (equivalent to removing a constant term from the original stride linear acceleration data) led to the following expression of the mean-subtracted k-th stride velocity data: It is noted that the mean subtraction in Equation (3) is implemented in recognition of practical difficulties in estimating the pelvis APV by time integration of noisy acceleration data [8].
One advantage of the Fourier regression (Equations (1) and (2)) is that, differently from, e.g., a polynomial regression, the coefficients computed up to a given order do not change when one is interested in adding further terms to the model (e.g., for achieving a better model fit to the original time function). Moreover, when dealing with raw data, effective smoothing can be achieved by using Fourier basis functions [26]. The Fourier analysis was performed using the methods of functional data analysis developed in [26], for the implementation of which a MATLAB toolbox is available.

Fourier Analysis for Estimating the APV
For the purpose of APV prediction, we propose to use the coefficients {c i (k)} M i=1 of the Fourier series in the phase-angle form, which were computed using the k-th stride signals from the shank IMU. Six measurement channels were available (three for the acceleration, three for the angular velocity). The feature vector built at the k-th stride included, in addition to the Fourier harmonic coefficients for each measurement channel, the k-th stride time T FS (k), and some additional temporal parameters of gait from the HMM: The feature vector had thus size d = 6M + 4. By preliminary testing, a good trade-off between fitting accuracy and APV predictive power was achieved by retaining the first two harmonic components of shank angular velocity and acceleration (i.e., M = 2 and d = 16).

Bayesian Linear Regression
The model was created from a training set of N observations is the feature vector and v(i) is the APV target paired to the i-th gait cycle. Henceforth, the feature vectors are denoted collectively by the N × D data matrix X, whose n-th row is x T n (n = 1, . . . , N), and the corresponding target values are given by the column vector For some parameter vector w, the targets v(i) are given by adding noise to a linear combination of the input variables: v where the residual ε(i) is modeled as a zero-mean Gaussian random variable with precision (inverse variance) α. When a bias weight w 0 or offset is included in the parameter (weight) vector w, as in Equation (6), the input vector is augmented with an additional element whose value is always one, which leads to the following expression of the N × (D + 1) design matrix Φ: In the BLR approach a zero-mean Gaussian prior with precision β i is assigned to each parameter Under the assumption of independence of the marginal distributions, the weight prior is written as follows: where B = diag(β 0 , . . . , β d ).
The likelihood function of the target data is given by: The posterior distribution p(w|X, B , α) for the weights is Gaussian, with mean and covariance given by: The values of the hyperparameters α, B are obtained by maximization of the log-marginal likelihood: where the N × N matrix C is written in the following form: The point estimates of B, α can be substituted back into Equation (11) to give an updated posterior distribution for the weights, from which re-compute the log-marginal likelihood. B and α can be thus iteratively updated for maximization [20]: where µ i is the i-th component of the posterior mean µ defined by Equation (11) and the quantity γ i is defined by in which Σ ii is the i-th diagonal component of the posterior covariance matrix Σ defined by Equation (11). The iterative process stops when a criterion of convergence is met, i.e., the largest change in the values of the hyperparameters is below a tolerance. During re-estimation, some β i 's can become very large, shrinking the corresponding posteriors p(w i |X, B, α) to zero; this implies that the corresponding i-th column in Φ can be "pruned" (Automatic Relevance Determination (ARD)). In our implementation, the β i 's that were 200 times larger than the data precision α were discarded, removing the corresponding feature variables from the input space [21].
The LOSO validation was conducted to evaluate the performance of each regression model. A difference exists in the way TW and OW datasets were managed in this regard: the LOSO validation of the HMM-based gait event detector and the LOSO validation of the ARD-BLR model were intertwined in the former case. Conversely, in the latter case the LOSO validation of the HMM was not performed, since the HMM structure emerging from analyzing the TW dataset was used.

Performance Assessment
Steady-state gait strides were considered as for the TW dataset; data from the first 30-s periods of recording were thus discarded; all gait strides were retained for analysis in the case of the OW dataset. The performance assessment for the cyclical component of velocity was based on the procedure devised for method comparison studies [27,28]. The Mean Difference (MD) and upper and lower Limit of Agreement (LA), namely MD ± 1.96 Standard Deviation (SD) of differences, were computed for the cyclical component of velocity in the ML, VT and AP directions, as it was obtained from IMU and OMCS data. Henceforth, the LA width is the difference between the upper and lower LAs. Scatter plots were produced to visualize differences between IMU and OMCS data against their mean.
The APV BLR values during steady-state gait cycles were compared with the treadmill speed (TW-APV REF ), yielding the estimation error e ln per speed condition. The index l run over the number of strides L n that were walked by the n-th participant in each walking trial. In the case of overground locomotion, the APV BLR values were compared with the corresponding target values (pelvis or shank OW-APV REF ).
The mean bias error (MBE) was the mean value of e ln where, before taking the mean over participants, the mean error per participant m n was computed [15]:

Results
The HMM-based gait event detector estimated the time occurrences of the FS, FF, HO and TO events. The mean, the SD and the Mean Absolute Value (MAV) of their difference from the OMCS reference values, averaged across subjects, are reported in Table 3. A total number of 3377 strides were analyzed; missed and additionally detected gait strides (deletions and insertions, respectively) were not observed. As for the OW dataset, mean = 7.0 ms, SD = 29.8 ms, and MAV = 37.1 ms were obtained by analyzing the difference between the FS time occurrences delivered by the HMM and the Zijlstra's method [29], which was applied to the pelvis accelerometer. A total number of 753 strides were analyzed, without observing any erroneous event. Scatter plots of the difference between IMU and OMCS estimates of the cyclical component of the pelvis instantaneous velocity over their mean are reported in Figure 3 (data from all subjects and speed conditions were collapsed in producing the scatter plots).
Slight tendencies are observed for differences being increasingly negative with increasing mean value of the velocity; hence, compared to OMCS, the Fourier-based integration method slightly underestimated and overestimated for respectively negative and positive values of the velocity resolved in the GGF. Moreover, slight tendencies are observed for the spread of the differences to vary over the measurement range, which may produce non-constant LAs. In the attempt to refine the LA computation, the regression approach for non-uniform differences was applied [27]; the correction equations needed to compute the LA widths are reported in Table 4.
For each direction, a representative value of the LA width was finally obtained by taking the average of the widths that were computed over a range of velocity covering 95% of the measured values; the following values were obtained: ±0.  Scatter plots of the difference between IMU and OMCS estimates of the cyclical component of the pelvis instantaneous velocity over their mean are reported in Figure 3 (data from all subjects and speed conditions were collapsed in producing the scatter plots).  Slight tendencies are observed for differences being increasingly negative with increasing mean value of the velocity; hence, compared to OMCS, the Fourier-based integration method slightly underestimated and overestimated for respectively negative and positive values of the velocity resolved in the GGF. Moreover, slight tendencies are observed for the spread of the differences to vary over the measurement range, which may produce non-constant LAs. In the attempt to refine the LA computation, the regression approach for non-uniform differences was applied [27]; the  The ARD process pruned some temporal features (i.e., T FF , T HO ) and all Fourier coefficients of the shank angular velocity; the stride time T FS and the duration of the stance phase T TO were retained, together with the Fourier series coefficients of the shank acceleration, yielding an input space of dimension d = 8, either in TW or OW conditions. The evolution of the hyperparameters in two representative runs of the ARD process is shown in Figure 4.  The performance metrics of the ARD-BLR are reported in Tables 5 and 6, for TW and OW datasets, respectively (cross-validation was task-specific, namely training and testing were done using data from the same walking condition). As for the OW dataset, two different targets were considered, namely the pelvis and the shank OW-APVREF. Especially when walking along a curved path, the velocity fields changed across the body: Table 6 also reports the statistics of the pelvis and The priors of the stride time and stance duration were given precisions β 1 and β 2 , respectively; the priors of the ML, AP and VT components of the shank acceleration were given precisions β 3β 4 , β 5β 6 and β 7β 8 , respectively; the Gaussian noise in the likelihood function was given precision α. The threshold for pruning the i-th column of the design matrix is depicted using the horizontal line in black.
The performance metrics of the ARD-BLR are reported in Tables 5 and 6, for TW and OW datasets, respectively (cross-validation was task-specific, namely training and testing were done using data from the same walking condition). As for the OW dataset, two different targets were considered, namely the pelvis and the shank OW-APV REF . Especially when walking along a curved path, the velocity fields changed across the body: Table 6 also reports the statistics of the pelvis and shank OW-APV REF (minimum and maximum value, designated min and max, respectively, and the 25%, 50%, and 75% percentile, designated Q 1 , Q 2 and Q 3 , respectively). The final step of analysis consisted of training with one dataset and testing with the other dataset. The performance metrics of the ARD-BLR are reported in Tables 7 and 8. Data in Table 7 were obtained when the regression model was trained with the OW dataset (target: shank OW-APV REF ) and were tested with the TW dataset. Data in Table 8 concern the case in which the regression model was trained with the TW dataset and were tested with the OW dataset. The generalization capabilities of the method across the two walking conditions improved when using the pelvis OW-APV REF as target for testing.  Figure 5 shows representative examples of pelvis instantaneous velocity, when cross-validation was task-specific; transient stride data were also inspected for the treadmill-related example.  Figure 5 shows representative examples of pelvis instantaneous velocity, when cross-validation was task-specific; transient stride data were also inspected for the treadmill-related example.

Discussion
The HMM-based gait phases segmentation method delivered accurate estimates of the gait events' time occurrences, especially as for the FS and TO. FF and HO were detected with greater difficulty [31] (Table 3). The algorithm was capable of generalizing well across different tested subjects and walking conditions, although overground walking is known to be different from treadmill walking [32,33]; moreover, the participants recruited for constructing the two datasets were different people, and even the experimental setups were different in the two cases (Tables 1  and 2). The Zijlstra's method, which is widely used to detect gait strides [34], matched closely the HMM-based method as for the estimates of the FS events in conditions of overground walking.
As for the estimation of the cyclic component of velocity, the regression approach for non-uniform differences allowed accounting for slight systematic differences between the IMU and OMCS data, and to establish average LA widths across the measurement range (Table 4). It is noted that the agreement appears to be slightly looser in OW conditions, as compared with TW conditions, particularly for the horizontal components (i.e., ML and AP), compared with the vertical component (VT). This fact can be explained by the much larger measurement volume of the OMCS and by the longer duration of the OW trials. Both these factors make the procedure of aligning the ULF to the GGF in the direction of travel less accurate [24].
The ARD automatically removed feature variables that were irrelevant to the regression task. The same input space was considered to predict the APV in either TW or OW conditions, with one important difference. As shown in Figure 3, the ARD gave less importance to the temporal feature

Discussion
The HMM-based gait phases segmentation method delivered accurate estimates of the gait events' time occurrences, especially as for the FS and TO. FF and HO were detected with greater difficulty [31] ( Table 3). The algorithm was capable of generalizing well across different tested subjects and walking conditions, although overground walking is known to be different from treadmill walking [32,33]; moreover, the participants recruited for constructing the two datasets were different people, and even the experimental setups were different in the two cases (Tables 1 and 2). The Zijlstra's method, which is widely used to detect gait strides [34], matched closely the HMM-based method as for the estimates of the FS events in conditions of overground walking.
As for the estimation of the cyclic component of velocity, the regression approach for non-uniform differences allowed accounting for slight systematic differences between the IMU and OMCS data, and to establish average LA widths across the measurement range (Table 4). It is noted that the agreement appears to be slightly looser in OW conditions, as compared with TW conditions, particularly for the horizontal components (i.e., ML and AP), compared with the vertical component (VT). This fact can be explained by the much larger measurement volume of the OMCS and by the longer duration of the OW trials. Both these factors make the procedure of aligning the ULF to the GGF in the direction of travel less accurate [24].
The ARD automatically removed feature variables that were irrelevant to the regression task. The same input space was considered to predict the APV in either TW or OW conditions, with one important difference. As shown in Figure 3, the ARD gave less importance to the temporal feature variables in OW conditions as compared with TW conditions (i.e., the corresponding priors were given higher precisions β 1 and β 2 ). Locomotion along a curved path is in fact a task that implies continuous deviation from straight-ahead locomotion, thereby requiring continuous adjustment of body movement [33]. For example, stride length is unchanged for the outer but decreases for the inner leg; however, the cadence is the same for both legs, in spite of the different length of the inner and outer strides. In other words, whilst the increase of the treadmill speed determined a marked reduction of stride time and duration of the stance phase [35], speed changes occurring during curved path locomotion were not accompanied by significant changes of either of them.
The predictive power of the input space slightly decreased in OW conditions compared with TW conditions when the validation was task-specific, Tables 5 and 6. ARMSEs of less than 7.5%-8% were reported at the speed of 3 km/h, which fell at slightly less than 4% at velocities greater than 4 km/h. These results compare favorably with the state-of-the-art reported in the literature (e.g., ARMSE = 4%, overground walking with velocities 4-6 km/h, using a 16-dimensional input space of time-and frequency-domain feature variables from a foot IMU, in combination with a Linear Least Squares model [15]). Since, in contrast with RMSE, the ARMSE accounts for the intra-subject variability [15], it is not surprising that the ARMSE was larger than the RMSE in OW conditions, compared with TW conditions, due to the high walking regularity enforced by the treadmill machine.
When the regression model trained in OW conditions was applied to TW data for testing purposes (Table 7), the predictive performances degraded slightly compared with those obtained by the task-specific validation of Table 5. Using the regression model trained in TW conditions directly on OW data performed poorly, especially when the target for testing was the pelvis OW-APV REF , Table 8. Using the TW dataset, the shank acceleration signal features were used to predict the treadmill velocity, and the shank APV was matched to that velocity, Table 5. Since the forward component of the pelvis IMU motion was small (ideally, null) when walking on the treadmill machine, no effort was spent to learn the relationship between the input space and the pelvis APV. The use of the shank OW-APV REF as target improved significantly the generalization abilities of the ARD-BLR across the two walking conditions. This was also true when the treadmill velocity was predicted using the regression model learnt from the OW data. As shown in Table 8, the best predictive performance were observed at 4 km/h, a value close to the mode of the distribution of the shank OW-APV REF ( Table 6); for the highest treadmill speeds, the ARD-BLR worked in extrapolation mode, which reflects in less predictive power. Under-estimation at 5 and 6 km/h treadmill speed were observed, which did not occur when validation was task-specific. ARMSEs of 9% were reported at the speed of 3 km/h, which fell at 5% at speeds greater than 4 km/h, implying a one percentage-point drop compared with when the validation was task-specific. Another factor, namely that the treadmill nominal speed might not perfectly match the actual walking speed, would also have contributed to this one percentage-point drop.
Finally, the good behavior of the proposed method is shown, qualitatively, in Figure 5; here the analysis started before the treadmill reached the preset velocity, although the ARD-BLR worked in extrapolation mode in these conditions. As for the Fourier-based method of integration, it is noted that the method could be successfully applied to data in the time interval between static upright posture and steady-state locomotion. The applicability of the proposed method during gait initiation is currently being studied in more detail, although we verified that the LA widths did not change appreciably when transient strides were accounted for in their calculation.
The model predictive power could be enhanced by subject-specific model calibration (i.e., biometric parameters, such as the height, are included in the input space) or personalization (i.e., few gait strides from a new subject for whom reference data are available are used to refine the model, thereby promoting a better match of estimates to the reference) [15]. Subject-specific calibration and personalization of the model were not considered in this paper and are left to our future work, which will be devoted to study mechanical energy changes of the BCOM in normal and pathologic walking using inertial motion sensors.

Conclusions
In this paper, an inertial sensor-based algorithm was developed with the aim of estimating the instantaneous velocity of an IMU attached to the pelvis during walking.
Under the assumption about the cyclical motion of human body parts during walking, the instantaneous velocity was modeled by the sum of two components, namely the cyclical component and the average progression velocity at each gait cycle. The algorithm made use of methods of Fourier harmonic analysis applied to pelvis and shank acceleration data to address two tasks: analytical time-integration of the linear acceleration, which enabled the estimation of the cyclical component, and regression using the ARD-BLR approach for estimating the average progression velocity. Analytical time-integration and regression required gait phases segmentation, which was efficiently done using an HMM-based gait event detector.
The inertial sensor-based algorithm was validated in conditions of treadmill and overground walking by healthy subjects. Analytical integration based on Fourier series coefficients was thus shown a useful approach to accurately estimate instantaneous velocity data from noisy acceleration measurements, whilst good generalizability of the ARD-BLR across different factors (namely, subjects, walking conditions, and IMU hardware) was demonstrated.