Estimation of Tri-Axial Walking Ground Reaction Forces of Left and Right Foot from Total Forces in Real-Life Environments

Continuous monitoring of natural human gait in real-life environments is essential in many applications including disease monitoring, rehabilitation, and professional sports. Wearable inertial measurement units are successfully used to measure body kinematics in real-life environments and to estimate total walking ground reaction forces GRF(t) using equations of motion. However, for inverse dynamics and clinical gait analysis, the GRF(t) of each foot is required separately. Using an experimental dataset of 1243 tri-axial separate-foot GRF(t) time histories measured by the authors across eight years, this study proposes the ‘Twin Polynomial Method’ (TPM) to estimate the tri-axial left and right foot GRF(t) signals from the total GRF(t) signals. For each gait cycle, TPM fits polynomials of degree five, eight, and nine to the known single-support part of the left and right foot vertical, anterior-posterior, and medial-lateral GRF(t) signals, respectively, to extrapolate the unknown double-support parts of the corresponding GRF(t) signals. Validation of the proposed method both with force plate measurements (gold standard) in the laboratory, and in real-life environment showed a peak-to-peak normalized root mean square error of less than 2.5%, 6.5% and 7.5% for the estimated GRF(t) signals in the vertical, anterior-posterior and medial-lateral directions, respectively. These values show considerable improvement compared with the currently available GRF(t) decomposition methods in the literature.


Introduction
The tri-axial walking ground reaction forces, moments and the trajectory of center of plantar pressure CoP(t) under each foot are critical inputs for gait analysis [1]. Regardless of their importance, developing practical technologies for long-term measurement of these parameters in real-life environment are still challenging.
Estimation of the total (superimposed left and right foot) walking GRF i (t), i ∈ {v : vertical, ap : anterior − posterior, ml : medial − lateral} using a limited number of inertial measurement units (IMUs) is a practical option to monitor the total GRF i (t) in real-life environments [2,3] due to the high durability, low power demand, low cost, and small size of the IMUs. In this method, the total GRF i (t) signals are often estimated based on the Second Newton law using the measured acceleration of the body segments. During the single-support (SS) phase of the gait, GRF i (t) of the stance foot is equal to the estimated total GRF i (t) signals. However, during the double-support (DS) phase, both feet are in contact with the ground, which creates an indeterminate closed kinematic 1.
The ground reaction forces and moments of the trailing foot reduce smoothly to zero during the DS phase. 2.
The ratios of the ground reactions to their values at contralateral heel strike (i.e., the non-dimensional ground reactions) can be expressed as functions of DS phase duration (termed transition functions).
The model proposed by Ren et al. [15] showed good results in the sagittal plane, but results were less promising in the frontal and transverse planes. They reported 5.6%, 10.9% and 20% relative root mean square error for estimation of separate feet GRF v (t), GRF ap (t), and GRF ml (t), respectively. Subsequently, several studies were carried out to increase the accuracy of the estimated left and right foot GRF v (t) [16][17][18][19]. However, the proposed methods in these studies require the trajectory of CoP(t) as input and, therefore, are mostly suitable for force plate application.
Recently, Karatsidis et al. [20] proposed a comprehensive methodology to estimate tri-axial GRF i (t) and moment signals for each foot using inertial motion capture data. The total external forces and moments were calculated directly from equation of motion. A distribution algorithm based on a smooth transition assumption was then proposed to solve the indeterminacy problem during the DS phase. They reported relative root mean square error of 5.3%, 9.4% and 13.1% for the estimated GRF v (t), GRF ap (t) and GRF ml (t), respectively, compared with the force plate measurements. An extensive review of the literature on estimation of walking ground reactions in real-life environments is presented in [21].
This study is aimed at developing a methodology to estimate the tri-axial left and right foot GRF i (t) from the corresponding total GRF i (t) signal with: (1) enhanced accuracy; (2) enhanced versatility (applicable to different walking speeds and step lengths); (3) extensively validated and (4) suitable for IMU application i.e., only total GRF i (t) signals are available as input. To achieve these objectives, a uniquely extensive experimental dataset of 1243 tri-axial left and right foot GRF i (t) time histories with different walking speeds were used to develop the 'Twin Polynomial Method' (TPM) to estimate the tri-axial left and right foot GRF i (t) from the corresponding total GRF i (t) signal. Section 2 of this paper explains the TPM methodology and specifications. TPM procedure is presented in Section 3.1 and is validated in both laboratory and outdoor environment in Section 3.2. Section 4 discusses the results and compares the performance of TPM with other methods in the literature. Finally, the conclusions are highlighted in Section 5.

