Modeling and Individualizing Continuous Joint Kinematics Using Gaussian Process Enhanced Fourier Series

Prosthetic discrete controller relies on finite state machines to switch between a set of predefined task-specific controllers. Therefore the prosthesis can only perform a limited number of discrete locomotion tasks and need hours to tune the parameters for each user. In contrast, the continuous controller treats a gait cycle in a unified way. Thus it is expected to better facilitate normative biomechanics by providing a gait predictive model to contribute a non-switching controller that supports a continuum of tasks. Furthermore, a better method is to train a personalized trajectory prediction model suitable for personal characteristics according to personal walking data. This paper proposes a Gaussian process enhanced Fourier series (GPEFS) method to construct a gait prediction model that represents the human locomotion as a continuous function of phase, speed and slope. Firstly the joint trajectories are transformed into the Fourier coefficient space by least square method. Then the relationship between each Fourier coefficient and task input can be learned by multiple Gaussian process regression (GPRs) model respectively. Compared with directly using GPR to fit the joint trajectory under multi task, our method greatly reduces the computational burden, so as to meet the real-time application scenario. In addition, in Fourier coefficient space, the difference in all tasks between the Fourier coefficient of personal data and the one of statistical data follows the same trend. Therefore, a personalized prediction model is built to predict an individual’s kinematics over a continuous range of slopes and speeds given only one personalized task at level ground and normal speed. The experimental results show that the gait prediction model and the personalized prediction model are feasible and effective.

can greatly improve the amputee's life quality by restoring normative biomechanics of many locomotion tasks [1]. There are two main types of control paradigms of powered prosthesis [2]: one is the discrete control paradigm in which a gait is divided into multiple stages [3], [4], [5], [6], [7], [8], and the other is the continuous control paradigm that treats a gait cycle in a unified way [9], [10], [11], [12], [13]. Although the former has achieved great success, its main disadvantages, such as there are a large number of parameters need to be adjusted manually, have hindered its large-scale clinical application [4]. On the contrary, the latter creates a unified controller for the whole gait cycle. Thus it has no inherent drawbacks of the discrete one as described above. Early continuous controllers, such as the echo controller [9], [10], were time-based preprogrammed and cannot reflect the movement intentions of leg in prosthesis side. Recently, Gregg et al. [11], [12], [13] developed a very attractive continuous control method for prosthetics by introducing virtual constraints controller (VCC) originally used for bipedal robot control. VCC in the prosthetic field has two main steps: the first is to obtain phase variables (PV) that can reflect wearers' movement intentions; and the second is to build the gait prediction models, which will accept the input of PV to generate the desired trajectories of the prosthetic joints under different modes of motion. With respect to the constructing of PV, Gregg et al. initially used the center of pressure (COP) of foot [11] and subsequently used residual thigh movement information [12] to obtain PV. It was reported that both of them all achieved promising results.
With regard to building the gait prediction models, various methods have been proposed to estimate joint trajectories under level walking task with fixed or variable speed [14], [15], [16], [17], [18]. However, it is a challenge to construct a uniform gait prediction models to predict lower limb joints trajectories under a task with continuously varying slopes and speeds [22]. If it is required to generate personalized joint trajectories according to different amputee wearers, which will further complicate the problem. In general, we will train a corresponding gait prediction model according to a specific wearer under different walking tasks, it will not only require a large amount of expensive data from individuals to train the model, but also places high demands on the trial site.

