Hostname: page-component-74656f9884-rdgjk Total loading time: 0 Render date: 2026-11-07T14:56:31.501Z Has data issue: false hasContentIssue false
  • English
  • Français

Optimal Wiener Filter for a Body Mounted Inertial Attitude Sensor

Published online by Cambridge University Press:  26 June 2008

Peter Rizun
Peter Rizun*
Affiliation:
(University of Calgary)
*
(Email: prizun@ucalgary.ca)
Rights & Permissions [Opens in a new window]

Abstract

An optimal attitude estimator is presented for a human body-mounted inertial measurement unit employing orthogonal triads of gyroscopes, accelerometers and magnetometers. The estimator continuously fuses gyroscope and accelerometer measurements together in a manner that minimizes the mean square error in the estimate of the gravity vector, based on known spectral characteristics for the gyroscope noise and the linear acceleration of points on the human body. The gyroscope noise is modelled as a white noise process of power spectral density δn2/2 while the linear acceleration is modelled as the derivative of a band-limited white noise process of power spectral density δv2/2. The estimator is robust to centripetal acceleration and guaranteed to have zero mean error regardless of the motion of the sensor. The mean square angular error in attitude is shown to be independent of the module's angular velocity and equal to 21/2g−1/2δn3/2δv1/2.

Keywords

Attitude estimationInertial/magnetic motion captureMultisensory fusionWiener filtering

Information

Type
Research Article
Information
The Journal of Navigation , Volume 61 , Issue 3 , July 2008 , pp. 455 - 472
Copyright
Copyright © The Royal Institute of Navigation 2008

1. INTRODUCTION

The emergence of microelectromechanical systems, or MEMS, has provided new opportunities in navigation and tracking technology. The small size and low power consumption of MEMS accelerometers and gyroscopes makes their use in human tracking completely feasible. Applications have evolved to include pedestrian tracking in environments with degraded GPS reception (Stirling, Reference Stirling, Fyfe and Lachapelle2005; Godha, Reference Godha, Lachapelle and Cannon2006) as well as complete body motion tracking for medicine in environments where traditional camera-based techniques are impractical (Mayagoitia, Reference Mayagoitia, Nene and Veltink2002; Biosyn, 2007). The dynamic movements of the human body, coupled with the noise characteristics of MEMS-based sensors, impose different requirements on tracking than do inertial navigation techniques on massive ships, airplanes or vehicles. Solutions which are optimal in the latter case are not, in general, optimal in the former case.

An important sub-problem in tracking is the sensing of up and down – or attitude. Earth's gravity presents an extremely robust signal of the vertical direction, deflecting only arc seconds in the presence of Earth's largest mountains (Analytic Sciences, Reference Gelb1974). In MEMS accelerometers, the effect of gravity is detected by observing the force it exerts on a small test mass. Complicating this is the fact that linear acceleration of the sensing structure also produces a force on the test mass, adding noise to the gravity signal and corrupting the measurement of attitude. Thus the best one can hope for is to gain knowledge of the low-frequency components of attitude through a filtering operation intended to attenuate the effect of linear acceleration (Luinge, Reference Luinge and Veltink2004; Bernmark, Reference Bernmark and Wiktorin2002). Fortunately, by employing gyroscopes, the rate of change of angular orientation can be measured directly, providing information on the high-frequency components of attitude. It is reasonable, and has been accomplished in the past (Foxlin, Reference Foxlin1993; Bachman, Reference Bachmann2000; Williamson, Reference Williamson and Andrews2001; Luinge, Reference Luinge2002), to combine the outputs of accelerometers and gyroscopes to produce an estimate of attitude that is more accurate than could be produced by either type of sensor alone.

The primary contribution of this manuscript is a signal processing algorithm for determining attitude from gyroscope and accelerometer measurements in a manner that is optimal. The algorithm is optimal in the sense that it minimizes the mean-square estimation error based on physically reasonable assumptions about the mechanical motion of the inertial measurement unit (IMU) and the noise characteristics of the MEMS sensors. Although similar estimators have been known to the literature since at least as early as 1993, the proposed estimator is novel for two distinct reasons. The first is that it is optimized for linear acceleration modelled as the derivative of a band limited white noise process; support for this model is given from both the theoretical and experimental level. The second reason the filter is novel is because it is shown to be robust to centripetal accelerations; this important filter property has been largely ignored in the existing literature. Secondary contributions include a discussion on how centripetal acceleration leads to bias error in many previous attitude estimators, and a discussion on how filters operating on Euler angles can also lead to biased estimates.

Section 2.1 formalizes the problem of inertial/magnetic orientation sensing and the sub-problem of attitude sensing with gyroscopes and accelerometers. Section 2.2 considers the problem of attitude determination based on the strap-down integration of rate gyroscopes measurements and shows that the presence of measurement noise results in a mean-square error that grows at a rate proportional to the product of the gyroscope noise power spectral density (NPSD) and the elapsed time. This result echoes the well-established fact that attitude estimation based on gyroscope measurements alone will always suffer from drift (Foxlin, Reference Foxlin1996). Section 2.3 considers the power spectrum of human acceleration and reproduces the important and well-known result that the average linear acceleration, expressed in the Earth frame and taken over sufficiently long time periods, tends to zero. This result is then extended to show that the power spectral density of the acceleration at zero frequency must be zero. Models of the acceleration as a white noise process, such as those used in previous attitude estimators (Luinge, Reference Luinge and Veltink2005), are therefore inconsistent with this result and not optimal. This section then shows, based on an energy argument and empirical data, that a more accurate but equally general model for the power spectrum of human acceleration is the derivative of a band-limited white noise process. Such a process has the desired trait that its zero-frequency component is zero and thus the process has, as it should, finite kinetic energy.

The remainder of Section 2 explores two subtle conditions that have not been fully appreciated in the literature under which attitude bias error can occur. Section 2.4 illustrates that although linear acceleration has zero mean value in the Earth frame, when expressed in the sensor frame, linear acceleration has a mean value proportional to the sensor's average centripetal acceleration. Since the accelerometer measurements are naturally expressed in the sensor frame, one must be careful with any type of accelerometer filtering, including “indirect” Kalman filtering discussed in Section 4. Section 2.5 then illustrates how filters operating on the angular output of inclinometers will, in general, lead to bias error when any type of acceleration is present.

Section 3.1 applies a differential coordinate transformation from the sensor frame to what we call the gyroscope-estimated Earth frame. We argue that in this frame, the accelerometer measurements represent the sum of the low-frequency attitude information lost in the gyroscope integration procedure and the corrupting linear acceleration. The problem of attitude estimation is then cast as one of Wiener filtering, separating a signal of known statistical characteristics (the low-frequency attitude information) from noise (the linear acceleration) with equally well known statistics. The optimal filter in the transformed coordinates is shown to be a second-order Butterworth filter acting on the transformed scalar accelerometer measurements. Section 3.3 applies an inverse transformation to recover the optimal estimator operating on the real measurements taken in the sensor coordinate frame. Section 3.4 evaluates the accuracy of the proposed attitude estimator and presents formulae for the mean-square estimation error as well as how the power in the estimation error is distributed over frequency. Section 4 discusses the results in relation to the existing literature. The manuscript concludes in Section 5.

2. FORMULATION OF THE PROBLEM

2.1. Estimator architecture

The estimator explored in this manuscript uses measurements from orthogonal triads of rate gyroscopes, yg, accelerometers, ya, and magnetometers, ym, to estimate its angular orientation relative to a reference frame fixed on earth. The gyroscopes measure the angular velocity vector, ω, distorted by random noise, ng, and a slowly changing offset, b:

(1)
{\bf y}_{g} \equals {\bmomega} \plus {\bf n}_{g} \plus {\bf b}.

The accelerometers measure the sum of the gravity vector, g, and the linear acceleration of the IMU, a:

(2)
{\bf y}_{a} \equals {\bf g} \plus {\bf a} \comma

where the gravity vector is defined to point vertically upwards. (The random noise due to shot effects inside the accelerometers is ignored because in practice it is dwarfed in magnitude by the linear acceleration.) The magnetometers measure the sum of Earth's magnetic field vector, m, a magnetic distortion vector due to nearby ferrous metal or stray magnetic fields, d, and a random noise component, nm:

(3)
{\bf y}_{m} \equals {\bf m} \plus {\bf d} \plus {\bf n}_{m}.

The estimator's purpose is to fuse the three vector measurements to extract the directional information contained in g and m. It is a well-known fact in classical mechanics that the orientation of a rigid body in space has three degrees of freedom and thus requires three scalar parameters to fully describe its angular orientation. The directions of the gravity vector and magnetic field vector provide a total of four parameters, and thus over-constrain the problem. To eliminate the extra constraint, only the angle of the magnetic field vector projected onto the plane normal to the gravity vector – the heading – is used.

Once the direction of the gravity vector and the heading are known in the sensor coordinate frame, an orthogonal matrix, , specifying the coordinate system of the sensor relative to the Earth can be constructed using well-established techniques. From this matrix, angular information such as Euler angles can be easily extracted. Refer to (Stirling, Reference Stirling, Fyfe and Lachapelle2005) for a graphical explanation of how the directions of the gravity vector and magnetic north are used to define orientation.

The task of orientation estimation can logically be divided into the sub-problems of estimating the direction of the gravity vector or attitude, and estimating a scalar that specifies the heading, as shown in Figure 1. This paper focuses on the former, leaving the equally challenging problem of the latter to be addressed in future work. We also proceed assuming that the gyroscope offset estimator shown in the figure is functioning, leaving out the details on how this can be accomplished in practice. For brevity, we will drop the subscripts on ya and ng for the remainder of the paper.

Figure 1.

Orientation estimator architecture. This manuscript addresses the attitude estimator boxed by the dashed lines.

2.2. Noise power spectrum of a gyroscope orientation estimate

An orthogonal triad of gyroscopes mounted to a rigid structure measures the structure's angular velocity vector. In this section, we consider integrating the angular velocity measurement as the sensor structure rotates to correctly adjust the estimate of the gravity vector. The basic kinematical law upon which the dynamical equations of motion for a rigid body are founded states that the rate of change of a vector in a stationary reference frame is related to its rate of change in a rotating frame by the operator equation (Goldstein, Reference Goldstein1980)

(4)
{}\vskip-12pt^{G}\left( {{d \over {dt}}} \right) \equals {}\vskip-12pt^{S}\left( {{d \over {dt}}} \right) \plus {\bmomega} \times

acting on the given vector. The G superscript denotes the stationary ground or Earth frame and the S superscript denotes the rotating sensor frame. Since the gravity vector is static in the ground frame, G(d g/dt)=0. Applying Equation (4) to the gravity vector furnishes its equation of motion in the rotating sensor frame,

(5)
{{d^{S} {\bf g}} \over {dt}} \equals \minus{}^{S}{\bmomega} \times {}^{S}{\bf g}.

The gravity vector, as observed from the sensor frame at any later time t, regardless of the sensor's motion, follows from the integration of Equation (5):

(6)
{{}^{S}{\bf g}\lpar t \rpar \equals \minus \int_{\setnum{0}}^{t} {}^{S}\bmomega \lpar t\rpar \times {}^{S}{\bf g}\lpar t\rpar dt \plus {}^{S}{\bf g}\lpar 0\rpar}.

Equation (6) can be implemented easily using analogue or digital means. In practice, this gyroscope integration technique fails for two distinct reasons. The first is that the constant of integration, Sg(0), is unknown and thus the estimate of Sg(t) will be biased. The second reason is that only an estimate of the angular velocity,

(7)
\hat{\bmomega }{\bf \equals }\bmomega \plus {\bf n}\comma

is available, which contains the true angular rate plus a random variable that represents measurement noise. As the noise component integrates, it produces an error that walks randomly on the surface of a sphere.

In the frequency domain, each scalar component of n is assumed to be Gaussian white noise of power spectral density

(8)
S_{n} \lpar \omega \rpar \equals \delta _{n}^{\setnum{2}} \sol 2.

S n(ω) is defined for both positive and negative frequencies but is written such that δn has units of radians per second per root Hertz for positive frequencies only (i.e. δn is the value commonly quoted in sensor data sheets). (Note the unfortunate notation where ω refers the angular velocity vector while ω refers to angular frequency.) It is a well-known result that integration in the time domain corresponds to division by iω in the frequency domain, and that power spectral density of a signal is equal to the square of its magnitude. If we denote the orientation error, ΔΦ, as a vector that specifies a small rotation about each of the coordinate axes, its power spectral density, for short time periods, becomes:

(9)
{S_{{\rmDelta}{\rmPhi}}} \lpar \omega \rpar \equals {{\delta_{n}^{\setnum{2}} \sol 2} \over {\omega ^{\setnum{2}} }}.

The double pole at zero frequency in Equation (9) indicates that the gyroscope-estimated coordinate system rotates slowly away from the true coordinate system, as expected. Using techniques well established in Norbert Wiener's famous treatise on time series (Wiener, Reference Wiener1949), one can show that the mean square value of each component of the orientation error grows at a rate proportional to the noise power spectral density of the gyroscopes and the elapsed time:

(10)
\left\langle {\rmDelta}{ \rmPhi_{i}}\lpar t\rpar ^{\setnum{2}} \right\rangle \equals {{\delta _{n}^{\setnum{2}} t} \over 2}\comma \quad t \ges 0.

Since only two scalar components of ΔΦ contribute to the attitude error, the resulting root mean square (RMS) attitude error, ΔΩRMS, again only valid for short periods of time, becomes

(11)
{\rmDelta}{\rmOmega}_{\rm RMS} \lpar t\rpar \equals \delta_{n} \sqrt t \comma \quad t \ges 0.

From a frequency domain perspective, gyroscope integration provides excellent knowledge of the high frequency part of attitude but lacks the low frequency part – the error – which has a power spectrum given by Equation (9). As previously outlined, we hope to recover the low frequency part from the accelerometer measurements. To proceed further, the power spectrum of human acceleration needs to be considered.

2.3. The power spectrum of human acceleration

The IMU's acceleration, as measured in the ground frame, can be thought of as noise that corrupts the desired gravity signal. The average acceleration over a time interval T is given by

(12)
\bar{\bf a} \equals {{{\bf v}\lpar o\rpar \minus {\bf v}\lpar t \minus T\rpar } \over T}.

Since the velocity of the IMU can be assumed finite, the average acceleration tends to zero as the time interval becomes infinitely long:

(13)
\mathop {\lim }\limits_{T \to \infty } \bar{\bf a} \equals\hskip-2pt \mathop {\lim }\limits_{T \to \infty }\hskip-2pt \left( {{{\bf v}\lpar o\rpar \minus {\bf v}\lpar t \minus T\rpar } \over T} \right) \hskip-2pt\equals 0.

Equation (13) is another way of saying that the zero-frequency component of the acceleration is zero.