TPM Methodology
The Twin Polynomial Method proposed in this study assumes that the total GRF i (t) signals (estimated using a wearable IMU system [21]) and the weight of the subject are the only known inputs to the model and the left and right foot GRF i (t) signals are the desired outputs.
As mentioned in Section 1, during the SS phase of the gait, only the stance foot is in contact with the ground. Therefore, the magnitude of the stance foot walking force (blue curves in Figure 1) is equal to the total GRF i (t) (dashed black curves in Figure 1) and the amplitude of the walking force of the swing foot is equal to zero. The challenge is, however, to estimate the left and right foot GRF i (t) during the DS phase (red curves in Figure 1). The core idea of TPM is to fit a polynomial curve of degree n to the known SS part of each of the left and right foot GRF i (t) signals (blue curve in Figure 1) to extrapolate the unknown DS part of the corresponding signals (red curves in Figure 1).
To simplify the estimation process, TPM divides each complete gait cycles into two halves and estimates the unknown DS parts of the left and right foot GRF i (t) signals for each half gait cycle, separately ( Figure 1). In the vertical direction, the window of the total GRF v (t) signal between the two consecutive SS local minima is taken equal to half of a complete gait cycle (Figure 1a). Similarly, in the AP direction, the window of the total GRF ap (t) signal between two consecutive zero crossing points (Figure 1b), and in ML direction the window of the total GRF ml (t) signal between the two consecutive SS local minima and maxima points ( Figure 1c) are taken equal to half of a complete gait cycle. The left and right foot GRF i (t) signals corresponding to half of a gait cycle are referred to as hGRF i,j (t), where i ∈ {v, ap, ml} indicates the axis and j ∈ {l, r} represents the left and right part of the half gait force signal, respectively.

TPM Methodology
The Twin Polynomial Method proposed in this study assumes that the total ( ) signals (estimated using a wearable IMU system [21]) and the weight of the subject are the only known inputs to the model and the left and right foot ( ) signals are the desired outputs.
As mentioned in Section 1, during the SS phase of the gait, only the stance foot is in contact with the ground. Therefore, the magnitude of the stance foot walking force (blue curves in Figure 1) is equal to the total ( ) (dashed black curves in Figure 1) and the amplitude of the walking force of the swing foot is equal to zero. The challenge is, however, to estimate the left and right foot ( ) during the DS phase (red curves in Figure 1). The core idea of TPM is to fit a polynomial curve of degree n to the known SS part of each of the left and right foot ( ) signals (blue curve in Figure 1) to extrapolate the unknown DS part of the corresponding signals (red curves in Figure 1). To simplify the estimation process, TPM divides each complete gait cycles into two halves and estimates the unknown DS parts of the left and right foot ( ) signals for each half gait cycle, separately ( Figure 1). In the vertical direction, the window of the total ( ) signal between the two consecutive SS local minima is taken equal to half of a complete gait cycle (Figure 1a). Similarly, in the AP direction, the window of the total ( ) signal between two consecutive zero crossing points (Figure 1b), and in ML direction the window of the total ( ) signal between the two consecutive SS local minima and maxima points ( Figure 1c) are taken equal to half of a complete gait cycle. The left and right foot ( ) signals corresponding to half of a gait cycle are referred to as ℎ , ( ), where ∈ { , , } indicates the axis and ∈ { , } represents the left and right part of the half gait force signal, respectively.

Experimental Data
As the walking force is characterized with high inter-and intra-subject variability, a large dataset of measured walking force was required for the analysis to represent these variations statistically. For each of the measured left and right foot ( ) signals, the heel-strike and toe-off points corresponding to each gait cycle were identified (Figure 1), and the signal corresponding to each half gait cycle was extracted. To make cross-comparison and statistical analysis possible, each ℎ , ( ) signal was then resampled to 100 points and normalized by the body weight of the test subject using MATLAB software [22]. This resulted in a dataset of over 100,000 of ℎ , ( ) signals. The overlaid plot of one foot ℎ , ( ) and ℎ , ( ) signals are presented in Figure 2a,b, respectively. As signals are resampled, the variable "t" in ℎ , ( ) signals does not represent the actual timing of the measured ℎ , ( ) signals any more. However, the instances of the resampled ℎ , ( ) signals are still referred to by the variable "t" in this study, indicating their nature.   Figure 1. A typical illustration of total (dashed black) and a single foot (Red and blue) GRF v (t) (a), GRF ap (t) (b) and GRF ml (t), (c) during a complete gait cycle.