A. Related Work
The related work of predicting or estimating the trajectory of lower limb joints can be divided into two categories: polynomial-based regression (PR) methods and neural network This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ (NN)-based methods. Works pertaining to both categories are presented following sequentially.
Fukuchi et al. [19] used a second-order polynomial to describe the relationship between a given walking speed and the knee gait trajectory. Their method is fast as the computation burden needed by a second-order polynomial for prediction is small. However, the prediction is less accurate due to the second-order polynomial is not enough to describe the complex curve of joint angle profiles. Moissenet et al. [20] proposed a multiple regression model that can estimate the relative contribution to lower limb joint profile of walking speed, gender, age and BMI. Yet, it is reported that a mismatch in predictors may lead to errors higher than 5 • on lower limb sagittal kinematics. Ren et al. [21] developed a personalized gait trajectory generation method based on Random Forest algorithm, the method modelling the relationship between fourteen subject-specific parameters (SSPs) and gait trajectory. Experimental results show that their method has a good prediction accuracy, nevertheless, the computational burden of that method is so heavy that it can not be used for real-time applications. Among all the work in the class of PR models, Embry et al. [22] proposed the most relevant work to this study. They developed a gait prediction model which is defined as the sum of the basis functions weighted by task function, and the task function accepts, phase variable, gait speed and slope as inputs. Experimental results demonstrated their method outperformed linear interpolation significantly when predicting the inter-subject mean joint kinematics. However, their method is prone to overfitting for it is not easy to set the form and the order of task function. In addition, their method cannot generate personalized joint trajectories depending on individual characteristics. Reznick et al. [23] proposed an approximate personalization method for joint trajectories by using the "difference" at the level ground task to approximate the "difference" at all other inclines, where the "difference" is computed by the same subject's experimental data and the inter-subject mean. They showed that the prediction accuracy of personal trajectory has been greatly improved after personalized. Nonetheless, because their method is personalized directly from the level of joint trajectory, thus it cannot be applied to online scenes. Also, their method is easy to mistakenly incorporate measurement noise into the individual's contribution.
Cunha et al. [24] used two neural networks (namely, the back propagation neural network (BNN) and the extreme learning machine (ELM)) to predict the personalized joint trajectory. The input of these NNs are four subject-specific parameters (SSPs) included the subject's age, weight, height and desired walking speed. However, four SSPs parameters selected base on priori may not represent all subject-specific gait properties. Liang et al. [25] proposed a method based on Long Short-Term Memory (LSTM) network to generate hip joint and knee joint gait profiles by exploiting inter-limb synergy. Whereas, due to the use of 21 SSPs, the proposed model required a large data and an ample amount of time and effort for angle profile prediction. He et al. [26] proposed an algorithm, which is based on extreme learning machine and AutoEncoder, to offer a natural and personalized gait trajectory for lower limb joint based on the body parameters. For the purpose of improving the efficiency and safety, a gait cell concept is proposed. Yet, like work in Liang et al. [25], He et al. [26] is also relatively slower methods and thus cannot be used for real-time applications. Bajpai and Joshi [27] proposed a deep neural network named as "MoveNet", which consists of three NN modules (encoder, mapper, and decoder), to predict joint profile across variable walking speeds and slopes. Like Reznick et al. [23], the method proposed by Bajpai and Joshi [27] can also predict subject-specific knee joint angle profiles for variable slopes from the data of knee joint angle profiles recorded at a flat surface. However, the input of MoveNet must be a complete joint profile, which limits its application in prosthetics. As the author points out, another limitation of MoveNet is that it can only be used to abled-bodied users. Yun et al. [28] proposed a method which model the relationship between fourteen body features and human gait trajectories by Gaussian process regression (GPR) algorithm. Although they showed that their method has a very good prediction performance, it will take more than one or two days to optimize the hyperparameters of GPR model.

