Estimation of Spatial-Temporal Gait Parameters Using a Low-Cost Ultrasonic Motion Analysis System

In this paper, a low-cost motion analysis system using a wireless ultrasonic sensor network is proposed and investigated. A methodology has been developed to extract spatial-temporal gait parameters including stride length, stride duration, stride velocity, stride cadence, and stride symmetry from 3D foot displacements estimated by the combination of spherical positioning technique and unscented Kalman filter. The performance of this system is validated against a camera-based system in the laboratory with 10 healthy volunteers. Numerical results show the feasibility of the proposed system with average error of 2.7% for all the estimated gait parameters. The influence of walking speed on the measurement accuracy of proposed system is also evaluated. Statistical analysis demonstrates its capability of being used as a gait assessment tool for some medical applications.


Introduction
The significance of spatial-temporal gait parameters measurement has been addressed in many research papers [1][2][3]. The quantitative analysis of such gait parameters can be helpful to diagnose impairments in balance control [4], monitor the progress in rehabilitation [5], and predict the risk of falling [6,7]. Such parameters include stride length, walking velocity, stride cadence, stride duration and asymmetry of stride. In particular, stride asymmetry has been shown to be more indicative the ultrasonic sensor node placed on human body is light and small. Additionally, the proposed motion analysis system is low-cost compared with the camera-based motion analysis system.

Related Work
The tracking techniques for locating a mobile device's position are studied using many approaches [28][29][30].
There are two major localization and tracking techniques, Receiver-Synchronization relative measurement (RS) and Global-Synchronization absolute range measurement (GS) [30]. RS range measurement only requires anchors to be synchronized and Time-Difference-of-Arrival (TDOA) technique is used for tracking and localization. In GS range measurement, both the mobile and the anchors are synchronized and the absolute distance can be estimated using Time-of-Arrival (TOA) technique. In our system, we prefer higher tracking and localization accuracy to accurately measure spatial-temporal gait parameters. Thus, we used the TOA-based tracking technique because the TDOA-based tracking technique has worse performance. RF signals are used in our system for synchronization between the anchors and the mobile. RF signal travels at the speed of light and the time it takes to reach mobile target is almost instantaneous and can be considered zero since the speed of ultrasound in air is much lower [31].
Under ideal range measurement case, an analytical localization method called trilateration, which uses only distance measurements, can be applied to identify the position of the mobile. For TOA-based localization technique, the target can be located at the intersecting point of several cycles that are formed by these anchors with known positions and distances to the mobile [31]. However, for a mobile target, it is not easy to track or localize because the range measurements are noisy and fluctuate, since the mobile can be located at anywhere in overlapped regions of such circles rather than being located at a single point at the intersection of the circles.
It is therefore desired to have accurate tracking and localization methods capable of filtering out the range measurement noises. One of the representative nonlinear state estimators is the least square (LS) method, which first transforms the nonlinear equations into linear ones and then solves the linear equations by LS-based estimator. Although the computation of this method is efficient, the tracking accuracy may not be sufficient [32]. Another typical method is proposed in [33], which begins with an initial guess and then applies least sum squared error to solve the nonlinear equation recursively. Although it provides better tracking performance, the initial guess should be carefully selected to guarantee the convergence of the iteration [34]. Therefore, many researchers proposed other methods to enhance the positioning performance. One representative implementation of indoor sensor network used to track a mobile is the Cricket of MIT [35,36]. It employs a hybrid approach involving the use of an Extended Kalman Filter (EKF) and Least Square Minimization to enhance the tracking and localization performance. EKF is the most commonly used nonlinear state estimator using the first or second order terms of the Taylor series expansion, which is most appropriate when the noise statistics is Gaussian distribution, to linearize the state and observation models [37]. Therefore, for some highly nonlinear dynamics, the linearization of EKF insufficiently characterizes the relationship. Therefore, we use Unscented Kalman Filter (UKF) to overcome such limitations of EKF, i.e., the requirement for the noises to be Gaussian and the poor linearization of first or second order approximation. We will explain the tracking algorithm in more detail in the following section.

