Open AccessArticle

Keeping a Good Attitude: A Quaternion-Based Orientation Filter for IMUs and MARGs

Roberto G. Valenti

¹,

Ivan Dryanovski

² and

Jizhong Xiao

^1,*

The City College of New York, The City University of New York, Convent Avenue and 140th Street, New York, NY 10031, USA

The Graduate Center, The City University of New York, 365 Fifth Avenue, New York, NY 10016, USA

Author to whom correspondence should be addressed.

Sensors 2015, 15(8), 19302-19330; https://doi.org/10.3390/s150819302

Submission received: 24 June 2015 / Accepted: 27 July 2015 / Published: 6 August 2015

(This article belongs to the Special Issue Inertial Sensors and Systems)

Download

Browse Figures

Figure 1
Simulated experiment showing the orientation output in quaternion form (left) and in Euler representation (right). (a) Orientation reference; (b) Orientation output of the presented method with noisy acceleration data and noise-free magnetometer data; (c) Orientation output of the presented method with noise-free acceleration data and magnetic readings affected by magnetic disturbances. "> Figure 2
Block diagram of the complementary filter for IMU implementation (no magnetometer data). "> Figure 3
Block diagram of the complementary filter for MARG implementation (with magnetometer data). "> Figure 4
2-D representation of two different interpolation methods between quaternion qa and qb resulting in quaternion qint. On the left is represented the linear interpolation (LERP), whereas on the right the spherical linear interpolation (SLERP) showing the correct result. "> Figure 5
The gain factor function used in the adaptive gain to reduce the attitude error deriving from the external linear acceleration because of vehicle motion. "> Figure 6
2D top view (left) and 3D side view (right) of the trajectory traveled by the quadrotor during the “1LoopDown” experiment. "> Figure 7
Comparison of the orientation Euler angles output of the three methods during the “1LoopDown” experiment. (a) Roll; (b) Pitch; (c) Yaw. "> Figure 8
Comparison of the different estimation methods in the presence of magnetic disturbances. (a) Norm of the measured magnetic field vector; (b) Roll; (c) Pitch; (d) Yaw. "> Figure 9
Pitch estimation under the condition of high nongravitational acceleration; (a) Non gravitational acceleration applied during the experiment; (b) Pitch angle output of the different algorithms ">

Versions Notes

Abstract

Orientation estimation using low cost sensors is an important task for Micro Aerial Vehicles (MAVs) in order to obtain a good feedback for the attitude controller. The challenges come from the low accuracy and noisy data of the MicroElectroMechanical System (MEMS) technology, which is the basis of modern, miniaturized inertial sensors. In this article, we describe a novel approach to obtain an estimation of the orientation in quaternion form from the observations of gravity and magnetic field. Our approach provides a quaternion estimation as the algebraic solution of a system from inertial/magnetic observations. We separate the problems of finding the “tilt” quaternion and the heading quaternion in two sub-parts of our system. This procedure is the key for avoiding the impact of the magnetic disturbances on the roll and pitch components of the orientation when the sensor is surrounded by unwanted magnetic flux. We demonstrate the validity of our method first analytically and then empirically using simulated data. We propose a novel complementary filter for MAVs that fuses together gyroscope data with accelerometer and magnetic field readings. The correction part of the filter is based on the method described above and works for both IMU (Inertial Measurement Unit) and MARG (Magnetic, Angular Rate, and Gravity) sensors. We evaluate the effectiveness of the filter and show that it significantly outperforms other common methods, using publicly available datasets with ground-truth data recorded during a real flight experiment of a micro quadrotor helicopter.

Keywords:

orientation estimation; inertial measurement unit; magnetic angular rate and gravity; quaternions, micro aerial vehicles

1. Introduction

The accurate estimation of the orientation of a rigid body, relative to an inertial frame, is required for a wide range of applications. For the purpose of navigation, such estimation has employed high precision inertial and magnetic sensors, but the recent development of low-cost and light-weight MicroElectroMechanical Systems (MEMS) has allowed smaller and cheaper inertial sensors to be adopted for a wider range of applications and even in the daily use of consumer electronics such as game consoles and mobile devices. In robotics, the evolution of Micro Aerial Vehicles (MAVs) has increased in the last decade, and research is moving toward their full autonomization. Navigation and control of MAVs are possible thanks to the MEMS-based Inertial Measurement Units (IMU) that meet the MAVs’ limited size and payload requirement.

Data provided by low-cost IMU is affected by high noise levels and time-varying biases. Therefore, a sensor-fusion algorithm must be used to process the data to obtain a smooth and bias-free estimation of the orientation maintaining a low computational cost for running on the onboard processor. The orientation can be generally represented in three principal forms: Euler angles, quaternion, and Direction Cosine Matrix (DCM). The orientation in Euler form is expressed using three angles; it is conceptually easy to understand, but may reach a singularity state commonly referred as “gimbal lock”. DCM and quaternion do not incur a singularity state but the DCM represents the orientation by a 3 × 3 matrix. Furthermore, the quaternion representation offers a linear formulation of the orientation dynamics.

In this article, we propose a deterministic approach to solve Whaba’s problem [1] given gravity and magnetic field observations provided by the MARG sensor. This method returns an estimation of the attitude in quaternion form without leading to ambiguous results and avoiding singularity problems. Moreover, it does not need a predefined knowledge of the magnetic field direction. We also propose a new approach to a complementary filter for fusing together gyroscope data with acceleration and magnetic field measurements from IMU and MARG sensors to obtain an orientation estimation in quaternion form. This algorithm is meant to be used onboard MAVs because of its low computational burden, robustness, and fast convergence. However, we believe our contribution is applicable beyond the field of MAV and aerial robotics.

This article is organized as follows. Section 2 explores the literature of attitude estimation methods previously proposed as solution to Whaba’s problem. Further, in the same section, we briefly survey several fusion algorithms that exploit the measurements given by gyroscope, accelerometer, and magnetometer of an IMU/MARG sensor, to obtain a more accurate evaluation of the orientation. Section 3 briefly provides an overview of some properties of unit quaternions and how they can be used to define rotations in the 3D space. Section 4 analyzes the novel approach taken to obtain the quaternion orientation of a rigid body from earth-field observations as a closed-form solution. Section 5 explains the novel quaternion-based complementary filter. Section 6 is dedicated to the experiments and results interpretations, and Section 7 concludes the article by summarizing the findings.

2. Previous Work

The problem of finding the optimal attitude estimation, given multiple unit vector correspondences across two frames, was formulated by Whaba in 1965 [1]. Whaba’s problem was meant to estimate the attitude of a satellite by finding the optimal orthogonal matrix (the attitude matrix) that minimizes a least-square loss function using n numbers of unit vectors pairs where

n \geq 2

. Many algorithmic solutions of Whaba’s problem have been proposed. They are generally classified into two categories: deterministic and optimal, according to the popular definition by Wertz [2]. Deterministic algorithms use a minimal set of data and derive the attitude by solving nonlinear equations, whereas optimal algorithms use more than a minimal set of measurements and compute the attitude by minimizing an appropriate cost function.

A very well known deterministic algorithm is the three axis attitude determination (TRIAD) [3]. It constructs two triads of orthonormal unit vectors by combining the normalized measurement of two nonparallel reference vectors and provides an estimation of the attitude matrix. QUaternion ESTimator (QUEST) [3] is a famous optimal algorithm that produces the attitude estimation in quaternion form given a set of 3D reference unit vectors in a fixed frame and their corresponding observations in the local frame. It is based on Davenport’s q-method [4] that finds the optimal quaternion by minimizing a quadratic gain, which is a transformation of Whaba’s loss function through the parametrization of the orientation matrix by the unit quaternion. Many other techniques have been proposed, such as Singular Value decomposition (SVD) [5], Polar Decomposition (PD) [6], Euler-n [7], Fast Optimal Matrix Algorithm (FOAM) [8], and Energy Approach Algorithm (EAA) [9]. All of them produce an optimal attitude estimation and they differ from each other in their computational speed. A complete survey and analysis of attitude estimation methods in quaternion form using vector observation is provided by Markley and Mortari in [10].

It is well known that in applications based on inertial/magnetic sensors, to estimate the attitude, the two local observations of the earth’s fields (gravity and local magnetic field) are compared to the reference vectors that are supposed to be fixed. This assumption can lead to problems when the local magnetic field is perturbed by ferromagnetic objects or electrical appliances. In these cases, when algorithms such as QUEST are employed, the attitude estimation would be subject to errors not only in the yaw component but also in the pitch and roll. To avoid this problem, Yun et al. [11] proposed the Factored Quaternion Algorithm (FQA), which computes the orientation of a rigid body based on earth gravity and magnetic field measurements. The quaternion orientation is estimated by analyzing a series of three sequential rotations. This approach allows to reduce the orientation error, caused by the presence of local magnetic disturbances, only into the error in the yaw component maintaining the accuracy of the QUEST algorithm. To obtain the quaternion orientation from gravity and magnetic field observations, the method presented in this article solves the problem described above by separating the quaternion in only two parts, one for the roll and pitch components and one for the yaw. Both quaternions are found as an algebraic solution of a system instead of the result of an optimization problem. As inertial sensors provide at most the observation of 2 pairs of vectors (given a priori the global-frame gravity and magnetic field vectors), the use of optimal algorithms does not produce any improvement but only reduces the convergence speed of the orientation.