B. Motivation and Contributions
After a comprehensive review of related work, although ADLs includes various combinations of walking speeds and slopes, we find that most of the previous studies neglected the demand of predicting joint trajectories for different slopes (except [22], [23], [28]). With regard to the subject-specific personalized trajectory prediction, most of these related work used handcrafted SSPs as the features to train their models. However, studies in [29] shown that the age, gender, BMI, walking speed, and other body parameters may not be enough to fully characterize the individual's gait. Besides, almost all related work based on a parametric model to build gait prediction model. In order to obtain good prediction accuracy and prevent overfitting, a large amount of training data is usually required. Nevertheless, recording gait data from amputee are expensive. Furthermore, for some work from the field of gait analysis, model complexity and computing time may not matter. But in the field of prosthetic control, the accuracy of model prediction and the rapidity of calculation must be balanced, so that the gait prediction model can be applied to real-time scenes.
Considering the limitations of previous methods, and the personalized and real-time requirements of prosthesis control for trajectory prediction under variable speeds and slopes tasks, This study presents a gaussian process enhanced Fourier series (GPEFS) model. Firstly, in order to obtain a low computational load and the convenience of control, gait trajectory is transformed into the finite Fourier coefficients (FFCs) vector space by the least-square method. Secondly, we use multiple GPRs models to fit the Fourier coefficients in FFCs vector space. As a nonparametric kernel-based probabilistic regression method, GPR can effectively avoid the overfitting problem. Hence, unlike work in [22], there no need to determine the form and order of the model according to the prior knowledge. Usually, the number of FFCs m is much smaller than the length of joint trajectory n (namely,m n), and fitting FFCs with multiple GPRs is equivalent to making multiple GPRs calculated in parallel. Therefore, the problem of high computational complexity of GPR (O N 3 ) [33] can be effectively avoided. Finally, in terms of trajectory personalization, we do not explicitly set SSPs, but believe that the difference between the level ground walking task data and statistical average (or inter-subject mean) data is enough to capture personal characteristics. And we operate that in FFCs vector space, thus it needn't to take a complete trajectory as input in an application. Thus, the proposed GPEFS model can be applicable to a real-time online scenes.
The rest of this paper is organized as following. The data set used to train the gait model are presented in Section II-A. The model format and the solving of Fourier series of the training data are presented in section II-B. Section II-C describes the model learning process with GPRs. And the Section II-D describes the method of individualizing. The numerical results are shown in Section III and are discussed in Section IV. Finally, the conclusion and future work is presented in Section V.

II. METHODS A. Data Set
Thanks for the open access dataset provided in [22] and [30], which greatly facilitates our research. The experimental protocol was approved by the Institutional Review Board at the University of Texas at Dallas (UTD). In this dataset, a total of 10 non-disabled subjects (5 female) participated in the experiment, and they are all provided written informed consent. The test method is that all subjects walking at a steady speed and incline on a Bertec instrumented split-belt treadmill for one minute, and the subjects' kinematics data are recorded by motion capture system; The specific test content for each subject including walking at a constant speed of 0.8m/s, 1.0m/s, and 1.2m/s with a constant ground slope ranging from −10 degrees to +10 degrees at 2.5 degree increments. There are 3 speeds and 9 slopes, resulting in 27 different tasks in the training dataset. We defined a task variable χ i ∈ R 2×1 made of speed v m ∈ R 1 with m = 1, 2, 3 and slope s n ∈ R 1 with n = 1, . . . , 9, and all observed tasks are aggregated into x ob , they are expressed as: where v m , as the subject's forward speed, is linearly mapped from the range of 0.6 m/s to 1.4 m/s to the range of 0 to 1; and s n , as the ground slope, is linearly mapped from the range of −10 degrees to 10 degrees to the range of 0 to 1. It is worth noting that χ i can be easily expanded to include other motion modes like a sit/stand or flat/stairs. In the data set, per stride was interpolated to contain N = 150 points, and all gait cycles begin at heel strike (corresponding to point 1), and end just before the next heel strike (corresponding to point 150).

B. Model Format and the Solving of Fourier Coefficient
In this study, gait model represents joint kinematics as a function of continuous phase variable ϕ ∈ R 1 and a continuous task variable χ ∈ R 2×1 . Where the phase variable ϕ represents the progress of gait cycle, and can be expressed as ϕ ∈ { R| 0 ϕ 1,φ > 0} [13]. The task variable χ represents the evolution of walking task. For humanity's gait kinematics has a continuous characteristic, χ not only includes the recorded task x ob but also includes tasks that beyond the fine granularity of record. In this study, we mainly consider the problem of interpolation estimation of joint profiles outside the trained data. So the continuous task variable χ can be expressed as |0 v 1, 0 s 1 , note that speed v and slope s are all continuous variable, and all have been normalized with the same range as v m and s n , respectively. The challenge is how to construct a unified model q(ϕ, χ) to predict the trajectories of lower limb joints under continuous changing tasks.
Previous study [22] provided a scheme that created a unified model based on the sum of K basis functions b k (ϕ) weighted by K task functions c k (χ): where E = 10 is usually selected to match the significant frequency content of human kinematics [32]. They defined c k (χ) as a form of Bernstein polynomials, and it's order γ and the form of function f (χ) needs to be determined based on domain knowledge or cross validation method. In contrast, this paper proposed a unified model that uses only one basis function. The central insight is that we can use the finite Fourier coefficients (FFCs) of the basis function to account for the trends of different tasks (instead of using task functions), and the relationship between FFCs and task variables under different tasks can be learned by machine learning algorithm. The prediction model is defined as: where E = 10 is also selected, and the Fourier basis vector b(ϕ) and FFCs vector A(χ) are defined as: One of the advantages of the proposed method is that, under the recorded task data x ob , the corresponding FFCs vector A(x ob ) can be solved easily by batch least-squares (BLS) or recursive least-squares (RLS) (for real-time online scenes) method. According to the solution process of BLS or RLS (it is not shown for the sake of brevity) [34], for a stride with N interpolation points in the task χ i , we can defined the Fourier basis matrix B(ϕ) ∈ R N×(2E+1 ) and the inter-subject mean joint kinematics vector Q(ϕ, χ i ) ∈ R N with the follow form: where q(ϕ k , χ i ) is the inter-subject mean joint kinematics recorded at discrete phase point ϕ k and discrete task vector χ i . Then FFCs corresponding to task χ i can be obtained by: For all recorded task x ob , the corresponding Fourier coefficient can be formed as a matrix: (6) So far, we are equivalent to transforming the joint trajectories into FFCs space, so the model learning can also be realized in this space.