Ultrasonic Sensor System
The acquisition system that we developed for wearable gait analysis uses the wireless sensor network concept, with all mobile nodes communicating wirelessly with the coordinator to enable patients to be monitored in unconstrained environment, as shown in Figure 1. The proposed measurement system consists of one ultrasonic transmitter (referred to as the mobile with form factor: 4 cm × 3 cm × 1.6 cm) and four ultrasonic receivers (referred to as the anchors with the same form factor) made by Embedream studio, China [38]. The foot displacements measured using the TOA-based tracking technique were expressed in a global coordinate system that described foot position relative to the ground, as shown in Figure 1a. The X-axis was defined as the direction of progression, i.e., anterior-posterior direction, and the Y-axis was defined as the vertical direction. The third axis of the coordinate system, i.e., the Z-axis, was determined in such a way to form a right-handed coordinate system. However, for healthy subjects, the 2-dimensional model is sufficient to obtain spatial-temporal gait parameters, because the sagittal plane is the plane where the majority of movements take place.   Figure 1b shows the configuration of the ultrasonic measurement system. A battery-powered ultrasonic transmitter node is attached to the heel of the subject's foot. The mobile sensor node establishes communication with the coordinator node through a low power 430 MHz RF transceiver RFM12B. The coordinator node is also wirelessly communicating with the computer through a wireless data transmission module. The wireless data transmission module forwards all collected information to a personal computer through RS232 cable for postprocessing.
In the system, ultrasonic range measurements are initiated by a periodic trigger input with a pulse duration of 10 µs. Then, the ultrasound transmitter is triggered to produce an ultrasonic burst consisting of 8 pulses with a frequency of 40 kHz. Meanwhile, the RF module on the mobile node is triggered synchronously, thus sending out a data package with a timer starter command (TSC) using broadcast address to notify the anchors that an ultrasound signal has been transmitted. Once the anchor receives TSC, it will start its 16-bit counter to record the propagation delay from the mobile to the anchor. The transmission time of the RF signal from the mobile is negligible, since the speed of light is much faster than the speed of ultrasound. The 16-bit counter will stop counting immediately after each of the transmitted burst is detected by the anchor. Then the counted steps will be converted to propagation delay by multiplying the time resolution (instruction cycle) of the microcontroller. From this delay, the distance between the mobile and the anchor can be calculated by: where d is the distance in meters, t is the propagation delay in seconds and v s is the speed of ultrasound in air. The ultrasound velocity can be approximated to [26]: where T c is the air temperature in degree Celsius. Together with the known positions of these anchors, the position of the mobile is located using the TOA-based tracking technique, which finds the intersection area of circles centered at each anchor with radius equal to the measured distances. The tracking algorithm is discussed in the following section.

Tracking Algorithm
In this section, we first explain how to establish a state space of nonlinear system to estimate the state of the moving target. Next we will apply UKF to enhance the performance of the tracking technique.

Motion and Measurement Model
The mobile target in 3-dimensional field is represented by its position and velocity in X-Y-Z plane: where P k = [x k y k z k ] T are the position coordinates along X-, Y-and Z-axes at time step k, anḋ P k = [ẋ kẏkżk ] T are the moving velocities with respect to these three axes at time step k. To formulate the dynamic transition process, the following state space equations are given where where ∆t k = t k+1 − t k is the sampling interval. w k = [w x w y w z ] T is the white Gaussian noise with zero mean and covariance matrix W = diag(σ 2 wx , σ 2 wy , σ 2 wz ). V = diag(e 2 1 , e 2 2 , e 2 3 , e 2 4 ) denotes the covariance matrix of the measurement noise v k .
Let d ik denote the measured distance at the ith anchor using the equation: where [x i y i z i ] is the known position of anchor i, e ik is the distance measurement error at anchor i,

Unscented Kalman Filtering
The aforementioned state space model is a nonlinear dynamical system to the measurement distances and the state of foot motion. The approximation of UKF is to find a transformation that captures the mean and covariance of state random variable of length n through a nonlinear function [39]. We summarize the algorithm as follows.
For each time step k, start from X k/k and P k/k ,

Generate sigma points
2. Compute the a-priori statistics 4. Compute the Kalman gain 5. Compute the a-posteriori statistics W is the associated weight matrix: Parameter λ is the scaling factor, which is defined as: where α and κ control the spread of the sigma points around the mean of the state (α is usually set to a small positive value, e.g., 10 −3 ) and κ is set to 0), β is related to the distribution of state variable (for Gaussian distribution, β = 2 is optimal).