In the class of deterministic methods, Euler-2 by Daniele Mortari [7] is a mathematical approach to compute the Euler axis whose application is restricted when only two unit-vector pairs are available. Overall, besides the TRIAD algorithm, deterministic attitude-determination algorithms are not very frequent and their study is not extensive. Nonetheless, in IMU applications, only two sets of vector observations are provided, thus even the optimal approaches would have the same level of accuracy as a deterministic method.

To obtain a better estimation of the orientation, acceleration and magnetic field data are fused together with angular rate readings from a gyroscope. Although many approaches have been adopted for filtering gyroscope data with inertial measurements, the most commonly used techniques are Extended Kalman filtering (EKF) and complementary filters. A survey of other nonlinear attitude estimation methods can be found in [12].

Kalman filtering based techniques adopt a probabilistic determination of the state modeled as a Gaussian distribution given the system’s model. They are widely used in aerospace applications [13,14], human motion analysis [15,16,17] and robotics [18,19].

A Complementary filter is a common alternative to the EKF because of its simplicity and effectiveness. It uses an analysis in the frequency domain to filter the signals and combine them together to obtain an estimation of the orientation without any statistical description. In UAV applications, the use of complementary filters [20,21,22] is often preferred to the EKF because EKFs can be complicated to implement and the convergence is slower because of the time required for the linear regression iterations.

Most of the recent sensor fusion algorithms for inertial/magnetic sensors provide orientation estimation in quaternion form. Quaternions are a useful mathematical tool that require less computation time because of their minimal number of parameters and do not result in singularity configurations as the Euler representation does. Further, rotations of vectors are simply performed by quaternion multiplications.

Marins et al. [23] propose two different Kalman filter approaches to estimate the orientation in quaternion form from a MARG sensor. Both methods have the same 7-state (angular velocities and quaternion) process model but different measurement models. The first measurement model uses each MARG output in a 9-element measurement vector resulting in a complicated EKF. The second approach uses an external Gauss-Newton algorithm to directly estimate the quaternion measurement that will be part, along with the angular velocities, of the measurement vector. In this case, the relation between the process and measurement model is linear allowing the use of a simpler Kalman Filter. The quaternion-based EKF presented by Sabatini [16] is similar to the first approach of Marins et al. [23], but the state is composed by the unit quaternion augmented with accelerometer and magnetometer bias vectors for on-line calibration. In [24] Sabatini presents a similar EKF where the state is augmented with magnetic distortions vector, modeled as a Gauss-Markov process, to reduce the heading drift in magnetically non-homogeneous environments. Choukroun et al. [14] present a novel linear pseudo-measurement model that combined with the linear measurement model, yields a linear Kalman Filter algorithm that eliminates the linearization procedure of an EKF. Bachman et al. [25] present an efficient quaternion-based complementary filter for human-body-motion tracking.

Euston et al. [21] present a nonlinear quaternion-based complementary filter to estimate the attitude of a UAV given measurement from a low-cost IMU. The filter is augmented by a first-order model of the vehicle dynamics to compensate for external centripetal acceleration. In the work by Madgwick et al. [26], a constant gain filter is adopted to estimate the attitude in quaternion form of a rigid body by using data from a MARG sensor. A first quaternion estimation is obtained by gyroscope output integration and then it is corrected by a quaternion from the accelerometer and magnetometer data computed through a gradient descent algorithm. Madgwick’s method ensures good attitude estimation at low computational cost. Further, it addresses the problem of the local magnetic disturbances that, when present, affect all the orientation components (roll, pitch, and yaw). By reducing the constraint of the magnetic field vector rotation, it is able to limit the effect of the magnetic disturbances to only affect the yaw component of the orientation. The last two constant gain filters, by Euston et al. [21], and by Madgwick et al. [26], are commonly used because they offer good performances at low computational cost and a comparative analysis of the two algorithms is presented in [27]. The adaptive-gain complementary filter proposed by Tian et al. [28] fuses a quaternion estimation from fast moving gyroscope signal with a quaternion, from slow moving accelerometer and magnetometer signals, computed through a Gauss-Newton algorithm. For more robust results, the gain is adaptively adjusted according to the convergence rate of the low-frequency estimation and the divergence rate of the high-frequency estimation. Fourati et al. [29] combine the output of a Levenberg-Marquardt algorithm, the inputs for which are acceleration and local magnetic field measurements, with the angular rate measurements in a complementary observer based on the multiplicative correction technique.

To improve the performance of rigid body orientation estimation from low cost IMU/MARG under dynamic motion, different approaches can be adopted. A simple switching method, as the one proposed in the complementary filter presented in [30], can be used to determine the gain value under varying dynamics sensed by the accelerometer. Another common method includes an adaptive measurement noise covariance matrix that, in a Kalman filter framework, can be tuned adaptively to yield optimal performance during the dynamic periods as proved in [31,32]. Alternatively, the acceleration can be modeled to directly estimate the external acceleration and thus used to reduce the attitude error, as in the Kalman filter proposed in [33].

The quaternion-based complementary filter proposed in this article can be used for both IMU and MARG sensors. It ensures fast convergence and robustness thanks to the analytical derivation of the correction quaternion. We address the problem of the magnetic disturbances by separating the quaternion corrections in two different correction sequences and adopting the method presented in [26] that also makes the algorithm independent on external parameters. Further, we adopt an adaptive-gain approach to reduce the estimation error during high dynamic motion.

3. Background Theory

Any arbitrary orientation in the 3D space of a frame A with respect to a frame B can be represented by a unit quaternion

_{A}^{B} q

defined as following:

_{A}^{B} q = {[\begin{matrix} q_{0} & q_{1} & q_{2} & q_{3} \end{matrix}]}^{T} = {[\begin{matrix} cos \frac{α}{2} & e_{x} sin \frac{α}{2} & e_{y} sin \frac{α}{2} & e_{z} sin \frac{α}{2} \end{matrix}]}^{T}

(1)

where α is the rotation angle and e is the unit vector that represents the rotation axis. The quaternion conjugate of

_{A}^{B} q

, given its unit norm, is equivalent to the inverse quaternion and describes the inverse rotation. Therefore, the conjugate quaternion can be used to represent the orientation of frame B relative to frame A, as defined below.

{_{A}^{B} q}^{*} =_{B}^{A} q = {[\begin{matrix} q_{0} & - q_{1} & - q_{2} & - q_{3} \end{matrix}]}^{T}

(2)

The orientation quaternion after a sequence of rotations can be easily found by quaternion multiplication where each quaternion represents the orientation of a frame with respect to the rotated one. For example, given three frames A, B and C, and given the quaternion

_{A}^{B} q

orientation of frame A expressed with respect to frame B and given the quaternion

_{B}^{C} q

orientation of frame B expressed with respect to frame C, the orientation of frame A with respect to frame C is characterized by:

_{A}^{C} q =_{B}^{C} q \otimes_{A}^{B} q

(3)

where quaternion multiplication, given two quaternions p and q, is defined as

p \otimes q = [\begin{matrix} p_{0} q_{0} - p_{1} q_{1} - p_{2} q_{2} - p_{3} q_{3} \\ p_{0} q_{1} + p_{1} q_{0} + p_{2} q_{3} - p_{3} q_{2} \\ p_{0} q_{2} - p_{1} q_{3} + p_{2} q_{0} + p_{3} q_{1} \\ p_{0} q_{3} + p_{1} q_{2} - p_{2} q_{1} + p_{3} q_{0} \end{matrix}]

(4)

Unit quaternions can be applied to operate rotations of 3D vectors. For example a vector ^Av, expressed with respect to the A frame, can be expressed with respect to the B frame by the following operation:

{_{}^{B} v}_{q} =_{A}^{B} q \otimes {_{}^{A} v}_{q} \otimes {_{A}^{B} q}^{*}

(5)

where the symbol ⊗ indicates the quaternion multiplication, and ^Av_q and ^Bv_q are the observations of the vector v, in the two reference frames, written as pure quaternions as shown in Equation (6).

v_{q} = {[\begin{matrix} 0 & v \end{matrix}]}^{T} = {[\begin{matrix} 0 & v_{x} & v_{y} & v_{z} \end{matrix}]}^{T}

(6)

The inverse rotation that describes the vector ^Bv relative to the frame A can be easily found by using the property of the conjugate quaternion and it is presented in Equation (7).

{_{}^{A} v}_{q} = {_{A}^{B} q}^{*} \otimes {_{}^{B} v}_{q} \otimes_{A}^{B} q =_{B}^{A} q \otimes {_{}^{B} v}_{q} \otimes {_{B}^{A} q}^{*}

(7)

The rotation defined in Equation (5) can be written in matrix form as in Equation (8)

_{}^{B} v = R (_{A}^{B} q)_{}^{A} v

(8)

where

R (_{A}^{B} q)

, which belongs to the 3D special orthogonal group

S O (3)

, is the direct cosine matrix (DCM) given in terms of the orientation quaternion

_{A}^{B} q

as shown below.

R (_{A}^{B} q) = [\begin{matrix} q_{0}^{2} + q_{1}^{2} - q_{2}^{2} - q_{3}^{2} & 2 (q_{1} q_{2} - q_{0} q_{3}) & 2 (q_{1} q_{3} + q_{0} q_{2}) \\ 2 (q_{1} q_{2} + q_{0} q_{3}) & q_{0}^{2} - q_{1}^{2} + q_{2}^{2} - q_{3}^{2} & 2 (q_{2} q_{3} - q_{0} q_{1}) \\ 2 (q_{1} q_{3} - q_{0} q_{2}) & 2 (q_{2} q_{3} + q_{0} q_{1}) & q_{0}^{2} - q_{1}^{2} - q_{2}^{2} + q_{3}^{2} \end{matrix}]

