Inertial Measurement Unit-Based Estimation of Foot Trajectory for Clinical Gait Analysis

Gait analysis is used widely in clinical practice to evaluate abnormal gait caused by disease. Conventionally, medical professionals use motion capture systems or make visual observations to evaluate a patient's gait. Recent biomedical engineering studies have proposed easy-to-use gait analysis methods employing wearable sensors with inertial measurement units (IMUs). IMUs placed on the shanks just above the ankles allow for long-term gait monitoring because the participant can walk with or without shoes during the analysis. To the knowledge of the authors, no IMU-based gait analysis method has been reported that estimates stride length, gait speed, stride duration, stance duration, and swing duration simultaneously. In the present study, we tested a proposed gait analysis method that uses IMUs attached on the shanks to estimate foot trajectory and temporal gait parameters. Our proposed method comprises two steps: stepwise dissociation of continuous gait data into multiple steps and three-dimensional trajectory estimation from data obtained from accelerometers and gyroscopes. We evaluated this proposed method by analyzing the gait of 19 able-bodied participants (mean age 23.9 years, 9 men and 10 women). Wearable sensors were attached on the participants' shanks, and we measured three-axis acceleration and three-axis angular velocity with the sensors to estimate foot trajectory during walking. We compared gait parameters estimated from the foot trajectory obtained with the proposed method and those measured with a motion capture system. Mean accuracy (± standard deviation) was 0.054 ± 0.031 m for stride length, 0.034 ± 0.039 m/s for gait speed, 0.002 ± 0.020 s for stride duration, 0.000 ± 0.017 s for stance duration, and 0.002 ± 0.024 s for swing duration. These results suggest that the proposed method is suitable for gait analysis, whereas there is a room for improvement of its accuracy and further development of this IMU-based gait analysis method will enable us to use such systems for clinical gait analysis.