Autocorrelation Procedure
The idea of analyzing gait data by autocorrelation procedure is first proposed by Barrey et al. [40] and Auvinet et al. [41]. Then, the difference between biased and unbiased autocorrelation procedure for gait data analysis has been discussed by Moe et al. [9]. Here we summarize the autocorrelation procedure as follows.
The autocorrelation coefficient shows the degree of similarity between the given observations a i (i = 1, 2, ..., N ) as a function of the time lag over successive time intervals, as given by: where m is the phase shift in the number of observations. The autocorrelation coefficients of a periodical signal will produce peak values for lag time equivalent to the cycle of the signal, which is the stride duration. Therefore, visual assessment of autocorrelation from the time series plot can be used to inspect the structure of a cyclic component. As discussed in [9], either biased or unbiased autocorrelation coefficient can be computed for gait data analysis, but biased autocorrelation is not suitable for comparing autocorrelation coefficient over different time lags. The biased autocorrelation is the result of the raw autocorrelation coefficient A divided by the number of the observations in Equation (14): In Equation (15), the denominator N is the number of samples in observation a i , which is independent of the time lag m. It means that the number of samples in the numerator will decrease as the time lag m increases, and then the autocorrelation coefficient will attenuate. However, this is not the case in unbiased autocorrelation estimator, expressed as: Since the number of terms in the numerator N − m is always equal to the value of the denominator, there is no noticeable attenuation in the unbiased estimator. Figure 2 shows the two different estimators for horizontal displacement during treadmill walking. The biased estimator shows clear periodicity but with attenuated amplitudes, while the unbiased estimator introduces no obvious attenuation except a deteriorated curve at the tails.   Figure 3 shows the normalized unbiased autocorrelation of horizontal and vertical foot displacement during treadmill walking. Since the first peak from the zero phase represents a phase shift of one stride duration, the autocorrelation coefficient at the periodic phase shift is defined as the regularity of the stride between neighboring strides, referred to as hR i for horizontal displacement and vR i for vertical displacement. Therefore, either for horizontal or vertical displacement, the closeness of hR i+1 /hR i or vR i+1 /vR i reflects the stride symmetry. Figure 4 demonstrates an example of asymmetric gait showing the unbiased autocorrelation sequence of the horizontal and vertical displacements.

Estimation of Gait Parameters
From the estimated foot displacements by the proposed algorithm, the following spatial-temporal gait parameters can be obtained. With respect to the jth gait cycle, the estimators of the spatial-temporal gait parameters are as follows: • Stride Length, S: where the functions M ax(x) and M in(x) return the maximum and minimum value of the variable x, and x j is the horizontal displacement in the jth gait cycle; where N S is defined as the stride length normalized by the number of strides n; • Stride Duration, T : where the function Index(max(x j )) returns the location of the maximum value in x j ; • Normalized Velocity, N V : where the normalized speed is the speed as percentage of the number of strides n; • Cadence, C: where the cadence is the number of strides in a second.

Experiment Setup
The proposed method was tested on 10 healthy subjects (age 25.7 ± 1.4 years, height 171.4 ± 6.5 cm, and weight 62.8 ± 5.6 kg) walking 5 min on a treadmill at slow, normal, and fast walking speeds, the results of which are presented in this paper. The subjects were recruited from students of Nanyang Technological University and none of them had a history of pathological gait disorders. To provide a more systematic validation, we conducted the experiments in a motion analysis lab with eight high speed cameras (Motion Analysis Eagle System, Santa Rosa, CA, USA) in the School of Mechanical and Aerospace Engineering at Nanyang Technological University. The Motion Analysis Eagle System consists of Eagle Digital Cameras and Cortex software, which captures complex 3D motion with extreme accuracy. System calibrations of the reference system should be done at both static (with 4-point calibration L-frame) and dynamic process (with 3-point calibration wand) to ensure an acceptable accuracy of the reference system. In our experiments, the accuracy of the reference system is 0.43 ± 0.18 mm (Average ± Standard deviation). Figure 1a shows the placement of ultrasonic sensor and reflective markers on the test subject's foot. There were four anchors used in our experiment with positions p 1 = [0 0 0] T , p 2 = [0.324m 0 0] T , p 3 = [0.324m 0.230m 0] T , p 4 = [0 0.230m 0] T . The ultrasonic transmitter was attached to the heel of the foot pointing towards the four anchors, using elastic straps. In our method, only one ultrasonic sensor (transmitter) is needed to attach to the foot, which minimizes user discomfort and avoids complex calibration procedures and synchronization issues. All data transmission between anchors, coordinator and transmitter are done wirelessly through the RF module. Therefore, it does not restrict the movement of subjects. The ultrasonic sensor data were acquired at 50 Hz. Data from the reference system were captured at 200 Hz. The difference between the sampling rate of these two systems was compensated by linear interpolation. All data were low-pass filtered by second order low-pass Butterworth filter at 10 Hz.