Experimental Data
As the walking force is characterized with high inter-and intra-subject variability, a large dataset of measured walking force was required for the analysis to represent these variations statistically. A set of 1243 measured right and left foot walking GRF i (t) time histories, each lasting between 60 and 240 s and measured from over 200 different test subjects (age: 26 ± 8 years, weight: 77 ± 26 kg, height: 1.74 ± 0.21 m, and walking speed: 1.22 ± 0.44 m/s) with healthy gait were used in the analysis. Measurements were carried out using a bespoke separate-belt instrumented treadmill at 1 kHz between 2008 and 2016 in the biomechanics laboratory at the University of Sheffield. The investigations were carried out following the rules of the Declaration of Helsinki of 1975 (https://www.wma.net/what-we-do/medical-ethics/declaration-of-helsinki/), revised in 2008. Ethics approval for the walking tests was issued by the ethics committee of the Faculty of Engineering, University of Sheffield, and all subjects provided an informed consent in accordance with the University of Sheffield ethical guidelines for research involving human participants.
For each of the measured left and right foot GRF i (t) signals, the heel-strike and toe-off points corresponding to each gait cycle were identified (Figure 1), and the signal corresponding to each half gait cycle was extracted. To make cross-comparison and statistical analysis possible, each hGRF i,j (t) signal was then resampled to 100 points and normalized by the body weight of the test subject using MATLAB software [22]. This resulted in a dataset of over 100,000 of hGRF i,j (t) signals. The overlaid plot of one foot hGRF v,l (t) and hGRF v,r (t) signals are presented in Figure 2a,b, respectively. As signals are resampled, the variable "t" in hGRF i,j (t) signals does not represent the actual timing of the measured hGRF i,j (t) signals any more. However, the instances of the resampled hGRF i,j (t) signals are still referred to by the variable "t" in this study, indicating their nature.

Experimental Data
As the walking force is characterized with high inter-and intra-subject variability, a large dataset of measured walking force was required for the analysis to represent these variations statistically. For each of the measured left and right foot ( ) signals, the heel-strike and toe-off points corresponding to each gait cycle were identified (Figure 1), and the signal corresponding to each half gait cycle was extracted. To make cross-comparison and statistical analysis possible, each ℎ , ( ) signal was then resampled to 100 points and normalized by the body weight of the test subject using MATLAB software [22]. This resulted in a dataset of over 100,000 of ℎ , ( ) signals. The overlaid plot of one foot ℎ , ( ) and ℎ , ( ) signals are presented in Figure 2a,b, respectively. As signals are resampled, the variable "t" in ℎ , ( ) signals does not represent the actual timing of the measured ℎ , ( ) signals any more. However, the instances of the resampled ℎ , ( ) signals are still referred to by the variable "t" in this study, indicating their nature.

TPM Specifications
In order to estimate the unknown DS part of the hGRF i,j (t) signals, TPM proposes to fit a polynomial curve of degree n to the known SS part of the hGRF i,j (t) signals (e.g., A-B in Figure 3a) and their corresponding zero point (e.g., point C in Figure 3a). In this study, the MATLAB [22] function 'polyfit' was used for curve fitting which uses Vandermonde matrix to estimate the coefficients of the best polynomial fit in a least-square sense.

TPM Specifications
In order to estimate the unknown DS part of the ℎ , ( ) signals, TPM proposes to fit a polynomial curve of degree n to the known SS part of the ℎ , ( ) signals (e.g., A-B in Figure 3a) and their corresponding zero point (e.g., point C in Figure 3a). In this study, the MATLAB [22] function 'polyfit' was used for curve fitting which uses Vandermonde matrix to estimate the coefficients of the best polynomial fit in a least-square sense.

Identification and Validation Experimental Datasets
The measured ℎ , ( ) signals corresponding to 20 randomly-selected subjects were initially extracted from the experimental dataset to be used only for model validation. The remainder of the dataset was used for model identification. As the ℎ , ( ) dataset was very large, for practicality, it was necessary to use a smaller sub-set of it for model identification which statistically represents the main dataset. To find the size of this subset, a Monte Carlo analysis was carried out with increasing

Identification and Validation Experimental Datasets
The measured hGRF i,j (t) signals corresponding to 20 randomly-selected subjects were initially extracted from the experimental dataset to be used only for model validation. The remainder of the dataset was used for model identification. As the hGRF i,j (t) dataset was very large, for practicality, it was necessary to use a smaller sub-set of it for model identification which statistically represents the main dataset. To find the size of this subset, a Monte Carlo analysis was carried out with increasing number of measured hGRF i,j (t) signals fitted with polynomials of degree three (cubic) to 12. Each randomly selected hGRF i,j (t) signal from the experimental dataset, was fitted with a Sensors 2018, 18,1966 6 of 18 polynomial P n (t) of degree n, and the corresponding peak-to-peak normalized root mean square error (NRMSE), defined by Equation (1), was calculated: In Equation (1), the values of P n (t) when the foot is in swing is set to zero. It was found that a subset of minimum 940 hGRF i,j (t) randomly-selected signals can represent statistically the whole dataset in all three axes with less than 1% NRMSE standard deviation σ NRMSE . Therefore, a set of 1000 hGRF i,j (t) randomly-selected signals was used in Section 2.2.2 to Section 2.2.4 to develop TPM.

Polynomial Degree Selection
To find the optimal polynomial order that can fit accurately the hGRF i,j (t) signals, in each direction, each of the N = 1000 hGRF i,j (t) randomly-selected measured signals was fitted with polynomials of degree three (cubic) to 12. Figure 3 compares the probability distribution function (pdf) of the NRMSE errors corresponding to each polynomial degree. There is a tradeoff between the fit accuracy and the computational efficiency of the fitting process: the higher the polynomial order of the fitting curve is, the more computationally demanding and challenging would be to find the best fit, but a more accurate fit can potentially be achieved.
Based on Figure 3, polynomials of degree five, eight, and nine were selected empirically to fit the hGRF i,j (t) signals in the vertical, anterior-posterior, and medial-lateral directions, respectively. The choice was made to maximize the accuracy of the fit while avoiding unnecessary computation. As it can be seen in Figure 3, in each direction, the reduction in NRMSE of polynomials with lower orders than V:5, AP:8 and ML:9 is significant and for the higher orders is subtle.

Added Constraints
As it can be seen in Figure 3a In another word, the SS part and zero point are not enough constraints to produce a unique polynomial fit of degree five, eight, and nine to hGRF v,j (t), hGRF ap,j (t), and hGRF ml,j (t) signals, respectively.
Using Monte Carlo analysis, it was found that one extra (guiding) point in the vertical direction (point D in Figure 3a), three points in the AP direction (points H, I and J in Figure 3c), and three points in the ML direction (points O, P and Q in Figure 3e) are enough to achieve polynomial fits with maximum 1% difference in µ NRMSE and σ NRMSE compared with the polynomials in Figure 3b,d,f. Polynomials of higher order failed to achieve comparable µ NRMSE and σ NRMSE values with the same number of guiding points (2-5% higher error).
The magnitude (y) and timing (t) of the guiding points of the best polynomial fits, were statically analyzed ( Figure 4) for the random sample of N = 1000 hGRF i,j (t) signals to find out the statistical distribution of t B , t C , y D , y H , t I , y J , y O , t P , and y Q parameters (and y D , y H , t I , y J , y O , t P and y Q for the counterpart hGRF i,j (t) signals) suitable for curve-fitting analysis (Table 1). It was empirically assumed that The range of µ ± 2σ (representing 95% of the hGRF i,j (t) signals) corresponding to each of these parameters ( Figure 4) are used in TPM curve-fitting procedure as explained in Section 3 ( Table 1).

Optimization Strategy
As it can be seen in Figure 3a,c,e, a family of polynomial fits can be generated using different guiding points that match the known SS part of the hGRF i,j (t) signal while they do not necessarily match the DS part of the hGRF i,j (t) signal. The strategy proposed by TPM to circumvent this problem is to use the fitted polynomials on both the hGRF i,l (t) (Figure 3a,c,e) and hGRF i,r (t) signals and compare their total curve with the known total hGRF i (t) for that half gait cycle. This ensures that the information available in the total hGRF i (t) signal during the DS phase is optimally used to find the pair of polynomial curves that best fit the hGRF i,l (t) and hGRF i,r (t) signals. Section 3 describes the TPM procedure in details.

TPM Procedure
The step-by-step procedure proposed by TPM to estimate hGRF i,l (t) and hGRF i,r (t) from the known total hGRF i (t) is: In the first step, the SS local minima points are identified from the total GRF v (t) (Figure 5a).
Sensors 2018, 18, x FOR PEER REVIEW 9 of 18 distributed between the estimated ℎ , ( ) (dashed pink curve in Figure 5e) and ℎ , ( ) (dashed green curve in Figure 5e) signals using Equations (3) and (4) when < < and = : In these equations, the linear scaling functions ( ) (pink line in Figure 5e) and ( ) (green line in Figure 5e) are empirically defined as: e. The estimated ℎ , ( ) and ℎ , ( ) signals (dashed curves in Figure 5f) are multiplied by the weight of the subject and resampled back to the actual length of the measured total ℎ ( ) signal for the current half gait cycle.
III. The next half gait cycle is selected and the estimation procedure described in Step A.II is repeated.

B. Anterior-posterior direction
I.
In the first step, the SS phase zero crossing points are identified from the total ( ) (Figure 6a).
II. For each segment of the total ( ) signal between two consecutive zero crossing points (total ℎ ( )): a. The total ℎ ( ) segment is resampled to 100 points (T = 100) and normalized to the weight of the subject (Figure 6b). b. The timing of the heel-strike ( ) and toe-off ( ) points of this half gait cycle are identified from the previously-estimated ℎ , ( ) and ℎ , ( ) signals corresponding to this half gait cycle. c. A set of polynomial curves of degree eight are generated using the SS part (F-G) and points H, I, J and K for the ℎ , ( ) (Figure 6b) where = + 0.05 , = ( + )/2, 0.02 < < 0.10, −0.02 < < 0 and 0.46 < < 0.65 (Table 1) The difference between the measured (blue curve in Figure 6d) and estimated total ℎ ( ) (dashed red curve in Figure 6d) curves during the DS phase ( ) , are distributed between the estimated ℎ , ( ) (dashed pink curve in Figure 6d) and ℎ , ( ) (dashed green curve in Figure 6d) signals using Equations (3) and (4) when < < and = (Figure 6e). g. The estimated ℎ , ( ) and ℎ , ( ) signals are multiplied by the weight of the subject and resampled back to the actual length of the measured total ℎ ( ) signal for the current half gait cycle.
III. The next half gait cycle is selected and the estimation procedure described in Step B.II is repeated.

II.
For each segment of the total GRF v (t) signal between two consecutive local minima (total hGRF v (t)): a. The total hGRF v (t) segment is resampled to 100 points (T = 100) and is normalized to the weight of the subject (Figure 5b). b.
For each pair of toe-off (t C ) and heel-strike (t B ) points selected from their initial ranges (27 < t B < 52 and 53 < t C < 85): i. A set of polynomial curves of degree five are generated using the SS part (A-B) and points C and D for the hGRF v,l (t) (Figure 5b), and E-C, B, and D for the hGRF v,r (t) (Figure 5c) where −1 < y D and y D < 3.5, t D = t C + 0.1T and t D = t B − 0.1T (Table 1). ii.
A pair of left (pink curve in Figure 5d) and right (green curve in Figure 5d) polynomial curves are found to represent this t C , t B combination that their total curve (red curve in Figure 5d) estimated the total hGRF v (t) with minimum NRMSE. iii.
The procedure is repeated for all possible combinations of t C and t B selected from their initial ranges (Table 1).
c. The pair of left-right polynomial curves that estimated the total hGRF v (t) with minimum NRMSE is found. These curves are considered accurate estimates of the hGRF v,l (t) and hGRF v,r (t) signals for the current half gait cycle. d.
The difference between the measured (blue curve in Figure 5e) and estimated total hGRF v (t) (dashed red curve in Figure 5e) curves during the DS phase Err v (t), are distributed between the estimated hGRF v,l (t) (dashed pink curve in Figure 5e) and hGRF v,r (t) (dashed green curve in Figure 5e) signals using Equations (3) and (4) when t B < t < t C and i = v: In these equations, the linear scaling functions S l (t) (pink line in Figure 5e) and S r (t) (green line in Figure 5e) are empirically defined as: e. The estimated hGRF v,l (t) and hGRF v,r (t) signals (dashed curves in Figure 5f) are multiplied by the weight of the subject and resampled back to the actual length of the measured total hGRF v (t) signal for the current half gait cycle.
III. The next half gait cycle is selected and the estimation procedure described in Step A.II is repeated.
B. Anterior-posterior direction I.
In the first step, the SS phase zero crossing points are identified from the total GRF ap (t) (Figure 6a).

II.
For each segment of the total GRF ap (t) signal between two consecutive zero crossing points (total hGRF ap (t)): a. The total hGRF ap (t) segment is resampled to 100 points (T = 100) and normalized to the weight of the subject (Figure 6b). b.
The timing of the heel-strike (t G ) and toe-off (t K ) points of this half gait cycle are identified from the previously-estimated hGRF v,r (t) and hGRF v,l (t) signals corresponding to this half gait cycle. c.
A set of polynomial curves of degree eight are generated using the SS part (F-G) and points H, I, J and K for the hGRF ap,l (t) (Figure 6b) where t H = t G + 0.05T, t J = (t I + t k )/2, 0.02 < y H < 0.10, −0.02 < y J < 0 and 0.46T < t I < 0.65T (Table 1). d.
A set of polynomial curves of degree eight are generated using the SS part L-K , H , I , J and G for the hGRF ap,r (t) (Figure 6c) t H = t K − 0.05T, t J = (t I + t G )/2, −0.10 < y H < 0.14, −0.01 < y J < 0.07 and 0.20T < t I < 0.46T (Table 1). e.
A pair of left (Figure 6d-pink curve) and right (Figure 6d-green curve) polynomial curves that estimated the known total hGRF ap (t) with minimum NRMSE is found. This pair of polynomial fits are considered to be an accurate estimation of the hGRF ap,r (t) and hGRF ap,l (t) signals for the current half gait cycle.
f. The difference between the measured (blue curve in Figure 6d) and estimated total hGRF ap (t) (dashed red curve in Figure 6d) curves during the DS phase Err ap (t), are distributed between the estimated hGRF ap,l (t) (dashed pink curve in Figure 6d) and hGRF ap,r (t) (dashed green curve in Figure 6d) signals using Equations (3) and (4) when t G < t < t K and i = ap (Figure 6e). g.
The estimated hGRF ap,l (t) and hGRF ap,r (t) signals are multiplied by the weight of the subject and resampled back to the actual length of the measured total hGRF ap (t) signal for the current half gait cycle.
III. The next half gait cycle is selected and the estimation procedure described in Step B.II is repeated.

I.
In the first step, the SS local minima and maxima points are identified from the total GRF ml (t) (Figure 7a). II.
For each segment of the total GRF ml (t) signal between two consecutive local minima and maxima (total hGRF ml (t)): a. The total hGRF ml (t) signal is resampled to 100 points (T = 100) and normalized to the weight of the subject (Figure 7b). b.
The timing of the heel-strike (t N ) and toe-off (t R ) points for this half gait cycle are identified from the previously-estimated hGRF v,r (t) and hGRF v,l (t) signals corresponding to this half gait cycle. c.
A set of polynomial curves of degree nine are generated using the SS part (M-N) and points O, P, Q, and R for the hGRF ml,l (t) (Figure 7b) where t O = t N + 0.05T, t Q = (t P + t R )/2, −0.021 < y O < 0.011, −0.032 < y Q < 0.032 and 0.42T < t P < 0.67T (Table 1). d.
A set of polynomial curves of degree nine are generated using the SS part S − R , O , P , Q and N for the hGRF ml,r (t) (Figure 7c) t O = t R − 0.05T, t Q = (t N + t P )/2, −0.018 < y O < 0.034, −0.032 < y Q < 0.052 and 0.26T < t P < 0.65T (Table 1). e.
A pair of left (pink curve in Figure 7d) and right (green curve in Figure 7d) polynomial curves that estimated the known total hGRF ml (t) with minimum NRMSE is found.
These polynomial fits are considered accurate estimates of the hGRF ml,r (t) and hGRF ml,l (t) signals for the current half gait cycle. f.
The difference between the measured (blue curve in Figure 7d) and estimated total hGRF ml (t) (dashed red curve in Figure 7d) curves during the DS phase Err ml (t), are distributed between the estimated hGRF ml,l (t) (dashed pink curve in Figure 7d) and hGRF ml,r (t) (dashed green curve in Figure 7d) signals using Equations (3) and (4) when t N < t < t R and i = ml (Figure 7e). g.
The estimated hGRF ml,l (t) and hGRF ml,r (t) signals are multiplied by the weight of the subject and resampled back to the actual length of the measured total hGRF ml (t) signal for the current half gait cycle.
III. The next half gait cycle is selected and the estimation procedure described in Step C.II is repeated.  (Table 1). Figure 7d) and right (green curve in Figure 7d) polynomial curves that estimated the known total ℎ ( ) with minimum NRMSE is found. These polynomial fits are considered accurate estimates of the ℎ , ( ) and ℎ , ( ) signals for the current half gait cycle. f.

