A Multi-Path Compensation Method for Ranging in Wearable Ultrasonic Sensor Networks for Human Gait Analysis

Gait analysis in unrestrained environments can be done with a single wearable ultrasonic sensor node on the lower limb and four fixed anchor nodes. The accuracy demanded by such systems is very high. Chirp signals can provide better ranging and localization performance in ultrasonic systems. However, we cannot neglect the multi-path effect in typical indoor environments for ultrasonic signals. The multi-path components closer to the line of sight component cannot be identified during correlation reception which leads to errors in the estimated range and which in turn affects the localization and tracking performance. We propose a novel method to reduce the multi-path effect in ultrasonic sensor networks in typical indoor environments. A gait analysis system with one mobile node attached to the lower limb was designed to test the performance of the proposed system during an indoor treadmill walking experiment. An optical motion capture system was used as a benchmark for the experiments. The proposed method gave better tracking accuracy compared to conventional coherent receivers. The static measurements gave 2.45 mm standard deviation compared to 10.45 mm using the classical approach. The RMSE between the ultrasonic gait analysis system and the reference system improved from 28.70 mm to 22.28 mm. The gait analysis system gave good performance for extraction of spatial and temporal parameters.


Introduction
Gait analysis is an important clinical tool to assess disorders due to neurovascular or musculoskeletal diseases and to study relations between gait and falls in the elderly population [1]. The gold standard system for human motion capture and gait analysis is an optical system with high-speed infrared cameras [2]. Optical motion tracking requires dedicated laboratories with controlled lighting conditions which makes this system unfit for unrestrained environments or a homely setup [3]. Non-traditional methods developed for human motion tracking include magnetic sensors [4], laser sensors [5,6], inertial measurement units [7,8], ultrawideband (UWB) ranging sensors [2,9], ultrasonic sensors [1,10] etc. Force sensing platforms or wearable force sensors can be used to the get ground reaction forces, the center of pressure and the joint moments [11]. Magnetic sensors are affected by the heterogeneity of the earth's magnetic field and other ferromagnetic materials [12]. Laser sensors take a long time to take one sample of the whole scenario due to their scanning pattern [13]. Inertial sensors are affected by the drift during integration and fluctuating offset values and require sensor fusion to obtain accurate results [14]. UWB sensors demand high accuracy clock synchronization and high-cost ADCs to provide accurate ranging measurements [1]. Ultrasonic sensors can provide high accuracy time of flight (ToF) measurements without using expensive hardware as the wave travels much slower The rest of the paper is organized as follows: Section 2 explains the proposed multi-path compensation method. Section 3 explains the simulations conducted and their results. Section 4 explains the setup used for the experiments with ultrasonic sensors and the gait analysis system. Results from the experiments, the performance of tracking and estimation of spatial and temporal parameters are mentioned in Section 5. Finally, conclusions from the work are explained in Section 6 with some suggestions for future work.