(9)

Given the properties of the elements of the

S O (3)

, the inverse rotation can be defined as:

_{}^{A} v = R (_{B}^{A} q)_{}^{B} v = R^{T} (_{A}^{B} q)_{}^{B} v

(10)

4. Quaternion from Earth-Field Observations

In this section, we analyze the algebraic derivation of a quaternion from the observation of the earth’s fields. For a clear understanding of the following derivation, let us first define a notation that will be used throughout this article. We refer to the local (sensor) frame as L, and the global (earth) frame as G. We can define the measured acceleration

_{}^{L} A

and the true earth gravitational acceleration

_{}^{G} g

as the unit vectors:

\begin{matrix} _{}^{L} a & = {[\begin{matrix} A_{x} & A_{y} & A_{z} \end{matrix}]}^{T}, & ∥ a ∥ = 1 \\ _{}^{G} g & = {[\begin{matrix} 0 & 0 & 1 \end{matrix}]}^{T} \end{matrix}

Similarly, we define the measured local magnetic field

_{}^{L} m

and the true magnetic field

_{}^{G} h

as the unit vectors:

\begin{matrix} _{}^{L} m & = {[\begin{matrix} m_{x} & m_{y} & m_{z} \end{matrix}]}^{T}, & ∥ m ∥ = 1 \\ _{}^{G} h & = {[\begin{matrix} h_{x} & h_{y} & h_{z} \end{matrix}]}^{T}, & ∥ h ∥ = 1 \end{matrix}

Finally, the gyroscopes measure the angular velocity

_{}^{L} ω

around the three sensor frame axes:

_{}^{L} ω = {[\begin{matrix} ω_{x} & ω_{y} & ω_{z} \end{matrix}]}^{T}

Note that most IMUs usually measure a non-normalized vector a and m. However, for the purposes of derivation in this article, we assume that the quantities have been normalized. The only relevant units are those of ω, which we assume are radians per second.

We present an algebraic derivation of the orientation quaternion

_{G}^{L} q

, of the global frame

(G)

relative to the local frame

(L)

, as a function of

_{}^{L} A

and

_{}^{L} m

. We have two independent sensors observing two independent fields; a straightforward way to formulate the quaternion is through the inverse orientation which rotates the measured quantities

_{}^{L} A

and

_{}^{L} m

into the reference quantities

_{}^{G} G

and

_{}^{G} h

