Abstract
Background
The present study aimed at devising a data processing procedure that provides an optimal estimation of the 3D instantaneous orientation of an inertial measurement unit (IMU). This result is usually obtained by fusing the data measured with accelerometers, gyroscopes, and magnetometers. Nevertheless, issues related to compensation of integration errors and high sensitivity of these devices to magnetic disturbances call for different solutions. In this study, a method based on the use of gyroscope data only is presented, which uses a Weighted Fourier Linear Combiner adaptive filter to perform a driftfree estimate of the 3D orientation of an IMU located on the lower trunk during walking.
Methods
A tuning of the algorithm parameters and a sensitivity analysis to its initial conditions was performed using treadmill walking data from 3 healthy subjects. The accuracy of the method was then assessed using data collected from 15 young healthy subjects during treadmill walking at variable speeds and comparing the pitch, roll, and yaw angles estimated from the gyroscopes data to those obtained with a stereophotogrammetric system. Root mean square (RMS) difference and correlation coefficients (r) were used for this purpose.
Results
An optimal set of values of the algorithm parameters was established. At all the observed speeds, and also in all the various subphases, the investigated angles were all estimated to within an average RMS difference of less than 1.2 deg and an average r greater than 0.90.
Conclusions
This study proved the effectiveness of the Weighted Fourier Linear Combiner method in accurately reconstructing the 3D orientation of an IMU located on the lower trunk of a subject during treadmill walking. This method is expected to also perform satisfactorily for overground walking data and to be applicable also to other “quasiperiodic” tasks, such as squatting, rowing, running, or swimming.
Background
Inertial measurement units (IMUs) have gained in popularity as a tool to quantify human motion [1], thanks to their easeofuse, robust design, and their small dimensions. In the same time, electronics massmarket companies provided low cost devices, such as cellphones, that contain embedded IMU devices along with recording and transmission capabilities. These advantages enable their use for extended periods outside the confines of a laboratory.
An IMU normally includes accelerometers and gyroscopes to measure three accelerations and three angular velocities, respectively. Theoretically, the determination of the position and orientation in space of such a device could be obtained by double and single integration of the above signals, respectively. Unfortunately, the IMU outputs are subject to drift over time, which jeopardizes the time integration of the raw signals [2]. The additional use of magnetometers has been proposed to compensate for the integration errors, but their effectiveness is limited by their high sensitivity to magnetic disturbances [3,4].
The present study focused on the estimate of orientation and aimed at devising a data processing procedure that would compensate for the abovementioned drift and provide an optimal estimation of the 3D instantaneous orientation of the IMU and, therefore, of the body segment it is attached to.
As reported in the literature, the use of recursive filters, such as the Kalman filter, allows the realtime accurate estimate of lower trunk 2D orientation (pitch and roll angles) from the measurement of three accelerations and three angular velocities [5]. A possible alternative is represented by modelbased adaptive filters, which can be used when it is possible to formulate reliable hypotheses about the shape and the time evolution of the signal, which is plausible when the type of motor task performed is known [2]. Among the most popular modelbased adaptive filters are the Fourier Linear Combiner (FLC) filters [6,7]. These filters, which model the measured signal by a Fourier series, are effective when dealing with periodic signals, which is hardly ever the case in biomechanics in general, and in human movement in particular. For this reason, while dealing with the realtime analysis and cancellation movements such as hand tremor, Riviere and his colleagues [8] proposed the use of Weighted Fourier Linear Combiner (WFLC) filters. These filters are an extension of the FLC to be used when investigating signals that display an oscillatory pattern but with a timevarying period. Human walking is a phenomenon that, although for analysis purposes is most of the time hypothesized to be periodic, it does exhibit minor features that vary from stride to stride and, as such, may be defined as quasiperiodic. WFLC filters could hence be considered to be suitable for analyzing gait related data.
This paper proposes the use of a WFLC adaptive filter to perform driftfree 3D orientation angle estimation starting from the measurement of three angular velocities as provided by three orthogonally mounted gyroscopes. The accuracy of the method was assessed using data concerning the orientation of the lower trunk and collected while volunteers walked on a treadmill at different speeds, both at quasi steadystate and in accelerating and decelerating conditions.
Methods
The proposed method for the estimate of the 3D sensor orientation is based on two main steps. First, a tracking of each of the three measured angular velocity components is performed by identifying the corresponding Fourier series coefficients using the WFLC. The identified Fourier series are then analytically integrated to estimate the three orientation angles.
The WFLC filter
As previously mentioned, the WFLC is an adaptive filter that allows the analytical tracking of a quasiperiodic signal. The architecture of the WFLC is presented in Figure 1. The input of the WFLC is the angular velocity signal as measured at the instant of time k (s_{k}). Depending on the instantaneous difference ε_{k} between the signal s_{k} and the output estimated by the WFLC, , the WFLC computes the Fourier series coefficients that will represent the measured signal at time k+1. This result is obtained by adjusting, at each iteration, the socalled filter weights, (the frequency weight, taking into account for the fundamental pulsation) and w_{k} (the vector containing the amplitude weights), using the least mean square (LMS) algorithm proposed by Widrow et al. [9]. Initial values of the weights, and w_{1,} need, of course, to be provided.
Figure 1. Block diagram of the weighted fourier linear combiner filter. Block diagram of the weighted fourier linear combiner filter.
The state vector used by the WFLC is composed of the sine and cosine functions computed using the frequency weight [68]:
where M is the order of the Fourier series representing the measured signal s.
At each instant of time, and are computed using the following equations:
with μ_{0} and μ being the so called frequency and amplitude adaptation gains, respectively. is introduced in the computation of ε_{k} to estimate the bias present in the signal [911], due to possible low frequency components and/or drift, and remove it. It is computed as:
where μ_{b} is the so called bias adaptation gain.
The three adaptation gains, μ_{0}, μ, and μ_{b} determine the convergence time, the accuracy of the algorithm in tracking the measured signal, and the algorithm stability at each sample of time. It has to be noted that the algorithm convergence time cannot be analytically computed [8] and that high values of the gains can improve tracking of the input signal, but can also cause the algorithm to diverge.
In order to allow the use of higher gains, the WFLC can be run twice [8]: the first time using high μ_{0} and low μ values to identify the frequency weight and the second time using the so identified and a higher μ (μ_{FLC}) to identify w_{k}. It has to be noted that the latter use of the WFLC algorithm, in which equation 1c is ignored, makes it equivalent to an FLC [68].
Analytical integration
The estimates of the IMU orientation require an integration of the gyroscopes’ data. Once the angular velocity is tracked with the above algorithms, its analytical and instantaneous representation is provided by the identified Fourier series. The analytical integration of this Fourier series removes drift issues, and can be computed at each sampled instant of time using the following method [2], [11,12]: first the vector is computed using the identified amplitude, and frequency weights:
where f_{s} is the sampling frequency. The instantaneous estimate of the integral of the measured angular velocity is then obtained as [2], [11,12]:
The above angle represents the result of the integration of each single angular velocity component along its corresponding sensor axis, which does not necessarily coincide with the actual rotation of the corresponding axis of the IMU local frame (Figure 2). To overcome this issue, an estimate of the actual rotation can be obtained through a rigid transformation that accounts for the fact that at every instant of time there is a rotation of the local system of reference (x, y, z) fixed to the IMU with respect to the one computed at the previous instant of time.
Figure 2. Results of the sensitivity analysis: initial frequency weight. The figure shows the J_{RMS }values plotted as a function of the observed range of initial frequency weight values (w_{0}).
Experimental session
A total of 18 volunteers (10 males, 8 females, age range: 24–64 years, stature: 1.76±0.09 m, mass: 78±11 kg) were included in the study after signing an informed consent. They were asked to perform three walking trials at natural walking speed (as measured over level ground), 80% and 120% of natural speed on a motorized treadmill. The subjects initially stood on the treadmill, which was then accelerated to the desired velocity. After 35 s of steady state walking, the treadmill was decelerated and stopped for 5 s, before being reaccelerated to the same velocity for additional 35 s and then stopped again, for a total recording time of 80 s. The transition and stopping phases were used to assess the ability of the proposed algorithm to provide accurate estimates also during nonperiodic motion over a short interval of time.
An IMU (Freesense, Sensorize srl) was located on the lower back of the subjects so that the unit local frame (ULF) axes were aligned with the anatomical axes of the lower trunk. In addition, three retroreflective markers were attached to the IMU sensor in order to define a markercluster local frame (MLF). Acceleration and angular velocity data were collected from the IMU (f_{s}=100 samples·s^{1}) while the marker trajectories were tracked by five infrared cameras (MX, Vicon, f_{s}=100 samples·s^{1}).
Pitch, roll and yaw angles, describing the orientation of the ULF, were estimated from the IMU data using the WFLC algorithm and those describing the orientation of the MLF were reconstructed using the stereophotogrammetric system. The time invariant offset of the MLF orientation relative to the ULF orientation was mathematically removed through a rigid transformation while the subject was standing still. In this way both instruments could be assumed to provide pitch, roll, and yaw angles of the same lower trunk anatomical frame. Pitch, roll and yaw angles, can hence be associated with lowertrunk frontal and lateral bending and axial rotation, respectively.
It should be noted that in this study the stereophotogrammetric errors [13] propagate to the angles of interest causing a maximal inaccuracy of 0.5 deg.
Tuning of the algorithm parameters
The order of the Fourier series M, representing the measured signal s, was set to M=1. This conservative choice was initially made following the literature [2,8] and was supported by a preliminary analysis, in which higher values of M caused convergence problems, leading to very inaccurate and unrealistic estimates for some of the investigated trials. Conversely, the choice of M=1 guaranteed the stability of the algorithm for all trials and all subjects.
In order to provide the best tradeoff between accuracy, stability and robustness of the above described method, data recorded from three randomly selected subjects at three different speeds were used to determine the optimal combination of the WFLC gains μ, μ_{0}, μ_{b}. To this purpose, the three combinations of gains that minimize the three cost functions defined for the pitch, roll, and yaw angles, J_{P}, J_{R,}J_{Y}, respectively, were searched for. Each of these cost functions was calculated as the average over the nine trials of the root mean square (RMS) differences between the angle calculated from the stereophotogrammetric data and the one calculated from the IMU.
Knowing that all the gains need to be positive, noting that previous implementations of the FLC algorithm suggest that in most cases a value of μ<0.5 is most likely to ensure convergence [6,7], and on the basis of the authors experience, the following range of values were considered in the gain identification process:
The WFLC algorithm requires an initial estimate of the value of the frequency weight w_{0} (equation 1c and Figure 1). A sensitivity analysis of the effects of varying w_{0} was performed to choose this value by searching for the frequency value that minimized the cost function J_{RMS} = (J_{P} + J_{Y} + J_{R})/3. The tested range of frequencies was 0.15 Hz, which is expected to include all possible walking related frequencies.
Assessment of estimate accuracy
The accuracy of the WFLC in the estimate of the pitch, roll, and yaw angles was assessed by comparing them to the corresponding angles obtained from the stereophotogrammetric data in the 45 trials that were not used in the gains identification process. The RMS, correlation coefficient r, and offset values between estimated and measured angles were calculated for the whole trials (80 s), for the steadystate subphase (25 s), and for the transient subphase, which included the deceleration, stopping and acceleration phases (15 s) (Figure 3).
Figure 3. Angular velocity tracking. Representative results of one walking trial showing the estimated (red line) and measured (black line) angular velocity components.
Results
Results of the gain identification procedure are shown in Figure 4. The red circles indicate the optimum weight combination for pitch (a); roll (b), and yaw (c) angles, respectively. The optimum gain values, which will be used on all the following computations, were: Pitch: μ=0.160, μ_{0}=1e^{5}, μ_{b}=5e^{4}; Roll: μ=0.160, μ_{0}=5e^{5}, μ_{b}=1e^{3}; Yaw: μ=0.135, μ_{0}=1e^{3}, μ_{b}=5e^{4}. The corresponding RMS values were J_{P}=0.8±0.2 deg, J_{R}=0.4±0.1 deg, and J_{Y}=1.0±0.5 deg.
Figure 4. Results of the sensitivity analysis: weight parameters. Results of the gains identification process: values of J_{P }for pitch (a), J_{R }for roll (b), and J_{Y }for yaw (c) angles are shown. Data have been plotted for the combinations of the 3 parameters μ, μ_{0}, μ_{b}, with the corresponding J values represented using a color scale. The red ellipses indicate the optimum weight combinations.
The results of the sensitivity analysis concerning the choice of w_{0} are shown in Figure 2. The values obtained for J_{RMS} ranged from 1.3 to 0.65 deg, with very small variations observed for w_{0} varying in the range of 1–3 Hz, indicating a low sensitivity of the outputs to this initial condition.
The minimum J_{RMS} value was found for w_{0}=2 Hz.
The ability of the WFLC algorithms to track a measured angular velocity is depicted in Figure 3, where the results are shown for one randomly selected entire trial. It should be noted that, after the stopping phase, the algorithm rapidly reconverges as soon as the subject starts walking again.
As shown in Figure 5 for one randomly selected trial, the proposed algorithm led to very satisfactory results in terms of accuracy for the pitch, roll, and yaw angles, not only in the steadystate walking phase (from 10 s to 35 s), but also during the decelerating (from 40 s to 42 s), stopping (from 42 s to 45 s), and accelerating (from 45 s to 48 s) phases (Table 1).
Figure 5. Representative trunk angles estimate. Representative results obtained for one randomly selected trial. The lower trunk orientation angles, as obtained from the stereophotogrammetric system (black curves) and using the proposed method (red curves) are shown. Angles are expressed in the ULF.
Table 1. Results of the accuracy analysis
The mean (standard deviation), RMS, r and offset values for the 45 trials that were not used for the gains identification are reported in Table 1. At all the observed speeds, all the investigated angles were estimated within an average of less than 1.2 deg and with average correlation coefficients greater than 0.90 (with the highest values found for the yaw angles). This applied to both the whole trials and to their subphases. An average offset of less than 3 deg was found, with lowest values observed for the roll angle.
Discussion
The aim of this study was to validate a method based on the use of a WFLC adaptive filter approach, to obtain a driftfree estimate of the 3D orientation of a sensor attached to the lower trunk for a prolonged period of time during treadmill walking, from angular velocities recorded using only one IMU.
A tuning of the WFLC was initially performed, to find optimal values of its gains. A sensitivity analysis was then performed to assess the effects of changes in algorithm frequency weight w_{0}, which is crucial for ensuring that equation 2, and hence the output of the proposed method, are not determinate. Results of this analysis showed that the outputs were always determinate for frequencies ranging between 0.1 Hz and 5 Hz, and that frequencies ranging from 1 Hz to 3 Hz led to very similar results. It has to be noted that these frequencies are actually those expected to be of interest when dealing with human locomotor tasks.
After the above tuning process, the method proved to be very accurate in estimating all the three angles, for all the observed speed conditions and also when the subjects were not walking at steady state. Interestingly, the convergence time of the algorithms, which generally depends on the signal properties, appeared to be negligible for the specific investigated application, as shown by the fact that the results obtained for the transition phases were almost identical to those obtained in the steady state phase (Table 1). This ability of the method to provide accurate angle estimations during nonperiodic motion (acceleration and deceleration phases) and during short intervals of almost no motion (stopping phase) opens the way to future applications, such as uncontrolled walking.
The accuracy of the estimates of lower trunk bending in the sagittal (pitch) and frontal (roll) planes is similar to that obtained in a previous study using a properly optimized Kalman filter [5]. A clear advantage of the proposed method is that, conversely from the previous approach, it uses only the angular velocity signals. Nevertheless, the Kalman filter approach is expected to be more robust for nonperiodic motions than the proposed method, since it does not require any apriori assumption about the signal characteristics.
It has been previously shown that when tracking signals that have a frequency content composed by many frequencies that are close to each other, the performance of the WFLC can be degraded [8]. A Bandlimited Multiple Fourier Linear Combiner (BMFLC) [12], [14] can be used to overcome this problem. However, the BMFLC filter requires an a priori determined set of frequencies, which is not always available when dealing with human movement analysis applications [15,16]. Numerical integration, as associated with WFLCBMFLC adaptive filters, has been recently used successfully for tremor cancelation [14]. This numerical approach, however, requires the use of a highpass filter, which allows easy separation of the tremor oscillations (high frequency) from the voluntary motion (low frequency). Unfortunately, this approach is not suitable for lower trunk angular velocity data recorded during walking, when the determination of highpass filter cutoff frequency is not straight forward due to the variability of walking speed and to the fact that most of the gyroscope signals power is within the low frequencies, which hinder the determination of a proper highpass cutoff frequency
In conclusion, this study proved the effectiveness of the WFLC method in accurately reconstructing the 3D orientation of an IMU located on the lower trunk of a subject during treadmill walking. This method is expected to also perform satisfactorily for overground walking data. The small differences in the values of the measured angular velocities which might be observed between treadmill and level walking data, might require an adjustment of the identified values of the algorithm gains. Further studies are needed to test the suitability of the method for the assessment of pathological gaits and to examine the generalizability of the method to other “quasiperiodic” tasks, such as squatting, rowing, running, or swimming.
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
VB implemented the algorithms, processed the data and contributed to the drafting of the manuscript. CM designed the study, supervised the experimental sessions and contributed to the drafting of the manuscript. JM carried out the experimental sessions and contributed to the drafting of the manuscript. AC contributed to the design and drafting of the manuscript. All authors read and approved the final manuscript.
Acknowledgments
This study is part of the SIVAM project cofunded by the Italian Ministry of Economic Development – ICE (framework agreement with CRUI) – call 2010, by Letsense srl, by ITOP Officine Ortopediche spa, and by the authors’ University.
References