Exactly how long the averaging time needs to be to eliminate a specific portion of the acceleration is linked to the power spectral density of the sensor's acceleration. An approximation to the spectrum can be made by reasoning that a human can generate kinetic energy equally well over a band-limited frequency interval between zero and an upper frequency limit of human motor control, ωc. Research on human motion also supports this assumption (Burdea, Reference Burdea1996). Since kinetic energy is proportional to the square of velocity, the velocity of the sensor can be modelled as band-limited white noise with power spectral density

(14)
S_{v} \lpar \omega \rpar \equals {{\delta _{v}^{\setnum{2}} \sol 2} \over {1 \plus \lpar \omega \sol \omega_{c} \rpar ^{\setnum{2}} }}.

The constant δv is proportional to the standard deviation of the sensor velocity but can equally well be thought of as a parameter that specifies the intensity of the motion, or the slosh as per the work of Eric Foxlin (Foxlin, Reference Foxlin1996). Finally, since acceleration is the derivative of velocity, the acceleration has a power spectral density given by the derivative of a band-limit white noise process,

(15)
S_{a} \lpar \omega \rpar \equals {{\lpar \delta _{v}^{\setnum{2}} \sol 2\rpar \omega ^{\setnum{2}} } \over {1 \plus \lpar \omega \sol \omega_{c} \rpar ^{\setnum{2}} }}.

At frequencies much less than ωc, Equation (15) can be approximated as

(16)
S_{a} \lpar \omega \rpar \equals \lpar \delta _{v}^{\setnum{2}} \sol 2\rpar \omega ^{\setnum{2}}.

Figure 2 shows the assumed power spectral density of human velocity and acceleration. If we imagine measuring the acceleration with body-mounted accelerometers, then these measurements, when expressed in the Earth frame, will contain acceleration signals with spectral characteristics described by Equation (15) added to a constant vertical gravity vector.

Figure 2.

Simplified velocity and acceleration power spectrum of points on a human body.

Before deriving the optimum attitude estimator in Section 3, we take a detour for the remainder of this section to examine two ways that bias can be introduced into our estimate.

2.4. Bias error from centripetal acceleration

Although the acceleration in the ground frame averaged over a long enough time tends to zero, this result does not apply to time averages of acceleration expressed in the sensor frame. Consider an IMU accelerating and rotating in an arbitrary manner in three-dimensional space. The sensor-frame acceleration is related to the earth-frame acceleration by Sa(t)=SRG(t)Ga(t), where SRG(t) is the orthonormal matrix that represents a coordinate tranformation from the ground to sensor frame. The time average of the sensor-frame acceleration is (dropping the explicit time dependence)

(17)
{1 \over T}\int_{ \minus T}^{\setnum{0}} {{}^{S}{\bf a}dt} \equals {1 \over T}\int_{ \minus T}^{\setnum{0}} {{}^{S}{\bf R}_{G} {}^{G}{\bf a}dt}.

The right hand side of Eq. (18) can be expanded using the technique of integration by parts: \int {^{S} {\bf R}_{G}\hskip0 ^{G} {\bf a}} dt \equals \int {^{S} {\bf R}_{G}\hskip0 ^{G} \dot{\bf {v}}} dt \equals ^{S} {\bf R}_{G}\hskip0 ^{G} {\bf v} \minus \int {^{S} \dot{\bf {R}}_{G}\hskip0 ^{G} {\bf v}} dt. Recognizing ^{S} \dot{\bf {R}}_{G} \equals \minus ^{S} {\bmomega } \times ^{S} {\bf R}_{G} and Sv=SRGGv, one obtains

(18)
{1 \over T}\int_{ \minus T}^{\setnum{0}} {{}^{S}{\bf a}dt} \equals {{{}^{S}{\bf v}\lpar 0\rpar \minus {}^{S}{\bf v}\lpar \minus T\rpar } \over T} \plus {1 \over T}\int_{ \minus T}^{\setnum{0}} {{}^{S}{\bmomega } \times {}^{S}{\bf v}dt}.

As the averaging time T tends to infinity, the boundary term vanishes because the velocity of the sensor must remain finite, as per our discussion in Section 2.3. The second term represents the centripetal acceleration of the sensor, which in general, has a non-zero mean value:

(19)
\eqalign{\tab \overline{{}^{S}{\bf a}} \equals \mathop {\lim }\limits_{T \to \infty } {1 \over T}\int_{ \minus T}^{\setnum{0}} {{}^{S}{\bf a}dt} \equals \mathop {\lim }\limits_{T \to \infty } \left\{ {{{{}^{S}{\bf v}\lpar 0\rpar \minus {}^{S}{\bf v}\lpar \minus T\rpar } \over T} \plus {1 \over T}\int_{ \minus T}^{\setnum{0}} {{}^{S}{\bmomega } \times {}^{S}{\bf v}dt} } \right\}. \cr \tab \equals \overline{^{S} {\bmomega } \times ^{S} {\bf v}} \cr}

What this reveals is that any direct averaging of the sensor-mounted accelerometer outputs will result in a steady-state error when centripetal acceleration is present. Further implications of this are explored in Section 4.

2.5. Filters operating on angles introduce bias

Since an estimate of angular orientation is ultimately sought, one is tempted to design a filter that operates on angles. However, in circumstances when the sensor experiences acceleration comparable to the acceleration due to gravity, the use of angles may lead to significantly biased results. Consider the case, as shown in Figure 3, where the sensor is subject to simple harmonic motion along a 45° inclined plane with peak accelerations of \sqrt 2 g. The mean acceleration of the sensor is of course zero and a suitable low-pass filter acting on the accelerometer output would attenuate the acceleration leaving behind an unbiased estimate of the gravity vector. If the angle θ were (linearly) filtered instead, the average angle estimate would clearly be biased towards the large (apparent) pitch angle near the block's highest point.

Figure 3.

Sinusoidal acceleration along an inclined plane. The mean acceleration is zero; however, the mean pitch angle found by averaging the apparent pitch angle is not equal to the pitch of the inclined plane.

3. OPTIMAL WIENER FILTER FOR ATTITUDE ESTIMATION

3.1. Optimal fusion of accelerometer and gyroscope measurements to infer attitude as a problem in Wiener filtering

The strap-down integration procedure outlined in Section 2.2 results in an angular error with a power spectral density given by Equation (9), while using all of the information available from the gyroscope measurements. It is logical to try to recover the missing information from the accelerometer measurements expressed in the gyroscope-estimated Earth frame. As we will see, the fusing of angular velocity becomes implicit in the coordinate transformation with this method. The problem can then be posed as one of Wiener filtering – separating a signal with known power spectrum from noise with known power spectrum. To begin, we envision a transformation from the sensor frame, where the measurements are available, into the gyroscope-estimated Earth frame via the differential transformation Equation (4). To continue, we need to make two subtle simplifications in order to formulate the Wiener filtering problem. The first is to notice that the error in the gyroscope's gravity vector estimate via Equation (9) is in the form of an angle, constraining the vector to a random walk along the surface of a sphere. Instead, we must express this error as a random vector added to the tip of the gravity vector, thereby allowing it to walk randomly in three-dimensional space. The length of this random vector is chosen such that its projection onto the sphere is equal to the radius over which the gravity vector could have moved along the sphere's surface. The gravity vector, expressed in the gyroscope-estimated Earth frame, can thus be considered composed of three scalar components each with power spectral density given by

(20)
S_{g} \lpar \omega \rpar \equals g^{\setnum{2}} {{\delta_{n}^{\setnum{2}} \sol 2} \over {\omega ^{\setnum{2}} }}