Processing of Measured Data
In order to compare the estimated spatial-temporal gait parameters at each recorded gait cycle, the foot trajectory estimate with proposed ultrasonic sensors was temporally delayed to match the trajectory estimated by the camera reference system, by finding the maximum values of cross-correlation between these two trajectories. To quantify the performance of the proposed system against the camera reference system, the mean and standard deviation (std) were calculated on the datasets of difference, as well as the Root Mean Square Error (RMSE). This is followed by using the analysis of variance (ANOVA) to test differences in the means of the ten subjects for statistical significance. Finally, walking speed was estimated using the proposed ultrasonic sensor configuration to check significant changes over different speeds. Two-sample t-tests were performed on the walking velocity and the extracted gait parameters to assess the significance of change in these gait parameters with speed, and thus investigate the effect of walking velocity on the difference between the proposed system and the reference system in gait parameters estimation.

Parameters Identification
As the system modelling we have adopted in Section 3.2.1., the process and measurement noise statistics should be estimated. A wooden pendulum was constructed using a uniaxial pivot so that it swung through an arc [42]. The ultrasonic transmitter was placed at the end of this pendulum, and a reflective marker was also located approximately in alignment with the ultrasonic transmitter head. The pendulum was raised up at an angle and allowed to drop freely until it came to a stable position. This action was repeated M times. The experiment helps to find suitable values of process noise W and measurement noise V . The measurements from camera system, r i , are referred to as the actual distance for test i, and there are N measurement samples m j i collected for each test, where j = 1, · · · , N .

Process Noise Statistics in Kalman Filter
As the process noise in UKF is an independent variable, it is difficult to get an exact value [31]. Here, we consider it as a velocity noise in X, Y and Z directions in mm/s. The process noise W was estimated using numerical methods. By varying the values of σ wx , σ wy and σ wz , we will get the corresponding trajectory of the mobile to compute the RMSE value. Typical values of σ wx , σ wy and σ wz will be selected when their corresponding RMSE is minimal. The typical values of W used in our experiments are σ wx = 30, σ wy = 25, σ wz = 10.

Measurement Noise Statistics in Kalman Filter
It is reasonable to assume that all anchors have independent distributed noise. Then, the mean and covariance of the measurement noise can be evaluated by the pendulum experiments. Using the data obtained from the specific experiments, straightforward calculations lead to the estimation of mean and variance of the measurement errors Typical value of V used in our experiments is V = diag (11, 9.3, 9, 9) with the units as mm 2 . In other words, the accuracy of distance measurement by each ultrasonic sensors is around 3 mm.
The results of pendulum experiment have been shown in Table 1. The Net RMSE is defined as The difference between the two systems was obtained with an RMSE value of 4.08 mm in horizontal direction (X), 0.72 mm in vertical direction (Y) and 1.08 mm in lateral direction. The Net RMSE value of 4.28 mm in 3D space of UKF estimator is achievable in the pendulum model.