MoeNilssen R, Helbostad JL: Estimation of gait cycle characteristics by trunk accelerometry.
J Biomech 2004, 37:121126. PubMed Abstract  Publisher Full Text

Tan UX, Veluvolu KC, Latt WT, Shee CY, Riviere CN, Ang WT: Estimating displacement of periodic motion with inertial sensors.

Caruso MJ: Applications of magnetic sensors for low cost compass systems. In Proceedings of IEEE Position, Location and Navigation Symp:13–16 March 2000. San Diego, CA; 2000:177184.

Kemp B, Janssen AJMW, van der Kamp B: Body position can be monitored in 3D using accelerometers and earthmagnetic field sensors.
Electroencephalogr ClinNeurophysiol 1998, 109:484488. Publisher Full Text

Mazzà C, Donati M, McCamley J, Picerno P, Cappozzo A: An optimized Kalman filter for the estimate of trunk orientation from inertial sensors data during treadmill walking.
Gait Posture 2011, 35:138142. PubMed Abstract  Publisher Full Text

Vaz CA, Kong X, Thakor NV: An adaptive estimation of periodic signals using a Fourier linear combiner.
IEEE Trans Signal Process 1994, 42:110. Publisher Full Text

Vaz CA, Thakor NV: Adaptive Fourier estimation of time varying evoked potentials.
IEEE Trans Biomed Eng 1989, 36:448455. PubMed Abstract  Publisher Full Text

