2.4. IMU Data Processing

This protocol is extracted from research article:

Orientation-Invariant Spatio-Temporal Gait Analysis Using Foot-Worn Inertial Sensors

**
Sensors (Basel)**,
Jun 7, 2021;
DOI:
10.3390/s21113940

Orientation-Invariant Spatio-Temporal Gait Analysis Using Foot-Worn Inertial Sensors

DOI:
10.3390/s21113940

Procedure

The following sections describe the methods for IMU data processing, inspired by typical gait analysis routines—including the stages of zero velocity detection, orientation estimation, double integration, events detection, and gait parameters estimation. The processing pipeline included evaluation of different approaches for orientation estimation and double integration. Every stage was formulated using an orientation-invariant approach.

Zero velocity intervals (ZVIs) were detected using the angular rate energy detector [28]. According to [28], angular rate provides rich information concerning the detection of ZVIs. Compared to other methods (e.g., the acceleration magnitude detector), the angular rate energy detector achieved the highest performance [28].

To calculate the energy of the angular rate magnitude, we used a sliding window with $0.15$ s. The size of the window was experimentally set to ensure an appropriate energy result, i.e., not too smoothed, nor too noisy. The window size of $0.15$ s resulted in a good compromise that could adapt to all walking speeds considered. To determine ZVIs, a threshold was applied to the energy of the angular rate magnitude. The threshold was calculated with basis on the average of the energy, to ensure that differences in walking speed (reflected as different signal amplitudes) would be considered. The threshold was experimentally set at $1/8$ of the average, to ensure that all strides in all velocities would be captured by the method.

After roughly detecting ZVIs, we proceeded with a dynamic adjustment of the intervals. For each ZVI, we did a threshold refinement search: starting with a very low threshold, we have progressively increased it (in steps of $1/20$ of the ZVI energy range), until a minimum interval size was obtained. The minimum interval size was set at $0.1$ s, unless an interval lower than $0.1$ s was registered in the trial. Using this process, we ensured a confident detection of all strides and, additionally, that all ZVIs were refined to potentially include only zero velocity instants. The remaining intervals, i.e., the moving intervals, were used to evaluate movement in the subsequent processing stages.

Orientation, or attitude, of the sensor relative to the global frame of reference was expressed using quaternions. The global frame of reference was defined by two perpendicular horizontal axes (x and y)—arbitrarily set—and a vertical axis, z, pointing to the sky. To obtain quaternions, we tested three methods: gyroscope integration and two complementary filters (CFs)—Madgwick and Euston.

The gyroscope integration method takes advantage of the knowledge of the moving and not-moving intervals. When the sensor is not moving (i.e., during a ZVI)—when all measured accelerations are due to the Earth’s gravity acceleration—the accelerometer signal is used to estimate sensor inclination. The vector ${\mathit{g}}_{s}=[{\overline{a}}_{x},{\overline{a}}_{y},{\overline{a}}_{x}]$, defined with the average acceleration values measured while the sensor is not moving, corresponds to the z-axis in the global frame of reference, i.e., the gravity axis. The horizontal axis—y-axis, perpendicular to gravity—is then defined arbitrarily. The initial quaternion is obtained from gravity and the horizontal vector, using the tri-axial attitude determination (TRIAD) algorithm [29].

The initial quaternion is updated each time the foot is in contact with the ground. During the intervals when the foot is moving, the quaternion $\mathit{q}\left(t\right)$ is updated resorting to the integration of angular rate $\mathbf{\omega}\left(t\right)$ measured by the gyroscope, as described in [8,13]. This is performed as defined by Equations (1) and (2):

where $\dot{\mathit{q}}\left(t\right)$ is the quaternion derivative and $\Delta t$ is the sampling interval. The function $\mathit{p}(.)$ denotes the quaternion representation of a vector and the ⊗ operator represents the quaternion product.

In Madgwick [30], the quaternion derivative $\dot{\mathit{q}}\left(t\right)$ used in Equation (2) is replaced by a corrected estimate $\dot{\widehat{\mathit{q}}}\left(t\right)$ that incorporates orientation information provided by the accelerometer. The method is based on the calculation of the gradient descent, as shown in Equation (3).

where *J* denotes the Jacobian, $\beta $ is the filter gain and $\mathbf{\epsilon}\left(t\right)$ denotes the error term, obtained by subtracting the accelerometer measurement in the sensor frame ${\mathit{a}}_{s}\left(t\right)$ to the theoretical gravity vector in sensor coordinates ${\mathit{g}}_{s}\left(t\right)$ (obtained by transforming gravity in global coordinates to sensor coordinates). $\beta $ can be defined as $\beta =\sqrt{3/4}\phantom{\rule{4pt}{0ex}}{\overline{\omega}}_{max}\pi /180$, where ${\overline{\omega}}_{max}$, expressed in degrees, represents the maximum gyroscope measurement error (mean zero gyroscope measurement error).