Performance Comparison
The mean and standard deviation in stride length, stride duration, and stride velocity estimation between the proposed system and the reference system together with RMSE value are reported in Tables 2-4 for all subjects walking at normal speed. On average, across all subjects, the estimates of stride length from the proposed method were 0.001 m less than the reference measurements. The overall RMSE value is about 0.027 m, which is 2.3% of the mean estimated stride length of the reference system. The mean and standard deviation of stride duration at normal walking speed is reported as 1.18 ± 0.02 s by the reference system and 1.18 ± 0.04 s by the proposed system, which shows no mean difference between the two systems. The average error across all subjects of RMSE of the estimated stride duration is 0.035 s with 3% percent error. The mean and standard deviation in the estimation of the stride velocity is reported in Table 4, which shows that the proposed method slightly overestimated the stride velocity by 0.001 m/s with an RMSE value of 0.036 m/s, occupying 3.6% of the proposed estimates of stride velocity.  Table 4. Mean and standard deviation (in meters per second) of the reference (Ref) and proposed (Pro) systems and RMSE in detecting stride velocity for each subject. Averaged values across the ten subjects are also reported. We have elaborated how gait cycle periodicity of foot displacement data can be used to extract stride regularity and symmetry by unbiased autocorrelation procedure in Section 3.3.2. . Tables 5 and 6 show the mean and standard deviation of the reference system and the proposed system together with RMSE values in detecting horizontal and vertical stride symmetry respectively for each subject. The mean and standard deviation data of horizontal stride symmetry are 1.001 ± 0.021 by the reference system and 0.999 ± 0.027 by the proposed system, which shows that the ultrasonic-based horizontal stride symmetry was underestimated by a negligible error of 0.002. An RMSE of 0.013 with 1.3% percent error is also reported for the estimates of horizontal stride symmetry across all subjects. In the contrary, the ultrasonic-based vertical stride symmetry was overestimated by 0.007, where the RMSE value is 0.034 with a percent error of 3.5%. In summary, all the numerical results show a clinically acceptable accuracy of the proposed system with an average percent error of 2.7% for all the estimated gait parameters. Table 5. Mean and standard deviation of the reference (Ref) and proposed (Pro) systems and RMSE in detecting horizontal stride symmetry (hS) for each subject. Averaged values across the ten subjects are also reported.  Table 6. Mean and standard deviation of the reference (Ref) and proposed (Pro) systems and RMSE in detecting vertical stride symmetry (vS) for each subject. Averaged values across the ten subjects are also reported.

Statistical Analysis
In this part, ANOVA has been performed to test differences in the means (for ten subjects) for statistical significance. We base this test on a comparison of the variance due to the between-groups variability (called Mean Square Effect, or M S ef f ect ) with the variance due to the within-group variability (called Mean Square Error, or M S error ). Before applying ANOVA, whether the distribution of the data is normal or not should be checked. Results are reported in Table 7 and Figure A1. In Table 7, H = 0 indicates that the null hypothesis ("mean is zero") cannot be rejected at the 5% significance level. The p-value is the probability of observing the given result by chance if the null hypothesis is true. Large value of p shows the validity of the null hypothesis. As in Table 7, not only all values of H are equal to zero and the values of p are equal to one, but also the means of estimates are located in the 95% confidence interval. Therefore, the estimated parameters are normally distributed. Under the null hypothesis (that there are no mean differences among subjects), we compare the M S ef f ect and M S error via the F-test, which tests whether the ratio of the two variance estimates is significantly greater than 1. Otherwise, we will accept the null hypothesis of no differences between the means, i.e., the means (in the population) are not statistically different from each other. Figure B1 shows the boxplots of stride length, stride duration, stride velocity, horizontal stride symmetry and vertical stride symmetry for each subject. The analysis of variance is summarized in Table C1. As shown in Table C1, for all estimated gait parameters, the small value of between-groups sum of squares likely indicates no differences among the subjects. Additionally, the values of F are less than 1, which indicates that the means of all gait parameters are not statistically different. Table 8 provides the numerical results of estimated gait parameters by the proposed system compared with those obtained from the reference system using the pair t-test. Significant difference was assumed when the null hypothesis can be rejected at p-value smaller than 0.05. The walking speed, on average, across all subjects was significantly different (p < 0.001 for the two measurement systems) among slow (0.54 ± 0.02 m/s), normal (0.99 ± 0.04 m/s)and fast (1.40 ± 0.04 m/s) speed. The influence of walking speed on all spatial-temporal gait parameters was tested by the mean and standard deviation values for the proposed and reference systems.