Riviere CN, Rader RS, Thakor NV: Adaptive cancelling of physiological tremor for improved precision in microsurgery.
IEEE Trans Biomed Eng 1998, 45:839846. PubMed Abstract  Publisher Full Text

Widrow B, Stearns SD: Adaptive Signal Processing. Englewood Cliffs: PrenticeHall; 1985.

Gallego JA, Rocon E, Roa JO, Moreno JC, Pons JL: Realtime estimation of pathological tremor parameters from gyroscope data.
Sensors J 2010, 10:21292149. Publisher Full Text

Veluvolu KC, Ang WT: Estimation of physiological tremor from accelerometers for realtime applications.
Sensors J 2011, 11:30203036. Publisher Full Text

Veluvolu KC, Latt WT, Ang WT: Double adaptive bandlimited multiple Fourier linear combiner for realtime estimation/filtering of physiological tremor.

Chiari L, Della Croce U, Leardini A, Cappozzo A: Human movement analysis using stereophotogrammetry. Part 2: instrumental errors.
Gait Posture 2005, 21:197211. PubMed Abstract  Publisher Full Text

Latt WT, Veluvolu KC, Ang WT: Driftfree position estimation of periodic or quasiperiodic motion using inertial sensors.
Sensors 2011, 11:5931595. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Pecoraro F, Mazzà C, Zok M, Cappozzo A: Assessment of levelwalking aperiodicity.
J Neuroeng Rehabil 2006, 3:28. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Balasubramanian CK, Neptune RR, Kautz SA: Variability in spatiotemporal step characteristics and its relationship to walking performance poststroke.