C. Inter-Subject Mean Gait Model Learning With Gaussian Process Regression Enhanced Fourier Series
Based on the insight that we can use FFCs, which can be denoted as a function A(χ), to account for the trends of continuously varying tasks. Therefore, now the core problem is how to define a regression function F(χ) that approximates as accurately as possible the unknown FFCs function A(χ), where F(χ) is defined as: There are a multitude of machine learning algorithms can be used to model the function F. Yet, considering the training data is small-sized (only 27 recorded tasks), and gaussian process regression (GPR), as a nonparametric kernel-based probabilistic models, can automatically achieve a balance between fitting accuracy and model complexity to prevent overfitting of a small-sized data set [33]. Hence we selected GPR to model F.
We used 2E + 1 independent GPRs to model each dimension of F. For f κ (χ) with κ = 1, 2, . . . , 2E + 1, the corresponding GPR model using the training dataset D κ = {(χ i , A ob (κ, i ), i = 1, 2, . . . , 27} can be defined as: where μf κ (χ) is the mean function, and kf κ χ, χ is the kernel function (covariance function). We can query the GP at a new task input χ * as: The mean and variance predictions of that GPR are computed using a kernel vector k f κ = k f κ (D κ , χ * ), and a kernel matrix In this paper, we used the exponential kernel with automatic relevance determination [33]: where σ f is the signal standard deviation and σ l is the characteristic length scale, and they are usually be denoted as = [σ f , σ l ], which as the hyperparameters of GPR. The training process of GPR is using the recorded training data to determine the hyperparameters , by maximizing the log marginal likelihood function as [33]: where n = 27 is the number of the recorded tasks.
After concatenated all trained f κ (χ) into F(χ), then the inter-subject mean joint kinematics model can be reconstructed as: This is the unified model learned by GPR enhanced Fourier series. It can address the question of predicting the inter-subject mean joint kinematics for continuously varying task.

D. Individualizing Personal Gait Model Based on the Task Data of Flat Ground and Normal Speed
Inter-subject mean gait model accounts for trends well over continuous varying task, however, it may not be enough to fit a personal characteristic kinematics. For individual subject p with p = 1, 2, . . . , 10, the most direct method to get the personalized model q p (ϕ, χ) is using a complete experimental data of all tasks with the same process of obtaining the inter-subject mean model. However, it is often difficult to measure a large number of tasks over varying speeds and slopes due to the limitations of experimental equipment and site, especially for clinical related tasks. Therefore, how to use as limited task data as possible to individualize personal gait model has strong practical significance. In this paper, we proposed a method that modifying the inter-subject mean model based on only a single personalized task of flat ground and normal speed, so as to obtain the individual gait model across all tasks.
We selected the normal speed (v = 1.0m/s) and flat ground (s = 0 degree) as the basis task (denoted as χ base ), which is the task an individual subject need to carry out. We defined the FFCs of the subject p corresponding the basis task as A p (χ base ), and the FFCs of subject p of task χ i is denoted as A p (χ i ) wiht i = 1, 2, . . . , 27. Noting that in practical application we do not need to calculate A p (χ i ) (except A p (χ base )) for they are only to facilitate analysis and description. The rationality behind the personalization method proposed in this paper is that, as showed in Fig. 4, we found that the difference FFCs between individual trajectory and mean trajectory follows the same trend in all tasks. Therefore, we can use the coefficient error between individual trajectory and mean trajectory in the base task χ base to approximate the coefficient error in other tasks, and they can be formed as: we considered the A p, (χ base ) as a measure of personal characteristics, and the A p, (χ j ) can be approximated by A p, (χ base ). Therefore, the FFCs of individual p in all tasks can be approximated as: After obtaining the Fourier coefficient A p (χ j ), we can fit the personalized model use GPR with the same process as obtaining the inter-subject mean model.

III. RESULTS
In this section, we will first model the inter-subject mean kinematics of hip, knee and ankle, and evaluate the ability of the unified mean model learned by GPR enhanced Fourier series to predict the kinematics of untrained tasks. Secondly, we will evaluate the feasibility and performance of the proposed personalization method. We used the percent gait to act as the phase variable, which are defined as:ϕ = t T , ϕ i = t i T , where T is the total time of a gait cycle, t is the continuous time, and t i is the discrete time.  Fig.1, we can see that the Fourier coefficients playing a major role in hip, knee and ankle are "w 0 ∼ w 4 ", "w 0 ∼ w 6 " and "w 0 ∼ w 10 " respectively, this is consistent with the complexity of the joint trajectory. Fig.2 shows the fitting effect of the GPR model on Fourier coefficients w 0 , w 1 and w 2 of hip, knee and ankle in all record tasks, where the "Target Value" is the Fourier coefficients solved by BLS. Due to the limited space of the paper, the fitting effect of other Fourier coefficients has not been shown. Table I shows the fitting-loss of GPR model in fitting Fourier coefficients "w 0 ∼ w 20 " of hip, knee and ankle, where the fitting-loss is computed by:

A. Performance of the GPR Enhanced Fourier Series in Modeling the Inter-Subject Mean Joint Kinematics
where subscripts i = 0, 1, . . . , 20 and j represent FFCs serial number and task serial number, respectively. From Fig.1 (as the target value solved by BLS) and Table I (as the fittingloss), we can see that FFCs are fitted well by GPRs in all tasks. Fig.3 shows the example of inter-subject mean model joint surfaces of hip, knee and ankle at normal speed (1.0m/s). Where the model outputs are represented as blue surface with a bright yellow gradient, and the inter-subject mean data used to train the model are shown in dotted cyan lines. We can note that our proposed method can not only fit the training data well, but also have good prediction ability and generalization ability in untrained areas (regions between dotted cyan lines denoting the training data). In addition, noting that the mean model learned by GPR enhanced Fourier series is C ∞ smooth, this is different from the surface obtained by linear interpolation, which is reported in [22] being prone to un-smooth when predicting untrained regions.

B. Feasibility Verification and Performance Evaluation of the Individualize Method
To make the paper be more concise, in this section, we chose subject 5 as the representative for discussion. Joint kinematics trajectories of hip, knee and ankle at normal speed (1.0m/s) of subject 5 are denoted by the solid red lines in Fig.3. From Fig.3 we can see that the mean model can account for trends well over continuous varying task, but it may not be enough to fit a personal characteristic kinematics.  in experiments with the solid red lines. By comparing Fig.5 and Fig.3, we can intuitively see that the individualized model can better fit the joint kinematics of individual. Noting that we have not made special adjustments except applied personal gait data at normal speed on flat ground. Table II shows the root mean squared error (RMSE) and max error across all tasks of inter-subject mean model and individualized model for all subjects and joints. The formula of RMSE is defined as: where N = 150 is the number of interpolated point in a gait period, and T = 27 is the number of the recorded tasks. Averaged across all subjects, after individualization, the reduced RMSE/Max-Error at the hip, knee, and ankle joint are 2.94/5.059, 1.822/5.82 and 1.156/4.51, respectively. Fig.6 shows the statistical information corresponding to Table II. We statistically compared the gait model before and after individualized to check whether the improvement is statistically significant. Since the number of recorded tasks is only 27, it may not satisfy the homogeneity of variance, so a two-tailed Mann Whitney test using an alpha of 0.05 was adopted. As showing in figure.6, the RMSE and Max-Error before and after individualized in all subjects for all joints  II  RMSE AND MAX ERROR OF INTER-SUBJECT MEAN MODEL AND INDIVIDUALIZED MODEL ACROSS ALL TASKS FOR ALL SUBJECTS AND JOINTS.  THE LAST TWO COLUMNS SHOW THE MEAN  are significantly improved. The p-value of RMSE/Max-Error of hip, knee and ankle are p = 0.0011/p = 0.0185, p = 0.0147/p = 0.0068 and p = 00186/p = 0012, respectively.

IV. DISCUSSION
Human beings are born to walk naturally in a continuous way. In addition to walking with the most comfortable speed on the flat ground, we often face the task of variable speed and variable slope in ADLs. However, when designing a continuous prosthetic controller, it is a challenge to build a gait prediction model that can accurately and continuously predict the trajectories of lower limb joints across various slopes and walking speeds. This study aims to contribute to this topic.

A. Advantages of GPR Enhanced Fourier Series Modeling
For highlighting the advantages of our proposed model, as shown in Fig.7, we compared the Mean/Maximum prediction error values of the proposed method (namely, GP enhanced Fourier series (abbreviated as GF)) with the ones of the state-of-the-art method proposed in [22] (namely, "basis function weighted by task function" (BM)) and a traditional linear interpolation (LN) method. The definition and calculation formula of prediction error can be referred to equation (10) of [22]. It is clearly shown that both the mean and maximum prediction errors in our proposed method (GF) are the smallest. Specifically, in hip, knee and ankle joint, the mean prediction errors in GF are only 46.5%/44.2%, 21.7%/17.5%, and 15.7%/13.5% of that in BM/LN, respectively; And the maximum prediction errors in GF are 17.6%/23.3%, 15.0%/8.8%, and 10.5%/8.8% of that in BM/LN, respectively.
In addition to good prediction accuracy, our method is also beneficial from the ease of implementation and a low computational burden. Previous work in [22] defined a task function c k (χ) (as shown in equation (2)) as a form of Bernstein polynomials to account the influence of different tasks on joint trajectory. However, how to set the order γ and the form of function f (χ) of the task function is not easy. It is particularly important to have an appropriate form of task function. If the task function c k (χ) is too simple, the finally prediction model may fail to fit the true characteristics of the data under variable speeds and variable slops; If the task function c k (χ) is too complex (e.g., the order (γ ) is too high or the function f (χ) is too complex), the finally prediction model may be over fitted and the generalization ability may be greatly reduced. In the proposed method, the only parameter that needs to be determined according to prior knowledge is the order of the Fourier series E. Fortunately, human joint profile has a relatively fixed range or shape in daily life, so when the order of Fourier coefficient is E = 10 [32], it can capture the main features of human kinematics. In addition, GPR is a fitting method that lets the data speak for itself, so it can effectively prevent over fitting on small data sets.
Previous work in [15] and [28] directly used GPR to construct gait prediction model. Their models have achieved good performance, but the computation burden is so heavy that they cannot be used for real-time scene. For example, it was reported in [28] that the time took to optimize the hyperparameters of their GPR based model more than one or two days. This was caused by the biggest problem of   Fourier series (abbreviated as GF ), the "basis function weighted by task function" (BM) method proposed in [22], and the linear interpolation (LN) method. Please refer to equation (10) of [22] for the specific formula and definition.
GPR that its computational complexity is O N 3 for the time series with length N. In this study, for the inter-subject mean data, there 27 task data and the length of each task data is 150. If we directly use GPR to build a gait prediction model, the computational complexity is N 4050 3 . That may be unbearable for real-time application. A novelty of this paper is that we first transformed the joint trajectory into Fourier coefficient space. Thus, we can use GPRs to learn the target Fourier coefficients in the Fourier coefficient space. Furthermore, we created multiple GPR models (can optionally be computed by parallel) and each of them are corresponding to each of the Fourier coefficients. Specifically, for the inter-subject mean model, the computational complexity is 21 * O 27 3 . It can be seen that the computational complexity has been greatly reduced. Besides, the target value of Fourier coefficients can be solved directly by the least square method framework rather than the optimization method, therefore, the solution process is greatly simplified. Moreover, the solution process is very robust owing to the condition number of B(ϕ), being 1.42, is very small (the closer the condition number to 1, the more robust of the solution of the least squares method). Where the the condition number is defined as Cond = λ max (B(ϕ)) λ min (B(ϕ)) , and λ max (B(ϕ)) and λ min (B(ϕ)) represent the maximum and minimum eigenvalues of matrix B(ϕ) respectively.

B. The Performance of the Personalized Gait Model
In the prosthetic control, the average model trained by inter-subject mean data can make the amputee walk normally, however, it is usually suboptimal for that the personal characteristics of the wearer are ignored. A better method is to train a personalized trajectory prediction model suitable for personal characteristics according to personal walking data.
As showed in Table III, we compared the proposed individualizing method with eight state-of-the-art methods. From the comparison results, we can see that the personalized trajectory prediction performance of [19], [20], and [23] is relatively poor due to the limited expressive ability of polynomials. Methods in [21], [25], [26], and [24] achieved excellent personalized trajectory prediction performance, however, they didn't involve the task of variable slopes. References [23] and [27] and this study are the only three methods involving variable speed and variable slope tasks, and more interestingly, all three methods follow the similar idea that believing the difference between the level ground walking task data and statistical average (or inter-subject mean) data is enough to capture personal characteristics, but the implementation strategies of the three methods are completely different. Method proposed in [27] regards a complete gait periodic trajectory as an onedimensional "image", and obtains excellent personalized trajectory prediction performance with the help of a convolutional neural network. However, for its input must be joint profiles, so it can only be utilized to off-line gait analysis rather than the real-time application scenario of prosthetic control. Method proposed in [23] directly obtain the so-called "personal contribution" by calculating the joint trajectories difference between all recorded tasks data with inter-subject mean data. Nevertheless, it is easy to mistakenly incorporate measurement noise into the individual's contribution. That may be the reason for its relatively poor prediction performance.
As showed in Fig.4, in this study we found that the Fourier coefficient errors between all recorded tasks data and statistical average (or inter-subject mean) data follow the same trend. Therefore, we can construct personalized methods directly in Fourier coefficient space (as shown in equation (14)). One of the advantages is that we can only perform this operation on the dominant Fourier coefficients (i.e. Fourier coefficients with small serial number), which can effectively prevent the negative impact of high-frequency measurement noise. In addition, our method does not need complete joint contour as input, so it can be applied to real-time application scenarios such as prosthetics.

V. CONCLUSION
In conclusion, this paper proposed a gait prediction model which models the human locomotion as a continuous function of phase, speed and slope. This model makes it feasible to realize a continuous controller for powered prosthesis in a wider range of tasks. The main contribution of this study is that we present a formal modeling method. The joint trajectories under different tasks are first transformed into Fourier coefficients space; then we can use GPR to fit the relationship between Fourier coefficient and joint trajectory. GPR can automatically balance the fitting performance and model complexity, so its use can avoid the model from over fitting on the sparse training data set. Compared with the complexity N 4050 3 of learning trajectory directly using GPR for 27 recorded gait data with a length of 150, the computational complexity of our method is only 21 * O 27 3 (where 21 is the selected order of Fourier series).
In addition, we implicitly capture the characteristics of individual gait trajectory in the Fourier coefficient space. The individualizing method can model individual kinematics across tasks based only on personalized data for level-ground walking at a normal speed, so it can greatly reduce the time it takes to tune multiple tasks. Applied to 10 human subjects, the individualization method reduced the RMSE between the model and subject's kinematics over all tasks by an average of 2.943 deg (max 5.059 deg) at the hip, 1.822 deg (max 5.82 deg) at the hip, and 1.156 deg (max 4.51 deg) at the ankle. The individualization method makes the powered prosthesis be more clinically viable over a variety of tasks without the need of an engineer.
A future work is to study the real-time estimation algorithm of walking speed and slope based on IMU, and combine them into the gait prediction model. Another work is to extend the exist personalized method to include tasks such as climbing stairs of different inclinations and sit-to-stand.