where g is the magnitude of the gravity vector. This simplification is of little consequence as small changes to the radial component of the estimated gravity vector have negligible effect on the estimation of attitude.

The second simplification involves assuming that the acceleration spectrum in the ground frame is equal to the acceleration spectrum in the equivalent frame estimated by the gyroscopes. Because the two frames differ only by a slow drift, only the very low-frequency components of the acceleration are affected by this simplification. Since the low-frequency components of the acceleration in the Earth frame have vanishing magnitude, this assumption also has negligible effect on the attitude estimate.

With these two assumptions, accelerometer measurements expressed in the gyroscope-estimated ground frame, S ag(ω), contain the desired signal. S g(ω), as well as the corrupting acceleration noise, S a(ω),

(21)
\eqalign{S_{ag} \lpar \omega \rpar \tab \equals g^{\setnum{2}} {{\lpar \delta_{n}^{\setnum{2}} \sol 2\rpar } \over {\omega ^{\setnum{2}} }} \plus {{\lpar \delta _{v}^{\setnum{2}} \sol 2\rpar \omega ^{\setnum{2}} } \over {1 \plus \lpar \omega \sol \omega_{c} \rpar ^{\setnum{2}} }} \equals {{\delta _{v}^{\setnum{2}} \omega _{c}^{\setnum{2}} \omega ^{\setnum{4}} \plus g^{\setnum{2}} \delta _{n}^{\setnum{2}} \omega ^{\setnum{2}} \plus g^{\setnum{2}} \delta _{n}^{\setnum{2}} \omega _{c}^{\setnum{2}} } \over {2\omega ^{\setnum{2}} \lpar \omega ^{\setnum{2}} \plus \omega _{c}^{\setnum{2}} \rpar ^{\setnum{2}} }} \cr \tab\equals {{\delta_{v}^{\setnum{2}} \omega_{c}^{\setnum{2}} } \over 2}{{\prod \nolimits_{n \equals \setnum{0}}^{\setnum{3}} {\lpar \omega \minus z_{n} \rpar } } \over {\prod\nolimits_{m \equals \setnum{0}}^{\setnum{3}} {\lpar \omega \minus p_{m} \rpar } }}.}

From the work of Norbert WienerFootnote 1 (Wiener, Reference Wiener1949) we know that the optimum filter for separating the desired signal, S g(ω), from the noise, S a(ω), is given by the solution to the Wiener-Hopf equation for the problem.

(22)
H\lpar \omega \rpar \equals {1 \over {S_{ag}^{ \plus } \lpar \omega \rpar }}\left[ {{\rm causal\ part\ of}{{S_{g} \lpar \omega \rpar } \over {S_{ag}^{ \minus } \lpar \omega \rpar }}}\right]

where S ag+(ω) is the part of S ag(ω) that has all its poles and zeros in the upper half plane (and hence corresponds to the causal part of S ag(ω)), and S ag−(ω) is the part of S ag(ω) that has all its poles and zeros in the lower half plane (and hence corresponds to the anti causal part of S ag(ω)).

The optimum filter, H(ω), can be found by first splitting S ag(ω) into causal and anti causal factors. Because the numerator of S ag(ω) is quadratic in ω2, it can be written as

(23)
{{\delta _{v}^{\setnum{2}} \omega _{c} ^{\setnum{2}} } \over 2}\left( {\omega ^{\setnum{2}} \minus Z_{a} } \right)\left( {\omega ^{\setnum{2}} \minus Z_{b} } \right)

Using the quadratic formula, one then obtains zeros of ω2 at

(24)
Z \equals \minus {{g^{\setnum{2}} \delta _{n}^{\setnum{2}} } \over {2\delta _{v}^{\setnum{2}} \omega _{c}^{\setnum{2}} }} \pm \sqrt {{{g^{\setnum{4}} \delta _{n}^{\setnum{4}} } \over {4\delta _{v}^{\setnum{4}} \omega _{c}^{\setnum{4}} }} \minus {{g^{\setnum{2}} \delta _{n}^{\setnum{2}} } \over {\delta _{v}^{\setnum{2}} }}}

With the definitions

(25)
\omega_{g} \equals \sqrt {{{g\delta _{n} } \over {\delta_{v} }}} \quad {\rm and} \quad \rmDelta \omega_{g}^{\setnum{2}} \equals {{g^{\setnum{2}} \delta _{n}^{\setnum{2}} } \over {2\delta_{v}^{\setnum{2}} \omega_{c}^{\setnum{2}} }} \equals {{\omega _{g}^{\setnum{4}} } \over {2\omega _{c}^{\setnum{2}} }}

and making the assumption that \omega _{c}} \gt \omega_{g}} \sol \sqrt 2, Equation (24) can be written as

(26)
Z_{a} \equals Z_{b} ^{\ast } \equals \minus \rmDelta \omega _{g}^{\setnum{2}} \plus i\sqrt {\omega _{g}^{\setnum{4}} \minus \rmDelta \omega _{g}^{\setnum{4}} } \equals \omega _{g}^{\setnum{2}} e^{ \plus i\left( {{\pi \sol {\setnum{2} \plus \rmDelta \phi }}} \right)}

where

(27)
\rmDelta \phi \equals \arcsin \left( {{{\rmDelta \omega _{g}^{\setnum{2}} } \over {\omega _{g}^{\setnum{2}} }}} \right)

The zeroes of ω follow easily by taking the square root of the polar form of Equation (26) defining Z0 and Z1 to be causal zeros:

(28)
z_{n} \equals \left\{ {\matrix{ \hskip-2pt{\omega _{g} e^{i\lpar \setnum{2}n \plus \setnum{1}\rpar \pi \sol \setnum{4}} e^{i\rmDelta \phi \sol \setnum{2}} \comma \quad \hskip7pt n \equals 0\comma 2} \cr {\hskip-2pt\omega _{g} e^{i\lpar \setnum{2}n \plus \setnum{1}\rpar \pi \sol \setnum{4}} e^{ \minus i\rmDelta \phi \sol \setnum{2}} \comma \quad n \equals 1\comma 3} \cr} } \right..

The poles of S ag(ω) include a double pole at the origin and symmetric poles at ±iωc:

(29)
p_{m} \equals \left\{ {\matrix{ 0 \comma\hfill\tab { m \equals 0\comma 2} \cr {\left( { \minus 1} \right)^{{{m \minus \setnum{1}} \over \setnum{2}}} i\omega _{c} \comma} \tab { m \equals 1\comma 3} \cr} } \right.

Splitting S ag(ω) into causal and anti causal factors yields

(30)
S_{ag} \lpar \omega \rpar \equals \underbrace {\left( {{{\delta _{v}^{\setnum{2}} \omega _{c} ^{\setnum{2}} } \over 2}{{\lpar \omega \minus z_{\setnum{0}} \rpar \lpar \omega \minus z_{\setnum{1}} \rpar } \over {\omega \lpar \omega \minus p_{\setnum{1}} \rpar }}} \right)}_{S_{{ag}}^{ \plus } \lpar \omega \rpar }\underbrace {\left( {{{\lpar \omega \minus z_{\setnum{2}} \rpar \lpar \omega \minus z_{\setnum{3}} \rpar } \over {\omega \lpar \omega \minus p_{\setnum{3}} \rpar }}} \right)}_{S_{{ag}}^{ \minus } \lpar \omega \rpar}

Dividing the power spectrum of the signal by the anti causal part of S ag(ω) yields

(31)
\eqalign{ {{S_{g} \lpar \omega \rpar } \over {S_{ag}^{ \minus } \lpar \omega \rpar }} \tab \equals{{g^{\setnum{2}} \delta _{n} ^{\setnum{2}} } \over {2\omega ^{\setnum{2}} }}{{\omega \lpar \omega \minus p_{\setnum{3}} \rpar } \over {\lpar \omega \minus z_{\setnum{1}} \rpar \lpar \omega \minus z_{\setnum{3}} \rpar }} \cr \tab \equals \underbrace { \minus {{{{g^{\setnum{2}} \delta _{n}^{\setnum{2}} p_{\setnum{3}} } \sol {z_{\setnum{2}} z_{\setnum{3}} }}} \over {2\omega }}}_{{\rm causal\ part}} \plus \underbrace {{{{{g^{\setnum{2}} \delta _{n}^{\setnum{2}} \lpar z_{\setnum{2}} \minus p_{\setnum{3}} \rpar } \sol {z_{\setnum{2}} \lpar z_{\setnum{2}} \minus z_{\setnum{3}} \rpar }}} \over {2\lpar \omega \minus z_{\setnum{1}} \rpar }} \plus {{{{g^{\setnum{2}} \delta _{n}^{\setnum{2}} \lpar z_{\setnum{3}} \minus p_{\setnum{3}} \rpar } \sol {z_{\setnum{3}} \lpar z_{\setnum{3}} \minus z_{\setnum{2}} \rpar }}} \over {2\lpar \omega \minus z_{\setnum{3}} \rpar }}}_{{\rm anticausal\ part}} \cr}

where the pole at the origin is taken, without consequence, to be causal. Applying Equation (22) gives

(32)
H\lpar \omega \rpar \equals \left( {{2 \over {\delta _{v}^{\setnum{2}} \omega _{c} ^{\setnum{2}} }}{{\omega \lpar \omega \minus p_{\setnum{1}} \rpar } \over {\lpar \omega \minus z_{\setnum{0}} \rpar \lpar \omega \minus z_{\setnum{1}} \rpar }}} \right)\left( { \minus {{{{g^{\setnum{2}} \delta _{n}^{\setnum{2}} p_{\setnum{3}} } \sol {z_{\setnum{2}} z_{\setnum{3}} }}} \over {2\omega }}} \right)

which after simplification yields

(33)
H\lpar \omega \rpar \equals {{\omega _{g}^{\setnum{2}} \lpar 1 \plus i\omega \sol \omega _{c} \rpar } \over { \minus \omega ^{\setnum{2}} \plus \sqrt 2 \omega _{g} \left( {\cos {{\rmDelta \phi } \over 2} \plus \sin {{\rmDelta \phi } \over 2}} \right)i\omega \plus \omega _{g}^{\setnum{2}}}}

In practice, however, ωg=|z n| is much smaller than ωc and thus the zero takes effect at frequencies where the gain of the filter is extremely small. If we ignore this zero by allowing the ratio ωc/ωg to tend to infinity then

(34)
H\lpar \omega \rpar \approx {{\omega _{g}^{\setnum{2}} } \over { \minus \omega ^{\setnum{2}} \plus \sqrt 2 \omega_{g} i\omega \plus \omega _{g}^{\setnum{2}} }}\comma

which is a second-order filter bearing the name of Stephen Butterworth who described its properties in a paper published in 1930 (Butterworth, Reference Butterworth1930). The Butterworth filter is optimally flat in the pass band, which has a certain intuitive appeal. What optimally flat means is that the rate of change of the filter gain with frequency is as small as possible near ω=0. For the present problem, we can take this to mean that the filter passes as completely as possible the desired gravity vector signal (centred around ω=0), and then falls off at 40 decibels per decade to attenuate the corrupting noise from the linear acceleration. The filter break point represents the optimal balance between losing too much of the desired gravity information and picking up too much of the corrupting linear acceleration.

3.2. The optimal Wiener filter expressed in the sensor coordinate frame

In the previous section, we used Equation (4) to transform to a coordinate system where the problem of attitude estimation could be formulated as one of Wiener filtering. Now that the optimal filter is known, we wish to transform this filter back to the sensor coordinate system where it can be applied directly to the gyroscope and accelerometer measurements. To do so, we apply the inverse transformation

(35)
{}\vskip-12pt^{S}\left( {{d \over {dt}}} \right) \equals {}\vskip-12pt^{G}\left( {{d \over {dt}}} \right) \minus {\hat{\bmomega}} \times

to our expression for the optimal filter in the time domain. The optimal Wiener filter implicit in Equation (34) can be cast into the Laplace domain by making the substitution iω→s, where s is the Laplace variable. In the Earth frame, the filter acts on (the scalar components of) the accelerometer measurement vector to produce the gravity vector estimate:

(36)
H \lpar s\rpar {\bf y}\lpar s\rpar \equals {\bf \hat{g}} \lpar s\rpar

where g(s) and y(s) represent the Laplace transforms of the corresponding signals in the time domain. Rearranging Eq. (36) yields

(37)
\lpar s^{\setnum{2}} \plus \sqrt 2 \omega_{g} s \plus \omega_{g}^{\setnum{2}} \rpar {\bf \hat{g}} \lpar s\rpar \equals \omega_{g}\hskip0pt^{\setnum{2}} {\bf y}\lpar s\rpar

from which the time domain filter equations can be written directly as

(38)
\overbrace {{1 \over {\omega_{g}^{\setnum{2}} }}{{d\hskip1.5pt^{\setnum{2}} {\bf \hat{g}}} \over {dt^{\setnum{2}} }} \plus {{\sqrt 2 } \over {\omega _{g} }}{{d{\bf \hat{g}}} \over {dt}} \plus {\bf \hat{g}} \equals {\bf y}}^{{\rm gyroscope} {\hbox{-}}{\rm estimated} \ {\rm ground} \ {\rm frame}}\hskip-2pt.

Inverse transformation by Eq. (35) produces the desired filter:

(39)
\overbrace {{1 \over {\omega _{g}^{\setnum{2}} }}{{d\hskip1.5pt^{\setnum{2}} {\bf \hat{g}}} \over {dt^{\setnum{2}} }} \plus {{\sqrt 2 } \over {\omega _{g} }}{{d{\bf \hat{g}}} \over {dt}} \plus {\bf \hat{g}} \equals {\bf y} \minus {2 \over {\omega _{g}^{\setnum{2}} }}\hat{\bmomega } \times {{d{\bf \hat{g}}} \over {dt}} \minus {1 \over {\omega _{g}^{\setnum{2}} }}{{d\hat{\bmomega }} \over {dt}} \times {\bf \hat{g}} \minus {1 \over {\omega _{g}^{\setnum{2}} }}\hat{\bmomega } \times \hat{\bmomega } \times {\bf \hat{g}} \minus {{\sqrt 2 } \over {\omega _{g} }}\hat{\bmomega } \times {\bf \hat{g}}}^{{\rm sensor}\ {\rm frame}}.

To implement this filter, it is convenient to express it as two first-order vector differential equations. It is easy to show that Equation (39) can be decomposed into

(40)
\eqalign{\tab {{d{\bf \hat{g}}_{\setnum{1}} } \over {dt}} \equals {{\omega _{g} } \over {\sqrt 2 }}\lpar 2{\bf y} \minus {\bf \hat{g}}_{\setnum{1}} \minus {\bf \hat{g}}\rpar \minus \hat{\bmomega } \times {\bf \hat{g}}_{\setnum{1}} \cr \tab {{d{\bf \hat{g}}} \over {dt}} \equals {{\omega _{g} } \over {\sqrt 2 }}\lpar {\bf \hat{g}}_{\setnum{1}} \minus {\bf \hat{g}}\rpar \minus \hat{\bmomega } \times {\bf \hat{g}} \cr}

where ĝ1 is an intermediate estimate of the gravity vector and ĝ is the desired optimal estimate. The flow diagram for the computation implicit in Equation (40) is shown in Figure 4.

Figure 4.

Flow diagram for the optimal Wiener filter in the sensor coordinate frame for attitude estimation based on the fusion of gyroscope and accelerometer measurements. The constant A=2−1/2ωg.

3.3. Estimator accuracy

The filter transformation possible via Equation (4) is also useful because it allows us to analyze estimator performance in the ground frame where the estimator equations are simple second-order filters acting on the scalar components of the transformed accelerometer output. The first interesting observation is that is absent from Equation (38); in fact, the only property of the gyroscopes that affects the filter accuracy is the noise power spectral density on the choice of estimator pole locations. This means that the performance of the attitude estimator is not affected by the angular velocity of the sensorFootnote 2.

The accuracy of the estimator can be determined by considering the power spectral density of the estimate error. The error spectrum (Simon, Reference Simon2006) for the scalar components of the gravity vector is given by

(41)
S_{\rmDelta g} \lpar \omega \rpar \equals \lsqb 1 \minus H\lpar \omega \rpar \rsqb \lsqb 1 \minus H\lpar \minus \omega \rpar \rsqb S_{g} \lpar \omega \rpar \plus H\lpar \omega \rpar H\lpar \hskip-1pt\minus \omega \rpar S_{a} \lpar \omega \rpar.

This quantity can be thought of as the sum of the signal power (gravity vector) that is rejected by the filter and the noise power (linear acceleration) that is accepted by the filter. Using Equations (9), (22) and (41), the error spectrum becomes

(42)
S_{\rmDelta g} \lpar \omega \rpar \equals g\delta _{n} \delta _{v} \left[ {{{1 \plus \lpar \omega \sol \omega _{g} \rpar ^{\setnum{2}} } \over {1 \plus \lpar \omega \sol \omega _{g} \rpar ^{\setnum{4}} }}} \right]\comma

which is plotted in Figure 5. The spectrum has a power density of gδnδv at zero frequency, a peak value of {1 \over 2}\lpar 1 \plus \sqrt 2 \rpar g\delta_{n} \delta _{v} at \omega \equals \sqrt {\sqrt 2 \minus 1} \omega _{g} and falls off to zero at 20 dB/decade at higher frequencies. Most of the error is therefore concentrated at frequencies between zero and a few times ωg, allowing clean estimates of derivatives to be made. This is advantageous to camera-based motion capture techniques where the angular error – arising from positional error in marker localization – is more or less equi-distributed over frequency, resulting in noisy first and second derivatives.

Figure 5.

Noise power spectral density of gravity vector estimate for positive frequencies.

The mean square estimation error is found by integrating the error spectrum over all frequencies,

(43)
\left\langle {\rmDelta g^{\setnum{2}} } \right\rangle \equals {1 \over {2\pi }}\int_{ \minus \infty }^{\infty } {S_{\rmDelta g} \lpar \omega \rpar d\omega }.

Making the substitution χ=ω/ωg and dω=ωgdχ we get

(44)
\left\langle {\rmDelta g^{\setnum{2}} } \right\rangle \equals {{g\delta _{n} \delta _{v} \omega _{g} } \over {2\pi }}\int_{ \minus \infty }^{\infty } {{{1 \plus \chi ^{\setnum{2}} } \over {1 \plus \chi ^{\setnum{4}} }}} d\chi.

The definite integral can be evaluated using residue theory

(45)
\int_{ \minus \infty }^{\infty } {{{1 \plus \chi ^{\setnum{2}} } \over {1 \plus \chi ^{\setnum{4}} }}} d\chi \equals \sqrt 2 \pi

and the mean square error of the scalar components becomes

(46)
\left\langle {\rmDelta g^{\setnum{2}} } \right\rangle \equals {1 \over {\sqrt 2 }}g\delta _{n} \delta _{v} \omega _{g} \equals {1 \over {\sqrt 2 }}g^{\setnum{3}\sol \setnum{2}} \delta _{n}^{\setnum{3}\sol \setnum{2}} \delta _{v}^{\setnum{1}\sol \setnum{2}}.

In practice, the gyroscope noise can be assumed constant whereas the noise caused by linear acceleration corrupting the gravity signal depends strongly on the motion of the sensor. The mean square estimation error can be determined again using Equation (41), however, this time the acceleration power spectrum is scaled by the positive quantity a 2=(δv′/δv)2

(47)
S_{a} \lpar \omega \rpar \equals {{a^{\setnum{2}} \delta _{v}^{\setnum{2}} \omega ^{\setnum{2}} } \over 2}.

After some simplification, the mean square error, now a function of the parameter a, works out to

(48)
\left\langle {\rmDelta g^{\setnum{2}} } \right\rangle \equals {{g\delta _{n} \delta _{v} \omega _{g} } \over {4\pi }}\int_{ \minus \infty }^{\infty } {{{2 \plus \lpar 1 \plus a^{\setnum{2}} \rpar \chi ^{\setnum{2}} } \over {1 \plus \chi ^{\setnum{4}} }}} d\chi.

The integral can again be evaluated using residue theory resulting in a mean square error of

(49)
\left\langle {\rmDelta g^{\setnum{2}} } \right\rangle \equals \left( {{{3 \plus a^{\setnum{2}} } \over {4\sqrt 2 }}} \right)g\delta _{n} \delta _{v} \omega _{g} \equals \left( {{{3 \plus a^{\setnum{2}} } \over {4\sqrt 2 }}} \right)g^{\setnum{3}\sol \setnum{2}} \delta _{n}^{\setnum{3}\sol \setnum{2}} \delta _{v}^{\setnum{1}\sol \setnum{2}}

If we assume that the error in the scalar components of the gravity vector estimate are small compared to g, then the error in attitude can be approximated as

(50)
\left\langle {\rmDelta \rmOmega ^{\setnum{2}} } \right\rangle \approx {{\rmDelta g_{ \bot }\vskip0pt^{\setnum{2}} } \over {\lpar g \plus \rmDelta g_{\parallel} \rpar ^{\setnum{2}} }} \approx {{2\rmDelta g^{\setnum{2}} } \over {g^{\setnum{2}} }} \equals \left( {{{3 \plus a^{\setnum{2}} } \over {2\sqrt 2 }}} \right){{\delta _{n}^{\setnum{3}\sol \setnum{2}} \delta _{v}^{\setnum{1}\sol \setnum{2}} } \over {g^{\setnum{1}\sol \setnum{2}} }}\comma

which is plotted in Figure 6.

Figure 6.

Mean square attitude error at various slosh levels for a filter optimized for slosh at the fixed level δv (solid) and for a filter optimized at the arbitrary slosh level δv′ (dashed).

We conclude this section by considering the root mean-square error for an attitude sensor that we have built using off-the-shelf components, and assuming slosh levels that correspond to vigorous human motion such as running with an IMU attached to the ankle. The results are presented in Table 1 and reveal that even under conditions of high linear acceleration, the achievable accuracy is better than 1° and likely suitable for most purposes.

Table 1.

Maximum achievable performance of typical MEMS attitude sensors.

4. DISCUSSION

4.1. Estimator architecture

The orientation sensor presented in Figure 1 divides the problem into the three natural sub-problems of estimating attitude, heading, and gyroscope offsets. Such a division makes understanding the inner workings of the filter easier without degrading filter performance. For instance, the heading estimate does not provide any additional information to help solve the attitude estimation problem beyond what was implicitly used to estimate the gyroscope offsets.

4.2. Linear acceleration power spectrum for human motion

Noise is often the result of random thermal motions of charge carriers in the electronic sensor chips. The energy associated with this random motion is commonly equi-distributed in frequency and manifests itself as a random voltage with a certain RMS amplitude per root Hertz. In the attitude estimation problem, the electrical noise on the accelerometer outputs is dominated by the “noise” caused by the mechanical motion of the sensor. We argued that, to a first approximation, a human tends to generate kinetic energy equi-distributed over a band-limited range of frequencies. From this, it was shown that velocity can then be modelled as a band-limited white noise process and linear acceleration as the derivative of such a process. In previous works, (for instance: Luinge, Reference Luinge and Veltink2005; Yun, Reference Yun, Aparicio, Bachmann and McGhee2005), the acceleration is modelled directly as a white noise process, implying unrealistically that the root mean square velocity of the IMU was infinite.

It is worthwhile to consider the accuracy of a filter designed for linear acceleration with a white spectrum when subject to the more realistic spectrum proposed in this manuscript. Such a filter is H′(ω)=ωg′/(iω+ωg′). With the proposed optimal second-order filter, H(ω), the mean square error of the gravity vector components, 〈Δg 2〉, was shown to be proportional to gδnδvωg. With the first-order filter, H′(ω), and setting ωg′=ωg, one can easily show that the mean square error will instead be proportional to gδnδvωc, where ωc is the frequency limit of human motion. The first order filter will thus produce a mean square error that is larger by a factor proportional to ωc/ωg. As Table 1 showed, ωg will often be on the order of 10−1 rad/s, whereas the frequency limit of human motion (ωc) is on the order of 101 rad/s such that ωc/ωg≈100. Note, that the first-order filter can be optimized with a significantly smaller cut-off frequency than the one previously assumed. Using the optimal cut-off frequency, the mean square error for the first-order filter will be larger by the reduced factor of (ωc/ωg)1/3. In practice one cannot achieve this performance because it requires a cut-off frequency so low that un-modelled drift in the gyroscope offsets become the dominant source of error. We thus expect considerably larger mean square errors with a filter optimized for acceleration with a white power spectrum.

4.3. Bias error from centripetal acceleration

In Section 2.4 we showed that the mean value of linear acceleration, when expressed in the sensor coordinate frame, is non-zero and depends on the centripetal acceleration of the sensor. This is an important result because it can be used to clarify some confusion in the literature regarding the role of accelerometers. The accelerometer measurements, expressed in the earth frame, indeed do provide low-frequency orientation information. The measurements contain the sum of the gravity signal plus the zero-mean linear acceleration. However, since the actual accelerometer measurements are made in the sensor frame, the mean value is shifted by an amount equal to the centripetal acceleration. Any averaging of the accelerometers then permanently fuses the gravity vector data (attitude) with the mean centripetal acceleration, causing uncorrectable bias error. The statement that when “averaged for a sufficiently long period, the output of an accelerometer triad can be used to measure the components of the gravity vector” (Bachmann, Reference Bachmann, Yun, McKinney, McGhee and Zyda2003), or that the accelerometers provide a good mean value (Vanagay, Reference Vaganay, Aldon and Fournier1993) must be used with care.

Filtering of the accelerometer outputs, or values derived from their outputs, is effectively an averaging operation. Such an operation will therefore cause bias error in the mean estimate of the gravity vector by an amount that depends on the level of centripetal acceleration and the strength of the filter employed. The result has implications to what have become known as indirect filters (Setoodah, Reference Gelb2004; Luinge, Reference Luinge and Veltink2005). The indirect Kalman filter operates on only the difference between the output of an inclinometer and a gyroscope-based attitude sensor. The point is that when centripetal acceleration is present, the mean value of the inclinometer will in general contain bias error and thus the output of the Kalman filter will likewise contain bias error. Only a filter operating directly and continuously on the accelerometers and gyroscopes can hope to remove centripetal acceleration in all cases.

Many indirect formulations have been discussed in the literature. Setoodah (Setoodah, Reference Setoodeh, Khayatian and Farjah2004) argues that the indirect filtering approach “has a superior performance from the viewpoint of computational complexity.” In (Yun, Reference Yun, Lizarraga, Bachmann and McGhee2003; Reference Yun, Aparicio, Bachmann and McGhee2005), the investigators use a quaternion-based Kalman filter that fuses an accelerometer-estimated orientation with a gyroscope-estimated orientation. The author is not aware of any analysis done that includes the filter's robustness to centripetal acceleration. Earlier work (Foxlin, Reference Foxlin1996) recommended that the drift from angular velocity integration only be corrected with accelerometers “during low acceleration periods” because of the problem of linear acceleration folding into the inclination estimate and causing bias error. However, we have shown that by fusing the gyroscope and accelerometer information in the method proposed in this manuscript, the filter will perform well and always extract useful information from the accelerometers regardless of the intensity of the motion.

4.4. Filter accuracy

An important result of this manuscript is an analytic description of the accuracy that can be achieved through accelerometer/gyroscope attitude measurements. Table 1 shows that with off-the-shelf electronics, the estimation error can quite easily be kept below 1 degree even during vigorous human motion. This makes one question the use of more advanced techniques such as covariance estimation through wavelet decomposition (Seetodah, Reference Setoodeh, Khayatian and Farjah2004). Our experience with this type of sensor is that a significant source of error is due to many subtle but predictable effects. These include non-orthogonality of the sensing elements, non-collocation of the accelerometer axes, non-simultaneous sampling of the individual sensing elements, and bias drift on the accelerometer outputs. Except for the last, each of these effects can be calibrated for and completely removed. The last effect can be significantly reduced by using accelerometers with small bias drift and good temperature stability.

5. CONCLUSION

This manuscript presented an optimal attitude estimator for a human body-mounted inertial measurement unit employing orthogonal triads of gyroscopes, accelerometers and magnetometers. It also presented formulae that describe the spectrum and root mean square value of the estimation error. From a historical perspective, the proposed attitude estimator is an extension of Wiener optimal filtering applied to inertial attitude sensing. A differential coordinate transformation was found that cast the problem of optimal attitude estimation through multi-sensory fusion as the solution to three non-coupled Wiener-Hopf equations. It was shown that linear filtering of the accelerometer measurements in the transformed coordinates eliminated bias error from centripetal acceleration, representing a new contribution to the literature. The two pieces of statistical information (human acceleration power spectrum and gyroscope noise characteristics) needed to apply Wiener filtering were defined. During this process, it was observed that, in a very general sense, the kinetic energy of the inertial measurement unit should be considered a band limited white noise process. The corollary of this is that the power spectrum of the linear acceleration should be modelled as the derivative of such a process, also representing a new contribution to the literature. The Wiener filtering problem was then solved to obtain the frequency domain representation for the optimal filter operating on the transformed accelerometer output. The optimal filter was shown to be second order – a consequence of modelling the acceleration as the derivative of a white noise process. The inverse Laplace transform was then used to determine the estimator's time-domain equations. Finally, the inverse differential coordinate transformation was employed to obtain the optimal multi-sensory fusion algorithm that can be directly applied to the accelerometer and gyroscope measurements. The accuracy of the proposed filter was then explored using conventional techniques. It was shown that root mean square attitude errors less than 1° were obtainable during conditions of vigorous human activity using off-the-shelf MEMS gyroscopes.

This manuscript did not discuss how the heading can best be inferred, nor how the gyroscope offset can be continually estimated. The former is a much more difficult problem than the latter, due to the distortion caused by nearby stray magnetic fields. Investigators (Roetenberg, Reference Roetenberg2006; Reference Gelb2007) have attacked this problem using a Kalman filter that estimates the “quality” of the magnetic field and then adjusts the blending between gyroscope and magnetometer to rely heavily on the gyroscopes when the magnetic field is judged to contain significant distortion. The limitation of this technique is that it is not useful for correcting static distortions of the magnetic field. In order for it to work effectively, the Kalman gain-weighted average of the magnetic field direction must be correct and the duration of the distortions must be sufficiently short such that the compounded error due to the integration of gyroscope noise does not become significant. If the IMU moves permanently into an area where the magnetic field is twisted, the heading error will eventually become equal to the magnitude of the twist. The problem of magnetic distortion will undoubtedly be an area of future endeavour for inertial/magnetic tracking technology.

Footnotes

1 Because the autocorrelation functions of S a(ω) or S g(ω) do not exist, this case is not strictly covered by the theory of Wiener filtering. However, both S a(ω) or S g(ω) can be considered limiting cases of processes that do possess autocorrelation functions. The former can be considered the output of a first order low pass filter of vanishing cutoff frequency suject to a white noise input, while, to the later, one can add poles at a frequency far greater than those that are important to the problem.

2 This assumes that the angular velocity estimates, except for the white noise corruption, are otherwise perfect. Specifically, it assumes that the gyroscope gain error is zero, that the gyroscope are not sensitive to linear acceleration, and that the gyroscope signals do not saturate for use outside their operating range. In practice, the performance does deteriorate slightly with increases in angular velocity.

References

REFERENCES

Analytic Sciences Corporation, Technical Staff (1974). Applied Optimal Estimation. Gelb, Arthur (Ed.). MIT Press.Google Scholar
Bachmann, E. (2000). Inertial and magnetic angle tracking of limb segments for inserting humans into synthetic environments. PhD thesis. Naval Postgraduate School. Monterey, California, USA.Google Scholar
Bachmann, E., Yun, X., McKinney, D., McGhee, R. and Zyda, M. (2003). Design and Implementation of MARG Sensors for 3-DOF Orientation Measurement of Rigid Bodies. Proceedings of the 2003 IEEE International Conference On Robotics & Automation. Taipei, Taiwan, September 14–19.Google Scholar
Bernmark, E., and Wiktorin, C. (2002). A triaxial accelerometer for measuring arm movements. Appl. Ergon., 33, 541–547.CrossRefGoogle ScholarPubMed
Biosyn Systems Inc. (2007). http://www.biosyn.ca.Google Scholar
Burdea, G. (1996). Force and Touch Feedback for Virtual Reality. Wiley: New York.Google Scholar
Butterworth, S. (1930). On the theory of filter amplifiers. Wireless Engineer, 7, 536–541.Google Scholar
Foxlin, E. (1993). Inertial Head-Tracking. M.Sc. Thesis, Department of Electrical Engineering and Computer Science, MIT, MA, USA.Google Scholar
Foxlin, E. (1996). Inertial head-tracker sensor fusion by a complementary separate-bias Kalman filter. Proceedings of VRAIS.Google Scholar
Godha, S., Lachapelle, G. and Cannon, M. E. (2006). Integrated GPS/INS system for pedestrian navigation in a signal degraded environment. Proceedings of ION GNSS, Fort Worth, Texas.Google Scholar
Goldstein, H. (1980). Classical Mechanics (second edition). Addison-Wesley.Google Scholar
Luinge, H. J., Veltink, P. H. (2005). Measuring orientation of human segments using miniature gyroscopes and accelerometers. Med. Biol. Eng. Comput., 43, 273–282.CrossRefGoogle ScholarPubMed
Luinge, H. J., and Veltink, P. H. (2004). Inclination measurement of human movement using a 3D accelerometer with autocalibration. IEEE Trans. Neural Syst. Rehabil. Eng., 12, 112–121.Google Scholar
Luinge, H. (2002). Inertial sensing of human motion. PhD thesis. University of Twente, Netherlands.Google Scholar
Mayagoitia, R. E., Nene, A. V., Veltink, P. H. (2002). Accelerometer and rate gyroscope measurement of kinematics: an inexpensive alternative to optical motion analysis systems. Journal of Biomechanics, 35, 537–542.Google Scholar
Roetenberg, D. (2006). Inertial and magnetic sensing of human motion. PhD thesis. University of Twente, Netherlands.Google Scholar
Roetenberg, D., Slycke, P. and Veltink, P. (2007). Ambulatory position and orientation tracking fusing magnetic and inertial sensing. IEEE Transactions on Biomedical Engineering. 54 (5) 883–890.Google Scholar
Setoodeh, P., Khayatian, A., Farjah, E. (2004). Attitude estimation by separate-bias kalman filter-based data fusion. The Journal of Navigation, 57, 261–273.CrossRefGoogle Scholar
Simon, D (2006). Optimal State Estimation: Kalman, H-infinity, and Nonlinear Approaches. John Wiley & Sons.CrossRefGoogle Scholar
Stirling, R., Fyfe, K. and Lachapelle, G. (2005). Evaluation of a new method of heading estimation for pedestrian dead reckoning using shoe mounted sensors. The Journal of Navigation, 58, 31–45.CrossRefGoogle Scholar
Vaganay, J., Aldon, M. J., Fournier, A. (1993). Mobile robot attitude estimation by fusion of inertial data. Proceedings of the International Conference on Robotics and Automation.Google Scholar
Wiener, N. (1949). Extrapolation, interpolation, and smoothing of stationary time series with engineering applications. MIT Press.CrossRefGoogle Scholar
Williamson, R. and Andrews, B. J. (2001). Detecting absolute human knee angle and angular velocity using accelerometers and rate gyroscopes. Med. Biol. Eng. Comput., 39, 294–302.Google Scholar
Yun, X., Lizarraga, M., Bachmann, E. R. and McGhee, R. B. (2003). An improved quaternion-based Kalman filter for real-time tracking of rigid body orientation. Proceedings of the International Conference on Intelligent Robots and Systems, Las Vegas, Nevada.Google Scholar
Yun, X., Aparicio, C., Bachmann, E. R. and McGhee, R. B. (2005). Implementation and experimental results of a quaternion-based Kalman filter for human body motion tracking. Proceedings of the International Conference on Robotics and Automation, Barcelona, Spain.Google Scholar
Figure 0

Figure 1. Orientation estimator architecture. This manuscript addresses the attitude estimator boxed by the dashed lines.

Figure 1

Figure 2. Simplified velocity and acceleration power spectrum of points on a human body.

Figure 2

Figure 3. Sinusoidal acceleration along an inclined plane. The mean acceleration is zero; however, the mean pitch angle found by averaging the apparent pitch angle is not equal to the pitch of the inclined plane.

Figure 3

Figure 4. Flow diagram for the optimal Wiener filter in the sensor coordinate frame for attitude estimation based on the fusion of gyroscope and accelerometer measurements. The constant A=2−1/2ωg.

Figure 4

Figure 5. Noise power spectral density of gravity vector estimate for positive frequencies.

Figure 5

Figure 6. Mean square attitude error at various slosh levels for a filter optimized for slosh at the fixed level δv (solid) and for a filter optimized at the arbitrary slosh level δv′ (dashed).

Figure 6

Table 1. Maximum achievable performance of typical MEMS attitude sensors.