Proposed Method
In the proposed method, linear chirp signals are used for ranging. Chirp pulse compression is the process of transforming a long duration frequency modulated pulse into a narrow pulse with much higher amplitude. The pulse compression process can be considered as an application of matched filter system. Let g(t) be the impulse response of the transmitter coding filter. Then, the response of the matched filter at the receiver should be g * (−t), where g * () represents the complex conjugate of g() and the output of the matched filter receiver can be given as shown in Equation (1).
In chirp compression process, consider g(t) as a linear chirp signal multiplied by a rectangular window function given by From Equation (1), the output of chirp compression process, y(t) can be derived as where, T is the chirp duration, B is the chirp bandwidth and f 0 is the centre frequency. If we transmit a linear chirp signal with signal duration T seconds, the output, y(t) will be a waveform with sync characteristics with signal duration 2T seconds. In traditional ultrasonic ranging methods, the time at which the peak of y(t) happens is taken as an estimate of the ToF. The low cost, low power piezoelectric transducers comes with a narrow bandwidth around 40 kHz. We used MA40S4S/R transducers from MURATA Electronics which has a bandwidth of about 2 kHz. Hence, we selected the frequency range from 39 to 41 kHz. We selected linear chirp signals as they are robust in terms of Doppler sensitivity and frequency selective fading [32]. Non-linear chirp signals give better results when the time-bandwidth product is higher, which was 14 in our experiments. Also, non-linear chirps give higher side lobs in the presence of Doppler effect [33].
In the chirp compression process, the range resolution of an ultrasonic ranging sensor is the minimum separation of two targets ranging using the same signal that can be identified as separate targets from the received signals. The range resolution of the ultrasonic ranging in ideal case with a chirp signal multiplied by rectangular window function can be given by v/δ f , where v is the speed of sound and δ f is the bandwidth of chirp. In this study, we consider linear chirp signal from 39-41 kHz, thus δ f = 2 kHz and the range resolution was found to be 172 mm when the speed of sound was 345.3 m/s. If the extra path length of the multi-path component lies within 172 mm, which is quite likely to happen in a multi-path environment, there might not be two separate peaks and the resultant peak can be shifted in the time axis. The earliest peak detection algorithm explained in [24] fails here, leading to errors in the measured value. Let T p be the position of the peak value of correlation signal mentioned in Equation (1) and P 0 be the peak value. In a close multi-path environment, the position of the peak will be shifted to T p + X. Where X represents the error. Now, instead of taking P 0 to calculate the ToF, we propose to take P 0 m crossing point to the left of earliest peak instead of the actual peak of the absolute value of the correlation waveform to calculate the ToF. The earliest 1 m th peak will be less affected by the close multi-path components as the multi-path components travel a higher distance than the LOS path and hence appears time-shifted in the cross-correlation of the received signal with the originally transmitted signal. Since the correlation output at the receiver is sync signal with the chances of corruption by the multi-path components higher at the descending part of the main peak as shown in Figure 1, the chances of shift in time domain at the earliest 1 m peak point (ascending part of the main peak) will be less. In this work, the value of m is chosen out of 2, 3 and 4 by conducting simulations with 10,000 iterations with different multi-path components each time. However, P 0 m should be kept at a safe margin above the channel noise floor. We selected the value of m as 2, 3 and 4 for simulation experiments keeping m = 1 for conventional method and the best value was chosen for implementation in the proposed gait analysis system. The pictorial representation of the proposed method is shown in Figure 1. Since we are removing any constant offset which is present in the range measurement, we are concerned only about the relative distance between two range measurements of the same transmitter. Even though the existence of a close multi-path shifts the time at which the correlation peak happens, the earliest 1 m ×peak point will be less affected compared to the actual peak. The chirp signals can be multiplied with various kinds of window functions to increase the sensitivity and to reduce the amplitude of side lobes in the signal [34]. We test the performance of three different window functions (rectangular, hamming and hann) along with the proposed method. Let S(t) be the transmitted chirp signal multiplied by the window function, W(t) given by Equation (4).
where, f s , f e and φ represents the starting frequency, the ending frequency and the initial phase respectively. The received signal without any multi-path components can be represented aŝ where, A 0 , h(t) and * represents the amplitude after path loss, channel impulse response and the convolution operator respectively. Considering M multi-path components, the correlation of received signal can be given by Equation (6).
where, , n(t), A m and τ m represents the correlation operator, channel noise, the amplitude and the extra path travelled by the multi-path components respectively. In Equation (6), the actual ToA information from LOS path is embedded in the first term. The second term in Equation (6) is the multi-path noise. Our hypothesis is that the earliest point at which the value of C(t) crosses 1 m (where, m > 1) of the peak of C(t) will be less affected by the multi-path components (the second term in Equation (6)) in time axis compared to the peak of C(t). This is because the multi-path components appear at a later time compared to the LOS path.

Doppler Compensation in the Received Signal
Moving sensor nodes emitting ultrasonic waves create Doppler effect and this alters the actual range measurements. A Doppler compensation method with linear up-and down-chirp signals was explained in [23]. Even though linear chirp signals show higher Doppler-tolerance, range-Doppler coupling leads to an error which is proportional to the moving velocity. For down-chirp signals, the error due to Doppler effect will be the same in magnitude with different polarity to that of up-chirp signals. In our experiments, in each ranging cycle of 40 ms, ranging with up-chirp was conducted in the initial 20 ms and with down-chirp in the next 20 ms duration. The chirp duration was selected as 7 ms as we allocated 20 ms for up-chirp and 20 ms for down-chirp and out of 20 ms, 13 ms time was provided between two signals so that the echo will die out completely. Let R 1 , R 2 and R 0 be the range estimated from the up-chirp, down-chirp and the actual range respectively. Then, whereṘ represents the Doppler velocity of the moving target

Simulation Results
A customized ranging environment with one transmitter and one receiver was simulated in Matlab to test the performance of the multi-path compensation technique. The amplitude of signal after attenuation in air was modeled as A = A 0 e −γ×δd , where A 0 and δd are the originally transmitted amplitude and the distance traveled through the medium respectively. The value of attenuation constant, γ was selected as 0.17 Np/m. The ranging was done with a chirp signal with linearly increasing frequency from 39 to 41 kHz in 7 ms. In order to test the performance of other non-linear chirp signals, a simulation was conducted to find the performance of linear, logarithmic and quadratic chirp signals to estimate the Doppler-velocity using the method explained in Section 2.1 with 1000 iterations. The results show that the performance of the linear chirp was better than the logarithmic and quadratic chirp signals. A close multi-path component (NLOS1) with a reflection coefficient of 0.9 was provided at random positions so that the difference between the LOS and the NLOS was less than the range resolution which is 172 mm in this case. The relative distance between two LOS values which were set to 1000 mm and 2000 mm respectively was calculated to find the improvement in ranging accuracy using our proposed method and the classical peak detection method. A Monte Carlo simulation with 10,000 iterations was done to estimate the distance at each position with different multi-path reflections. We selected the time index at which the cross-correlation of the received signal with the original signal crosses 1 m of its earliest peak value for the first time as an estimate for the LOS path. The plots of the correlation waveform from the LOS, NLOS and the combined signal with m = 2, 3 and 4 points marked are shown in Figure 1. The results obtained when the transmitted signals are multiplied by three different kinds of window functions are as shown in Table 1, where, I av , I max , D av and D max represent the average improvement, maximum improvement, average degradation, and maximum degradation in absolute error compared to peak detection respectively. Out of 10,000 iterations, N(I) and N(D) are the number of cases with improvement and degradation respectively and the ratio of average improvement to average degradation was calculated as R = I av /D av . In order to find the changes in the ratio, R as the value of m changes, we repeated the same simulation experiment with 10,000 iterations and a rectangular window function with 10 dB channel SNR for each value of m increasing from 1 to 4 in steps of 0.1 and the plot of R against m was shown in Figure 2. It can be observed that R slightly decreases after m = 1.9 as the value of m increases. We can also observe a sudden increase in R as m changes from 1 to 1.1. We found that the "hann" and "hamming" windows showed better results in terms of the ratio, R for m = 3 and 4. This can be explained by lower side lobs obtained for the correlation output y(t). With the attenuated side lobs, the time shift by the multi-path component in the 1 m peak also decreases thus resulting in better R values. However, it should be noted that as the value of m increases, the estimate is getting closer to the noise floor and chances of errors and D max also increase. Hence, optimum value of m depends on the peak value of cross correlation and noise power. It can also be observed that the maximum degradation of the accuracy is higher for "hann" and "hamming" windows. The plot of average improvement and average degradation in the absolute error while ranging with linear chirp signals multiplied with rectangular window and a channel with 10 dB SNR is shown in Figure 3.
The Empirical Cumulative Distribution Function (ECDF) for the errors in range estimation from the proposed method and the conventional peak detection method was calculated for rectangular window function and m = 2 and plotted in Figure 4. The plot shows that the proposed method outperforms the conventional peak detection method in a high multi-path environment. Thus, we selected m = 2 and a rectangular window function at the transmitter side for our tests using ultrasonic sensors for further analysis.