\{\begin{matrix} R^{T} (_{G}^{L} q)_{}^{L} a =_{}^{G} g \\ R^{T} (_{G}^{L} q)_{}^{L} m =_{}^{G} h \end{matrix}

(11)

This system, however, is overdetermined - each of the two equations provides two independent constraints on the orientation

_{G}^{L} q

, whereas the orientation only has three degrees of freedom. In the case when there is a disagreement between the gravitational and magnetometer readings, the system will not have a solution. The disagreement could arise from random sensor noise or unmodeled field disturbances (nongravitational accelerations or magnetic field variations). A possible solution would be to define an error metric and find the quaternion which minimizes this error. However, this could still result in disturbances in the magnetic field affecting the roll and pitch, which we are trying to avoid.

To address this problem, we present a modified system of equations. The definition of the system (but not its solution) is based on the approach presented in [26]. First, we redefine the global coordinate frame G to be aligned with the magnetic north. Specifically, the global frame’s x-axis points in the same direction as the local magnetic field (the z-axis remains vertical). Obviously, this global frame is only “fixed” in the case when the local magnetic field does not change its heading.

Next, we modify the system in Equation (11) so that the second equation provides only one constraint. Let

_{}^{G} Π_{z x^{+}}

be the half-plane which contains all points that lie in the global

x z

-plane such that x is non-negative. We require that the magnetic reading, when rotated into the global frame, lies on the half-plane

_{}^{G} Π_{z x^{+}}

. Thus, we guarantee that the heading will be measured with respect to magnetic north, but do not enforce a constraint on the magnetic inclination.

\{\begin{matrix} R^{T} (_{G}^{L} q)_{}^{L} a =_{}^{G} g \\ R^{T} (_{G}^{L} q)_{}^{L} m \in_{}^{G} Π_{z x^{+}} \end{matrix}

(12)

Note that when defined in this manner, the system no longer needs a priori knowledge of the direction of the earth’s magnetic field

_{}^{G} h

In the remainder of the section, we present a novel algebraic solution to obtain

_{G}^{L} q

as a function of

_{}^{L} a

and

_{}^{L} m

. We begin by decomposing

_{G}^{L} q

into two auxiliary quaternions,

q_{ACC}

and

q_{m AG}

, such that:

_{G}^{L} q = q_{ACC} \otimes q_{m AG}

(13)

and:

R (_{G}^{L} q) = R (q_{ACC}) R (q_{m AG})

(14)

We further define

q_{m AG}

to have only a single degree of freedom, by setting it to:

q_{m AG} = {[\begin{matrix} {q_{0}}_{m AG} & 0 & 0 & {q_{3}}_{m AG} \end{matrix}]}^{T}

(15)

It follows from the quaternion definition in Equation (1) that

q_{m AG}

represents a rotation around the z-axis only. Informally,

q_{ACC}

rotates a vector from the sensor frame to the horizontal plane of the global frame, and

q_{m AG}

rotates it around the z axis to point North. In the following subsections, we present an algebraic derivation of

q_{ACC}

and

q_{m AG}

4.1. Quaternion from Accelerometer Readings

In this subsection, we present a derivation for the auxiliary quaternion

q_{ACC}

as a function of

_{}^{L} a

. The observations of the gravity vector in the two reference frames allows us to find the quaternion that performs the transformation between the two representations. The rotation in the first equation of system Equation (12) can be re-written as:

R (_{G}^{L} q)_{}^{G} g =_{}^{L} a

(16)

and decomposed by using Equation(14) obtaining:

R (q_{ACC}) R (q_{m AG}) [\begin{matrix} 0 \\ 0 \\ 1 \end{matrix}] = [\begin{matrix} A_{x} \\ A_{y} \\ A_{z} \end{matrix}]

(17)

The representation of the gravity vector in the global frame only has a component on the z-axis; therefore any rotation about this axis does not produce any change on it. Consequently, Equation (17) is equivalent to:

R (q_{ACC}) [\begin{matrix} 0 \\ 0 \\ 1 \end{matrix}] = [\begin{matrix} A_{x} \\ A_{y} \\ A_{z} \end{matrix}]

(18)

expanding the multiplication we obtain the following system:

\{\begin{matrix} 2 (q_{1_{ACC}} q_{3_{ACC}} + q_{0_{ACC}} q_{2_{ACC}}) = A_{x} \\ 2 (q_{2_{ACC}} q_{3_{ACC}} - q_{0_{ACC}} q_{1_{ACC}}) = A_{y} \\ q_{0_{ACC}}^{2} - q_{1_{ACC}}^{2} - q_{2_{ACC}}^{2} + q_{3_{ACC}}^{2} = A_{z} \end{matrix}

(19)

It is clear that the above system is underdetermined and has an infinite number of solutions. This is not an unexpected result because the alignment of the gravity vector from its representation in the global frame into the local frame does not give any information about the rotation around the z-axis (yaw). Thus, such alignment can be achieved by infinite rotations with definite roll and pitch angles and arbitrary yaw. To restrict the solutions to a finite number we choose

q_{3_{ACC}} = 0

simplifying the system in Equation (19) to:

\{\begin{matrix} 2 q_{0_{ACC}} q_{2_{ACC}} = A_{x} \\ - 2 q_{0_{ACC}} q_{1_{ACC}} = A_{y} \\ q_{0_{ACC}}^{2} - q_{1_{ACC}}^{2} - q_{2_{ACC}}^{2} = A_{z} \end{matrix}

(20)

The above system is fully determined; solving it results in four solutions for

q_{ACC}

. Two can be discarded since they have a negative norm. The remaining two are equivalent, with all the quaternion elements switching signs between one solution and the other. For convenience, we choose the solution with positive quaternion scalar (

q_{0}

), which corresponds to the shortest-path quaternion formulation [34]. Thus, we get:

{q_{ACC} = [\begin{matrix} \sqrt{\frac{A_{z} + 1}{2}} & - \frac{A_{y}}{\sqrt{2 (A_{z} + 1)}} & \frac{A_{x}}{\sqrt{2 (A_{z} + 1)}} & 0 \end{matrix}]}^{T}

(21)

The formulation in Equation (21) is valid for all values of

A_{z}

except

A_{z} = - 1

in which it has a singularity. Further, it may arise numerical instability when in proximity of the singularity point. To address this issue we provide an alternative solution to the system in Equation (19). By simply setting

q_{2_{ACC}} = 0

instead of

q_{3_{ACC}} = 0

in Equation (19) we obtain the reduced system:

\{\begin{matrix} 2 q_{1_{ACC}} q_{3_{ACC}} = A_{x} \\ 2 q_{0_{ACC}} q_{1_{ACC}} = A_{y} \\ q_{0_{ACC}}^{2} - q_{1_{ACC}}^{2} - q_{3_{ACC}}^{2} = A_{z} \end{matrix}

(22)

which admits the following two real solutions:

{q_{ACC 1} = [\begin{matrix} - \frac{A_{y}}{\sqrt{2 (1 - A_{z})}} & \sqrt{\frac{1 - A_{z}}{2}} & 0 & \frac{A_{x}}{\sqrt{2 (1 - A_{z})}} \end{matrix}]}^{T}

(23)

and

{q_{ACC 2} = [\begin{matrix} \frac{A_{y}}{\sqrt{2 (1 - A_{z})}} & - \sqrt{\frac{1 - A_{z}}{2}} & 0 & - \frac{A_{x}}{\sqrt{2 (1 - A_{z})}} \end{matrix}]}^{T}

(24)

and we choose the solution in Equation (23). This formulation for

q_{ACC}

has a singularity at

A_{z} = 1

. Therefore, the final formulation of

q_{ACC}

that avoids the singularity problem can be obtained by combining Equations (21) and (23):

q_{ACC} = \{\begin{matrix} \begin{matrix} {[\begin{matrix} \sqrt{\frac{A_{z} + 1}{2}} & - \frac{A_{y}}{\sqrt{2 (A_{z} + 1)}} & \frac{A_{x}}{\sqrt{2 (A_{z} + 1)}} & 0 \end{matrix}]}^{T}, & A_{z} \geq 0 \\ {[\begin{matrix} - \frac{A_{y}}{\sqrt{2 (1 - A_{z})}} & \sqrt{\frac{1 - A_{z}}{2}} & 0 & \frac{A_{x}}{\sqrt{2 (1 - A_{z})}} \end{matrix}]}^{T}, & A_{z} < 0 \end{matrix} \end{matrix}

(25)

Effectively, we solve the singularity problem by having two separate formulations for

q_{ACC}

depending on the hemisphere in which a is pointing. Note that, defined in this manner,

q_{ACC}

is not continuous at the

A_{z} = 0

point. However, we will demonstrate that this problem is resolved with the formulation of

q_{m AG}

in the following subsection.

4.2. Quaternion from Magnetometer Readings

In this subsection, we present a derivation for the auxiliary quaternion

q_{m AG}

as a function of

_{}^{L} m

and

q_{ACC}

. First we use the quaternion

q_{ACC}

to rotate the body frame magnetic field vector

_{}^{L} m

into an intermediate frame whose z-axis is the same as the global coordinate frame with orthogonal x-, y- axes pointing in unknown directions due to the unknown yaw of

q_{ACC}

R^{T} (q_{ACC})_{}^{L} m = l

(26)

where l is the rotated magnetic field vector. Next, we find the quaternion (

q_{m AG}

) that rotates the vector l into the vector that lies on the

_{}^{G} Π_{z x^{+}}

of Equation (12) using the following system:

R^{T} (q_{m AG}) [\begin{matrix} L_{x} \\ L_{y} \\ L_{z} \end{matrix}] = [\begin{matrix} \sqrt{Γ} \\ 0 \\ L_{z} \end{matrix}]

(27)

where

Γ = L_{x}^{2} + L_{y}^{2}

(28)

This quaternion performs a rotation only about the global z-axis by aligning the x-axis of the intermediate frame into the positive direction of the magnetic north pole, which coincides with the x direction of our global frame. This rotation will only change the heading component of the orientation without affecting the roll and pitch components. Therefore, in presence of magnetic disturbances, their influence is only limited on affecting the heading angle. The quaternion

q_{m AG}

has the following form:

q_{m AG} = {[\begin{matrix} q_{0_{m AG}} & 0 & 0 & q_{3_{m AG}} \end{matrix}]}^{T}

(29)

By reordering the system in Equation (27) and substituting

q_{m AG}

with its components we find the following simplified system:

\{\begin{matrix} (q_{0_{m AG}}^{2} - q_{0_{m AG}}^{2}) \sqrt{Γ} = L_{x} \\ 2 q_{0_{m AG}} q_{3_{m AG}} \sqrt{Γ} = L_{y} \\ (q_{0_{m AG}}^{2} + q_{3_{m AG}}^{2}) L_{z} = L_{z} \end{matrix}

(30)

the solution of the above system which ensures the shortest rotation is the following:

{q_{m AG} = [\begin{matrix} \frac{\sqrt{Γ + L_{x} \sqrt{Γ}}}{\sqrt{2 Γ}} & 0 & 0 & \frac{L_{y}}{\sqrt{2} \sqrt{Γ + L_{x} \sqrt{Γ}}} \end{matrix}]}^{T}

(31)

It is clear from the formulation of

q_{m AG}

that the latter quaternion incurs in a singularity state for negative

L_{x}

and zero

L_{y}

. To avoid the singularity of

q_{m AG}

we prevent the l vector from having negative x- component using the following procedure. If

L_{x} < 0

we rotate l of 180° around the world z- frame, applying the quaternion

q_{π} = {[\begin{matrix} 0 & 0 & 0 & 1 \end{matrix}]}^{T}

. Finally the rotated vector is used to find

q_{m A G^{+}}

, which has the same form of Equation (31) and aligns l with the magnetic north. The sequences of rotations are summarized in the quaternion multiplications below:

q_{m A G^{+}}^{*} \otimes q_{π}^{*} \otimes q_{ACC}^{*} \otimes m_{q} \otimes q_{ACC} \otimes q_{π} \otimes q_{m A G^{+}}

(32)

where

m_{q}

is the local magnetic field vector written as pure quaternion and

q_{ACC}^{*} \otimes m_{q} \otimes q_{ACC} = l_{q}

. For the sake of simplicity we consider the quaternion product between

q_{π}

and

q_{m A G^{+}}

as the alternative formulation of

q_{m AG}

in the case of

L_{x} < 0

as following.

q_{m AG} = q_{π} \otimes q_{m A G^{+}}

(33)

The result of the above multiplication is shown in Equation (34):

{q_{m AG} = [\begin{matrix} \frac{L_{y}}{\sqrt{2} \sqrt{Γ - L_{x} \sqrt{Γ}}} & 0 & 0 & \frac{\sqrt{Γ - L_{x} \sqrt{Γ}}}{\sqrt{2 Γ}} \end{matrix}]}^{T}

(34)

The complete formulation of

q_{m AG}

that avoids the singularity problem discussed above, is eventually obtained by combining Equation (31) with Equation (34):

q_{m AG} = \{\begin{matrix} \begin{matrix} {[\begin{matrix} \frac{\sqrt{Γ + L_{x} \sqrt{Γ}}}{\sqrt{2 Γ}} & 0 & 0 & \frac{L_{y}}{\sqrt{2} \sqrt{Γ + L_{x} \sqrt{Γ}}} \end{matrix}]}^{T}, & L_{x} \geq 0 \\ {[\begin{matrix} \frac{L_{y}}{\sqrt{2} \sqrt{Γ - L_{x} \sqrt{Γ}}} & 0 & 0 & \frac{\sqrt{Γ - L_{x} \sqrt{Γ}}}{\sqrt{2 Γ}} \end{matrix}]}^{T}, & L_{x} < 0 \end{matrix} \end{matrix}

(35)

Finally, we can generalize the quaternion orientation of the global frame relative to the local frame as the multiplication of the two quaternions

q_{ACC}

and

q_{m AG}

as below.

_{G}^{L} q = q_{ACC} \otimes q_{m AG}

(36)

Note that the quaternion

_{G}^{L} q

does not suffer from the discontinuity problem of the yaw angle given by the switching formulation of

q_{ACC}

of Equation (25) thanks to the multiplication with

q_{m AG}

, which performs the alignment of the intermediate frame into the global frame as previously discussed.

Figure 1. Simulated experiment showing the orientation output in quaternion form (left) and in Euler representation (right). (a) Orientation reference; (b) Orientation output of the presented method with noisy acceleration data and noise-free magnetometer data; (c) Orientation output of the presented method with noise-free acceleration data and magnetic readings affected by magnetic disturbances.

To test the effectiveness of our method, we implemented a MatLab simulation. The results of which are shown in Figure 1. We simulate three consecutive constant velocity (1 rad/s) rotations of 360° about each axis. First, we test the result of the quaternion solution of our proposed method by using perfect accelerometer and magnetometer data obtained by rotating the global frame vectors (gravity and magnetic field) into the local frame by using the orientation reference computed through the integration of the angular velocity vector, obtaining the same output of Figure 1a. Note, on the right side of Figure 1, that during the rotation about the y-axis, besides the variation of the pitch angle, we also observe an instantaneous jump of the roll and yaw angles because of the singularity configuration that affects the Euler representation. However, as shown on the left side of the figure, the quaternion representation does not incur in a singularity state; moreover, our formulation ensures the continuity of the quaternion throughout the rotation. In Figure 1b we prove that our method works with noisy acceleration data affected by Gaussian noise. When the sensor is close to some configuration where the yaw angle’s value is either π or −π, because of the noise, the yaw flips between the two values, which are two alternate representation of the same rotation. Finally, in the simulation of Figure 1c, we show the result of our orientation estimation from noise-free acceleration data and magnetometer data affected by magnetic disturbances. It is clear, from the Euler representation of the orientation, that the magnetic disturbances affect only the yaw angle.

5. Quaternion-Based Complementary Filter

A complementary filter uses an analysis in the frequency domain of the signals to combine them to obtain a better estimation of a particular quantity. If the signals are affected by noises with different frequency, two filters, with an appropriate bandwidth, can be applied such that the sum of the two filtered signals cover the full range of useful frequency. For attitude estimation from IMU readings, a complementary filter performs high-pass filtering on the orientation estimated from gyroscope data affected by low-frequency noise, and low-pass filter on accelerometer data affected by high-frequency noise. The fusion between the two filtered estimations will ideally obtain an all-pass and noise-free attitude estimation. The term “complementary” derives from the cut-off frequency value, which is the same for both filters. Thus, its correct value is found as a trade-off between the preserved bandwidth of each single signal.

The complementary filter we propose in this article can be used for both IMU and MARG sensors, and is described in the diagrams of Figure 2 and Figure 3. It fuses attitude estimation in quaternion form from gyroscope data, with accelerometer data in the form of a delta quaternion, which serves as correction only for roll and pitch components of the attitude maintaining the yaw estimation from the gyroscope. If magnetometer data is provided, a second step is added to the algorithm where a delta quaternion, from magnetic field readings, is derived to correct the heading of the previous estimation by performing a small rotation about the global z-axis in order to align the current frame with the magnetic field.

Figure 2. Block diagram of the complementary filter for IMU implementation (no magnetometer data).

Figure 3. Block diagram of the complementary filter for MARG implementation (with magnetometer data).

5.1. Prediction

In the prediction step, the angular velocity vector, measured by the tri-axis gyroscope, is used to compute a first estimation of the orientation in quaternion form. Assuming known the initial conditions, we first calculate the quaternion describing the rate of change of the orientation as a quaternion derivative, by multiplying the previous state with the angular velocity vector arranged as a pure quaternion as in Equation (6). In the literature, the quaternion derivative from an angular rate measurement is usually calculated for the forward quaternion, that is, the one representing the orientation of the local frame with respect to the global frame. Because in this article we use the inverse orientation, for the sake of clarity, we compute the quaternion derivative for our convention, starting from the formula for the forward quaternion.

It is well known that the angular rate

_{}^{L} ω

and the “forward” quaternion derivative are related by the following identity:

{_{L}^{G} \dot{q}}_{ω, t_{k}} = \frac{1}{2} {_{L}^{G} q}_{t_{k - 1}} \otimes {_{}^{L} ω}_{q, t_{k}}

(37)

where

{_{}^{L} ω}_{q, t_{k}}

is the angular velocity vector arranged as pure quaternion at time instant

t_{k}

and

{_{L}^{G} q}_{t_{k - 1}}

is the previous estimate of the orientation. In our case, we want to determine the derivative of the inverse unit quaternion, which is simply the conjugate of Equation (37), therefore:

{_{G}^{L} \dot{q}}_{ω, t_{k}} = {_{L}^{G} \dot{q}}_{ω, t_{k}}^{*} = \frac{1}{2} {({_{L}^{G} q}_{t_{k - 1}} \otimes {_{}^{L} ω}_{q, t_{k}})}^{*} = \frac{1}{2} {_{}^{L} ω}_{q, t_{k}}^{*} \otimes {_{L}^{G} q}_{t_{k - 1}}^{*}

as the conjugate of the pure quaternion

{_{}^{L} ω}_{q}

is just its negative (zero scalar component), we can finally write:

{_{G}^{L} \dot{q}}_{ω, t_{k}} = - \frac{1}{2} {_{}^{L} ω}_{q, t_{k}} \otimes {_{G}^{L} q}_{t_{k - 1}}

(38)

The above equation can be rewritten in matrix form as shown in Equation (39):

{_{G}^{L} \dot{q}}_{ω, t_{k}} = Ω ({_{}^{L} ω}_{t_{k}}) {_{G}^{L} q}_{t_{k - 1}}

(39)

where

Ω ({_{}^{L} ω}_{t_{k}}) = [\begin{matrix} 0 & {_{}^{L} ω}_{t_{k}}^{T} \\ - {_{}^{L} ω}_{t_{k}} & - [{_{}^{L} ω}_{t_{k}} \times] \end{matrix}]

(40)

and, regardless of the time dependence, the term

[_{}^{L} ω \times]

denotes the cross-product matrix that is associated with

_{}^{L} ω

and is equal to:

[_{}^{L} ω \times] = [\begin{matrix} 0 & - ω_{z} & ω_{y} \\ ω_{z} & 0 & - ω_{x} \\ - ω_{y} & ω_{x} & 0 \end{matrix}]

(41)

The orientation of the global frame relative to local frame at time

t_{k}

can be finally computed by numerically integrating the quaternion derivative using the sampling period

Δ t = t_{k} - t_{k - 1}

{_{G}^{L} q}_{ω, t_{k}} = {_{G}^{L} q}_{t_{k - 1}} + {_{G}^{L} \dot{q}}_{ω, t_{k}} Δ t

(42)

5.2. Correction

The correction adopted is based on a multiplicative technique. The predicted quaternion

{_{G}^{L} q}_{ω}

(for clarity we omit the time t) is corrected by means of two delta quaternions as in Equation (43):

_{G}^{L} q = {_{G}^{L} q}_{ω} \otimes {\hat{Δ q}}_{ACC} \otimes {\hat{Δ q}}_{m AG}

(43)

Each delta quaternion is computed and filtered by the high-frequency noise independently; thus we divide the correction into two separated steps. The first step corrects the predicted quaternion only in the roll and pitch components by the application of

{\hat{Δ q}}_{ACC}

computed with data from the accelerometer. The second step is used only when the magnetic field readings are provided and corrects the yaw component of the quaternion orientation by applying

{\hat{Δ q}}_{m AG}

5.2.1. Accelerometer-Based Correction

We use the inverse predicted quaternion

{_{L}^{G} q}_{ω}

to rotate the normalized body frame gravity vector

_{}^{L} a

(measured by the accelerometer) into the global frame as shown below.

R ({_{L}^{G} q}_{ω})_{}^{L} a = {_{}^{G} g}_{p}

(44)

obtaining the vector

{_{}^{G} g}_{p}

that we call “predicted gravity”. The predicted gravity will have a small deviation from the real gravity vector; therefore, we compute the delta quaternion,

Δ q_{ACC}

, which rotates

_{}^{G} g

into

{_{}^{G} g}_{p}

(because we represent the orientation of the global frame relative to the local frame), by using the following system:

R (Δ q_{ACC}) [\begin{matrix} 0 \\ 0 \\ 1 \end{matrix}] = [\begin{matrix} G_{x} \\ G_{y} \\ G_{z} \end{matrix}]

(45)

where

_{}^{G} g = [\begin{matrix} 0 \\ 0 \\ 1 \end{matrix}], {_{}^{G} g}_{p} = [\begin{matrix} G_{x} \\ G_{y} \\ G_{z} \end{matrix}]

Note that the system in Equation (45) is similar to the system in Equation (18); therefore, we proceed in the same manner to find a closed-form solution. We define the component

Δ q_{3_{ACC}} = 0

thus obtaining the following simplified system.

\{\begin{matrix} 2 Δ q_{0_{ACC}} Δ q_{2_{ACC}} = G_{x} \\ - 2 Δ q_{0_{ACC}} Δ q_{1_{ACC}} = G_{y} \\ Δ q_{0_{ACC}}^{2} - Δ q_{1_{ACC}}^{2} - Δ q_{2_{ACC}}^{2} = G_{z} \end{matrix}

(46)

Proceeding as for system in Equation (18) we find the following solution:

{Δ q_{ACC} = [\begin{matrix} \sqrt{\frac{G_{z} + 1}{2}} & - \frac{G_{y}}{\sqrt{2 (G_{z} + 1)}} & \frac{G_{x}}{\sqrt{2 (G_{z} + 1)}} & 0 \end{matrix}]}^{T}

(47)

which is the one that ensures the shortest rotation. Although Equation (47) has a singularity for

G_{z} = - 1

, we do not resolve it because the delta quaternion is applied at each step to correct the small deviation between the predicted gravity and the real gravity. Therefore, the value of

G_{z}

will always be very close to 1. However, for a more rigorous formulation, the approach used in Equation (22) can be adopted.

As the delta quaternion is affected by the accelerometer’s high frequency noise, before applying it to the predicted quaternion, we scale it down by using an interpolation with the identity quaternion

q_{I}

. We use two different interpolation approaches based on the angle between

q_{I}

and

Δ q_{ACC}

. Given two quaternions p and q, the cosine of the angle Ω subtended by the arc between them is equal to the dot product of the two quaternions as shown below:

cos Ω = p \cdot q = p_{0} q_{0} + p_{1} q_{1} + p_{2} q_{2} + p_{3} q_{3}

(48)

In our case, as the identity quaternion has the following components

{q_{I} = [\begin{matrix} 1 & 0 & 0 & 0 \end{matrix}]}^{T}

(49)

the dot product is equal to the

Δ q_{0_{ACC}}

component of

Δ q_{ACC}

. If

Δ q_{0_{ACC}} > ϵ

, where ϵ is a threshold value (in our case

ϵ = 0.9

), we use a simple Linear intERPolation (LERP) as in the equation below.

{\bar{Δ q}}_{ACC} = (1 - α) q_{I} + α Δ q_{ACC}

(50)

where

α \in [0, 1]

is the gain that characterizes the cut-off frequency of the filter [35]. LERP does not maintain the unit norm of the delta quaternion

{\bar{Δ q}}_{ACC}

that is then normalized immediately after:

{\hat{Δ q}}_{ACC} = \frac{{\bar{Δ q}}_{ACC}}{∥ {\bar{Δ q}}_{ACC} ∥}

(51)

Δ q_{0_{ACC}} \leq ϵ

, we use the Spherical Linear intERPolation (SLERP) [36]. The SLERP algorithm gives a correct evaluation of the weighted average of two points lying on a curve. In the case of quaternions, the points lie on the surface of the 4D sphere (hypersphere). A simple 2D visualization of LERP and SLERP mapping processes is shown in Figure 4. The SLERP formula we used is the following:

{\hat{Δ q}}_{ACC} = \frac{sin ([1 - α] Ω)}{sin Ω} q_{I} + \frac{sin (α Ω)}{sin Ω} Δ q_{ACC}

(52)

The

Δ q_{ACC}

may have a large value when the predicted gravity has a significant deviation from the real gravity. If that condition does not occur, the delta quaternion is very small; thus, we prefer using the LERP formula because it is computationally more efficient. Finally, the quaternion estimated from gyroscope data is multiplied by the filtered delta quaternion, providing correction in the roll and pitch component, while the heading is preserved from the predicted orientation.

_{G}^{L} q^{'} = {_{G}^{L} q}_{ω} \otimes {\hat{Δ q}}_{ACC}

(53)

Figure 4. 2-D representation of two different interpolation methods between quaternion q_a and q_b resulting in quaternion q_int. On the left is represented the linear interpolation (LERP), whereas on the right the spherical linear interpolation (SLERP) showing the correct result.

5.2.2. Magnetometer-Based Correction

If magnetic field data is provided, a second step performs the correction on the heading component. We proceed in the same way as in the first step by computing a delta quaternion. First, we use the quaternion inverse of

_{G}^{L} q^{'}

, calculated in Equation (53), to rotate the body frame magnetic field vector

_{}^{L} m

into the world frame magnetic vector.

R (_{L}^{G} q^{'})_{}^{L} m = l

(54)

Next, we find the delta quaternion (

Δ q_{m AG}

), which rotates the vector l into the vector that lies on the

x z

-semiplane defined in Section 4 such that:

R^{T} (Δ q_{m AG}) [\begin{matrix} L_{x} \\ L_{y} \\ L_{z} \end{matrix}] = [\begin{matrix} \sqrt{L_{x}^{2} + L_{y}^{2}} \\ 0 \\ L_{z} \end{matrix}]

(55)

This delta quaternion performs a rotation only about the global z-axis by aligning the global x-axis into the positive direction of the magnetic north pole. Further, with this formulation, such rotation does not affect the roll and pitch components even in the presence of magnetic disturbances, limiting their influence only on the heading angle. Thus, the delta quaternion

Δ q_{m AG}

has the following shape:

Δ q_{m AG} = {[\begin{matrix} Δ q_{0_{m AG}} & 0 & 0 & Δ q_{3_{m AG}} \end{matrix}]}^{T}

(56)

By rearranging the system in Equation (55) and substituting

Δ q_{m AG}

with its component, we find the following simplified system:

\{\begin{matrix} (Δ q_{0_{m AG}}^{2} - Δ q_{3_{m AG}}^{2}) \sqrt{L_{x}^{2} + L_{y}^{2}} = L_{x} \\ 2 Δ q_{0_{m AG}} Δ q_{3_{m AG}} \sqrt{L_{x}^{2} + L_{y}^{2}} = L_{y} \\ (Δ q_{0_{m AG}}^{2} + Δ q_{3_{m AG}}^{2}) L_{z} = L_{z} \end{matrix}

(57)

the solution of the above system that ensures the shortest rotation is the following:

Δ q_{m AG} = {[\begin{matrix} \frac{\sqrt{Γ + L_{x} \sqrt{Γ}}}{\sqrt{2 Γ}} & 0 & 0 & \frac{L_{y}}{\sqrt{2 (Γ + L_{x} \sqrt{Γ})}} \end{matrix}]}^{T}

(58)

with Γ as per Equation (28). The delta quaternion

Δ q_{m AG}

is affected by the noise of the magnetometer, which is filtered by using the procedure we adopted for

Δ q_{ACC}

, switching between LERP and SLERP according to the same criterion. Another advantage of using two correction steps is the possibility to apply two different filtering gains. As the two delta quaternions are completely independent to each other, and each of them is related to a particular sensor (accelerometer and magnetometer), they might be affected by different frequency noise. Thus, The LERP and SLERP formulas applied to

Δ q_{m AG}

are the same as the ones in Equations (50) and (52), respectively, where the gain α is replaced with β, obtaining

{\hat{Δ q}}_{m AG}

. Eventually, the delta quaternion

{\hat{Δ q}}_{m AG}

is applied to the quaternion in Equation (53) obtaining the final quaternion, which expresses the orientation of the global frame with respect to the local frame.

_{G}^{L} q =_{G}^{L} q^{'} \otimes {\hat{Δ q}}_{m AG}

(59)

5.3. Adaptive Gain

A drawback of a typical implementation of the complementary filter is the constant gain that often causes inaccurate orientation estimation during highly dynamic motion. When the vehicle moves with high acceleration, the magnitude and direction of the total measured acceleration vector are different from gravity; therefore the attitude is evaluated using a false reference, resulting in significant, possibly critical errors. However, the gyroscope readings are not affected by linear acceleration, thus they can still be used to compute a relatively accurate orientation estimation that, under this condition, should be treated as the main source of the estimate. A constant gain fusion algorithm cannot overcome the aforementioned problem if the optimal gain has been evaluated for static conditions. In this paper we address this issue adopting an adaptive gain whose strategy is slightly different from the switching approach proposed in [30].

First we define the magnitude error

e_{m}

as in the following equation:

e_{m} = \frac{| ∥_{}^{L} \tilde{a} ∥ - G |}{G}

(60)

where

∥_{}^{L} \tilde{a} ∥

is the norm of the measured local frame acceleration vector before normalization and g = 9.81 m/s

^{2}

. Given Equations (50) and (52), we make the filtering gain α function of the magnitude error

e_{m}

through the gain factor f, meaning:

α = α (e_{m}) = \bar{α} f (e_{m})

(61)

where

\bar{α}

is the constant value that gives the best filtering result in static conditions and

f (e_{m})

is what we call the “gain factor”, which is a piecewise continuous function of the magnitude error as shown in Figure 5. The gain factor is constant and equal to 1 when the magnitude of the nongravitational acceleration is not high enough to overcome the acceleration of gravity and the value of the error magnitude does not reach the first threshold. If the nongravitational acceleration rises and the error magnitude exceeds that first threshold, the gain factor decreases linearly with the increase of the magnitude error until reaching zero for error magnitude equal to the second threshold and over. Empirically, we found that the two thresholds’ values that give the best results are respectively 0.1 and 0.2.

Figure 5. The gain factor function used in the adaptive gain to reduce the attitude error deriving from the external linear acceleration because of vehicle motion.

In this article we do not address the problem of the centripetal acceleration; however, if conditions warrant, one can easily add the centripetal force model used in [21].

5.4. Filter Initialization and Bias Estimation

The complementary filter proposed in this article is initialized with the orientation quaternion calculated using the procedure explained in Section 4. The values of the current body-frame acceleration and magnetic field vectors are used to produce the quaternion representing the initial orientation of the rigid body in any configuration. Therefore, for the initialization, the filter does not need any assumption and it is performed in one single step.

Before using the angular velocity vector for the quaternion prediction, we correct it for the bias. The bias of a sensor’s reading is a slow-varying signal considered as low frequency noise. Therefore, we adopt a low-pass filter to separate the bias from the actual angular velocity reading. To avoid filtering useful information, the low-pass filtering is applied only when the sensor is in a steady-state condition that is previously checked. If the sensor is in the steady-state condition, the bias is updated, otherwise it is assumed to be equal to the previous step value. The estimated bias is then subtracted from the gyroscope reading obtaining a bias-free angular velocity measurement. The low-pass filter cut-off frequency can be selected by the user as well as parameters for the steady-state condition, to apply the bias-estimation to different sensors.

6. Experiments

In this section, we evaluate the performances of the proposed complementary filter under different conditions and we compare it against other common estimation methods. At first we evaluate the overall performances during MAV flights, then we evaluate and compare the behavior of the different methods under, respectively, magnetic disturbances and high nongravitational acceleration.

6.1. MAV Flight Experiment

In the first experiment, we evaluated the accuracy of the orientation estimation algorithm using publicly-available MAV datasets [37] with ground-truth orientation in Euler-angles form from a motion-capture system.

The datasets are recorded using an AscTec “Pelican” quadrotor, flying in an indoor environment of size 10 m×10 m×10 m, equipped with eight Vicon cameras, performing 1, 2, and 3 loops, respectively. Figure 6 shows the trajectory traveled by the quadrotor during the “1loop” experiment, tracked by the motion capture system. The datasets include acceleration and angular velocity readings from the IMU, and no magnetic field data. We compare our method against the orientation estimation proposed by Madgwick et al. in [26], the quaternion-based EKF proposed by Sabatini [24], and the algorithm embedded in the low-level processor of the AscTec quadrotor, whose output, in Euler angles form, is provided in the datasets. Table 1, Table 2 and Table 3 show the Root Mean Square (RMS) error for roll, pitch, and yaw angles, respectively, in the three datasets, proving that our algorithm outperforms the other methods for all the sets. Note that the proposed method and one of the algorithms against which we compare it, by Madgwick et al. are constant gain filters. Therefore, their performance will vary based on the value of the chosen gain. To obtain a fair qualitative comparison, we chose the gains that minimize the RMS error of the estimation for each method. Figure 7 shows the orientation in Euler angles for the four different algorithms using the data from the “1LoopDown” dataset. Note that as the magnetic field readings are not provided, no correction is available for the yaw component. However, enabling our bias estimation allows a significant reduction of the drift error.

Figure 6. 2D top view (left) and 3D side view (right) of the trajectory traveled by the quadrotor during the “1LoopDown” experiment.

Table 1. RMS roll angle error [radians].

**Table 1.** RMS roll angle error [radians].
Dataset	Proposed	Madgwick	AscTec	EKF
1LoopDown	0.0233	0.0370	0.0464	0.0287
2LoopsDown	0.0292	0.0470	0.0338	0.0314
3LoopsDown	0.0277	0.0405	0.0315	0.0331

Table 2. RMS pitch angle error [radians].

**Table 2.** RMS pitch angle error [radians].
Dataset	Proposed	Madgwick	AscTec	EKF
1LoopDown	0.0209	0.0336	0.0369	0.0284
2LoopsDown	0.0223	0.0369	0.0313	0.0384
3LoopsDown	0.0202	0.0360	0.0329	0.0392

Table 3. RMS yaw angle error [radians].

**Table 3.** RMS yaw angle error [radians].
Dataset	Proposed	Madgwick	AscTec	EKF
1LoopDown	0.1429	0.2543	0.3388	0.1888
2LoopsDown	0.1309	0.9229	0.3182	0.3345
3LoopsDown	0.2890	1.3327	0.3255	0.3545

Figure 7. Comparison of the orientation Euler angles output of the three methods during the “1LoopDown” experiment. (a) Roll; (b) Pitch; (c) Yaw.

6.2. Magnetic Disturbances

In the second experiment, we used a Phidgets Spatial 3/3/3 sensor. It provides acceleration, angular rate, and magnetic field strength measurements in three axes. In this case, we computed the orientation using acceleration as well as magnetic field data. We held the sensor in a steady state while reading the output of our complementary filter, the EKF, and the filter by Madgwick et al. We induced magnetic disturbances by approaching a magnet to the sensor and repeated this action thrice. The first two times we applied the magnetic disturbance only for a couple of seconds, whereas for the third time, we retained the disturbance until the end of the experiment. Figure 8 shows the results of the experiment. The roll and pitch angles from our complementary filter are completely immune to magnetic disturbances. The roll and pitch from the Madgwick filter are affected by the magnetic disturbances when we approach the magnet to the sensor, converging back to the right value when the disturbance is removed. When we induce the disturbance the third time, the roll and pitch of the Madgwick filter present a similar behavior, but with a slower convergence, showing that the roll and pitch angles are only immune to constant magnetic flux. The EKF has a random behavior because it is already affected by the magnetic flux present in the room caused by electric devices, as a consequence of the coupled nature of the acceleration and magnetic field correction. In fact, although the norm of the measured magnetic field is constant at the beginning of the experiment, its value is different from the norm of the reference magnetic field vector (about 0.54 Gauss in this area), vector that is needed as input to the EKF. At the bottom of Figure 8 we show the yaw angles output of the three algorithms, during the same experiment, compared against the reference yaw angle extracted from the magnetic field vector. Unlike the Madgwick filter, the proposed complementary filter has the advantage of having a dedicated gain (β) to filter the magnetic field measurement noise. This allows us to change the yaw sensitivity to the reference variations. For this experiment, we chose the gain value that gives the best trade-off between the error because of the magnetic disturbances and the rate of the convergence. Note that although the Madgwick filter can vary its magnetic sensitivity by changing the filtering gain, this would also cause a performance degradation in the estimation of the roll and pitch components because only one gain controls the cut-off frequency of the signals coming from two different sensors. In the case of the EKF, when the measurement covariance comes directly from the sensor’s standard deviation, the yaw output is very sensitive to the magnetic disturbances as can be seen in the figure. However, it is well known that the measurement covariance can be changed to make the algorithm less sensitive to the selected measurement.

6.3. High Nongravitational Acceleration

In this experiment, we analyze the efficiency of our complementary filter adopting the adaptive gain introduced in Section 5.3 under the condition of high nongravitational acceleration. We used the previously described Phidgets Spatial 3/3/3 sensor, attached to a slider, free to translate over a camera rail track, such that the sensor’s x-axis is pointing along the direction of the translation. This setup allows us to create accelerated motion even as maintaining the IMU on the same plane throughout the experiment with constant zero attitude. We induced a high linear acceleration on the

x -

body axis by moving the slider abruptly back and forth. Note that because the sensor moves in the horizontal plane, the inertial forces measured by the IMU are completely decoupled: gravitational acceleration on the z-axis, and total body motion acceleration on the x-axis.

Figure 8. Comparison of the different estimation methods in the presence of magnetic disturbances. (a) Norm of the measured magnetic field vector; (b) Roll; (c) Pitch; (d) Yaw.

We evaluate the improvement of our method from the addition of the adaptive gain. Figure 9 shows the pitch angle output of the proposed algorithm with constant and adaptive gain, as well as the pitch angle output of the EKF and the Madgwick filter. The acceleration because of the rigid-body motion has a varying value reaching a maximum of 50 m/s

^{2}

. This causes an error of the output of the constant gain complementary filter to reach 0.3 radians, whereas the output of the adaptive gain filter has an error which never exceeds 0.02 radians. Moreover, the proposed adaptive gain complementary filter clearly outperforms the other two methods that do not have a solution to reduce the attitude error during accelerated motion.

Figure 9. Pitch estimation under the condition of high nongravitational acceleration; (a) Non gravitational acceleration applied during the experiment; (b) Pitch angle output of the different algorithms

Ultimately, we analyzed the computational time of the three methods. We profiled a reasonably optimized C++ implementation of all the algorithms running on an Intel core I7, 3.6 GHz processor. Table 4 summarizes the results for the average execution time of a prediction-correction update cycle. The results demonstrate that the proposed complementary filter and the Madgwick algorithm have very similar processing times suitable for small embedded processors, whereas the EKF has a runtime almost six times higher. Note that given the higher standard deviation of the average time measurement relative to the Madgwick filter, we can consider negligible the small difference of 0.14

μ s

between the two faster approaches.

Table 4. Computational time of estimation algorithms.

**Table 4.** Computational time of estimation algorithms.
Algorithm	Average Time ( $μ s$ )	Standard Deviation ( $μ s$ )
Proposed	1.4243	0.4761
Madgwick	1.2839	0.7101
EKF	7.0408	0.2342

7. Conclusions

In this article we presented a novel deterministic solution to the Whaba’s problem to find the orientation in quaternion form given two pairs of vector observations from inertial sensors. We are able to find the orientation of the global frame with respect to the local frame without leading to ambiguous results or singularity problems. The subdivision of the quaternion in two parts makes the roll and pitch components immune to magnetic distortions. Further, the procedure adopted to find the heading component of the orientation eliminates the need for a direction magnetic field to be predefined. We demonstrated analytically and empirically the validity of the proposed method.

We also presented a novel quaternion-based complementary filter, which fuses together gyroscope data with accelerometer and magnetic field data. The predicted quaternion is calculated by using the well-known linear formulation of the quaternion rate of change from the gyroscope’s angular rate. The correction is applied by means of two delta quaternions that are computed using an approach based on the method described in Section 4. The presented algorithm offers several advantages over other implementations:

Fast initialization in quaternion form allowing any starting configuration.
Fast convergence of the orientation quaternion because of the algebraic formulation of the delta quaternions.
Two different gains to separately filter acceleration and magnetic field noises.
Magnetic distortion compensation that involves a two-fold advantage: it avoids the impact of the magnetic disturbances on the roll and pitch components of the orientation when the sensor is surrounded by unwanted magnetic flux and eliminates the need of a predefined magnetic field direction.

Moreover, our algorithm does not compute any matrix inversion or matrix multiplication, maintaining a low computational cost making it suitable for onboard implementation.

We evaluated our complementary filter using publicly available data of a micro quadrotor helicopter flying in an indoor environment. We compared our results against other algorithms, showing better performances of our method.

Acknowledgments

This study has been supported in part by U.S. Army Research Office under grant No. W911NF-09-1-0565, U.S. National Science Foundation under grants No. IIS-0644127 and No. CBET-1160046.

Author Contributions

Roberto G. Valenti and Ivan Dryanovski conceived the mathematical derivation of the presented approach and wrote the implementation of the algorithm; Roberto G. Valenti performed the experiments, analyzed the data and wrote the paper; Jizhong Xiao contributed material/analysis tools and mentorship.

Conflicts of Interest

The authors declare no conflict of interest.

References

Wahba, G. A Least Square Estimate of Spacecraft Attitude. Soc. Ind. Appl. Math. (SIAM) Rev. 1965, 7, 409–1965. [Google Scholar]
Wertz, J.R. Spacecraft Attitude Determination and Control; Wertz, J.R., Ed.; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1978; pp. 427–428. [Google Scholar]
Shuster, M.; Oh, S. Three-axis attitude determination from vector observations. J. Guid. Control 1981, 4, 70–77. [Google Scholar] [CrossRef]
Keat, J. Analysis of Least-Squares Attitude Determination Routine DOAOP; Report CSC/TM-77/6034; Computer Sciences Corporation: Greenbelt, MD, USA, February 1977. [Google Scholar]
Markley, F.L. Attitude Determination Using Vector Observations and the Singular Value Decomposition. J. Astronaut. Sci. 1988, 36, 245–258. [Google Scholar]
Bar-Itzhack, I.Y. Polar Decomposition for Attitude Determination from Vector Observations. In Proceedings of the AIAA Guidance, Navigation, and Control, Hilton Head Island, SC, USA, 10–12 August 1992.
Mortari, D. Euler-2 and Euler-n Algorithms for Attitude Determination from Vector Observations. In Proceedings of the IFAC Space Technology International Conference on Intelligent Autonomous Control in Aerospace, Beijing, China, 14–16 August 1995.
Markley, F.L. Attitude Determination Using Vector Observations: A Fast Optimal Matrix Algorithm. J. Astronaut. Sci. 1993, 41, 261–280. [Google Scholar]
Mortari, D. Energy Approach Algorithm for Attitude Determination from Vector Observations. Adv. Austronaut. Sci. 1995, 89, 773–784. [Google Scholar]
Markley, F.L.; Mortari, D. Quaternion Attitude Estimation Using Vector Observations. J. Austronaut. Sci. 2000, 48, 359–380. [Google Scholar]
Yun, X.; Bachmann, E.R.; Mcghee, R.B. A Simplified Quaternion-Based Algorithm for Orientation Estimation From Earth Gravity and Magnetic Field Measurements. IEEE Trans. Instrum. Meas. 2008, 57, 638–650. [Google Scholar]
Crassidis, J.L.; Markley, F.L.; Cheng, Y. Survey of nonlinear attitude estimation methods. J. Guid. Control Dyn. 2007, 30, 12–28. [Google Scholar] [CrossRef]
Gebre-Egziabher, D.; Hayward, R.C.; Powell, J.D. Design of Multi-sensor Attitude Determination Systems. IEEE Trans. Aerosp. Electron. Syst. 2004, 40, 627–649. [Google Scholar] [CrossRef]
Choukroun, D.; Bar-Itzhack, I.Y.; Oshman, Y. Novel Quaternion Kalman Filter. IEEE Trans. Aerosp. Electron. Syst. 2006, 42, 174–190. [Google Scholar] [CrossRef]
Tao, W.; Zheng, R.; Feng, H. Gait Analysis Using Wearable Sensors. Sensors 2012, 12, 2255–2283. [Google Scholar] [CrossRef] [PubMed]
Sabatini, A.M. Quaternion-based extended Kalman filter for determining orientation by inertial and magnetic sensing. IEEE Trans. Biomed. Eng. 2006, 53, 1346–1356. [Google Scholar] [CrossRef] [PubMed]
Hung, T.N.; Suh, Y.S. Inertial Sensor-Based Two Feet Motion Tracking for Gait Analysis. Sensors 2013, 13, 5614–5629. [Google Scholar] [CrossRef] [PubMed]
Barshan, B.; Durrant-Whyte, H.F. Inertial navigation systems for mobile robots. IEEE Trans. Robot. Autom. 1995, 11, 328–342. [Google Scholar] [CrossRef] [Green Version]
Jun, M.; Roumeliotis, S.I.; Sukhatme, G.S. State estimation of an autonomous helicopter using Kalman filtering. In Proceedings of the IEEE/RSJ International Conference on Intelligent Robots and Systems, Kyongju, Korea, 17–21 October 1999; Volume 3, pp. 1346–1353.
Baerveldt, A.J.; Klang, R. A low-cost and low-weight attitude estimation system for an autonomous helicopter. In Proceedings of the IEEE International Conference on Intelligent Engineering Systems, Budapest, Hungary, 15–17 September 1997; pp. 391–395.
Euston, M.; Coote, P.; Mahony, R.; Hamel, T. A complementary filter for attitude estimation of a fixed-wing UAV. In Proceedings of the IEEE/RSJ International Conference on Intelligent Robots and Systems, Nice, France, 22–26 September 2008; pp. 340–345.
Mahony, R.; Hamel, T.; Pflimlin, J. Nonlinear Complementary Filters on the Special Orthogonal Group. IEEE Trans. Autom. Control 2008, 53, 1203–1218. [Google Scholar] [CrossRef] [Green Version]
Marins, J.L.; Yun, X.; Bachmann, E.R.; McGhee, R.B.; Zyda, M.J. An extended Kalman filter for quaternion-based orientation estimation using MARG sensors. In Proceedings of the IEEE/RSJ International Conference on Intelligent Robots and Systems, Maui, HI, USA, 29 October–3 November 2001; pp. 2003–2011.
Sabatini, A.M. Kalman-Filter-Based Orientation Determination Using Inertial/Magnetic Sensors: Observability Analysis and Performance Evaluation. Sensors 2011, 11, 9182–9206. [Google Scholar] [CrossRef] [PubMed]
Bachmann, E.R.; McGhee, R.B.; Yun, X.; Zyda, M.J. Inertial and magnetic posture tracking for inserting humans into networked virtual environments. In Proceedings of the ACM Symposium on Virtual Reality Software and Technology VRST, Baniff, AB, Canada, 15–17 November 2001; pp. 9–16.
Madgwick, S.O.H.; Harrison, A.J.L.; Vaidyanathan, A. Estimation of IMU and MARG orientation using a gradient descent algorithm. In Proceedings of the IEEE International Conference on Rehabilitation Robotics, Zurich, Switzerland, 29 June–1 July 2011; pp. 1–7.
Alam, F.; ZhaiHe, Z.; Jia, H. A Comparative Analysis of Orientation Estimation Filters using MEMS based IMU. In Proceedings of the International Conference on Research in Science, Engineering and Technology, Dubai, UAE, 21–22 March 2014.
Tian, Y.; Wei, H.; Tan, J. An adaptive-gain complementary filter for real-time human motion tracking with MARG sensors in free-living environments. IEEE Trans. Neural Syst. Rehabil. Eng. 2013, 21, 254–264. [Google Scholar] [CrossRef] [PubMed]
Fourati, H.; Manamanni, N.; Afilal, L.; Handrich, Y. Complementary observer for body segments motion capturing by inertial and magnetic sensors. IEEE/ASME Trans. Mechatronics 2014, 19, 149–157. [Google Scholar] [CrossRef]
Yoo, T.S.; Hong, S.K.; Yoon, H.M.; Park, S. Gain-Scheduled Complementary Filter Design for a MEMS Based Attitude and Heading Reference System. Sensors 2011, 11, 3816–3830. [Google Scholar] [CrossRef] [PubMed]
Li, W.; Wang, J. Effective adaptive kalman filter for mems-imu/magnetometers integrated attitude and heading reference systems. J. Navig. 2013, 66, 99–113. [Google Scholar] [CrossRef]
Makni, A.; Fourati, H.; Kibangou, A.Y. Adaptive kalman filter for MEMS-IMU based attitude estimation under external acceleration and parsimonious use of gyroscopes. In Proceedings of the European Control Conference (ECC), Strasbourg, France, 24–27 June 2014; pp. 1379–1384.
Keun Lee, J.; Park, E.J.; Robinovitch, S.N. Estimation of attitude and external acceleration using inertial sensor measurement during various dynamic conditions. IEEE Trans. Instrum. Meas. 2012, 61, 2262–2273. [Google Scholar]
Breckenridge, W.G. Quaternions—Proposed Standard Conventions; NASA Jet Propulsion Laboratory, Technical Report Interoffice Memorandum IOM 343-79-1199; 1999. [Google Scholar]
De Franceschi, M.; Zardi, D. Evaluation of Cut-Off Frequency and Correction of Filter-Induced Phase Lag and Attenuation in Eddy Covariance Analysis of Turbulence Data. Bound.-Layer Meteorol. 2003, 108, 289–303. [Google Scholar] [CrossRef]
Shoemaker, K. Animating Rotation with Quaternion Curves. In Proceedings of the Special Interest Group on Graphics and Interactive Techniques (SIGGRAPH), 22–26 July 1985; pp. 245–254.
Lee, G.H.; Achtelik, M.; Fraundorfer, F.; Pollefeys, M.; Siegwart, R. A Benchmarking Tool for MAV Visual Pose Estimation. In Proceedings of the Conference on Control, Automation, Robotics and Vision (ICARCV), Singapore, 7–10 December 2010; pp. 1541–1546.

© 2015 by the authors; licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution license (http://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Valenti, R.G.; Dryanovski, I.; Xiao, J. Keeping a Good Attitude: A Quaternion-Based Orientation Filter for IMUs and MARGs. Sensors 2015, 15, 19302-19330. https://doi.org/10.3390/s150819302

AMA Style

Valenti RG, Dryanovski I, Xiao J. Keeping a Good Attitude: A Quaternion-Based Orientation Filter for IMUs and MARGs. Sensors. 2015; 15(8):19302-19330. https://doi.org/10.3390/s150819302

Chicago/Turabian Style

Valenti, Roberto G., Ivan Dryanovski, and Jizhong Xiao. 2015. "Keeping a Good Attitude: A Quaternion-Based Orientation Filter for IMUs and MARGs" Sensors 15, no. 8: 19302-19330. https://doi.org/10.3390/s150819302

APA Style

Valenti, R. G., Dryanovski, I., & Xiao, J. (2015). Keeping a Good Attitude: A Quaternion-Based Orientation Filter for IMUs and MARGs. Sensors, 15(8), 19302-19330. https://doi.org/10.3390/s150819302

Article Menu

Keeping a Good Attitude: A Quaternion-Based Orientation Filter for IMUs and MARGs

Abstract

1. Introduction

2. Previous Work

3. Background Theory

4. Quaternion from Earth-Field Observations

4.1. Quaternion from Accelerometer Readings

4.2. Quaternion from Magnetometer Readings

5. Quaternion-Based Complementary Filter

5.1. Prediction

5.2. Correction

5.2.1. Accelerometer-Based Correction

5.2.2. Magnetometer-Based Correction

5.3. Adaptive Gain

5.4. Filter Initialization and Bias Estimation

6. Experiments

6.1. MAV Flight Experiment

6.2. Magnetic Disturbances

6.3. High Nongravitational Acceleration

7. Conclusions

Acknowledgments

Author Contributions

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI