—Objective: This paper describes the development and preliminary ofﬂine validation of an algorithm facilitating automatic, self-contained learning of ground terrain transitions in a lower limb prosthesis. This method allows for continuous, in-ﬁeld convergence on an optimal terrain prediction accuracy for a given walking condition, and is thus not limited by the speciﬁc conditions and limited sample size of an in-lab training scheme. Methods: We asked one subject with a below- knee amputation to traverse level ground, stairs, and ramps using a high-range-of-motion powered prosthesis while internal sensor data were remotely logged. We then used these data to develop a dynamic classiﬁcation algorithm which predicts the terrain of each stride and then continuously updates the predictor using both data from the previous stride and an accurate terrain back-estimation algorithm. Results: Across 100 simulations randomizing stride order, our method attained a mean next-stride prediction accuracy of 96% . This value was ﬁrst reached after 200 strides, or about 5minutes of walking. Conclusion and signiﬁcance: These results demon- strate a method for automatically learning the gait patterns preceding terrain transitions in a prosthesis without relying on any external devices. By virtue of its dynamic learning scheme, application of this method in real-time would allow for continuous, in-ﬁeld optimization of prediction accuracy across a variety of walking variables including physiological conditions, variable terrain geometries, control methodologies, and users. I. INTRODUCTION With the advent of powered lower limb prostheses, in re- cent years there has been considerable interest in developing suitable control algorithms facilitating efﬁcient, comfortable, and biomimetic gait for people with lower limb amputations. In particular, a large body of work has focused on developing algorithms for the anticipation and adaptation to different walking terrains such as level ground, ramps, and stairs [1]– [15]. Studies have explored a multitude of machine learning algorithms and sensing modalities, many reporting a high degree of success in a laboratory environment, and almost always on a small cohort of subjects. Training procedures for developing these machine learning algorithms are often time consuming. For a given set of sensors and a powered prosthesis platform, studies involve logging data, often over multiple days, from multiple sub- jects traversing a given terrain. This is followed by manual ofﬂine terrain labeling and a pattern recognition analysis. For real-time control tasks, ofﬂine algorithms must then be translated into embedded languages, with no guarantee This work was supported by the National Defense Science and Engineer- ing Graduate Fellowship, Defense Advanced Research Projects Agency, and the MIT Media Lab Consortia. 1MIT Media Lab - Biomechatronics Group, Massachusetts Institute of Technology, Cambridge, MA USA These authors contributed equally to this work.that the predictive models will perform as well as they did in simulation when subjected to new terrains, users, control methodologies, prosthesis platforms, or physiological conditions – it is impractical to include all such possible conditions in one training set. Some attempts have been made to automate aspects of the machine learning process. Zhang et. al. used an external system during data collection to automatically identify terrain with high accuracy and train a machine learning model online, without the overhead of ofﬂine labeling and analysis [17]. However, the main drawback of this approach was the requirement for an external system for labeling, making re- training in alternate conditions burdensome and impractical. Spanias et. al. used an adaptive EMG-based machine learning model to compensate for EMG disturbances [10], but this system did not include a method for labeling the terrain and thus did not directly optimize for prediction accuracy. Additionally, the part of the model employing intrinsic sensor signals was static. In this work, we address the problem of non-generalized models and burdensome training routines by developing a method to automatically and continuously train a pattern recognition algorithm using only the intrinsic sensors on- board a below-knee prosthesis. This model relies on an accurate back-estimation step, which uses a heuristic to label strides after they have been taken, and thus continually updates the predictor with new training data. The back- estimation step was enabled largely by a novel high range- of-motion (ROM) prosthesis which allows the use of ankle angle during stance to distinguish inclined terrains from level terrains. Currently, the only commercially available powered below-knee prosthesis, the EmPower by Ottobock, has a lim- ited 22 degree ROM with zero degrees of dorsiﬂexion [19], [20], which is insufﬁcient to span even the biological range of level-ground walking (10 degrees dorsiﬂexion to 18 degrees plantarﬂexion [18]), let alone alternative terrains. Similarly, most powered prostheses in the research environment are designed to operate on level-ground [21]–[24]. The novel mechanical system used for this study enables 115 degrees ROM, spanning the entire mean biological ROM and conse- quently allowing for signiﬁcant biomechanical differentiation between ground terrains. By employing an incremental learning algorithm leverag- ing the backward estimation of terrain labels and a high- ROM prosthesis, we achieve a ﬁeld-usable automatic training method that requires no manual processing steps or external devices. This method would enable a powered prosthesis to automatically update a customized terrain predictor that continuously converges on the optimal prediction accuracy for a given walking condition. This method is computation- ally efﬁcient, employing a physics-based heuristic for back estimation and an order O(f2)incremental learning step once per stride, where fis the number of model features. Finally, while the method was developed on training data obtained from a transtibial prosthesis, the method is also applicable in the transfemoral case. II. METHODS A. Data collection 1) Overview: We asked a subject with a unilateral transtibial amputation to don a novel powered lower limb prosthesis (described below) and traverse various terrains including level ground, stairs, and ramps. While the subject was walking we remotely logged data from internal pros- thesis sensors, including ankle angle, ankle torque, and raw inertial measurements from a six degree-of-freedom (three accelerometers, three gyroscopes) inertial measurement unit (IMU), and ﬁlmed the subject to facilitate manual labeling of terrains to provide a ground truth terrain identity. We developed and tested our method in an ofﬂine simulation using the data collected during these trials. 2) Prosthesis: The prosthesis, referred to as TF8 and ﬁrst described in [16], is a new system built to achieve biological kinetics and kinematics that enable operation over a range of terrain conditions. The TF8 is a torque controlled powered prosthesis designed around a series elastic actuator that can provide peak torques up to 180 Nm across a 115 degree total operational ROM. The system architecture, shown in Figure 1, consists of a large gap radius motor (manufactured by T- Motor) modiﬁed to integrate a ballscrew into the rotor. The ballscrew applies a linear force to an output moment arm that generates a torque about the ankle joint. An axial load cell directly measures the force in the screw. This force signal is evaluated along with the joint encoder measurements to determine the effective joint torque with an accuracy of 0:5Nm. The joint encoder is a 14-bit absolute encoder, AS5048 (manufactured by Austria Microsystems). Inertial measurements are performed by the motion tracking MPU- 9250 (InvenSense) included in the control unit printed circuit board assembly. The control unit consists of a customized embedded system platform based on the FlexSEA system designed by Dephy, Inc. The system includes a motor driver and a separate mid-level control system that is based on the STM32F427 32-bit Cortex M-4 microcontroller operating at 180 Mhz. A control loop running at 1 kHz reads sensor values, and runs a state-machine that deﬁnes the desired operating condition based on evaluated system parameters. A state deﬁnes a desired joint torque utilizing the impedance control parameter =K(md) +B_deﬁned by Hogan [25]. A closed-loop torque controller then asserts a torque on the joint by converting the torque command to a desired motor current that the FlexSEA motor drivers internal current controller commands at the motor. Bidirectional data is communicated to a high-level controller across a wireless Joint EncoderControl Unit Force SensorIMU (a) Mid-Lev el Controller Impedance ControllerSensors Control Unit Mechanical SystemJoint Encoder Force SensorExternal Peripherals High-Lev el Controller IMU Torque Controller Physical Actuator Motor Driver (b) (c) Fig. 1: (a) The TF8 mechatronic system architecture is a reaction force series elastic actuator with an on-board em- bedded control system. (b) The Control Unit is a derivative of the FlexSEA embedded system from Dephy, Inc. that includes a motor driver unit and mid-level controller. A state machine runs on the mid-level controller that deﬁnes a desired impedance command and runs a closed-loop torque controller to deﬁne the behavior of the physical actuator. (c) The TF8 has 115 degree total ROM with 35 degrees of dorsiﬂexion. t=tstance(n)t=tswing(n)t=tswing+tcutoff(n)t=tstance(n+1)t=tswing(n+1)Fig. 2: Stride time point snapshots. Stride objects include both initial and target stance periods to allow back-estimation of stride terrain using data from the target stance period. Initialize Predictor Predict Back-estimateUpdate Predictorstride 0datapredictor label predictor features back-estimated labelPredict Back-estimateUpdate Predictorstride 1datapredictor label predictor features back-estimated labelpredictor parameters predictor parameters Predict Back-estimateUpdate Predictorstride n datapredictor label predictor features back-estimated labelpredictor parameters Fig. 3: Overall architecture of the incremental learning algorithm. An LDA classiﬁer is initialized with zero means and an identify covariance matrix. Subsequently, each new stride undergoes a prediction and back-estimation step. In prediction, the current LDA classiﬁer is applied to predict the next terrain. In back-estimation, the stride is labeled post completion and the LDA classiﬁer parameters are updated with the new training data. Bluetooth protocol. These experiments utilized a laptop com- puter to log system state data through the wireless channel. To enable terrain agnostic data collection a single state control scheme was implemented with a saturation threshold on joint torque. An impedance controller was set to a joint stiffness of 3 Nm/deg and 0.2 Nm s/deg damping. When the joint anglemreached beyond a threshold of 10 degrees in either dorsiﬂexion or plantarﬂexion direction the mvalue was set to this threshold value, th. This saturation provided stability for the subject while also enabling the joint angle to adapt to the terrain. B. Ofﬂine processing 1) Initial processing: We divided all data into individual strides using thresholds on ankle torque and timers for swing and stance phases. Individual strides were deﬁned by the period from the beginning of one stance period to the end of the subsequent stance period, so as to include enough information for the back-estimation of stride terrains. The resulting stride list thus contained all strides n= 1 :Nwith time points t(n) stance ,t(n) swing ,t(n+1) stance , andt(n+1) swing (illustrated in Figure 2), ankle angle, ankle torque, and six IMU signals in the window [t(n) stance;t(n+1) swing ]. Finally, we manually labeled all strides using trial videos as ground truth evaluation of the back-estimator. 2) Signal extraction: We were then interested in ex- panding the set of signals upon which to perform pattern recognition. This was shown to be advantageous in our previous work [13]. In particular, we extracted ﬁrst-orderintegrals and derivatives of all IMU signals for each stride, producing 12 additional signals. All derivatives were ﬁltered using a ﬁrst-order recursive low-pass ﬁlter. Next, we used the IMU signals to extract inverse kinematics of the ankle joint, using a method similar to that employed in [13]. This step generated signals for shank pitch and its two-dimensional sagittal-plane accelerations, velocities, and positions. 3) Feature extraction: For each stride, we extracted the mean, range, maximum, minimum, time of maximum, time of minimum, and ﬁnal value of all IMU-derived signals (6 raw, 12 integrals and derivatives, 7 inverse kinemat- ics) as well as ankle angle and torque in the window [t(n) static;t(n) swing +tcutoff ], wheret(n) static was a zero velocity time point identiﬁed in the window [t(n) stance;t(n) swing ]using the inverse kinematics algorithm and tcutoff was a deﬁned time aftert(n) swing chosen to cutoff data collection for the current gait cycle in order to allow enough time for prosthesis actuation in the case of real-time adaptive control. 4) Incremental learning simulation: We performed 100 incremental learning simulations, each time randomly order- ing the stride list so as to simulate receiving training data in an arbitrary order. For each simulation, we employed an algorithm for incrementally training a linear discriminant analysis (LDA) classiﬁer by updating intra-class and inter- class feature means and an interclass feature covariance matrix at every stride. All means were initialized at zero, and the covariance matrix was initialized as an identity matrix. The overall architecture of the incremental learning algo- US DS UR DR F/uni0394 >/uni0394 /uni0394 </uni0394 /uni0394 (+1)>,(+1)<, (+1)(a) Δ𝑧𝑛𝜃 𝜃 t (b) Fig. 4: (a) Heuristic back-estimation of stride terrain based on vertical ankle joint position z(n)and mean joint posi- tion(n+1) in[t(n+1) stance;t(n+1) swing ]. (b) A visualization of these parameters for a representative terrain transition. rithm is illustrated in Figure 3. In particular, the following steps were performed sequentially on every stride: 1) Stride classiﬁcation using the extracted feature set for prediction and LDA classiﬁer with current class means and feature covariance matrix. This algorithm can be described as: (n) k(~ xn) =~ xT n1~ k1 2~ T k1~ k+log(k)(1) 8k2CwhereCis the ﬁve-element list of all candidate terrains, ~ xnrepresents the feature vector calculated on stride n,1represents the feature covariance matrix inverse, ~ krepresents feature means for each of the possible terrains k2C, andkis the prior probability of terrain k. For each stride, the predicted terrain ^k(n)was calculated using: ^k(n)(~ xn) = argmax k2C(n) k(~ xn) (2) 2) Back-estimation of the stride terrain k=~k(n)used a heuristic, physics-based algorithm following the simple decision tree depicted in Figure 4. This algorithm relied predominantly on two signals, namely the anklevertical position at t(n+1) stance and the mean ankle angle in the window [t(n+1) stance;t(n+1) swing ]. 3) Update of class population pkand~ kusing~ xnand the back estimated stride label k: p(n) k=p(n1) k+ 1 (3) ~ (n) k=(p(n) k1)~ (n1) k+~ xn p(n) k(4) 4) Update of within-class covariance matrix using~ xn. This method for incrementally updating a covariance matrix is known as Welford’s algorithm: ~ (n)=(n1)~ (n1)+~ xn n(5) (n)=(n1)(n1)+ (~ xn~ (n))(~ xn~ (n1))T n (6) 5) Analysis: Across simulations, we calculated the mean next-stride prediction accuracy. Additionally, for each indi- vidual simulation we calculated the accuracy of the last 100 strides. III. RESULTS A. Initial stride list The subject took a total of 551 strides, including 329 level ground, 29 ascending a ramp, 35 descending a ramp, 82 ascending stairs, and 76 descending stairs. The back- estimation heuristic correctly identiﬁed >99% of strides. Mean prediction accuracy for the nth stride across 100 simulations is shown in Figure 5. The simulation includes only 482 strides because 69 strides (all of which were made on level ground) were automatically labeled as level because their swing periods were too short (that is, t(n+1) stancet(n) swing< Tcutoff ). Mean accuracy for the past 100 strides across ﬁve representative simulations are shown in Figure 6. Finally, an LDA model trained on all strides achieves a back-substitution accuracy of 98:1%. IV. DISCUSSION Our results indicate that accurate back-estimation and incremental learning can be used to gradually improve the accuracy of a prosthetic terrain prediction system, without the need for an ofﬂine training process. If employed on prosthetic ﬁrmware or on a device that is frequently connected to a prosthesis, the algorithm can be used to continuously update a dynamic machine learning model such that it will provide maximum accuracy for a given walking condition. Figure 5 demonstrates that the incremental LDA learner approaches an optimum accuracy over time, across arbitrary stride orders. The time constant to optimum accuracy appears to be approximately 200 strides. Given an average of 1.5 seconds per stride in our data set, this equates to about 300 seconds or 5 minutes to reach a high accuracy. Beyond this point, the accuracy continues to grow but more slowly, gaining about 2% on average between 200 strides and 450 strides. Fig. 5: Mean prediction accuracy of nth stride among 100 simulations. Standard deviations of accuracy for all time points were below 0.5%. Fig. 6: Mean prediction accuracy of last 100 strides for ﬁve representative simulations. The learning does not necessarily occur monotonically; Figure 5 shows a local minimum in accuracy at approxi- mately 50 strides, after a peak of about 80%. The reasons for this dip are unclear, but it may have to do with the preponderance of level ground strides in the data set. This would lead to situations in which substantial learning has been completed on level ground strides before many strides of other classes are encountered. This result would be veriﬁed with additional subject testing. Of note in Figure 6 is the variable rates of learning seen among simulations. We see that learning in any particular simulation is almost never monotonic, exhibiting potentially several substantial local minima on the way to the ﬁnal accuracy and in some cases (green) not improving over 200 strides before quickly catching up. Beyond the variability, this ﬁgure also suggests that we are training data limited. While they exhibit local minima, we see that all simulations do trend upward in their accuracy; it is likely that, given enough training examples, all simulations would reach the optimum achieved in Figure 5 and local ﬂuctuations would become less severe. Finally, the high back-substitution accuracy ( 98:1%) of anLDA model trained on all strides provides further evidence that it is possible to attain a signiﬁcantly higher accuracy than that in Figure 5, given a greater amount of training data. Perhaps the main reason we do not achieve the higher accuracy is the covariance update (Welford”s) algorithm, which only approximates the next covariance matrix rather than measuring it exactly. We approximated this error by ﬁnding the difference between the norms of the approximated covariance matrix and the true covariance matrix at different time points in a simulation, and found that the difference appears to grow approximately linearly over time. Therefore, one aspect of future work would likely involve correcting the covariance matrix. This could be done by periodically directly measuring the covariance for the last nstrides using suitable processing and data storage systems. A. Future work Avenues for future work include the addition of new intrinsic sensing modalities to further improve prediction accuracy, and incorporation of this algorithm into an em- bedded microcontroller for in-ﬁeld learning and analysis. Additionally, beyond automatic updates of means and co- variances, it is useful to determine methods of automatically updating the feature set. In particular, it is possible that providing the ability to select different features as input in a predictor depending on the situation could further improve capability. Finally, it would generally be useful to increase the amount of training data and incorporate intentionally modiﬁed walking conditions so as to study the behavior of the algorithm over time. ACKNOWLEDGMENTS The authors would like to acknowledge Dephy for con- tributing to the design of the TF8 embedded systems, Tony Shu for his low level software contributions, and the subject who participated in the study.