The Madgwick filter is applied at all instants of the signal (moving and not moving intervals) considering as basis the initial quaternion determined using TRIAD, as described previously.

The explicit complementary filter, also known as Euston filter, was implemented as described in [13,18]. In Euston, instead of replacing the value of the quaternion, the measured angular rate $\mathbf{\omega}\left(t\right)$ is replaced by a corrected angular rate signal, resulting in the following filter dynamics:

in which the error term $\delta $ is obtained following Equations (5) and (6), where the term $\mathit{e}\left(t\right)$ describes the angular mismatch between theoretical (${\mathit{g}}_{s}\left(t\right)$) and measured (${\mathit{a}}_{s}\left(t\right)$) direction of gravity [18].

The Euston filter has two adjustable parameters, the proportional gain ${k}_{P}$—to separate low- and high-frequency estimates of orientation—and the integrator gain ${k}_{I}$—to compensate for gyroscope bias [13,18]. Similarly to the Madgwick filter, we apply Euston to all instants of the signal, considering as basis the initial quaternion determined using TRIAD [29].

After obtaining orientation quaternions, $\mathit{q}\left(t\right)$, we calculate linear acceleration in global coordinates, hereinafter represented as $\mathit{a}\left(t\right)$ for simplicity. To that purpose, we first estimate the value of the gravity vector in global coordinates, ${\mathit{g}}_{w}=[0,0,{\overline{a}}_{zv}]$, where ${\overline{a}}_{zv}$ is the average acceleration magnitude measured during ZVIs. Linear acceleration is then obtained as shown in Equation (7).

where ${\mathit{a}}_{s}\left(t\right)$ represents raw acceleration, as measured by the sensor.

To obtain displacements, we integrate linear acceleration two times. On the first integration, an estimate of velocity, $\widehat{\mathit{v}}\left(t\right)$, is obtained. Integrals are computed using the Trapezoidal Rule, as shown in Equation (8).

To bound the errors, two different methods—linear dedrifting and direct and reverse integration—are employed, as we explain next. After obtaining trajectories, a novel approach to correct the final vertical position (assuming walking on a flat surface) is tested. The method rotates the trajectories so that the final height of each stride is zero. To that purpose, a rotation quaternion is calculated, using as basis the angle with the horizontal plane at the end of the stride and the rotation vector calculated as the cross product between the vertical axis and the stride displacement vector. The rotation quaternion is used to rotate trajectories within each moving interval.

Double integration is performed between ZVIs, on a stride-by-stride basis. To fulfil the zero-velocity assumption—on which moving intervals are bounded by zero velocity instants—a linear drift function (${\mathit{d}}_{v}\left(t\right)$) is estimated and subtracted from the estimated velocity, as described in [14] and shown in Equation (9).

Trajectory $\mathit{s}\left(t\right)$ is obtained by integrating again velocity $\mathit{v}\left(t\right)$.

The direct and reverse integration method fuses the regular integral with a time-reversed integral so that the boundary conditions, in this case, the zero-velocity conditions, are satisfied in the initial and final values of the integral [15,16]. The result of direct (${v}_{\to}\left(t\right)$) and reverse (${v}_{\leftarrow}\left(t\right)$) integration is combined using a sigmoid weighting function ($w\left(t\right)$), as shown in Equation (10).

The sigmoid $w\left(t\right)$, specified in Equation (11), is shaped using the steepness parameter, $\eta $, and the inflection point, ${t}_{i}$, defined between the temporal bounds ${t}_{n}$ and ${t}_{n+1}$ of each moving interval. To define ${t}_{i}$, a proportion ${\alpha}_{i}$, between 0 and 1, of the moving interval is used.

Position is estimated by integrating velocity.

To detect gait events avoiding the need of determining angular rate or acceleration in body coordinates (where the alignment of the sensor on the body would need to be known [8,15,16]), we used acceleration magnitude and the vertical component of acceleration (in global coordinates).

We observed that FO events can generally be found in an acceleration magnitude perturbation before swing. To approximate this instant, we filtered acceleration magnitude using a zero-lag bidirectional 2nd order Butterworth low-pass filter with a cutoff of 7 Hz. The cutoff frequency was experimentally chosen to ensure that (i) the perturbations were smoothed resulting in a single peak value, and (ii) the resultant peak approximates the acceleration magnitude perturbation where annotated FO events are observed. The FO event was considered the first maximum peak above the average of the moving interval appearing between two ZVIs (i.e., the first value above the average surrounded by two values with lower magnitude), as illustrated in Figure 3.