Materials and Methods
Initially, we conduct ranging experiments with one static transmitter and receiver in an indoor environment followed by the gait analysis experiment with four fixed receivers and one mobile node. Murata MA40S4S and MA40S4R ultrasonic transducers are used for ranging with chirp signals. These sensors have a narrow bandwidth with 40 kHz center frequency. The transmitter node was designed on an STM32F4-discovery development board which has a form factor, 66 mm × 97 mm. However, in the future, this system can be printed on a flexible wearable band to make it easier and unobtrusive to use. The mobile node was powered by an 11.1v Li-Po battery and a 5v voltage regulator. So the driving voltage of the ultrasonic transmitter was 11.1v. No wired connections were limiting the movement of the subject as we used wireless clock synchronization between the mobile node and the anchor nodes using 2.4 GHz wireless modules. The attachment of the mobile node to the lower limb is as shown in the Figure 5. The chirp signals were stored in the memory of the development board and are sent periodically to the transmitters through digital to analog converters and a driver IC (SN754410) using direct memory access. At the receiver side, we used a bandpass filter, an amplifier and a data acquisition board. ADCs with the same resolution (which was 125 kilo-samples/s in our experiments) can replace the data acquisition board in the future.
The power consumption of the mobile node is important in wireless sensors nodes as the mobile nodes are battery powered. Here the main part of power consumption comes from the memory access and the driving of ultrasonic transducers. We used STM32F4-discovery board to develop the sensor node and SN754410 as the driver for MA40S4S transmitter. Murata MA40S4S transmitters have 120 dB Sound Pressure Level (SPL) at 40 kHz (where, 0 dB = 0.02 mPa). An 11.1v battery and IC-SN754410 drove these transmitters. At ambient temperature, the absolute maximum continuous power dissipation in the driver IC is 2.075 W from the datasheet. The power dissipation in the development board assuming dynamic run mode at 180 MHz clock frequency is 0.22 W. The total power consumption is 2.295 W. Hence, the total energy required for one ranging cycle is 0.0918 J.
Here we neglected the power consumption in the voltage regulators. In the future, custom designed circuits can greatly reduce the power consumption.
An optical motion capture system with six high-speed infrared cameras were used as the reference for the experiments. The reference system was calibrated for static and dynamic conditions and the accuracy was found to be 0.36 ± 17 mm. Reflective markers were used for this system.

Ultrasonic Ranging Experiment
For ranging experiment we used both up-chirp and down-chirp and the average value was calculated using the proposed half peak method (m = 2) and classical peak detection method. The transmitted signals were multiplied by rectangular window functions. Measurements were taken with ten different positions of the transducers each for a low multi-path environment and a high multi-path environment. In the low multi-path environment, the transducers were placed about 750 mm above the ground with a line of sight and without any obstacles in the vicinity. In the high multi-path environment, three obstacles with flat surface were introduced near the transmission path to reflect the ultrasonic waves.

Gait Analysis System
A wireless gait analysis system was designed to test the proposed method. An ultrasonic mobile node was attached to the left lower limb of the subject as shown in Figure 5. For spherical localization, four static anchor nodes at the corners of a 250 × 200 mm rectangle were mounted on a vertical board and kept parallel to the sagittal plane at about 1000 mm from the treadmill as shown in Figure 6. The reflective markers of the motion capture system were attached beside the ultrasonic transmitters and the receivers and the offsets between the markers and ultrasonic transducers were adjusted during the post-processing. Five healthy subjects (four males and one female) with age in between 24 to 31 were recruited for the gait analysis study. Each subject was asked to do a treadmill walk for one minute for each walking speed of 1, 2 and 3 kph with both the ultrasonic mobile node and the camera marker attached to one lower limb. All participants gave their informed consent for inclusion before they participated in the study.  Doppler shift affects the performance of the system. We took the average of range estimated from linearly increasing and decreasing chirp signals to offset these errors as explained in Section 2.1 [23]. Each ranging cycle consists of two ranging each with 20 ms duration, first one with an up-chirp from 39 to 41 kHz in 7 ms and second one with a down-chirp from 41 to 39 kHz in 7 ms. A separation of 13 ms was thus provided between any transmitted signals. Thus the total update rate of the system was about 25 Hz. Both half peak detection method m = 2 and classical peak detection method were implemented for all the range measured from all the four ultrasonic anchor nodes in the gait analysis system. The state space equation for calculating the x, y and z coordinates from the measured distances is shown in the Equation (8).
where, A and g(X i ) are the state transition matrix and the non-linear output function given by Equation (9).
Here, I is a 3 × 3 diagonal matrix and x 1,2,3,4 , y 1,2,3,4 and z 1,2,3,4 represents the coordinates of the anchor nodes. The process noise and the measurement noise are given by q i and r i respectively which are designed as Gaussian function with noise co-variances Q and R respectively. Here, R = diag(e 2 , e 2 , e 2 , e 2 ) and e was set to 5 mm in our system. The state matrix, X i = [x i y i z iẋiẏiżi ] T and measurement matrix, Y i = [d 1,i d 2,i d 3,i d 4,i ] T where x i , y i , z i ,ẋ i , y i andż i represents the three dimensional coordinates and the velocities along the directions and d 1,i , d 2,i , d 3,i and d 4,i represents the measured distances from each anchor nodes at ith instance. An Unscented Kalman Filter (UKF) [35] along with an Unscented Raunch-Tung Striebel smoother [36] was implemented to obtain the three dimensional trajectory from the measured distances and the state space equation explained in Equation (8). All the trajectories were low pass filtered using a butter-worth filter with cut off frequency 10 Hz. The reader can refer [23] for further details of the filters applied.