INTRODUCTION
Analysis of abnormal gait can provide important information about diseases and injuries. For example, patients with Parkinson's disease (PD) often exhibit shuffling, festinating, and freezing of gait. The most widely used clinical rating scale for PD, the Unified Parkinson's Disease Rating Scale, includes observation of gait (Goetz et al., 2008). Patients with cerebellar disorders sometimes have a wide-based (atactic) gait, and those with cerebral vascular disease sometimes exhibit a hemiplegic gait. Recent articles have reported changes in gait, such as reduced gait velocity and stride length, in diseases with gait disorders and in other conditions, such as Alzheimer's disease (Mielke et al., 2013) and depression (Lemke et al., 2000).
Clinical gait analysis is performed mostly by health-care providers using visual observation (Krebs et al., 1985). Although this method is the most readily accessible means of gait analysis available to health-care providers (Barker et al., 2006), it is a subjective and qualitative method that is inadequate for assessing changes in gait features during ongoing treatment interventions. It is also difficult for clinicians to share this information with health-care providers and patients. Motion capture systems are used in clinical research for gait analysis (Mcginley et al., 2009) and scientific research (Lieberman et al., 2010). Because they provide well-quantified and accurate results, these systems are currently considered to be the criterion standard for clinical gait analysis (Cameron and Wagner, 2011). However, because the special equipment needed for motion capture is expensive and requires a large space, few medical institutions can use these systems for clinical gait analysis (Barker et al., 2006).
Several studies have proposed gait analysis methods using inertial measurement units (IMUs) to solve the problems described above (Sabatini et al., 2005;Moore et al., 2007;Mariani et al., 2010;Rebula et al., 2013;Kitagawa and Ogihara, 2016). IMUs used in these methods are inexpensive and wearable (Fujiki et al., 2009). In particular, we focused on methods that estimate trajectories of a foot because such methods can be used to obtain several spatial gait parameters. Sabatini et al. (2005) proposed an IMU-based gait analysis method that estimates a two-dimensional trajectory in the sagittal plane of a foot during walking. Other studies have proposed gait analysis methods that estimate the three-dimensional foot trajectory during walking in a stepwise manner to obtain values of foot clearance (Mariani et al., 2010(Mariani et al., , 2012Kitagawa and Ogihara, 2016). The trajectory estimation methods reported in several studies (Sabatini et al., 2005;Mariani et al., 2012;Kitagawa and Ogihara, 2016) use an IMU attached on the dorsum of the foot and are better for obtaining this gait feature. As described above, and to the knowledge of the authors, there is no report of a method that estimates three-dimensional foot trajectory from an IMU attached on the shank during waking in a stepwise manner to calculate simultaneously spatial and temporal clinical gait parameters, including stride length, gait speed, stride duration, stance duration, and swing duration.
In this study, we propose a novel gait analysis method for clinical purposes that uses IMUs attached on the shanks to estimate foot trajectory and to obtain estimated clinical gait parameters. We conducted an experimental evaluation of the proposed method. Here, we report an example of the application of the proposed method to PD patients.

PROPOSED METHOD Sensors Used and Wearing Method
Our proposed gait analysis system is illustrated in Figure 1A. For gait analysis, we used two IMUs (TSND121, ATR-Promotions, Kyoto, Japan; Figure 1B) with a triaxial accelerometer (±8 G range), triaxial gyroscope (±1,000 • per range), and Android OS tablet (ZenPad10, ASUSTeK Computer Inc., Taipei, Taiwan; Figure 1C). Raw accelerometer and gyroscope signals were sampled at 100 Hz (16 bits per sample). The size of the IMU is 37 mm × 46 mm × 12 mm and its weight is about 22 g. IMUs are attached on the shanks (just above the ankles) with bands ( Figure 1B). The inertial coordinate system used to represent foot orientation and position is shown in Figure 1B. Acceleration and angular velocity data of both the shanks measured during walking are transmitted to the tablet through Bluetooth.

Algorithm for Trajectory Estimation
Our proposed method comprises two steps: dissociation of continuous gait data into multiple steps and three-dimensional trajectory estimation in a stepwise manner. Each process is described as follows.

Stepwise Dissociation From Angular Velocity Signals
We dissociated the signals into four steps: (1) smoothing and finding the local maximums; (2) finding the heel-strike (HS) and toe-off (TO) points; (3) quadratic regression; and (4) calculating the split point. One walking cycle is defined as a single step, and the starting point of each cycle is defined as the split point between HS and immediate TO points. We identify the split point in each cycle based on raw angular velocity data in the z-axis ω z of the gyroscope (right-handed; Figure 1B, left).

Smoothing and Finding the Local Maximums
The raw data contained much noise, and a median filter (window length: 5) was first used to smooth the data. We found the local maximum l k that is larger than the threshold (200 • /s; Figure 2A).

Finding the Local Minimums Near the HS and TO Points
We then found the k-th local minimums near HS mhs k and TO mto k points between the l k and l k+1 . mhs k is the local minimum closest to the right of l k , and mto k is the local minimum closest to the left of l k+1 ( Figure 2B).

Quadratic Regression
We assume ω z between mhs k and mto k points to represent a quadratic curve ( Figure 2C): whereω z,t is the best-fitting result for angular velocity in the z-axis. A quadratic regression is then used to fit the raw data and to calculate the parameters a k , b k , c k : where M k is the number of points between mhs k and mto k .

Calculating the Split Point
Finally, the segmentation point sp k was defined as the maximum point of the quadratic fitting result ( Figure 2D): The k-th cycle was defined by the data between sp k and sp k+1 . In each cycle, we estimate the trajectory.

Estimating the Trajectory of a Step
The foot trajectory in each cycle can be calculated by integrating the acceleration between each segmentation point. Two coordinate systems are applied: a laboratory coordinate system (e) and a sensor coordinate system (s). Because the raw data from the sensor are represented by the time-variant sensor coordinates, we need to transpose them into timeinvariant laboratory coordinates. Acceleration in the laboratory coordinates can be converted from measured sensor acceleration using a rotational matrix: where θ , φ, and ψ are the Euler angles around the x-, y-, and z-axes. We divide the estimation process of the foot trajectory into five steps: (1) calculate the initial Euler angles; (2) calculate the time derivative of the Euler angles; (3) calculate the Euler angles using the integral; (4) transpose the accelerations using a rotational matrix; and (5) calculate the trajectory using the double integral.

Calculating the Initial Euler Angles
At the beginning of each cycle, we can assume that the foot is in full contact with the floor and is momentarily stationary. The accelerometer is assumed to detect only gravitational acceleration g at the outset: where a s,0 is the initial acceleration vector in the sensor coordinate, and φ 0 and ψ 0 are the initial Euler angles around the y-and z-axes. Therefore, the initial Euler angles vector θ 0 can be calculated as: Calculating the Time Derivative of the Euler Angles The relation between the ith angular velocity ω i from the gyroscope and the time derivation of Euler anglesθ i can be calculated as: (6) where φ i and ψ i are the ith Euler angles around the y-and z-axes.

Calculating the Euler Angles by Integral
Euler angles are derived by the integration of Euler angle derivation: where t is the sampling rate.

Transposing Accelerations Using a Rotational Matrix
The rotational matrix calculated by Equation (3) is used to estimate the acceleration in laboratory coordinates R s→e a s , which contains the gravitational acceleration. Linear acceleration a e in laboratory coordinates can then be calculated by simply subtracting gravity: Calculating Trajectory Using the Double Integral The i-th velocity v e,i in the laboratory coordinates is estimated by integration of linear acceleration a e,i : and the i-th foot trajectory p e,i is estimated by integration of v e,i :

Reducing of Brownian Noise
Integration may drift because of IMU sensor error, and the calculations of velocity and trajectory are corrected by the constraint condition. In each cycle, both the initial value and the end value of the velocity in three FIGURE 2 | Dissociation of the continuous gait signal (see the section Proposed method). The raw data contains much noise, and a median filter was first used to smooth the data. (A) We found local maxima that are larger than a threshold and (B) then the heel-strike and toe-off points. (C) We assume a quadratic curve between the heel-strike and toe-off points. (D) Finally, the split point was defined as the maximum point of the quadratic fitting result. The gait cycle was defined by the data from the split point to the next split point.
directions and the trajectory in the vertical direction can be assumed as 0. The algorithm below shows how to estimate velocity.
First, the forward integral and back integral are calculated separately from linear acceleration a e,i : where the superscripts f and b mean forward and backward, respectively. The correction result can then be derived by the weighted average of the forward and backward integral: where w i is the weight and w i ∈ [0, 1].
Because v f e,i is more accurate near the starting point and v b e,i is more accurate near the end point, the function for calculating w i should increase monotonically. Here, we choose the sigmoid function to calculate w i : where N is the number of points in the current cycle and m is a hyperparameter calculated from experiment, which we choose as m = 0.1. The trajectory in the vertical direction can also be calculated using the algorithm above.

Estimating of Spatial and Temporal Parameters for Clinical Gait Analysis
The gait events included the HS and TO points, which were extracted first. The HS and TO events detection algorithm was based on the peak detection of the raw angular velocity in sagittal plane ω z . At the end of the swing period, several of negative peaks can be observed in ω z and the first one is associated with the HS instant (Salarian et al., 2004). Before the swing period, a negative peak is associated with the TO instant (Salarian et al., 2004). For the definition of each gait event search region, the pattern of shank tilt angle, which was inspired from an instep-based previous study (Tunca et al., 2017) was introduced. At the end of the swing period, the shank will rotate around the knee axis caused by knee extension and reach the maximum forward. Then, the rotation of the shank around the ankle axis will start and the foot contact with the ground will produce the HS instant. In this process, θ z will appear to increase first and then decrease. Thus a positive peak will appear in θ z before the HS instant where θ z was computed via the integration of ω z with the sampling interval t. For the convenience of the description, we refer to this instant where peak occurrence as shank-max-forward (SMF). Similarly, after the TO instant, the ankle will lift slightly and rotate until reaching a certain height. Then it will start to rotate and present a negative peak in θ z . We refer to this instant as the shank-maxbackward (SMB) for convenience. As a result, SMF and SMB of θ z can be used to define a proper search interval of the HS and TO. We found the SMF and SMB via a peak search algorithm signal.find_peaks in SciPy (version 1.2.0) which can find proper peaks via the prominence (define intrinsic height of a peak) and the distance (define the distance between peaks) properties. Then, SMFs are the positive peaks and SMBs are the negative peaks whose prominence is larger than 0.2 [rad] and the distance is at least 0.4 [s]. Finally, HS is defined as the first peak appears after each SMF instant, and the TO is the minimum of angular velocity in the interval (SMB − 0.3 [s], SMB).
Estimated stride length SL is calculated by the trajectory in the y and z direction: where p e,i = (p x,i , p y,i , p z,i ) T . Estimated stride duration is calculated as the time from one HS to the next ipsilateral HS. Estimated gait speed is calculated as the value obtained by dividing stride length by stride duration. Estimated stance duration is calculated as the time from HS to ipsilateral TO. Estimated swing duration time is calculated as the time from TO to the next ipsilateral HS.

Overview of the Evaluation
We conducted an experimental evaluation of our proposed method to validate the accuracy of the trajectory estimation of the shanks (just above the ankles) to verify whether it allows the analysis of gait for clinical purposes. Twenty healthy people participated in the experiment, and we used a motion capture system as the criterion standard for the evaluating gait in a clinical setting. We evaluated the accuracy of our proposed method for calculating the estimated trajectory and clinical gait parameters. We used an IMU attached on the shanks for the experimental evaluation. An optical motion capture system (Nobby Tech. Ltd., Tokyo, Japan) was used as the reference system. We used 12 cameras, and the motion capture volume was about 2 m × 7 m × 1 m (width, length, and height; Figure 3A). The position error of the markers of the motion capture system during the calibration was <1 mm. The dimensions of the floor of the room were 18 m × 7 m (length and width; Figure 3A). Three optical markers were attached on each foot as shown in Figure 3B. Two of the three markers were attached on the heel and the toe (metatarsal head II) to assess gait parameters. The third marker was attached on the IMU to evaluate the trajectory estimation. To synchronize the IMU and the motion capture system, the participants hit their heel to the floor before the gait measurement. The peaks of the spike waveforms, that were caused by the heel hits of the participants, in both the IMU signals and the motion capture system data were used to define time 0.

Participants
We recruited 20 able-bodied volunteer participants, with no history of gait abnormalities, from the Tokyo Institute of Technology for the experimental evaluation. The Ethics Committee of Tokyo Institute of Technology approved the protocols for the evaluation, and all participants provided their written informed consent. Because of technical problems with the motion capture system, the data for one of the 20 volunteer participants were excluded from the analyses. Ultimately, we used data from 19 participants (mean age 23.9 years, 9 men and 10 women, mean height 1.66 ± 0.07 m and mean body mass index 20.2 ± 2.7).
As an example of the application of the proposed method to PD patients, we used the method to analyze gait in one healthy elderly participant and four patients with PD. The healthy elderly participant was recruited from a public interest incorporated association in Machida City, Tokyo that provides human resource services for elderly people. Patients with PD were recruited from Kanto Central Hospital, Tokyo. PD had been diagnosed by a physician. The exclusion criteria for this study were past history of other neurological or orthopedic disorders that can affect gait or posture (excluding PD). The healthy elderly participant and the four PD participants provided written informed consent in accordance with requirements of the Ethics Committee of Tokyo Institute of Technology. The Kanto Central Hospital Ethics Committee and the Ethics Committee of Tokyo Institute of Technology approved the protocol for the application of the proposed method.

Experimental Task
To construct the spike waveforms for synchronization between the IMU and the motion capture system (for the synchronization method, see "Overview of the evaluation"), the participants hit their heel to the floor before the gait measurement. In two trials, the participants walked on a flat floor at their own self-selected natural pace and a slow pace. In each trial, the participants walk straight in the motion capture volume and turned outside of the motion capture volume. Thus, each trial comprised four straight walks (two round trips; Figure 3A, left). We used the gait data obtained as the participants walked in the motion capture volume and removed the gait data as the participants turned outside of the motion capture volume.

Validation of the Location of IMUs
To consider the validity of the estimation of foot trajectory from IMUs attached on the shanks, we calculated correlations between stride length estimated with the proposed method, measured with a motion capture marker attached on the IMU, FIGURE 4 | Errors between the estimated foot trajectory identified with our proposed method and reference data from the motion capture system projected in the sagittal plane.

RESULTS
To evaluate the proposed method, the shank (just above the ankle) trajectory and clinical gait parameters calculated by our proposed method were compared with those collected by the motion capture system. The comparisons of trajectory information were conducted in the sagittal plane. Five clinical gait parameters (stride length, gait speed, stride duration, stance duration, and swing duration) were compared.

Comparison of the Trajectory Estimated With the Proposed Method and the Motion Capture System
The trajectories of our proposed method and the reference data are shown in Figure 4. The R value between displacement in the direction of forward movement calculated with the proposed method and measured with a marker attached on the IMU was 0.978 ( Figure 5A, left). The R value between the maximum vertical displacement calculated with the proposed method and measured with a marker attached on the IMU was 0.925 ( Figure 5A, right).

Validation of the Location of IMUs
The R value between displacement in the direction of forward movement calculated with the marker attached on the IMU and that measured with the marker attached on the heel was 0.994 ( Figure 5B).

Estimation of Clinical Gait Parameters
The means and standard deviations (SD) of the gait parameters compared between the proposed method and the motion capture system are summarized in Table 1. Figure 6 shows the agreement between the proposed method and the motion capture system in Bland-Altman plots. The mean ±1 SD accuracy of stride length was 0.054 ± 0.031 m ( Figure 6A). The R value between displacement in the direction of forward movement calculated with the proposed method and measured with a marker attached on the heel was 0.978 ( Figure 6A). The mean ±1 SD accuracy were as follows: 0.034 ± 0.039 m/s for gait speed ( Figure 6B); 0.002 ± 0.020 s for stride duration ( Figure 6C); 0.000 ± 0.017 s for stance duration ( Figure 6D); and 0.002 ± 0.024 s for swing duration (Figure 6E).

Application of the Proposed Method to Patients With a Gait Disorder
The shank trajectory over 15 steps for each participant is shown in Figure 7. The mean clinical gait parameters of the PD patients are summarized in Table 2.

DISCUSSION
We have proposed a new method for gait analysis that uses IMUs attached on the shanks to estimate foot trajectory and then to The proposed gait analysis method comprises two IMUs with a triaxial accelerometer, triaxial gyroscope, and tablet computer. This method can be applied in a variety of locations outside of the gait laboratory and is less expensive than conventional gait analysis methods such as motion capture systems. The clinical advantage is that the patient burden is low because of the light weight (about 24 g) and easy attachment of the IMUs. We therefore anticipate that the proposed method would be suitable for clinical gait analysis. As for the location of the IMUs, the R value between displacement in the direction of forward movement as measured with the marker attached on the IMU and as measured with the marker attached on the heel (0.994) indicates that the location of the IMUs is valid at least for estimating the stride length. The R value between displacement in the direction of forward movement as estimated by the proposed method and as measured with the marker of the motion capture system attached on the IMU indicates that displacement in the direction of forward movement estimated by the proposed method explained 96% of the variation in displacement in the direction of forward movement as measured with the motion capture system.
The mean error of stride length estimated with the proposed method was 0.054 ± 0.031 m ( Table 1). This result suggests that the proposed method can estimate clinical gait parameters such as stride length. A previous method in which the location of a IMU is on the dorsum of a foot found that the accuracy mean accuracy ± precision was 0.015 ± 0.068 m (Mariani et al., 2010). The IMU location on the shank may cause bias of this order of accuracy. We expected that further development of the method will overcome this limitation of performance. Several studies (Stolze et al., 2001;Curtze et al., 2015) have found that stride length is shorter in patients with PD than it is in healthy controls as observed in the example of the application of the proposed method ( Table 2). For example, Morris et al. (2005) reported that stride length in PD patients in the off state was   0.96 ± 0.19 m, which was shorter than the stride length of 1.46 ± 0.08 m as measured in healthy age-matched controls.
The R value between displacement in the direction of forward movement estimated by the proposed method and measured with the marker of the motion capture system attached on the heel indicates that stride length estimated by the proposed method explained 96% of the variation measured with the motion capture system. This result suggests that IMUs are potentially useful in clinical gait analysis. We expect that further development of the proposed method to evaluate the gait in people with PD. We expect that further development of this method or other methods will enable us to evaluate quantitatively the effects of drugs and interventions such as rehabilitation in patients with gait disorders. In the future, we plan to assess patients with gait abnormalities, such as those caused by PD. We will validate the proposed method to determine whether it can identify abnormal gait patterns, including shuffle, short-steppage, and hemiplegic gaits. The sampling frequency that we used in the present study was 100 Hz. We plan to investigate the effect of the sampling frequency on the estimation of the gait parameter in the next study.

CONCLUSION
Our results suggest that the proposed method is suitable for gait analysis whereas there is a room for improvement of its accuracy. Unlike methods that use motion capture systems, this method can be used in a variety of locations, such as in the corridor of a medical center. Further development of our proposed method is expected to enable clinicians to share objective information about gait features with health-care providers and patients.

DATA AVAILABILITY STATEMENT
The datasets for this article are not made publicly available because there is an agreement for data exchange between Tokyo Institute of Technology and Kanto Central Hospital.

ETHICS STATEMENT
The Ethics Committee of Tokyo Institute of Technology approved the protocols for this study, and all participants provided written informed consent.

AUTHOR CONTRIBUTIONS
KH, YO, YH, YMa, HO, SO, and YMi designed the research. YO, KH, YH, HS, and AI performed the experiments. YMa, YO, HO, KH, and YMi analyzed the data. HO, YO, and YMi wrote the paper.

FUNDING
This study was partly supported by the Core Research for Evolutional Science and Technology (CREST) Program Nano inertia detection device and system of the Japan Science and Technology Agency (JST), a Ministry of Education, Culture, Sports, Science and Technology/Japan Society for the Promotion of Science (MEXT/JSPS) grant (16K12950), and a COI-JST grant to the Research Center for the Earth Inclusive Sensing Empathizing with Silent Voices.