e. A pair of left (pink curve in
The difference between the measured (blue curve in Figure 7d) and estimated total ℎ ( ) (dashed red curve in Figure 7d) curves during the DS phase ( ) , are distributed between the estimated ℎ , ( ) (dashed pink curve in Figure 7d) and ℎ , ( ) (dashed green curve in Figure 7d) signals using Equations (3) and (4) when < < and = (Figure 7e). g. The estimated ℎ , ( ) and ℎ , ( ) signals are multiplied by the weight of the subject and resampled back to the actual length of the measured total ℎ ( ) signal for the current half gait cycle.
III. The next half gait cycle is selected and the estimation procedure described in Step C.II is repeated.

Model Validation
The ultimate goal of TPM is to be used to estimate left and right foot tri-axial GRF i (t) signals from corresponding total GRF i (t) signals in real-life environment. Therefore, performance of TPM is assessed for both laboratory and outdoor environment.

Performance of TPM in Laboratory Environment
TPM was used to estimate tri-axial hGRF i,l (t) and hGRF i,r (t) signals corresponding to one thousand half gait cycles, randomly selected from the validation dataset of measured hGRF i,j (t) signals. The NRMSE errors between the measured and estimated hGRF i,j (t) signals were then calculated and is shown in Figure 8a-c for the vertical, anterior-posterior, and medial-lateral directions, respectively. It was found that TPM estimated the left and right foot GRF v (t), GRF ap (t) and GRF ml (t) signals with mean NRMSE value of µ = 2.29%, µ = 6.27%, and µ = 7.22%, respectively.

Model Validation
The ultimate goal of TPM is to be used to estimate left and right foot tri-axial ( ) signals from corresponding total ( ) signals in real-life environment. Therefore, performance of TPM is assessed for both laboratory and outdoor environment.

Performance of TPM in Laboratory Environment
TPM was used to estimate tri-axial ℎ , ( ) and ℎ , ( ) signals corresponding to one thousand half gait cycles, randomly selected from the validation dataset of measured ℎ , ( ) signals. The NRMSE errors between the measured and estimated ℎ , ( ) signals were then calculated and is shown in Figure 8a-c for the vertical, anterior-posterior, and medial-lateral directions, respectively. It was found that TPM estimated the left and right foot ( ), ( ) and ( ) signals with mean NRMSE value of μ = 2.29%, μ = 6.27%, and μ = 7.22%, respectively.

Performance of TPM in Real-Life Environment
The walking gait in real-life environment is characterized with high variability in both magnitude and timing compared with the treadmill-measured ( ) signals. A set of tests was carried out where 10 subjects walked around the University of Sheffield campus building (in paved urban environment), while wearing a pair of Tekscan F-Scan in-shoe pressure insoles [23] to measure the benchmark ( ) signals for comparison. The walking pathway was characterized with flat parts as well as uphills and downhills. Subjects were asked to walk normally and no further instructions were given to keep the experiments as realistic as possible. The normal plantar pressures measured under each foot were used to calculate the left and right foot ( ) signals, and these signals were then summed up to calculate the measured total ( ). Before and after each trial, subjects walked with their normal speed on the instrumented treadmill while wearing pressure insoles, and the ( ) signals measured by the treadmill were used to calibrate the pressure insole measurements for each outdoor walking test. This helped to calibrate more accurately the F-scan pressure data, resulting in more accurate estimation of ( ) signals.

Performance of TPM in Real-Life Environment
The walking gait in real-life environment is characterized with high variability in both magnitude and timing compared with the treadmill-measured GRF i (t) signals. A set of tests was carried out where 10 subjects walked around the University of Sheffield campus building (in paved urban environment), while wearing a pair of Tekscan F-Scan in-shoe pressure insoles [23] to measure the benchmark GRF i (t) signals for comparison. The walking pathway was characterized with flat parts as well as uphills and downhills. Subjects were asked to walk normally and no further instructions were given to keep the experiments as realistic as possible. The normal plantar pressures measured under each foot were used to calculate the left and right foot GRF v (t) signals, and these signals were then summed up to calculate the measured total GRF v (t). Before and after each trial, subjects walked with their normal speed on the instrumented treadmill while wearing pressure insoles, and the GRF v (t) signals measured by the treadmill were used to calibrate the pressure insole measurements for each outdoor walking test.
This helped to calibrate more accurately the F-scan pressure data, resulting in more accurate estimation of GRM v (t) signals.
One hundred half-gait cycles were randomly selected from the measured total GRF v (t) signal and used as input to the TPM method to estimate the corresponding hGRF v,r (t) and hGRF v,l (t) signals. Figure 8d shows the NRMSE errors between the measured and estimated hGRF v,r (t) and hGRF v,l (t) signals. As it can be seen in this figure, the errors are comparable with the laboratory results with the mean NRMSE value of µ = 2.01%. It must be emphasized that, as these NRMSE values are calculated against the GRM v (t) signals measured by F-Scan as reference, they only represent the errors due to the GRF decomposition process (TPM) and do not include the errors due to the estimation of GRM v (t) signals using F-scan pressure insoles. Furthermore, the counter-intuitive fact that the real-life mean NRMSE value is marginally lower than the corresponding laboratory value must be interpreted in light of the fact that the benchmark hGRF v,r (t) signals used to calculate NRMSE values in real-life environment are measured using pressure insoles that has lower high-frequency sensitivity and higher error values compared with the instrumented treadmill measurements.

Application of TPM on IMU Data
Six healthy male subjects (age: 21 ± 1 years, weight: 77 ± 16 kg and height: 1.82 ± 0.08 m) participated in a set of walking gait measurement, where for each subject the treadmill speed was set to 60%, 70%, 80%, 90%, 100% and 110% of their normal walking speed. The tri-axial walking GRF(t) signals pertinent to each foot were recorded using the instrumented treadmill. A set of 12 Opal IMUs [24] were used to measure the tri-axial acceleration and orientation signals at the seventh cervical vertebrae (C7), fifth lumbar vertebrae (L5), upper arms, fore arms, thighs, shanks, and fourth metatarsals with 128 Hz sampling rate. The detailed explanation of the test protocol is presented in [25].
Model 2 proposed by Shahabpoor et al. [25] with subject specific training was initially used to estimate the tri-axial total walking GRF(t) signals, only from IMU measurements at C7, L5, and one of the thighs. TPM was then used to estimate tri-axial hGRF i,l (t) and hGRF i,r (t) signals corresponding to one thousand half gait cycles, randomly selected from these IMU-estimated total GRF(t) signals (Figure 9a,c,e). The mean NRMSE errors between the treadmill-measured and IMU-estimated hGRF i,j (t) signals were found equal to µ v = 7.12%, µ ap = 16.24%, and µ ml = 16.08%, for the vertical (Figure 9b One hundred half-gait cycles were randomly selected from the measured total ( ) signal and used as input to the TPM method to estimate the corresponding ℎ , ( ) and ℎ , ( ) signals. Figure 8d shows the NRMSE errors between the measured and estimated ℎ , ( ) and ℎ , ( ) signals. As it can be seen in this figure, the errors are comparable with the laboratory results with the mean NRMSE value of μ = 2.01%. It must be emphasized that, as these NRMSE values are calculated against the ( ) signals measured by F-Scan as reference, they only represent the errors due to the GRF decomposition process (TPM) and do not include the errors due to the estimation of ( ) signals using F-scan pressure insoles. Furthermore, the counter-intuitive fact that the real-life mean NRMSE value is marginally lower than the corresponding laboratory value must be interpreted in light of the fact that the benchmark ℎ , ( ) signals used to calculate NRMSE values in real-life environment are measured using pressure insoles that has lower highfrequency sensitivity and higher error values compared with the instrumented treadmill measurements.

Application of TPM on IMU Data
Six healthy male subjects (age: 21 ± 1 years, weight: 77 ± 16 kg and height: 1.82 ± 0.08 m) participated in a set of walking gait measurement, where for each subject the treadmill speed was set to 60%, 70%, 80%, 90%, 100% and 110% of their normal walking speed. The tri-axial walking ( ) signals pertinent to each foot were recorded using the instrumented treadmill. A set of 12 Opal IMUs [24] were used to measure the tri-axial acceleration and orientation signals at the seventh cervical vertebrae (C7), fifth lumbar vertebrae (L5), upper arms, fore arms, thighs, shanks, and fourth metatarsals with 128 Hz sampling rate. The detailed explanation of the test protocol is presented in [25]. Model 2 proposed by Shahabpoor et al. [25] with subject specific training was initially used to estimate the tri-axial total walking ( ) signals, only from IMU measurements at C7, L5, and one of the thighs. TPM was then used to estimate tri-axial ℎ , ( ) and ℎ , ( ) signals corresponding to one thousand half gait cycles, randomly selected from these IMU-estimated total ( ) signals (Figure 9a, c, and e). The mean NRMSE errors between the treadmill-measured and IMU-estimated ℎ , ( ) signals were found equal to = 7.12%, = 16.24%, and = 16.08%, for the vertical (Figure 9b), anterior-posterior (Figure 9d), and medial-lateral (Figure 9f)

Discussion
The accuracy of the results of the TPM method cannot be directly cross-compared with the literature due to fundamental differences in methodology. The NRMSE values reported in [15,20] correspond to the estimation of separate-feet ( ) signals directly from full kinematic measurements, and the methods suggested in [16][17][18][19] use measured ( ) as input.
However, the TPM methodology is developed using the most extensive experimental dataset presented so far in the literature, and the model is validated both in indoor and outdoor environments. The similarity of NRMSE values for indoor and outdoor environments confirms that the model is versatile and can estimate separate-feet ( ) signals for time-varying walking speeds and stride length in outdoor environments. The method, furthermore, only needs the total ( ) signals as input and therefore can be used with wearable IMU systems for monitoring separate-feet ( ) signals in real-life environment. The TPM however is developed and validated based on the data from young and healthy subjects. The measured data were also limited to smooth laboratory and outdoor walking surfaces.

Conclusions
In the present study, a method called Twin Polynomial Method was proposed based on a uniquely extensive dataset of measured walking ground reaction forces, to estimate the tri-axial left and right foot walking ground reaction forces from their corresponding total ( ) . The polynomials of degree five, eight, and nine were found to be the best candidates to fit right and left foot ( ) signals in the vertical, anterior-posterior, and medial-lateral directions, respectively.
Results of the TPM method both in the laboratory and in real-life environment showed an average normalized RMSE error of less than 2.5%, 6.5% and 7.5% in the vertical, anterior-posterior, and medial-lateral directions, respectively, compared with treadmill/pressure insole measurements. Further investigation is required for different pathological gaits and for walking on rough terrains to identify the required adjustments to the method.  These NRMSE values represent the errors due to both estimation of total GRF(t) signal using Model 2 and estimation of each foot GRF(t) using TPM method and the NRMSE pertinent to each method cannot be computed independently. Using the µ v = 7%, µ ap = 13% and µ ml = 13% NRMSE values reported in [25] for Model 2, decomposing the total GRF(t) into left and right foot GRF(t) using TPM method has increased the total NRMSE values by only µ v = 0%, µ ap = 3%, and µ ml = 3% which are favorably comparable with the NRMSE values reported in Sections 3.2.1 and 3.2.2 for TPM.

Discussion
The accuracy of the results of the TPM method cannot be directly cross-compared with the literature due to fundamental differences in methodology. The NRMSE values reported in [15,20] correspond to the estimation of separate-feet GRF i (t) signals directly from full kinematic measurements, and the methods suggested in [16][17][18][19] use measured CoP(t) as input.
However, the TPM methodology is developed using the most extensive experimental dataset presented so far in the literature, and the model is validated both in indoor and outdoor environments. The similarity of NRMSE values for indoor and outdoor environments confirms that the model is versatile and can estimate separate-feet GRF i (t) signals for time-varying walking speeds and stride length in outdoor environments. The method, furthermore, only needs the total GRF i (t) signals as input and therefore can be used with wearable IMU systems for monitoring separate-feet GRF i (t) signals in real-life environment. The TPM however is developed and validated based on the data from young and healthy subjects. The measured data were also limited to smooth laboratory and outdoor walking surfaces.

Conclusions
In the present study, a method called Twin Polynomial Method was proposed based on a uniquely extensive dataset of measured walking ground reaction forces, to estimate the tri-axial left and right foot walking ground reaction forces from their corresponding total GRF i (t). The polynomials of degree five, eight, and nine were found to be the best candidates to fit right and left foot GRF i (t) signals in the vertical, anterior-posterior, and medial-lateral directions, respectively. Results of the TPM method both in the laboratory and in real-life environment showed an average normalized RMSE error of less than 2.5%, 6.5% and 7.5% in the vertical, anterior-posterior, and medial-lateral directions, respectively, compared with treadmill/pressure insole measurements. Further investigation is required for different pathological gaits and for walking on rough terrains to identify the required adjustments to the method.

Declarations
Ethics approval and consent to participate: The investigations were carried out following the rules of the Declaration of Helsinki of 1975 (https://www.wma.net/what-we-do/medical-ethics/ declaration-of-helsinki/), revised in 2008. Ethics approval for the walking tests was issued by the ethics committee of the Faculty of Engineering, University of Sheffield, and all subjects provided an informed consent in accordance with the University of Sheffield ethical guidelines for research involving human participants.
Consent for publication: All subjects provided an informed consent in accordance with the University of Sheffield ethical guidelines for publication of the anonymized test data.
Availability of data and material: The GRF dataset analyzed during the current study are not publicly available as most of the measurements were done as part of confidential industrial projects and authors are not the owners of the dataset. However, a small subset of the data could be available from the corresponding author on reasonable request.
Author Contributions: Both authors have been involved in all stages of the study including measurements, data analysis, model development, verification and preparing the manuscript.

Funding:
The authors acknowledge the financial support provided by the UK Engineering and Physical Sciences Research Council (EPSRC) for the following research grants: Frontier Engineering Grant EP/K03877X/1 (Modelling complex and partially identified engineering problems: Application to the individualized multiscale simulation of the musculoskeletal system); and Platform Grant EP/G061130/2 (Dynamic performance of large civil engineering structures: an integrated approach to management, design and assessment).

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: A, B, C, D, E, and A', B', C', D', E' Guiding points related to hGRF v,j (t) (Figure 3a) ap Anterior-posterior direction CoP(t) Trajectory of center of plantar pressure DS Double-support phase of the walking gait Err i (t) Error in estimation of the total hGRF i (t) i ∈ {v, ap, ml} (Equation (2)  Polynomial of degree n pdf Probability Distribution Function S j (t) Error scaling function j ∈ {l, r} (Equations (5) and (6)) SS Single-support phase of the walking gait T Length of the resampled hGRF i,j (t) signals (T = 100)