Mahony Orientation Filter¶
This estimator proposed by Robert Mahony et al. [Mahony2008] is formulated as a deterministic kinematic observer on the Special Orthogonal group SO(3) driven by an instantaneous attitude and angular velocity measurements.
By exploiting the geometry of the special orthogonal group a related observer, the passive complementary filter, is derived that decouples the gyro measurements from the reconstructed attitude in the observer inputs.
Direct and passive filters are extended to estimate gyro bias online. This leads to an observer on SO(3), termed the Explicit Complementary Filter, that requires only accelerometer and gyro outputs, suitable for hardware implementation, and providing good estimates as well as online gyro bias computation.
Sensor Models¶
The gyroscopes measure angular velocity in the bodyfixed frame, whose error model is:
where \(\Omega\) is the true angular velocity, \(b\) is a constant (or slow timevarying) bias, and \(\mu\) is an additive measurement noise.
An ideal accelerometer measures the instantaneous linear acceleration \(\dot{\mathbf{v}}\) of the body frame minus the gravitational acceleration field \(\mathbf{g}_0\). In practice, the output \(\mathbf{a}\) of a triaxial accelerometer has also an added bias \(\mathbf{b}_a\) and noise \(\mu_a\).
where \(\mathbf{g}_0\approx 9.8\). Under quasistatic conditions it is common to normalize this vector so that:
where \(e_3=\frac{\mathbf{g}_0}{\mathbf{g}_0}=\begin{bmatrix}0 & 0 & 1\end{bmatrix}^T\).
A magnetometer provides measurements for the magnetic field
where \(\mathbf{h}\) is Earth’s magnetic field as measured at the inertial frame, \(\mathbf{B}_m\) is the local magnetic disturbance and \(\mu_b\) is the measurement noise.
The magnetic intensity is irrelevant for the estimation of the attitude, and only the direction of the geomagnetic field will be used. Normally, this measurement is also normalized:
The measurement vectors \(\mathbf{v}_a\) and \(\mathbf{v}_m\) are used to build an instantaneous algebraic rotation \(\mathbf{R}\):
where \(\mathbf{v}_m^*\) is the referential direction of the local magnetic field.
The corresponding weights \(\lambda_1\) and \(\lambda_2\) are chosen depending on the relative confidence in the sensor outputs.
Two degrees of freedom in the rotation matrix are resolved using the acceleration readings (tilt) and the final degree of freedom is resolved using the magnetometer (heading.)
The system considered is the kinematics:
where \(\lfloor\Omega\rfloor_\times\) denotes the skewsymmetric matrix of \(\Omega=\begin{bmatrix}\Omega_X & \Omega_Y & \Omega_Z\end{bmatrix}^T\):
The inverse operation taking the skewsymmetric matrix into its associated vector is \(\Omega=\mathrm{vex}(\lfloor\Omega\rfloor_\times)\).
The kinematics can also be written in terms of the quaternion representation in SO(3):
where \(\mathbf{q}=\begin{pmatrix}q_w & \mathbf{q}_v\end{pmatrix}=\begin{pmatrix}q_w & q_x & q_y & q_z\end{pmatrix}\) represents a unit quaternion, and \(\mathbf{p}(\Omega)\) represents the unitary pure quaternion associated to the angular velocity \(\mathbf{p}(\Omega)=\begin{pmatrix}0 & \Omega_X & \Omega_Y & \Omega_Z\end{pmatrix}\)
Warning
The product of two quaternions \(\mathbf{p}\) and \(\mathbf{q}\) is the Hamilton product defined as:
Error Criteria¶
We denote \(\hat{\mathbf{R}}\) as the estimation of the bodyfixed rotation matrix \(\mathbf{R}\). The used estimation error is the relative rotation from the bodyfixed frame to the estimator frame:
Mahony’s proposed observer, based on Lyapunov stability analysis, yields the cost function:
The goal of attitude estimation is to provide a set of dynamics for an estimate \(\hat{\mathbf{R}}(t)\in SO(3)\) to drive the error rotation \(\tilde{\mathbf{R}}(t)\to \mathbf{I}_3\), which in turn would drive \(\hat{\mathbf{R}}\to\mathbf{R}\).
The general form of the observer is termed as a Complementary Filter on SO(3):
where \(k_P>0\) and the term \(\mathbf{R}\Omega + k_P\hat{\mathbf{R}}\omega\) is expressed in the inertial frame.
The innovation or correction term \(\omega\), derived from the error \(\tilde{\mathbf{R}}\), can be thought of as a nonlinear approximation of the error between \(\mathbf{R}\) and \(\hat{\mathbf{R}}\) as measured from the frame of reference associated with \(\hat{\mathbf{R}}\).
When no correction term is used (\(\omega=0\)) the error rotation \(\tilde{\mathbf{R}}\) will be constant.
If \(\tilde{\mathbf{q}}=\begin{pmatrix}\tilde{q}_w & \tilde{\mathbf{q}}_v\end{pmatrix}\) is the quaternion related to \(\tilde{\mathbf{R}}\), then:
Explicit Complementary Filter¶
Let \(\mathbf{v}_{0i}\in\mathbb{R}^3\) denote a set of \(n\geq 2\) known directions in the inertial (fixed) frame of reference, where the directions are not collinear, and \(\mathbf{v}_{i}\in\mathbb{R}^3\) are their associated measurements. The measurements are bodyfixed frame observations of the fixed inertial directions:
where \(\mu_i\) is a process noise. We assume that \(\mathbf{v}_{0i}=1\) and normalize all measurements to force \(\mathbf{v}_i=1\).
For \(n\) measures, the global cost becomes:
where \(\mathbf{M}>0\) is a positive definite matrix if \(n\geq 3\), or positive semidefinite if \(n=2\):
with:
The weights \(k_i>0\) are chosen depending on the relative confidence in the measured directions.
Lowcost IMUs typically measure gravitational, \(\mathbf{a}\), and magnetic, \(\mathbf{m}\), vector fields.
In this case, the cost function \(E_\mathrm{mes}\) becomes:
Tip
When the IMU is subject to high magnitude accelerations (takeoff, landing manoeuvres, etc.) it may be wise to reduce the relative weighing of the accelerometer data (\(k_1 \ll k_2\)) compared to the magnetometer data. Conversely, when the IMU is mounted in the proximity of powerful electric motors leading to low confidence in the magnetometer readings choose \(k_1 \gg k_2\).
Expressing the kinematics of the Explicit Complementary Filter as quaternions we get:
The estimated attitude rate of change \(\dot{\hat{\mathbf{q}}}\) is multiplied with the samplerate \(\Delta t\) to integrate the angular displacement, which is finally added to the previous attitude, obtaining the newest estimated attitude.
References
[Mahony2008]  Robert Mahony, Tarek Hamel, and JeanMichel Pflimlin. Nonlinear Complementary Filters on the Special Orthogonal Group. IEEE Transactions on Automatic Control, Institute of Electrical and Electronics Engineers, 2008, 53 (5), pp.12031217. (https://hal.archivesouvertes.fr/hal00488376/document) 
[Mahony2005]  Robert Mahony, Tarek Hamel, and JeanMichel Pflimlin. Complementary filter design on the special orthogonal group SO(3). Proceedings of the 44th IEEE Conference on Decision and Control, and the European Control Conference 2005. Seville, Spain. December 1215, 2005. (https://folk.ntnu.no/skoge/prost/proceedings/cdcecc05/pdffiles/papers/1889.pdf) 
[Euston]  Mark Euston, Paul W. Coote, Robert E. Mahony, Jonghyuk Kim, and Tarek Hamel. A complementary filter for attitude estimation of a fixedwing UAV. IEEE/RSJ International Conference on Intelligent Robots and Systems, 340345. 2008. (http://users.cecs.anu.edu.au/~Jonghyuk.Kim/pdf/2008_Euston_iros_v1.04.pdf) 
[Hamel]  Tarek Hamel and Robert Mahony. Attitude estimation on SO(3) based on direct inertial measurements. IEEE International Conference on Robotics and Automation. ICRA 2006. pp. 21702175 (http://users.cecs.anu.edu.au/~Robert.Mahony/Mahony_Robert/2006_MahHamPflC68.pdf) 

class
ahrs.filters.mahony.
Mahony
(gyr: numpy.ndarray = None, acc: numpy.ndarray = None, mag: numpy.ndarray = None, frequency: float = 100.0, k_P: float = 1.0, k_I: float = 0.3, q0: numpy.ndarray = None, b0: numpy.ndarray = None, **kwargs)¶ Mahony’s Nonlinear Complementary Filter on SO(3)
If
acc
andgyr
are given as parameters, the orientations will be immediately computed with methodupdateIMU
.If
acc
,gyr
andmag
are given as parameters, the orientations will be immediately computed with methodupdateMARG
.Parameters:  acc (numpy.ndarray, default: None) – Nby3 array with measurements of acceleration in m/s^2
 gyr (numpy.ndarray, default: None) – Nby3 array with measurements of angular velocity in rad/s
 mag (numpy.ndarray, default: None) – Nby3 array with measurements of magnetic field in mT
 frequency (float, default: 100.0) – Sampling frequency in Herz
 Dt (float, default: 0.01) – Sampling step in seconds. Inverse of sampling frequency. Not required if frequency value is given
 k_P (float, default: 1.0) – Proportional filter gain
 k_I (float, default: 0.3) – Integral filter gain
 q0 (numpy.ndarray, default: None) – Initial orientation, as a versor (normalized quaternion).
Variables:  gyr (numpy.ndarray) – Nby3 array with N gyroscope samples.
 acc (numpy.ndarray) – Nby3 array with N accelerometer samples.
 mag (numpy.ndarray) – Nby3 array with N magnetometer samples.
 frequency (float) – Sampling frequency in Herz.
 Dt (float) – Sampling step in seconds. Inverse of sampling frequency.
 k_P (float) – Proportional filter gain.
 k_I (float) – Integral filter gain.
 q0 (numpy.ndarray) – Initial orientation, as a versor (normalized quaternion)
Raises: ValueError
– When dimension of input array(s)acc
,gyr
, ormag
are not equal.Examples
Assuming we have 3axis sensor data in Nby3 arrays, we can simply give these samples to their corresponding type. The Mahony algorithm can work solely with gyroscope samples, although the use of accelerometer samples is much encouraged.
The easiest way is to directly give the full array of samples to their matching parameters.
>>> from ahrs.filters import Mahony >>> orientation = Mahony(gyr=gyro_data, acc=acc_data) # Using IMU
The estimated quaternions are saved in the attribute
Q
.>>> type(orientation.Q), orientation.Q.shape (<class 'numpy.ndarray'>, (1000, 4))
If we desire to estimate each sample independently, we call the corresponding method.
orientation = Mahony() Q = np.tile([1., 0., 0., 0.], (num_samples, 1)) # Allocate for quaternions for t in range(1, num_samples): Q[t] = orientation.updateIMU(Q[t1], gyr=gyro_data[t], acc=acc_data[t])
Further on, we can also use magnetometer data.
>>> orientation = Mahony(gyr=gyro_data, acc=acc_data, mag=mag_data) # Using MARG
This algorithm is dynamically adding the orientation, instead of estimating it from static observations. Thus, it requires an initial attitude to build on top of it. This can be set with the parameter
q0
:>>> orientation = Mahony(gyr=gyro_data, acc=acc_data, q0=[0.7071, 0.0, 0.7071, 0.0])
If no initial orientation is given, then an attitude using the first samples is estimated. This attitude is computed assuming the sensors are straped to a system in a quasistatic state.
A constant sampling frequency equal to 100 Hz is used by default. To change this value we set it in its parameter
frequency
. Here we set it, for example, to 150 Hz.>>> orientation = Mahony(gyr=gyro_data, acc=acc_data, frequency=150.0)
Or, alternatively, setting the sampling step (\(\Delta t = \frac{1}{f}\)):
>>> orientation = Mahony(gyr=gyro_data, acc=acc_data, Dt=1/150)
This is specially useful for situations where the sampling rate is variable:
orientation = Mahony() Q = np.tile([1., 0., 0., 0.], (num_samples, 1)) # Allocate for quaternions for t in range(1, num_samples): orientation.Dt = new_sample_rate Q[t] = orientation.updateIMU(Q[t1], gyr=gyro_data[t], acc=acc_data[t])
Mahony’s algorithm uses an explicit complementary filter with two gains \(k_P\) and \(k_I\) to correct the estimation of the attitude. These gains can be set in the parameters too:
>>> orientation = Mahony(gyr=gyro_data, acc=acc_data, kp=0.5, ki=0.1)
Following the experimental settings of the original article, the gains are, by default, \(k_P=1\) and \(k_I=0.3\).

updateIMU
(q: numpy.ndarray, gyr: numpy.ndarray, acc: numpy.ndarray, dt: float = None) → numpy.ndarray¶ Attitude Estimation with a IMU architecture.
Parameters:  q (numpy.ndarray) – Apriori quaternion.
 gyr (numpy.ndarray) – Sample of triaxial Gyroscope in rad/s.
 acc (numpy.ndarray) – Sample of triaxial Accelerometer in m/s^2.
 dt (float, default: None) – Time step, in seconds, between consecutive Quaternions.
Returns: q – Estimated quaternion.
Return type: numpy.ndarray
Examples
>>> orientation = Mahony() >>> Q = np.tile([1., 0., 0., 0.], (num_samples, 1)) # Allocate for quaternions >>> for t in range(1, num_samples): ... Q[t] = orientation.updateIMU(Q[t1], gyr=gyro_data[t], acc=acc_data[t])

updateMARG
(q: numpy.ndarray, gyr: numpy.ndarray, acc: numpy.ndarray, mag: numpy.ndarray, dt: float = None) → numpy.ndarray¶ Attitude Estimation with a MARG architecture.
Parameters:  q (numpy.ndarray) – Apriori quaternion.
 gyr (numpy.ndarray) – Sample of triaxial Gyroscope in rad/s.
 acc (numpy.ndarray) – Sample of triaxial Accelerometer in m/s^2.
 mag (numpy.ndarray) – Sample of triaxail Magnetometer in uT.
 dt (float, default: None) – Time step, in seconds, between consecutive Quaternions.
Returns: q – Estimated quaternion.
Return type: numpy.ndarray
Examples
>>> orientation = Mahony() >>> Q = np.tile([1., 0., 0., 0.], (num_samples, 1)) # Allocate for quaternions >>> for t in range(1, num_samples): ... Q[t] = orientation.updateMARG(Q[t1], gyr=gyro_data[t], acc=acc_data[t], mag=mag_data[t])