Events detection from low pass filtered acceleration magnitude and low pass filtered vertical acceleration (*FO*—Foot off; *FC*—Initial foot contact; *MS*—Mid-stance; *Ann-FO*—Annotated FO; *Ann-FC*—Annotated FC).

FCs are detected between FOs and the beginnings of the next ZVIs. FCs were considered the absolute minimum of vertical acceleration measured between these two instants (Figure 3). Before detecting FC events, vertical acceleration was low-pass filtered using a zero-lag bidirectional 1st order Butterworth filter with a cutoff of 30 Hz. The cutoff frequency was experimentally chosen to ensure the attenuation of high frequency noise that could hinder the detection of FC events.

After calculating orientation, position and determining FO and FC events, temporal and spatial parameters were estimated for each gait cycle *n*. Temporal parameters—stride, swing and stance duration—were determined as defined in [14]. Cadence was obtained as the inverse of stride duration, expressed in steps per minute. Spatial parameters (illustrated in Figure 2) were calculated using information of moving intervals, defined by the temporal bounds of ${t}_{n}$ and ${t}_{n+1}$. To estimate SL and SW, we used trajectories on the horizontal plane, ${\mathit{s}}_{xy}$, as determined by Equations (12) and (13), where ${\overrightarrow{\mathit{s}}}_{n}\left(t\right)$ represents a displacement vector relative to the final stride position at ${t}_{n+1}$, obtained as ${\mathit{s}}_{xy}\left({t}_{n+1}\right)-{\mathit{s}}_{xy}\left(t\right)$). In Equation (13), the symbol *∠* denotes the angle between two vectors.

Gait speed was obtained by dividing SL by its corresponding stride duration.

To calculate MTC, we used a method inspired on the work by Kanzler, C. [9]. To estimate toe trajectory, we have first estimated the distance between the sensor and the toe, *r*, using as a basis the angle produced by the foot at FO, $\alpha \left(F{O}_{n}\right)$, in each gait cycle *n*, as illustrated in Figure 4.

Variables involved in the calculation of toe trajectory.

To obtain the angle $\alpha \left(t\right)$, we did a series of vector transformations. First, we have converted the vertical vector $[0,0,1]$ to sensor coordinates, using the quaternion at the beginning of the moving interval, i.e., at the foot flat at ${t}_{n}$. Then, we transformed the resultant vector back to global coordinates using the quaternion at FO. This vector, ${\overrightarrow{v}}_{w}$, was used to estimate the medio-lateral vector, ${\overrightarrow{l}}_{w}$, using the cross product between ${\overrightarrow{v}}_{w}$ and $[0,0,1]$. The forward vector, ${\overrightarrow{f}}_{w}$, was calculated using the cross product between $[0,0,1]$ and ${\overrightarrow{l}}_{w}$, which was then converted back to sensor coordinates using the quaternion at ${t}_{n}$. This vector, ${\overrightarrow{f}}_{s}$, parallel to the ground at foot flat and pointing forward towards the toes, was used to estimate the angle $\alpha \left(t\right)$, as depicted in Equation (14).

The distance between the sensor and the toe (*r*) was obtained by the average of the values determined in each stride *n*, as shown in Equation (15).

where *N* is the total number of strides and ${\mathit{s}}_{z}\left(t\right)$ represents the z component of the trajectory of the sensor (i.e., its vertical displacement). MTC was considered the minimum peak vertical toe displacement measured during the swing phase of walking. This vertical displacement, $m\left(t\right)$, was estimated as shown in Equation (16), where $r\times sin\left(\alpha \right(t\left)\right)$ represents the vertical distance between the sensor and the toe (shown as $j\left(t\right)$ in Figure 4).

To calculate turning angles, we converted an arbitrary horizontal vector (e.g., the vector $[0,1,0]$) to sensor coordinates, using as basis the quaternions estimated at ${t}_{n}$ and at ${t}_{n+1}$. The resultant vectors represented the orientation of the sensor in the horizontal plane, so that the angle between these two vectors corresponded to the turning angle.

Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Note: The content above has been extracted from a research article, so it may not display correctly.

Q&A

Your question will be posted on the Bio-101 website. We will send your questions to the authors of this protocol and Bio-protocol community members who are experienced with this method. you will be informed using the email address associated with your Bio-protocol account.