Estimation of Gait Parameters
We estimated the spatial and temporal gait parameters from the 3D coordinates of one lower limb. The estimation of some gait parameters from the x and z coordinates of a single gait cycle is as shown in the Figure 7. The peaks in the x-coordinate can be assumed as the heel-strikes and the troughs in the x-coordinate can be assumed as the toe-off points [16]. The following parameters were extracted:

Stride Time (ST)
Stride Time is defined as the time taken from one heel-strike to the next heel-strike of the same foot and is given by the time between two minimum points of the x-coordinate.

Stride Length (SL)
Stride length is defined as the distance covered in one stride and can be estimated as twice the difference between the maximum and the minimum point of x-coordinate.

Swing Time (SW)
Time taken from the toe off to next heel strike is defined as swing time.

Stance Time (STT)
Time taken from the heel strike to next toe off is defined as stance time.

Maximum Foot Clearance (MFC)
The foot elevation during swing phase is termed as foot clearance. Maximum foot clearance is the difference between the maximum and the minimum points of z-coordinate during one cycle.

Results and Discussion
Static measurements were taken with one transmitter and one receiver placed at 10 different positions each for a high multi-path and relatively low multi-path environments. The first measurement was subtracted from the remaining nine measurements to obtain the differential distances. The differential distance calculated from the proposed and the classical peak detection method was compared to the camera distance measurements to obtain the mean error. The standard deviation of distance calculated at each static position was also calculated. The results are provided in Table 2. Table 2. Comparison of mean error in differential distance (ME) calculated from the proposed and the classical methods and the standard deviation in the static measurement using up-chirp (S u ) and down-chirp (S d ).

Method
High It can be observed that the we were able to get consistent and considerable improvement in the tracking accuracy using the proposed method. In both high and low multi-path conditions, the average standard deviation for static measurements decreased considerably from 10.44 mm in classical method to 2.45 mm in the proposed method. Figure 8 shows the static measurement reading obtained from the classical peak detection and the proposed half-peak detection methods after removing the offset between the measurements. It is evident that the proposed method provides stable measurements readings compared to the classical method. The mean error and standard deviation of error (STD Error) and root mean square error (RMSE) between the ultrasonic system and the camera system are provided in Table 3. Any temporal delay between the foot trajectory from the ultrasonic system and the the camera system was removed before analysis. The net RMSE from the proposed method was found to be 22.28 mm compared to 28.70 mm from the classical peak detection method

Comparison of Estimated Gait Parameters
To quantify the performance of the proposed system with respect to the reference system, the mean, the standard deviation and the RMSE between the parameters from both the systems during treadmill walk with walking speed 3 kph are calculated and listed in Table 4. The proposed system gave good results for estimation of stance time, swing time and maximum foot clearance. However, our system slightly over-estimated the stride length consistently for all the subjects. This can be solved by removing the constant offset value from all the measurements. A Wilcoxon rank sum test was conducted to test the hypothesis that the gait parameters estimated from the ultrasonic and camera system came from a continuous distribution with equal medians for each subject. The tests failed to reject the null hypothesis with p > 0.05 for all subjects. In the case of stride length, the constant trend was removed before testing the hypothesis.

Conclusions
A new method of multi-path tolerant ultrasonic ranging and localization was proposed. The simulation experiments with earliest 1 m value detection shows that m = 3 and 4 performed better in terms of ratio, R, when the transmitted wave is multiplied by a "hann" or a "hamming" window. However, the maximum value of degradation increased for these window functions and m values. We implemented the half peak detection method with m = 2 and a rectangular window function on chirp signals and ultrasonic transducers. The proposed half-peak detection was found to perform better than the classical peak detection for coherent receiver in ultrasonic ranging system with static transducers as well as the localization system for gait analysis with one moving transmitter. The gait analysis system with one ultrasonic marker attached to one lower limb was found to give better results with net RMSE 22.28 mm for the proposed method compared to 28.70 mm for the classical method. The system gave good results for the estimation of spatial and temporal parameters from the proposed system. In future, the effects of SNR and different values of m on the ranging performance can be studied. The use of multiple sensors to increase the capture volume or use of cone-reflector to make the ultrasonic signals omni-directional [37] are also potential research opportunities in the future.