The Effect of Walking Speed on the Measurement of Gait parameters
The measurement errors of estimated S, NS, NV, C, and vS were not affected significantly by the changes in walking velocity (p > 0.05). Particularly, there is no difference in cadence estimation between the proposed and reference systems. The influence of speed on the measurement errors of stride duration T was found to be significantly higher (p < 0.05) at fast speed, but it was not significant for V and hS. This can be interpreted as the lower temporal resolution at higher walking speed. Figure 5 shows significant changes in T and C, but there is no significant change in other parameters. Although the means of both horizontal stride symmetry and vertical stride symmetry are not statistically significant, the largest variations at slow speeds were observed. Therefore, the stride symmetry can be used as warning sign of walking disorders.

Discussion and Conclusions
In this paper, a low-cost ultrasonic motion analysis system using an ultrasonic transmitter and four receivers to track the foot displacement in 3D space is developed. The proposed motion analysis system has been validated against camera-based system with 10 healthy subjects, and shown to produce accurate estimates of some spatial-temporal gait parameters including stride length with RMSE value of 0.027 m (2.3%), stride duration with RMSE value of 0.035 s (3%), stride velocity with RMSE value of 0.036 m/s (3.6%), horizontal stride symmetry with RMSE value of 0.013 (1.3%) and vertical stride symmetry with RMSE value of 0.034 m (3.5%). We have further evaluated the influence of walking speed on these gait parameters by paired t-test.
The proposed system includes some ultrasonic sensors and micro-controllers, estimated today at about a cost of $100, which is inexpensive compared with current commercial camera-based system. With the rapid development of technology, the performance of these sensors will continue to improve while becoming available at even lower price. Therefore, low-cost in-home monitoring system for clinical applications can be possible.
As the work stated here is a first step to evaluate the feasibility of the proposed ultrasonic system, only ten healthy subjects participated in the experiments and were instructed to walk 5 min on treadmill at different speeds. The walking experiments were chosen on treadmill due to the limited measurement volume of the reference camera-based system. In addition, we can get a cyclic signal on horizontal displacement to analyze the stride symmetry. Although the proposed ultrasonic motion analysis system also has such limitations, the maximum propagation distance of the ultrasonic signal used in our system is 20 m, which is large enough for indoor applications.
Although the positive results showed the feasibility of applying such a system for in-home monitoring, there is an issue to be addressed in further research, i.e., how to deal with the multipath propagation. All the experiments in this study are conducted under line-of-sight condition, where the ultrasonic transmitter faces all the receivers without any obstacles between them. Therefore, for 3D displacements, according to spherical positioning technique, a minimum of 4 anchors with known positions are required. The method used in our experiment to mitigate the multipath propagation is by setting an inhibit time, i.e., the ultrasound detector will be disabled within the inhibit time to detect an ultrasound signal, and will be enabled again after the inhibit time has passed. Another possible solution is that we can use more receivers, which can not only account for multipath propagation, but also increase the measurement volume and accuracy of the proposed system [34].
Long-term monitoring is expected to be more challenging as demonstrated in some studies [18,43]. In [18,43], foot clearance measurement using inertial sensors is proposed and investigated. The displacement estimation requires double integration of measured accelerations from inertial sensors, which involves error accumulation over long time monitoring. Even though the growth uncertainty that arises from the integration of acceleration error can be mitigated by periodic corrections like ZUPT, the prerequisite is that the initial and/or terminal contact should be detected correctly, but it may be difficult in some type of abnormal gait. Although not specifically studied under long term monitoring, the proposed system does not have significant error accumulation for a 5 min walk.
In summary, we used a low-cost ultrasonic motion analysis system to extract spatial-temporal gait parameters, and tested the feasibility of the system against a reference camera-based system. The positive results demonstrated a great potential in using this low-cost system for clinical applications such as rehabilitation, gait analysis, and sports. For further work, experiments conducted with patients in collaboration with a hospital are being planned using our system. Figure A1. Histograms of all estimates of stride length, stride duration, stride velocity, horizontal stride symmetry and vertical stride symmetry.   Table C1. Summary of the analysis of variance of stride length, stride duration, stride velocity, horizontal stride symmetry, and vertical stride symmetry. ANOVA