An efficient orientation filter for inertial and

inertial/magnetic sensor arrays

Sebastian O.H. Madgwick
April 30, 2010

Abstract
This report presents a novel orientation filter applicable to IMUs consisting of
tri-axis gyroscopes and accelerometers, and MARG sensor arrays that also include
tri-axis magnetometers. The MARG implementation incorporates magnetic distortion
and gyroscope bias drift compensation. The filter uses a quaternion representation,
allowing accelerometer and magnetometer data to be used in an analytically derived
and optimised gradient-descent algorithm to compute the direction of the gyroscope
measurement error as a quaternion derivative. The benefits of the filter include: (1)
computationally inexpensive; requiring 109 (IMU) or 277 (MARG) scalar arithmetic
operations each filter update, (2) e↵ective at low sampling rates; e.g. 10 Hz, and (3)
contains 1 (IMU) or 2 (MARG) adjustable parameters defined by observable system
characteristics. Performance was evaluated empirically using a commercially available
orientation sensor and reference measurements of orientation obtained using an optical
measurement system. A simple calibration method is presented for the use of the
optical measurement equipment in this application. Performance was also benchmarked
against the propriety Kalman-based algorithm of orientation sensor. Results indicate
the filter achieves levels of accuracy exceeding that of the Kalman-based algorithm;
< 0.6 static RMS error, < 0.8 dynamic RMS error. The implications of the low
computational load and ability to operate at low sampling rates open new opportunities
for the use of IMU and MARG sensor arrays in real-time applications of limited power
or processing resources or applications that demand extremely high sampling rates.

1
Contents
1 Introduction 3

2 Quaternion representation 4

3 Filter derivation 6
3.1 Orientation from angular rate . . . . . . . . . . . . . . . . . . . . . . . . . . 6
3.2 Orientation from vector observations . . . . . . . . . . . . . . . . . . . . . . 6
3.3 Filter fusion algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
3.4 Magnetic distortion compensation . . . . . . . . . . . . . . . . . . . . . . . . 11
3.5 Gyroscope bias drift compensation . . . . . . . . . . . . . . . . . . . . . . . 12
3.6 Filter gains . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

4 Experimentation 14
4.1 Equipment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
4.2 Orientation from optical measurements . . . . . . . . . . . . . . . . . . . . . 14
4.3 Calibration of frame alignments . . . . . . . . . . . . . . . . . . . . . . . . . 16
4.4 Experimental proceedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

5 Results 19
5.1 Typical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
5.2 Static and dynamic performance . . . . . . . . . . . . . . . . . . . . . . . . 21
5.3 Filter gain vs. performance . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
5.4 Sampling rate vs. performance . . . . . . . . . . . . . . . . . . . . . . . . . . 23
5.5 Gyroscope bias drift . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

6 Discussion 24

7 Conclusions 25

A IMU filter implementation optimised in C 29

B MARG filter implementation optimised in C 30

2
1 Introduction
The accurate measurement of orientation plays a critical role in a range of fields includ-
ing: aerospace [1, 2, 3], robotics [4, 5], navigation [6, 7] and human motion analysis [8, 9]
and machine interaction [10]. Whilst a variety of technologies enable the measurement of
orientation, inertial based sensory systems have the advantage of being completely self con-
tained such that the measurement entity is constrained neither in motion nor to any specific
environment or location. An IMU (Inertial Measurement Unit) consists of gyroscopes and
accelerometers enabling the tracking of rotational and translational movements. In order to
measure in three dimensions, tri-axis sensors consisting of 3 mutually orthogonal sensitive
axes are required. A MARG (Magnetic, Angular Rate, and Gravity) sensor is a hybrid IMU
which incorporates a tri-axis magnetometer. An IMU alone can only measure an attitude
relative to the direction of gravity which is sufficient for many applications [4, 2, 8, 1]. MARG
systems, also known as AHRS (Attitude and Heading Reference Systems) are able to provide
a complete measurement of orientation relative to the direction of gravity and the earth’s
magnetic field.
A gyroscope measures angular velocity which, if initial conditions are known, may be inte-
grated over time to compute the sensor’s orientation [11, 12]. Precision gyroscopes, ring laser
for example, are too expensive and bulky for most applications and so less accurate MEMS
(Micro Electrical Mechanical System) devices are used in a majority of applications [13]. The
integration of gyroscope measurement errors will lead to an accumulating error in the calcu-
lated orientation. Therefore, gyroscopes alone cannot provide an absolute measurement of
orientation. An accelerometer and magnetometer will measure the earth’s gravitational and
magnetic fields respectively and so provide an absolute reference of orientation. However,
they are likely to be subject to high levels of noise; for example, accelerations due to motion
will corrupt measured direction of gravity. The task of an orientation filter is to compute
a single estimate of orientation through the optimal fusion of gyroscope, accelerometer and
magnetometer measurements.
The Kalman filter [14] has become the accepted basis for the majority of orientation filter
algorithms [4, 15, 16, 17] and commercial inertial orientation sensors; xsens [18], micro-strain
[19], VectorNav [20], Intersense [21], PNI [22] and Crossbow [23] all produce systems founded
on its use. The widespread use of Kalman-based solutions are a testament to their accuracy
and e↵ectiveness, however, they have a number of disadvantages. They can be complicated to
implement which is reflected by the numerous solutions seen in the subject literature [3, 4, 15,
16, 17, 24, 25, 26, 27, 28, 29, 30, 31, 32]. The linear regression iterations, fundamental to the
Kalman process, demand sampling rates far exceeding the subject bandwidth; for example, a
sampling rate between 512 Hz [18] and 30 kHz [19] may be used for a human motion caption
application. The state relationships describing rotational kinematics in three-dimensions
typically require large state vectors and an extended Kalman filter implementation [4, 17, 24]
to linearise the problem.
These challenges demand a large computational load for implementation of Kalman-
based solutions and provide a clear motivation for alternative approaches. Many previous
approaches to address these issues have implemented either fuzzy processing [2, 5] or fixed
filters [33] to favour accelerometer measurements of orientation at low angular velocities and
the integrated gyroscope measurements at high angular velocities. Such an approach is simple

3
but may only be e↵ective under limited operating conditions. Bachman et al [34] proposed
an alternative approach where the filter achieves an optimal fusion of measurements data at
all angular velocities. However, the process requires a least squares regression, which also
brings in an associated computational load. Mahony et al [35] developed the complementary
filter which is shown to be an efficient and e↵ective solution; however, performance is only
validated for an IMU.
This report introduces novel orientation filter that is applicable to both IMUs and MARG
sensor arrays addressing issues of computational load and parameter tuning associated with
Kalman-based approaches. The filter employs a quaternion representation of orientation (as
in: [34, 17, 24, 30, 32]) to describe the coupled nature of orientations in three-dimensions and
is not subject to the problematic singularities associated with an Euler angle representation1 .
A complete derivation and empirical evaluation of the new filter is presented. Its performance
is benchmarked against an existing commercial filter and verified with optical measurement
system. Innovative aspects of the proposed filter include: a single adjustable parameter
defined by observable systems characteristics; an analytically derived and optimised gradient-
descent algorithm enabling performance at low sampling rates; an on-line magnetic distortion
compensation algorithm; and gyroscope bias drift compensation.

2 Quaternion representation
A quaternion is a four-dimensional complex number that can be used to represent the ori-
entation of a ridged body or coordinate frame in three-dimensional space. An arbitrary
orientation of frame B relative to frame A can be achieved through a rotation of angle ✓
around an axis A r̂ defined in frame A. This is represented graphically in figure 1 where the
mutually orthogonal unit vectors x̂A , ŷA and ẑA , and x̂B , ŷB and ẑB define the principle axis
of coordinate frames A and B respectively. The quaternion describing this orientation, A B q̂,
is defined by equation (1) where rx , ry and rz define the components of the unit vector A r̂ in
the x, y and z axes of frame A respectively. A notation system of leading super-scripts and
sub-scripts adopted from Craig [37] is used to denote the relative frames of orientations and
vectors. A leading sub-script denotes the frame being described and a leading super-script
denotes the frame this is with reference to. For example, A B q̂ describes the orientation of
A
frame B relative to frame A and r̂ is a vector described in frame A. Quaternion arithmetic
often requires that a quaternion describing an orientation is first normalised. It is therefore
conventional for all quaternions describing an orientation to be of unit length.

A
⇥ ⇤ ⇥ ⇤
B q̂ = q1 q2 q3 q4 = cos 2✓ rx sin 2✓ ry sin 2✓ rz sin 2✓ (1)
The quaternion conjugate, denoted by ⇤ , can be used to swap the relative frames described
by an orientation. For example, B A
A q̂ is the conjugate of B q̂ and describes the orientation of
A
frame A relative to frame B. The conjugate of B q̂ is defined by equation (2).

A ⇤
⇥ ⇤
B q̂ =B
A q̂ = q1 q2 q3 q4 (2)
1
Kuipers [36] o↵ers a comprehensive introduction to this use of quaternions.

4
 ẑ A
ẑ B

A
r̂
✓ ŷ B

ŷ A

x̂A x̂B

Figure 1: The orientation of frame B is achieved by a rotation, from alignment with frame
A, of angle ✓ around the axis A r.

The quaternion product, denoted by ⌦, can be used to define compound orientations.

For example, for two orientations described by A B A
B q̂ and C q̂, the compounded orientation C q̂
can be defined by equation (3).

A
C q̂ =B A
C q̂ ⌦ B q̂ (3)
For two quaternions, a and b, the quaternion product can be determined using the
Hamilton rule and defined as equation (4). A quaternion product is not commutative; that
is, a ⌦ b 6= b ⌦ a.

⇥ ⇤ ⇥ ⇤
a ⌦ b = a1 a2 a3 a4 ⌦ b 1 b 2 b 3 b 4
2 3T
a1 b 1 a2 b 2 a3 b 3 a4 b 4
6a1 b2 + a2 b1 + a3 b4 a4 b3 7 (4)
=6 7
4 a1 b 3 a2 b 4 + a3 b 1 + a4 b 2 5
a1 b 4 + a2 b 3 a3 b 2 + a4 b 1
A three dimensional vector can be rotated by a quaternion using the relationship de-
scribed in equation (5) [36]. A v and B v are the same vector described in frame A and frame
B respectively where each vector contains a 0 inserted as the first element to make them 4
element row vectors.

B A ⇤
v=A A
B q̂ ⌦ v ⌦ B q̂ (5)
The orientation described by A B q̂ can be represented as the rotation matrix A
B R defined
by equation (6) [36].
2 2 3
2q1 1 + 2q22 2(q2 q3 + q1 q4 ) 2(q2 q4 q1 q3 )
A 4
B R = 2(q2 q3 q1 q4 ) 2q12 1 + 2q32 2(q3 q4 + q1 q2 )5 (6)
2(q2 q4 + q1 q3 ) 2(q3 q4 q1 q2 ) 2q12 1 + 2q42
Euler angles , ✓ and in the so called aerospace sequence [36] describe an orientation
of frame B achieved by the sequential rotations, from alignment with frame A, of around

5
 A
ẑB , ✓ around ŷB , and around x̂B . This Euler angle representation of B q̂ is defined by
equations (7), (8) and (9).

= Atan2 2q2 q3 2q1 q4 , 2q12 + 2q22 1 (7)

1
✓= sin (2q2 q4 + 2q1 q3 ) (8)

= Atan2 2q3 q4 2q1 q2 , 2q12 + 2q42 1 (9)

3 Filter derivation
3.1 Orientation from angular rate
A tri-axis gyroscope will measure the angular rate about the x, y and z axes of the senor
frame, termed !x , !y and !z respectively. If these parameters (in rads 1 ) are arranged into
the vector S ! defined by equation (10), the quaternion derivative describing the rate of
change of orientation of the earth frame relative to the sensor frame SE q̇ can be calculated
[38] as equation (11).

S
⇥ ⇤
! = 0 !x !y !z (10)

1S
S
q̂ ⌦ S !
E q̇ = (11)
2E
The orientation of the earth frame relative to the sensor frame at time t, E S q!,t , can
S
be computed by numerically integrating the quaternion derivative E q̇!,t as described by
equations (12) and (13) provided that initial conditions are known. In these equations, S !t
is the angular rate measured at time t, t is the sampling period and SE q̂est,t-1 is the previous
estimate of orientation. The sub-script ! indicates that the quaternion is calculated from
angular rates.

S 1S
E q̇ !,t = q̂est,t-1 ⌦ S !t (12)
2E
S
E q!,t = SE q̂est,t-1 + SE q̇ !,t t (13)

3.2 Orientation from vector observations

A tri-axis accelerometer will measure the magnitude and direction of the field of gravity in the
sensor frame compounded with linear accelerations due to motion of the sensor. Similarly, a
tri-axis magnetometer will measure the magnitude and direction of the earth’s magnetic field
in the sensor frame compounded with local magnetic flux and distortions. In the context
of an orientation filter, it will initially be assumed that an accelerometer will measure only
gravity, a magnetometer will measure only the earth’s magnetic field.

6
 If the direction of an earth’s field is known in the earth frame, a measurement of the
field’s direction within the sensor frame will allow an orientation of the sensor frame relative
to the earth frame to be calculated. However, for any given measurement there will not be
a unique sensor orientation solution, instead there will infinite solutions represented by all
those orientations achieved by the rotation the true orientation around an axis parallel with
the field. In some applications it may be acceptable to use an Euler angle representation
allowing an incomplete solution to be found as two known Euler angles and one unknown [5];
the unknown angle being the rotation around an axis parallel with the direction of the field.
A quaternion representation requires a complete solution to be found. This may be achieved
through the formulation of an optimisation problem where an orientation of the sensor, SE q̂,
is that which aligns a predefined reference direction of the field in the earth frame, E d, ˆ
S
with the measured direction of the field in the sensor frame, ŝ, using the rotation operation
described by equation (5). Therefore SE q̂ may be found as the solution to (14) where equation
(15) defines the objective function. The components of each vector are defined in equations
(16) to (18).

ˆ S ŝ)
min f (SE q̂, E d, (14)
S q̂2<4
E

ˆ S ŝ) = S q̂ ⇤ ⌦ E dˆ ⌦ S q̂
f (SE q̂, E d, S
ŝ (15)
E E

S
⇥ ⇤
E q̂ = q1 q2 q3 q4 (16)
⇥ ⇤
E
dˆ = 0 dx dy dz (17)

S
⇥ ⇤
ŝ = 0 sx sy sz (18)
Many optimisation algorithms exist but the gradient descent algorithm is one of the
simplest to both implement and compute. Equations (19) describes the gradient descent
algorithm for n iterations resulting in an orientation estimation of SE q̂n+1 based on an ‘initial
guess’ orientation SE q̂0 and a step-size µ. Equation (20) computes the gradient of the solution
surface defined by the objective function and its Jacobian; simplified to the 3 row vectors
defined by equations (21) and (22) respectively.
ˆ S ŝ)
rf (SE q̂ k , E d,
S
E qk+1 = SE q̂ k µ , k = 0, 1, 2...n (19)
S E
rf ( q̂ , d, ŝ)ˆ S
E k

ˆ S ŝ) = J T (S q̂ , E d)f
rf (SE q̂ k , E d, ˆ (S q̂ , E d,
ˆ S ŝ) (20)
E k E k

2
2dx ( 12 q32 q42 ) + 2dy (q1 q4 + q2 q3 )+
ˆ S ŝ)
f (SE q̂ k , E d, = 42dx (q2 q3 q1 q4 ) + 2dy ( 12 q22 q42 )+
2dx (q1 q3 + q2 q4 ) + 2dy (q3 q4 q1 q2 )+
3 (21)
2dz (q2 q4 q1 q3 ) sx
2dz (q1 q2 + q3 q4 ) sy 5
2dz ( 12 q22 q32 ) sz

7
 2
2dy q4 2dz q3 2dy q3 + 2dz q4
ˆ
J (SE q̂k , E d) = 4 2dx q4 + 2dz q2 2dx q3 4dy q2 + 2dz q1
2dx q3 2dy q2 2dx q4 2dy q1 4dz q2
3 (22)
4dx q3 + 2dy q2 2dz q1 4dx q4 + 2dy q1 + 2dz q2
2dx q2 + 2dz q4 2dx q1 4dy q4 + 2dz q3 5
2dx q1 + 2dy q4 4dz q3 2dx q2 + 2dy q3

Equations (19) to (22) describe the general form of the algorithm applicable to a field
predefined in any direction. However, if the direction of the field can be assumed to only
have components within 1 or 2 of the principle axis of the global coordinate frame then
the equations simplify. An appropriate convention would be to assume that the direction of
gravity defines the vertical, z axis as shown in equation (23). Substituting E ĝ and normalised
accelerometer measurement S â for E dˆ and S ŝ respectively, in equations (21) and (22) yields
equations (25) and (26).

E
⇥ ⇤
ĝ = 0 0 0 1 (23)

S
⇥ ⇤
â = 0 ax ay az (24)
2 3
2(q2 q4 q1 q3 ) ax
fg (SE q̂, S â) = 4 2(q1 q2 + q3 q4 ) ay 5 (25)
2( 12 q22 q32 ) az
2 3
2q3 2q4 2q1 2q2
Jg (SE q̂) = 4 2q2 2q1 2q4 2q3 5 (26)
0 4q2 4q3 0
The earth’s magnetic field can be considered to have components in one horizontal axis
and the vertical axis; the vertical component due to the inclination of the field which is
between 65 and 70 to the horizontal in the UK [39]. This can be represented by equa-
tion (27). Substituting E b̂ and normalised magnetometer measurement S m̂ for E dˆ and S ŝ
respectively, in equations (21) and (22) yields equations (29) and (30).

E
⇥ ⇤
b̂ = 0 bx 0 bz (27)

S
⇥ ⇤
m̂ = 0 mx my mz (28)
2 3
2bx (0.5 q32 q42 ) + 2bz (q2 q4 q1 q3 ) mx
fb (SE q̂, E b̂, S m̂) = 4 2bx (q2 q3 q1 q4 ) + 2bz (q1 q2 + q3 q4 ) my 5 (29)
2bx (q1 q3 + q2 q4 ) + 2bz (0.5 q22 q32 ) mz

8
 2
2bz q3 2bz q4 4bx q3 2bz q1
Jb (SE q̂, E b̂) = 4 2bx q4 + 2bz q2 2bx q3 + 2bz q1 2bx q2 + 2bz q4
2bx q3 2bx q4 4bz q2 2bx q1 4bz q3
3 (30)
4bx q4 + 2bz q2
2bx q1 + 2bz q3 5
2bx q2

As has already been discussed, the measurement of gravity or the earth’s magnetic field
alone will not provide a unique orientation of the sensor. To do so, the measurements and
reference directions of both fields may be combined as described by equations (31) and
(32). Whereas the solution surface created by the objective functions in equations (25) and
(29) have a minimum defined by a line, the solution surface define by equation (31) has a
minimum define by a single point, provided that bx 6= 0.

S S E S fg (SE q̂, S â)
fg,b (E q̂, â, b̂, m̂) = (31)
fb (SE q̂, E b̂, S m̂)

S E JgT (SE q̂)
Jg,b (E q̂, b̂) = (32)
JbT (SE q̂, E b̂)
A conventional approach to optimisation would require multiple iterations of equation
(19) to be computed for each new orientation and corresponding senor measurements. Ef-
ficient algorithms would also require the step-size µ to be adjusted each iteration to an
optimal value; usually obtained based on the second derivative of the objective function, the
Hessian. However, these requirements considerably increase the computational load of the
algorithm and are not necessary in this application. It is acceptable to compute one iteration
per time sample provided that the convergence rate governed by µt is equal or greater than
the physical rate of change of orientation. Equation (33) calculates the estimated orienta-
tion SE q r,t computed at time t based on a previous estimate of orientation SE q̂ est,t-1 and the
objective function gradient rf defined by sensor measurements S ât and S m̂t sampled at
time t. The form of rf is chosen according to the sensors in use, as shown in equation
(34). The sub-script r indicates that the quaternion is calculated using the gradient descent
algorithm.

S rf
E q r,t = SE q̂ est,t-1 µt (33)
krf k
⇢
JgT (SE q̂est,t-1 )fg (SE q̂est,t-1 , S ât )
rf = T S (34)
Jg,b (E q̂est,t-1 , E b̂)fg,b (SE q̂est,t-1 , S â, E b̂, S m̂)
An optimal value of µt can be defined as that which ensures the convergence rate of SE q r,t
is limited to the physical orientation rate as this avoids overshooting due an unnecessarily
large step size. Therefore µt can be calculated as equation (35) where t is the sampling
period and SE q̇!,t is the physical orientation rate measured by gyroscopes and ↵ is an aug-
mentation of µ to account for noise in accelerometer and magnetometer measurements.

9
 S
µt = ↵ E q̇!,t t, ↵ > 1 (35)

3.3 Filter fusion algorithm

An estimated orientation of the sensor frame relative to the earth frame, SE qest,t , is obtained
through the fusion of the orientation calculations, SE q!,t and SE qr,t ; calculated using equations
(13) and (33) respectively. The fusion of SE q̂!,t and SE qr,t is described by equation (36) where
t and (1 t ) are weights applied to each orientation calculation.

S S S
E qest,t = t E qr,t + (1 t )E q!,t , 0 t 1 (36)
An optimal value of t can be defined as that which ensures the weighted divergence
of SE q! is equal to the weighted convergence of SE qr . This is represented by equation (37)
where µtt is the convergence rate of SE qr and is the divergence rate of SE q! expressed as
the magnitude of a quaternion derivative corresponding to the gyroscope measurement error.
Equation (37) can be rearranged to define t as equation (38).
µt
(1 t) = t (37)
t

t = µt (38)
t
+
Equations (36) and (38) ensure the optimal fusion of SE q!,t and SE qr,t assuming that the
convergence rate of SE q r governed by ↵ is equal or greater than the physical rate of change
of orientation. Therefore ↵ has no upper bound. If ↵ is assumed to be very large then
µt , defined by equation (35), also becomes very large and the orientation filter equations
simplify. A large value of µt used in equation (33) means that SE q̂est,t-1 becomes negligible
and the equation can be re-written as equation (39).

S rf
E q r,t ⇡ µt (39)
krf k
The definition of t in equation (38) also simplifies as the term in the denominator
becomes negligible and the equation can be rewritten as equation (40). It is possible from
equation (40) to also assume that t ⇡ 0.
t
t ⇡ (40)
µt
Substituting equations (13), (39) and (40) into equation (36) directly yields equation
(41). It is important to note that in equation (41), t has been substituted as both as
equation (39) and 0.
✓ ◆ ⇣ ⌘
S t rf S S
q
E est,t = µ t + (1 0) q̂
E est,t-1 + q̇
E !,t t (41)
µt krf k

10
 Equation (41) can be simplified to equation (42) where SE q̇est,t is the estimated rate of
change of orientation defined by equation (43) and SE q̂˙ ✏,t is the direction of the error of SE q̇est,t
defined by equation (44).

S
E qest,t = SE q̂est,t-1 + SE q̇est,t t (42)

S S ˙
E q̇est,t = SE q̇ !,t E q̂✏,t (43)

S ˙ rf
E q̂✏,t = (44)
krf k
It can be seen from equations (42) to (44) that the filter calculates the orientation SE qest
by numerically integrating the estimated orientation rate SE q̇ est . The filter computes SE q̇est
as the rate of change of orientation measured by the gyroscopes, SE q̇ ! , with the magnitude
of the gyroscope measurement error, , removed in the direction of the estimated error, SE q̂˙ ✏ ,
computed from accelerometer and magnetometer measurements. Figure 2 shows a block
diagram representation of the complete orientation filter implementation for an IMU.

Accelerometer S ât J Tg (SE q̂ est,t S S

1 )f g (E q̂ est,t 1 , ât )

rf
krf k 1
z

Z
1S S
q S
Gyroscope S ! t q̂ 1 ⌦ !t .dt
kqk E q̂ est,t
2 E est,t
S
E q̇ est,t 1
z

Figure 2: Block diagram representation of the complete orientation filter for an IMU imple-
mentation

3.4 Magnetic distortion compensation

Measurements of the earth’s magnetic field will be distorted by the presence of ferromagnetic
elements in the vicinity of the magnetometer. Investigations into the e↵ect of magnetic
distortions on an orientation sensor’s performance have shown that substantial errors may
be introduced by sources including electrical appliances, metal furniture and metal structures
within a buildings construction [40, 41]. Sources of interference fixed in the sensor frame,
termed hard iron biases, can be removed through calibration [42, 43, 44, 45]. Sources of
interference in the earth frame, termed soft iron, cause errors in the measured direction of

11
the earth’s magnetic field. Declination errors, those in the horizontal plane relative to the
earth’s surface, cannot be corrected without an additional reference of heading. Inclination
errors, those in the vertical plane relative to the earth’s surface, may be compensated for as
the accelerometer provides an additional measurement of the sensor’s attitude.
The measured direction of the earth’s magnetic field in the earth frame at time t, E ĥt , can
be computed as the normalised magnetometer measurement, S m̂t , rotated by the estimated
orientation of the sensor provided by the filter, SE q̂est,t-1 ; as described by equation (45). The
e↵ect of an erroneous inclination of the measured direction earth’s magnetic field, E ĥt , can
be corrected if the filter’s reference direction of the earth’s magnetic field, E b̂t , is of the same
inclination. This is achieved by computing E b̂t as E ĥt normalised to have only components
in the earth frame x and z axes; as described by equation (46).

E
⇥ ⇤ ⇤
ĥt = 0 hx hy hz = SE q̂est,t-1 ⌦ S m̂t ⌦ SE q̂est,t -1 (45)

E
⇥ p ⇤
b̂t = 0 h2x + h2y 0 hz (46)
Compensating for magnetic distortions in this way ensures that magnetic disturbances
are limited to only a↵ect the estimated heading component of orientation. The approach also
eliminates the need for the reference direction of the earth’s magnetic field to be predefined;
a potential disadvantage of other orientation filter designs [17, 24].

3.5 Gyroscope bias drift compensation

The gyroscope zero bias will drift over time, with temperature and with motion. Any
practical implementation of an IMU or MARG sensor array must account for this. An
advantage of Kalman-based approaches is that they are able to estimate the gyroscope bias
as an additional state within the system model [26, 30, 15, 24]. However, Mahony et al [35]
showed that gyroscope bias drift may also be compensated for by simpler orientation filters
through the integral feedback of the error in the rate of change of orientation. A similar
approach will be used here.
The normalised direction of the estimated error in the rate of change of orientation, SE q̂˙ ✏ ,
may be expressed as the angular error in each gyroscope axis using equation (47); derived
as the inverse to the relationship defined in equation (11). The gyroscope bias, S !b , is
represented by the DC component of S !✏ and so may removed as the integral of S !✏ weighted
by an appropriate gain, ⇣. This would yield the compensated gyroscope measurements S !c ,
as shown in equations (48) and (49). The first element of S !c is always assumed to be 0.

S ⇤ S ˙
!✏,t = 2SE q̂est,t -1 ⌦ E q̂✏,t (47)
X
S S
!b,t = ⇣ !✏,t t (48)
t

S
!c,t = S !t S
!b,t (49)
The compensated gyroscope measurements, S !c , may then be used in place of the of
the gyroscope measurements, S !, in equation (11). The magnitude of the angular error in

12
each axis, S !✏ is equal to a quaternion derivative of unit length. Therefore the integral gain
⇣ directly defines the rate of convergence of the estimated gyroscope bias, S !b , expressed
as the magnitude of a quaternion derivative. As this process requires the use of the filter
estimate of a complete orientation, SE q̂est , it is only applicable to a MARG implementation
of the filter. Figure 3 shows a block diagram representation of the complete filter implemen-
tation for a MARG sensor array, including the magnetic distortion and gyroscope bias drift
compensation.

S
E q̂ est,t-1 ⌦ S m̂t ⌦ SE q̂ ⇤est,t-1

h q i
0 h2x + h2y 0 hz
Group 1

Magnetometer S m̂t
J Tg,b (SE q̂ est,t-1, E b̂t )f g,b (SE q̂ est,t-1, S â, E b̂t , S m̂)
Accelerometer S ât
Z
rf
.dt 2SE q̂ ⇤est,t-1 ⌦ SE q̂˙ ✏,t
krf k
z -1

⇣
Group 2
Z
1S q S
S
Gyroscope ! t q̂ ⌦ S ! c,t .dt E q̂ est,t
2 E est,t-1 kqk
S
E q̇ est,t
z -1

Figure 3: Block diagram representation of the complete orientation filter for an MARG im-
plementation including magnetic distortion (Group 1) and gyroscope drift (Group 2) com-
pensation

3.6 Filter gains

The filter gain represents all mean zero gyroscope measurement errors, expressed as the
magnitude of a quaternion derivative. The sources of error include: sensor noise, signal
aliasing, quantisation errors, calibration errors, sensor miss-alignment, sensor axis non-
orthogonality and frequency response characteristics. The filter gain ⇣ represents the rate
of convergence to remove gyroscope measurement errors which are not mean zero, also ex-
pressed as the magnitude of a quaternion derivative. These errors represent the gyroscope
bias. It is convenient to define and ⇣ using the angular quantities ! and !˙ ⇣ respectively,
where ! ˜ represents the estimated mean zero gyroscope measurement error of each axis and
!˙ ⇣ represents the estimated rate of gyroscope bias drift in each axis. Using the relation-
ship described by equation (11), may be defined by equation (50) where q̂ is any unit
quaternion. Similarly, ⇣ may be described by equation (51).

13
 r
1 ⇥ ⇤ 3
= q̂ ⌦ 0 !
˜ !
˜ !
˜ = !
˜ (50)
2 4
r
3˜
⇣= !˙ ⇣ (51)
4

4 Experimentation
4.1 Equipment
The filter was tested using the xsens MTx orientation sensor [18] containing 16 bit resolu-
tion tri-axis gyroscopes, accelerometers and magnetometers. The device and accompanying
software o↵er a mode of operations where raw sensor data may be logged at a rate of 512 Hz
and then post-processed to provide calibrated sensor measurements. The calibrated sensor
measurements could then be processed by the proposed filter to provide the estimated ori-
entation of the sensor. The software also incorporates a propriety Kalman-based orientation
filter to provide an additional estimate orientation. As both the Kalman-based algorithm
and proposed filter’s outputs could be computed using identical sensor data, the perfor-
mance of each algorithm could be evaluated relative to one-another, independent of sensor
performance.
A Vicon system, consisting of 8 MX3+ cameras connected to an MXultranet server
[46] and Nexus [47] software was used to provide reference measurement of the orientation
sensor’s actual orientation. The system is an array of IR (Infrared) sensitive cameras with
incorporated IR flood lights. The cameras are fixed at calibrated positions and orientations so
that the measurement subject is within the field of view of multiple cameras. The Cartesian
positions of IR reflective optical markers fixed to the measurement subject may then be
computed in the coordinate frame of the camera array. Cameras were fixed at a height of
approximately 2.5 m, evenly distributed around the perimeter of a 4 m by 4 m enclosure.
Each camera was orientated to face toward the centre of the room, approximately 30 to
60 to horizontal. Experiments were conducted with the measurement subject in the centre
of the room at a height of approximately 1 m. To measure the orientation of the sensor,
it was fixed to an optical orientation measurement platform specifically designed for this
application. The system was used to log the positions of optical markers at a rate of 120 Hz.

4.2 Orientation from optical measurements

The orientation measurement platform is comprised of 3 500 mm, mutually orthogonal rods,
rigidly connected at the central position along each length. Optical markers were positioned
at both ends of each rod and the orientation sensor fixed to a platform at the point where the
rods coincide. The platform was constructed from an aluminium central hub, carbon fibre
rods and assembled using adhesives to ensure that it had no magnetic properties that may
interfere with the orientation sensor’s magnetometer. Additional optical markers were placed
at arbitrary but dissimilar positions along the lengths of the rods to break the rotational
symmetry and aid the identification of each rod within the measurement data. Figure 4

14
shows an annotated photograph of the orientation measurement platform where C istart , C iend ,
C
jstart , C jend , C kstart and C kend define the measured position of each marker within the
camera frame. These positions may be used to define 3 mutually orthogonal unit vectors,
C
x̂M , C ŷM , and C ẑM within the camera frame representing the direction the x, y and z axes
of the orientation measurement platform coordinate frame; as described by equations (52),
(53) and (54). These vectors define the rotation matrix describing the orientation of the
measurement platform in the camera frame; as shown in equation (55).

C
z end
C
j end

Sensor
C
xstart

C
xend
C
j start

C
z start

Figure 4: Photograph of the orientation measurement platform

C C
C iend istart
x̂M = C Ci
(52)
k iend start k

C C
C jend jstart
ŷM = C Cj
(53)
k jend start k

C C
C kend kstart
ẑM = C Ck
(54)
k kend start k

C
⇥C C C
⇤
MR = x̂M ŷM ẑM (55)
Due to measurement errors and tolerances in the marker frame construction, the rotation
matrix defined by equation (55) cannot be considered orthogonal and so does not represent a
pure rotation. Bar-Itzhack provides a method [48] where by an optimal ‘best fit’ quaternion
may be extracted from an imprecise and non-orthogonal rotation matrix. The method
requires the construction of the symmetric 4 by 4 matrix, K, defined by equation (56),
where rmn corresponds to the element of the mth row and nth column of C M R. The optimal
C
quaternion M q̂ is found as the normalised Eigen vector corresponding to the maximum
Eigen value of K. Equation (57) defines the optimal quaternion accounting the alternative

15
quaternion element order convention assumed by the method, where v1 , v2 , v3 and v4 define
elements of the normalised Eigen vector.
2 3
r11 r22 r33 r21 + r12 r31 + r13 r23 r32
1 6 r21 + r12 r22 r11 r33 r32 + r23 r31 r13 7
K= 6 7 (56)
3 4 r31 + r13 r32 + r23 r33 r11 r22 r12 r21 5
r23 r32 r31 r13 r12 r21 r11 + r22 + r33
C
⇥ ⇤
M q̂ = v 4 v 1 v 2 v 3 (57)

4.3 Calibration of frame alignments

In order to compare the optical measurement of the platform orientation in the camera frame,
C S
M q̂, and the orientation filter’s estimated orientation of the earth in the sensor frame, E q̂est ,
it is necessary to know the alignment of the earth frame relative to the camera frame, C E q̂,
S
and the alignment of the measurement platform relative to the sensor frame, M q̂. Once these
quantities are found the optical measurement of the sensor frame orientation in the earth
frame, SE q̂meas , may be defined by equation (58). Although the use of optical measurement
equipment in this application is documentaed [26, 24, 41], little discussion is o↵ered on the
calibration of these two quantities.

S
E q̂meas =C M S
E q̂ ⌦ C q̂ ⌦ M q̂ (58)
The earth frame x and z axes are defined by the earth’s magnetic and gravitational fields
respectively. Measurements of these fields in the camera frame can be used to define the
alignment C E q̂. The direction of gravity was measured using a pendulum constructed from a
1 m length of cotton thread with a small weight fixed to one end. Optical markers where
fixed at either end of the length of cotton and the pendulum left to come to rest. Additional
optical markers were required at arbitrary fixed positions relative to the static pendulum
to break the rotational symmetry of the optical marker constellation. Figure 5 shows an
annotated photograph of the pendulum where C pstart and C pend define the position of the
optical markers in the camera frame. The mean position of each marker over a period of
time defines the direction of the pendulum in the camera frame, C p̂. This directly defines
the earth frame z axis in the camera frame, C ẑE ; as shown in equation (59).
2 3
C C z1
p̄end p̄start
C
p̂ = C = ẑE = z2 5
C 4 (59)
k p̄end C p̄start k
z3
The direction of earth’s magnetic field was measured using a magnetic compass con-
structed from a 1 m carbon-fibre rod with neodymium magnets fixed to each end; polarising
each end to either magnetic north or south. The compass was hung from a 1 m length
of cotton thread and left to come to rest. Optical markers were fixed to either end of the
rod as well as to an arbitrary position along the length, but o↵-set form the rod’s axis, to
beak the rotational symmetry optical marker constellation. Figure 60 shows an annotated
photograph of the compass where C cstart and C cend define the position of the optical markers

16
 C
pend

C
pstart

Figure 5: Photograph of the pendulum and optical markers used to measure the direction
of gravity in the camera frame

in the camera frame. The mean position of each marker over a period of time define the
direction of the compass in the camera frame, C ĉ; as shown in equation (60).

C
cstart C
cend

Figure 6: Photograph of the magnetic compass and optical markers used to measure the
direction of the earth’s magnetic field in the camera frame

C C
C c̄end c̄start
ĉ = C C c̄
(60)
k c̄end start k

Due to measurement errors and the imbalance of the suspended magnetic compass, C ĉ
cannot be assumed to be orthogonal to the direction of gravity defined by C p̂ and so cannot
be used to directly define the earth x axis. The non-orthogonal component of C ĉ can be
computed as the vector projection of C ĉ on C p̂. This can then be removed from C ĉ to define

17
the the earth x axis direction in the camera frame, C xE ; as shown in equation (61). Once
normalised, this defines the earth x axis; C x̂E , as shown in equation (62).
C
C ĉ · C p̂ C
xE = C ĉ p̂ (61)
kC p̂k2
2 3
C x1
xE
C
x̂E = C = x2 5
4 (62)
k xE k
x3
The earth frame y axis in the camera frame, C ŷE , can be calculated as the vector that is
orthogonal to both C x̂E and C ẑE and so be defined by equation (63) where signs are chosen
to satisfy the relative axis direction direction convention. The alignment of the earth frame
may be defined as the rotation matrix C E R, constructed from
C
x̂E , C ŷE , and C ẑE . The
quaternion representation, C E q̂, can then be extracted form this using Bar-Itzhack’s method
[48].
2 p 3
±p1 x21 z12
C
ŷE = 4±p1 x22 z22 5 , C x̂E · C ŷE = 0, C ẑE · C ŷE = 0 (63)
2 2
± 1 x3 z 3
C
⇥C ⇤
ER = x̂E C ŷE C ẑE (64)
In order to find the alignment SM q̂ it is assumed that the static error of the orientation
filter’s Kalman-based algorithm is mean zero. The mean algorithm’s output, SE q̄ˆKalman , was
computed for the measurement platform held stationary for a period of approximate 10
seconds. This was used with the alignment E C
C q̂ and optical measurement M q̂ to define the
alignment of the measurement platform in the sensor frame, SM q̂, as equation (65).

S
M q̂ =C E Sˆ
M q̂ ⌦ C q̂ ⌦ E q̄Kalman (65)

4.4 Experimental proceedure

The optical measurement data and raw orientation sensor data were logged simultaneously.
The raw orientation sensor data was then processed by the accompanying software to provide
the calibrated sensor data and Kalman-based algorithm output. This data then synchronised
the optical measurement data, with the optical measurement data interpolated to match the
higher sampling rate of the orientation sensor data. The calibrated sensor data was then
processed through both the IMU and MARG implementations of the proposed orientation
filter and calibrated orientation measurements were extracted from the optical measurement
data using the methods described in sections 4.2 and 4.3.
The proposed filter’s gain was set to 0.033 for the IMU implementation and 0.041 for
the MARG implementation. Trials summarised in section 5.3, found these values to provide
optimal performance. However, an initial value of 2.5 was used for the first 10 seconds of
any experiment to ensure the convergence of algorithm states from initial conditions. The
gain ⇣, applicable to the MARG implementation of the proposed filer, was set to 0 as the
calibrated orientation sensor data was not subject to gyroscope bias drift.

18
 Data was obtained for a sequence of rotations preformed by hand. The measurement
platform was initially held stationary for 20 to 30 seconds to allow time for algorithm states
to converge to steady-state values. The platform was then rotated 90 around its x axis,
then 180 in the opposite direction, and then 90 to bring the platform back to the starting
position. The platform was held stationary for 3 to 5 seconds between each rotation. This
sequence was then repeated around the y and then z axes. The peak angular rate measured
during each rotation was between 110 /s and 190 /s. The experiment was repeated 8 times
to compile a dataset representative of system performance.

5 Results
It is common [24, 26, 18, 19, 20, 21] to quantify orientation sensor performance as the static
and dynamic RMS (Root-Mean-Square) errors in the decoupled Euler parameters describing
the pitch, roll, and heading components of an orientation. Pitch, , roll, ✓ and heading,
correspond to rotations around the sensor frame x, y, and z axis respectively. An Euler angle
representation has the advantage that the decoupled angles may be more easily interpreted
or visualised. The disadvantage of an Euler representation is that it fails to described the
coupling between each of the parameters and will subject to large and erratic errors if the
Euler angle sequence reaches a singularity.
Euler parameters were computed directly from quaternion data using equations (7), (8)
and (9). A total of 4 sets of of Euler parameters were computed, corresponding to the
calibrated optical measurements of orientation, the Kalman-based algorithm estimated ori-
entation and the proposed filter estimates orientation for both the MARG and IMU imple-
mentations. The errors of estimated Euler parameters, ✏ , ✓✏ and ✏ , were computed as the
di↵erence between the Euler parameters of the calibrated optical measurements and those
of each of the corresponding estimated values.

5.1 Typical results

Figures 7, 8 and 9 show results typical of the 8 experiments for both the Kalman-based
algorithm and the proposed filter MARG implementation. In each figure, the 3 traces of the
upper plot represent the optically measured angle, the Kalman-based algorithm estimated
angle, and the proposed filter estimated angle. The 2 traces of the lower plot represent the
calculated error in each of the estimated angles.

19
 Measured and estimated angle φ
100
Measured
50 Kalman

φ (degrees)
Proposed filter

−50

−100
22 24 26 28 30 32 34 36 38 40 42 44

Error
3
Kalman
2
Proposed filter
φε (degrees)

1
0
−1
−2
−3
22 24 26 28 30 32 34 36 38 40 42 44
time (seconds)

Figure 7: Typical results for measured and estimated angle (top) and error (bottom)

Measured and estimated angle θ

100
Measured
50 Kalman
θ (degrees)

Proposed filter

−50

−100
42 44 46 48 50 52 54 56 58 60 62

Error
5
Kalman
Proposed filter
θε (degrees)

−5
42 44 46 48 50 52 54 56 58 60 62
time (seconds)

Figure 8: Typical results for measured and estimated angle ✓ (top) and error (bottom)

20
 Measured and estimated angle ψ
100
Measured
50 Kalman

ψ (degrees)
Proposed filter

−50

−100
60 65 70 75 80

Error
5
Kalman
Proposed filter
ψε (degrees)

−5
60 65 70 75 80
time (seconds)

Figure 9: Typical results for measured and estimated angle (top) and error (bottom)

5.2 Static and dynamic performance

The static and dynamic RMS values of ✏ , ✓✏ , and ✏ were calculated where a static state
was assumed when the measured corresponding angular rate was < 5 /s, and a dynamic
state when 5 /s. This threshold was chosen to be suitably high enough above the noise
floor of the data. Each RMS value was calculated for the period of time framing only the
rotation sequence of the corresponding Euler parameter; as indicated in figures 7, 8 and
9. This was to prevent errors due to initial convergence or singularities in the Euler angle
representation from corrupting results; that is, when ✓ = ±90. The results are summarised
in table 1. Each value, represents the mean of all 8 experiments. Values of RMS[ ✏ ] were
not computed for the IMU implementation of the proposed filter as the filter cannot, nor is
intended to, compensate for an accumulating error in this parameter.

Euler parameter Kalman-based Proposed filter Proposed filter

algorithm (MARG) (IMU)
RMS[ ✏ ] static 0.789 0.581 0.594
RMS[ ✏ ] dynamic 0.769 0.625 0.623
RMS[✓✏ ] static 0.819 0.502 0.497
RMS[✓✏ ] dynamic 0.847 0.668 0.668
RMS[ ✏ ] static 1.150 1.073 N/A
RMS[ ✏ ] dynamic 1.344 1.110 N/A

Table 1: Static and dynamic RMS error of Kalman-based algorithm and proposed filter IMU
and MARG implementations

Results indicate that the proposed filter achieves higher levels of accuracy than the
Kalman-based algorithm. The manufacturer of orientation sensor specify the typical per-
formance of the Kalman-based algorithm as having a static RMS error of < 0.5 in and

21
✓, and < 1 in ; and a dynamic RMS error of < 2 in , ✓, and [18]. These values
do not conform with those listed in table 1. Other studies [49] have shown that accuracy
may far less than that quoted by the manufacture and that quoted levels of accuracy are
only achieved once the devices is recalibrated. The lower levels of accuracy in the heading,
✏ , are to due characteristics of the sensor’s measurements of the earth’s magnetic field.
The inclination of the earth’s magnetic field during testing was between 65 and 70 to the
horizontal [39]. As a consequence the component of the magnetic flux vector available as
a reference of heading is relatively small. The larger component of the vector serves as an
additional reference for pitch, , and roll, ✓, alongside the reference measurement of gravity;
hence errors in pitch, ✏ , and roll, ✓✏ , may be expected to be less than those in heading, ✏ .
The magnetometer is specified [18] as having a bandwidth of 10 Hz which, relative to the
30 Hz and 40 Hz bandwidth of the accelerometer and gyroscope respectively, suggests an
increased heading error, ✏ , in dynamic conditions.

5.3 Filter gain vs. performance

The results of an investigation into the e↵ect of the filter adjustable parameter, , on the
proposed filter performance are summarised in Figure 10. The static and dynamic perfor-
mance is quantified as the mean of the corresponding RMS values of ✏ and ✓✏ , for the IMU
implementation and ✏ , ✓✏ and ✏ for the MARG implementation. The experimental data
was processed though the separate proposed filter IMU and MARG implantations, using
fixed values of between 0 to 0.5. There is a clear optimal value of for each filter imple-
mentation; high enough to minimises errors due to integral drift but sufficiently low enough
that unnecessary noise is not introduced by large steps of gradient descent iterations.

Effect of β on proposed filter performance (IMU) Effect of β on proposed filter performance (MARG)
2 2
Static performance Static performance
mean[ RMS[φε], RMS[θε], RMS[ψε] ] (degrees)

1.8 Dynamic performance 1.8 Dynamic performance

mean[ RMS[φε], RMS[θε] ] (degrees)

1.6 1.6

1.4 β = 0.033 1.4 β = 0.041

1.2 1.2

1 1

0.8 0.8

0.6 0.6

0.4 0.4

0.2 0.2

0 0
0 0.1 0.2 0.3 0.4 0 0.1 0.2 0.3 0.4
β β

Figure 10: The e↵ect of the adjustable parameter, , on the performance of the proposed
filter IMU (left) and MARG (right) implementations

22
5.4 Sampling rate vs. performance
The results of an investigation into the e↵ect of sampling rate on filter performance is sum-
marised in Figure11. The experimental data was processed though the separate proposed
filter IMU and MARG implantations, using the previously defined, optimal values where
the experimental data was decimated to simulate sampling rates between 1Hz and 512 Hz.
It can be seen from Figure11 that the propsoed filter acheives similar levels of performance
at 50 Hz as at 512 Hz. Both filter implemntation are able to acheive a static error < 2
and dynamic error < 7 while sampling at 10 Hz. This level of accracuy may be sufficent be
e↵ective for applications such as human motion capture.

Effect of sampling rate on proposed filter performance (IMU) Effect of sampling rate on proposed filter performance (MARG)
30 30
Static performance Static performance

mean[ RMS[φε], RMS[θε], RMS[ψε] ] (degrees)

Dynamic performance Dynamic performance
25 25
mean[ RMS[φε], RMS[θε] ] (degrees)

20 20

15 15

10 10

5 5

0 0
0 1 2 0 1 2
10 10 10 10 10 10
Sampling rate (Hz) Sampling rate (Hz)

Figure 11: The e↵ect of sampling rate on the performance of the proposed filter IMU (left)
and MARG (right) implementations

5.5 Gyroscope bias drift

The calibrated gyrscope data used by the propsoed filter did not contain any bias errors. To
investiagate the propsoed filter’s ability to compensate for bias drift, errors were artifically
introduced to the 8 experimental datasets. A constant drift of 0.2 /s/s was introduced to
the gyroscope x axis measuremnts, !x and a constant bias error of 0.2 /s as added to the
gyroscope x axis measuremnts, !y . The filter gain ⇣ was set to 0 for the first 10 seconds of
each experiment whilst the filter states converge from intial conditions. After this a value of
0.015 was used corrisponding a maximum convergene rate of the esitmated gyroscope bias
of 1 /s/s.
Figure 12 shows resutls typical of teh 8 experimetns, showing the gyroscope x and y axis
bias estimated by the filter, plotted against the actual. From these results, the fitler can
be seen to scuse↵uly esitmate teh gyroscope biases witht eh rate limtied rate of convergese.
The filter performance under these conditions was constent with that described in Table 1.

23
 Gyroscope bias drift compensation
12
ωx,b (actual)
10 ωx,b (estimated)

8 ωy,b (actual)
ωy,b (estimated)
6
degrees/s

−2

−4
10 15 20 25 30 35 40 45 50 55 60
seconds

Figure 12: Filter tracking of gyroscope bias drift

6 Discussion
The derivation of the filter initially assumed that the accelerometer and magnetometer would
only measure gravity and the earth’s magnetic field. In practise, accelerations due to motion
will result in an erroneous observed direction of gravity and so a potentially corrupt the
estimated attitude and local magnetic distortions may still corrupt the estimated heading.
In many applications it can be assumed that accelerations due to motion and local magnetic
distortions are present for only short periods of time. Therefore the magnitude of the filter
gain may be chosen low enough that the divergence caused by the erroneous gravitational
and magnetic field observations is reduced to an acceptable level over the period. The
minimum acceptable value of is limited by the gyroscope measurement error. In many
applications it may be beneficial to use dynamic values of gains and ⇣. This will the
influence accelerometers and magnetometers have on the estimated orientation to be reduced
during potential problematic periods; for example, when large accelerations are detected.
The use of large gains during the filter initialisation may also improve the filter’s convergence
from initial conditions. For example, it was found that a and gain of 10 enabled the
filter states to converge within 5 seconds when initiated with a gyroscope bias error of 1000
deg/s in each axis.
The structure of the filter implantation for a MARG sensor array is similar to that
proposed by Bachman et al [34]. Both filters estimate the gyroscope measurement error as
the gradient of a error surface created by the magnetometer and accelerometer measurements.
Bachman’s filter computes this using a Gauss-Newton approach which requires numerical
di↵erentiation and a matrix inversion. The filter proposed in this report uses an analytical
derivation of the Jacobian and operates on a normalised gradient of the error surface. As
a result, the filter proposed in this paper provides a substantial reduction in computational
load and enables the derivation of an optimal filter gain based on system characteristics.
Appendix A and B describe the filter IMU and MARG implementations respectively,
each implemented in C where the code has been optimised to reduce the number of arith-
metic operations at the expense of data memory. The IMU implementation requires 108
arithmetic operations each filter update, and the MARG implementation requires 277 arith-

24
metic operations per filter update; which includes the magnetic distortion and gyroscope
drift compensation. These low computational requirements make the algorithm accessible
for low power embedded systems systems enabling the use of low cost, low power hardware.
The experimental procedure used to evaluated the filter performances has a number
of limitations; the filter performance was not evaluated for simultaneous rotations around
more than one rotational axis and rotational velocities were limited to short periods and
in magnitude. These limitations were necessary so that repeatable and quantifiable and
practical.

7 Conclusions
In this work, we have introduced a novel orientation filter, applicable to both IMUs and
MARG sensor arrays, that significantly ameliorates the computational load and parameter
tuning burdens associated with conventional Kalman-based approaches. The filter is based
on a Newton optimization using an analytic formulation of the gradient that is derived from
a quaternion representation of motion. Novel aspects of the filter include:

• Analytic derivation of the Jacobian matrix, which eliminates significant computational

load, allowing implementation at lower sampling frequencies and onto lower power,
smaller platforms.

• The need to tune only one or two filter gains ( and ⇣)), defined by the gyroscope
measurement error. Least squares curve fitting and complex tuning processes are
therefore eliminated.

The filter derivation, magnetic distortion anf gyroscope bias drift compensation, and
experimental testing have been detailed. Empirical testing and benchmarking has shown
that the filter performs as well as a high quality commercial Kalman-based system, even
with a full order of magnitude in reduction of sampling rate. The filter is both simple to
implement and simple to tune. The implications of the low computational load and ability
to operate at low sampling rates open a very wide range of new opportunities for the use of
IMU and MARG sensor arrays in real-time applications. Applications where limited power
or processing resources may be available are particularly well suited for the new filter. The
filter also has great potential to alleviate computational load for applications that demand
extremely high sampling rates.

References
[1] Mark Euston, Paul Coote, Robert Mahony, Jonghyuk Kim, and Tarek Hamel. A com-
plementary filter for attitude estimation of a fixed-wing uav with a low-cost imu. In 6th
International Conference on Field and Service Robotics, July 2007.

[2] Sung Kyung Hong. Fuzzy logic based closed-loop strapdown attitude system for un-
manned aerial vehicle (uav). Sensors and Actuators A: Physical, 107(2):109 – 118,
2003.

25
 [3] A. Kallapur, I. Petersen, and S. Anavatti. A robust gyroless attitude estimation scheme
for a small fixed-wing unmanned aerial vehicle. pages 666 –671, aug. 2009.

[4] B. Barshan and H. F. Durrant-Whyte. Inertial navigation systems for mobile robots.
11(3):328–342, June 1995.

[5] L. Ojeda and J. Borenstein. Flexnav: fuzzy logic expert rule-based position estimation
for mobile robots on rugged terrain. In Proc. IEEE International Conference on Robotics
and Automation ICRA ’02, volume 1, pages 317–322, May 11–15, 2002.

[6] David H. Titterton and John L. Weston. Strapdown inertial navigation technology. The
Institution of Electrical Engineers, 2004.

[7] S. Beauregard. Omnidirectional pedestrian navigation for first responders. In Proc. 4th
Workshop on Positioning, Navigation and Communication WPNC ’07, pages 33–36,
March 22–22, 2007.

[8] H. J. Luinge and P. H. Veltink. Inclination measurement of human movement using a

3-d accelerometer with autocalibration. 12(1):112–121, March 2004.

[9] Huiyu Zhou and Huosheng Hu. Human motion tracking for rehabilitation–a survey.
volume 3, pages 1 – 18, 2008.

[10] E. A. Heinz, K. S. Kunze, M. Gruber, D. Bannach, and P. Lukowicz. Using wearable

sensors for real-time recognition tasks in games of martial arts - an initial experiment.
In Proc. IEEE Symposium on Computational Intelligence and Games, pages 98–102,
May 22–24, 2006.

[11] J. E. Bortz. A new mathematical formulation for strapdown inertial navigation. (1):61–
66, January 1971.

[12] M B Ignagni. Optimal strapdown attitude integration algorithms. In Guidance, Control,

and Dynamics, volume 13, pages 363–369, 1990.

[13] N. Yazdi, F. Ayazi, and K. Najafi. Micromachined inertial sensors. 86(8):1640–1659,

August 1998.

[14] R. E. Kalman. A new approach to linear filtering and prediction problems. Journal of
Basic Engineering, 82:35–45, 1960.

[15] E. Foxlin. Inertial head-tracker sensor fusion by a complementary separate-bias kalman

filter. In Proc. Virtual Reality Annual International Symposium the IEEE 1996, pages
185–194,267, March 30–April 3, 1996.

[16] H. J. Luinge, P. H. Veltink, and C. T. M. Baten. Estimation of orientation with gyro-

scopes and accelerometers. In Proc. First Joint [Engineering in Medicine and Biology
21st Annual Conf. and the 1999 Annual Fall Meeting of the Biomedical Engineering
Soc.] BMES/EMBS Conference, volume 2, page 844, October 13–16, 1999.

26
[17] J. L. Marins, Xiaoping Yun, E. R. Bachmann, R. B. McGhee, and M. J. Zyda. An ex-
tended kalman filter for quaternion-based orientation estimation using marg sensors. In
Proc. IEEE/RSJ International Conference on Intelligent Robots and Systems, volume 4,
pages 2003–2011, October 29–November 3, 2001.

[18] Xsens Technologies B.V. MTi and MTx User Manual and Technical Documentation.
Pantheon 6a, 7521 PR Enschede, The Netherlands, May 2009.

[19] MicroStrain Inc. 3DM-GX3 -25 Miniature Attutude Heading Reference Sensor. 459
Hurricane Lane, Suite 102, Williston, VT 05495 USA, 1.04 edition, 2009.

[20] VectorNav Technologies, LLC. VN -100 User Manual. College Station, TX 77840 USA,
preliminary edition, 2009.

[21] InterSense, Inc. InertiaCube2+ Manual. 36 Crosby Drive, Suite 150, Bedford, MA
01730, USA, 1.0 edition, 2008.

[22] PNI sensor corporation. Spacepoint Fusion. 133 Aviation Blvd, Suite 101, Santa Rosa,
CA 95403-1084 USA.

[23] Crossbow Technology, Inc. AHRS400 Series Users Manual. 4145 N. First Street, San
Jose, CA 95134, rev. c edition, February 2007.

[24] A. M. Sabatini. Quaternion-based extended kalman filter for determining orientation

by inertial and magnetic sensing. 53(7):1346–1356, July 2006.

[25] H. J. Luinge and P. H. Veltink. Measuring orientation of human body segments us-
ing miniature gyroscopes and accelerometers. Medical and Biological Engineering and
Computing, 43(2):273–282, April 2006.

[26] David Jurman, Marko Jankovec, Roman Kamnik, and Marko Topic. Calibration and
data fusion solution for the miniature attitude and heading reference system. Sensors
and Actuators A: Physical, 138(2):411–420, August 2007.

[27] Markus Haid and Jan Breitenbach. Low cost inertial orientation tracking with kalman
filter. Applied Mathematics and Computation, 153(4):567–575, June 2004.

[28] D. Roetenberg, H. J. Luinge, C. T. M. Baten, and P. H. Veltink. Compensation of

magnetic disturbances improves inertial and magnetic sensing of human body segment
orientation. 13(3):395–405, September 2005.

[29] John L. Crassidis and Landis F. Markley. Unscented filtering for spacecraft attitude
estimation. Journal of guidance, control, and dynamics, 26(4):536–542, 2003.

[30] D. Gebre-Egziabher, R.C. Hayward, and J.D. Powell. Design of multi-sensor atti-
tude determination systems. Aerospace and Electronic Systems, IEEE Transactions
on, 40(2):627 – 649, april 2004.

27
[31] N.H.Q. Phuong, H.-J. Kang, Y.-S. Suh, and Y.-S. Ro. A DCM based orientation esti-
mation algorithm with an inertial measurement unit and a magnetic compass. Journal
of Universal Computer Science, 15(4):859–876, 2009.

[32] Daniel Choukroun. Novel methods for attitude determination using vector observations.
PhD thesis, Israel Institute of Technology, May 2003.

[33] R. A. Hyde, L. P. Ketteringham, S. A. Neild, and R. J. S. Jones. Estimation of upper-

limb orientation based on accelerometer and gyroscope measurements. 55(2):746–754,
February 2008.

[34] E. R. Bachmann, I. Duman, U. Y. Usta, R. B. McGhee, X. P. Yun, and M. J. Zyda.

Orientation tracking for humans and robots using inertial sensors. In Proc. IEEE Inter-
national Symposium on Computational Intelligence in Robotics and Automation CIRA
’99, pages 187–194, November 8–9, 1999.

[35] R. Mahony, T. Hamel, and J.-M. Pflimlin. Nonlinear complementary filters on the
special orthogonal group. Automatic Control, IEEE Transactions on, 53(5):1203 –1218,
june 2008.

[36] J. B. Kuipers. Quaternions and Rotation Sequences: A Primer with Applications to

Orbits, Aerospace and Virtual Reality. Princeton University Press, 1999.

[37] John J. Craig. Introduction to Robotics Mechanics and Control. Pearson Education
International, 2005.

[38] David R. Pratt Robert B. McGhee Joseph M. Cooke, Michael J. Zyda. Npsnet: Flight
simulation dynamic modelling using quaternions. Presence, 1:404–420, 1994.

[39] John Arthur Jacobs. The earth’s core, volume 37 of International geophysics series.
Academic Press, 2 edition, 1987.

[40] E. R. Bachmann, Xiaoping Yun, and C. W. Peterson. An investigation of the e↵ects of

magnetic variations on inertial/magnetic orientation sensors. In Proc. IEEE Interna-
tional Conference on Robotics and Automation ICRA ’04, volume 2, pages 1115–1122,
April 2004.

[41] C.T.M. Baten F.C.T. van der Helm W.H.K. de Vries, H.E.J. Veeger. Magnetic distortion
in motion labs, implications for validating inertial magnetic sensors. Gait & Posture,
29(4):535–541, 2009.

[42] Speake & Co Limited. Autocalibration algorithms for FGM type sensors. Application
note.

[43] Michael J. Caruso. Applications of Magnetoresistive Sensors in Navigation Systems.

Honeywell Inc., Solid State Electronics Center, Honeywell Inc. 12001 State Highway 55,
Plymouth, MN 55441.

28
[44] J. F. Vasconcelos, G. Elkaim, C. Silvestre, P. Oliveira, and B. Cardeira. A geomet-
ric approach to strapdown magnetometer calibration in sensor frame. In Navigation,
Guidance and Control of Underwater Vehicles, volume 2, 2008.

[45] Demoz Gebre-Egziabher, Gabriel H. Elkaim, J. David Powell, and Bradford W. Parkin-
son. Calibration of strapdown magnetometers in magnetic field domain. Journal of
Aerospace Engineering, 19(2):87–102, 2006.

[46] Vicon Motion Systems Limited. Vicon MX Hardware. 5419 McConnell Avenue, Los
Angeles, CA 90066, USA, 1.6 edition, 2004.

[47] Vicon Motion Systems Limited. Vicon Nexus Product Guide - Foundation Notes. 5419
McConnell Avenue, Los Angeles, CA 90066, USA, 1.2 edition, November 2007.

[48] Itzhack Y Bar-Itzhack. New method for extracting the quaternion from a rotation
matrix. AIAA Journal of Guidance, Control and Dynamics, 23(6):10851087, Nov.Dec
2000. (Engineering Note).

[49] M. A. Brodie, A. Walmsley, and W. Page. The static accuracy and calibration of
inertial measurement units for 3D orientation. Computer Methods in Biomechanics and
Biomedical Engineering, 11(6):641–648, December 2008.

A IMU filter implementation optimised in C

The following source code is an implementation of the orientation filter for an IMU, in C.
The code has been optimised minimise the required number of arithmetic operations at the
expense of data memory. Each filter update requires 109 scalar arithmetic operations (18
additions, 20 subtracts, 57 multiplications, 11 divisions and 3 square roots). The implemen-
tation requires 40 bytes of data memory for global variables and 100 bytes of data memory
for local variables during each filter update function call.

// Math library required for ‘sqrt’

#include <math.h>

// System constants
#define deltat 0.001f // sampling period in seconds (shown as 1 ms)
#define gyroMeasError 3.14159265358979f * (5.0f / 180.0f) // gyroscope measurement error in rad/s (shown as 5 deg/s)
#define beta sqrt(3.0f / 4.0f) * gyroMeasError // compute beta

// Global system variables

float a_x, a_y, a_z; // accelerometer measurements
float w_x, w_y, w_z; // gyroscope measurements in rad/s
float SEq_1 = 1.0f, SEq_2 = 0.0f, SEq_3 = 0.0f, SEq_4 = 0.0f; // estimated orientation quaternion elements with initial conditions

void filterUpdate(float w_x, float w_y, float w_z, float a_x, float a_y, float a_z)
{
// Local system variables
float norm; // vector norm
float SEqDot_omega_1, SEqDot_omega_2, SEqDot_omega_3, SEqDot_omega_4; // quaternion derrivative from gyroscopes elements
float f_1, f_2, f_3; // objective function elements
float J_11or24, J_12or23, J_13or22, J_14or21, J_32, J_33; // objective function Jacobian elements
float SEqHatDot_1, SEqHatDot_2, SEqHatDot_3, SEqHatDot_4; // estimated direction of the gyroscope error

// Axulirary variables to avoid reapeated calcualtions

float halfSEq_1 = 0.5f * SEq_1;
float halfSEq_2 = 0.5f * SEq_2;
float halfSEq_3 = 0.5f * SEq_3;
float halfSEq_4 = 0.5f * SEq_4;
float twoSEq_1 = 2.0f * SEq_1;
float twoSEq_2 = 2.0f * SEq_2;
float twoSEq_3 = 2.0f * SEq_3;

29
 // Normalise the accelerometer measurement
norm = sqrt(a_x * a_x + a_y * a_y + a_z * a_z);
a_x /= norm;
a_y /= norm;
a_z /= norm;

// Compute the objective function and Jacobian

f_1 = twoSEq_2 * SEq_4 - twoSEq_1 * SEq_3 - a_x;
f_2 = twoSEq_1 * SEq_2 + twoSEq_3 * SEq_4 - a_y;
f_3 = 1.0f - twoSEq_2 * SEq_2 - twoSEq_3 * SEq_3 - a_z;
J_11or24 = twoSEq_3; // J_11 negated in matrix multiplication
J_12or23 = 2.0f * SEq_4;
J_13or22 = twoSEq_1; // J_12 negated in matrix multiplication
J_14or21 = twoSEq_2;
J_32 = 2.0f * J_14or21; // negated in matrix multiplication
J_33 = 2.0f * J_11or24; // negated in matrix multiplication

// Compute the gradient (matrix multiplication)

SEqHatDot_1 = J_14or21 * f_2 - J_11or24 * f_1;
SEqHatDot_2 = J_12or23 * f_1 + J_13or22 * f_2 - J_32 * f_3;
SEqHatDot_3 = J_12or23 * f_2 - J_33 * f_3 - J_13or22 * f_1;
SEqHatDot_4 = J_14or21 * f_1 + J_11or24 * f_2;

// Normalise the gradient

norm = sqrt(SEqHatDot_1 * SEqHatDot_1 + SEqHatDot_2 * SEqHatDot_2 + SEqHatDot_3 * SEqHatDot_3 + SEqHatDot_4 * SEqHatDot_4);
SEqHatDot_1 /= norm;
SEqHatDot_2 /= norm;
SEqHatDot_3 /= norm;
SEqHatDot_4 /= norm;

// Compute the quaternion derrivative measured by gyroscopes

SEqDot_omega_1 = -halfSEq_2 * w_x - halfSEq_3 * w_y - halfSEq_4 * w_z;
SEqDot_omega_2 = halfSEq_1 * w_x + halfSEq_3 * w_z - halfSEq_4 * w_y;
SEqDot_omega_3 = halfSEq_1 * w_y - halfSEq_2 * w_z + halfSEq_4 * w_x;
SEqDot_omega_4 = halfSEq_1 * w_z + halfSEq_2 * w_y - halfSEq_3 * w_x;

// Compute then integrate the estimated quaternion derrivative

SEq_1 += (SEqDot_omega_1 - (beta * SEqHatDot_1)) * deltat;
SEq_2 += (SEqDot_omega_2 - (beta * SEqHatDot_2)) * deltat;
SEq_3 += (SEqDot_omega_3 - (beta * SEqHatDot_3)) * deltat;
SEq_4 += (SEqDot_omega_4 - (beta * SEqHatDot_4)) * deltat;

// Normalise quaternion
norm = sqrt(SEq_1 * SEq_1 + SEq_2 * SEq_2 + SEq_3 * SEq_3 + SEq_4 * SEq_4);
SEq_1 /= norm;
SEq_2 /= norm;
SEq_3 /= norm;
SEq_4 /= norm;
}

B MARG filter implementation optimised in C

The following source code is an implementation of the orientation filter for a MARG sensor
array including magnetic distortion and gyroscope drift compensation, in C. The code has
been optimised minimise the required number of arithmetic operations at the expense of
data memory. Each filter update requires 277 scalar arithmetic operations (51 additions,
57 subtracts, 155 multiplications, 14 divisions and 5 square roots). The implementation
requires 72 bytes of data memory for global variables and 260 bytes of data memory for local
variables during each filter update function call.

// Math library required for ‘sqrt’

#include <math.h>

// System constants
#define deltat 0.001f // sampling period in seconds (shown as 1 ms)
#define gyroMeasError 3.14159265358979 * (5.0f / 180.0f) // gyroscope measurement error in rad/s (shown as 5 deg/s)
#define gyroMeasDrift 3.14159265358979 * (0.2f / 180.0f) // gyroscope measurement error in rad/s/s (shown as 0.2f deg/s/s)
#define beta sqrt(3.0f / 4.0f) * gyroMeasError // compute beta
#define zeta sqrt(3.0f / 4.0f) * gyroMeasDrift // compute zeta

// Global system variables

float a_x, a_y, a_z; // accelerometer measurements
float w_x, w_y, w_z; // gyroscope measurements in rad/s
float m_x, m_y, m_z; // magnetometer measurements
float SEq_1 = 1, SEq_2 = 0, SEq_3 = 0, SEq_4 = 0; // estimated orientation quaternion elements with initial conditions
float b_x = 1, b_z = 0; // reference direction of flux in earth frame

30
float w_bx = 0, w_by = 0, w_bz = 0; // estimate gyroscope biases error

// Function to compute one filter iteration

void filterUpdate(float w_x, float w_y, float w_z, float a_x, float a_y, float a_z, float m_x, float m_y, float m_z)
{
// local system variables
float norm; // vector norm
float SEqDot_omega_1, SEqDot_omega_2, SEqDot_omega_3, SEqDot_omega_4; // quaternion rate from gyroscopes elements
float f_1, f_2, f_3, f_4, f_5, f_6; // objective function elements
float J_11or24, J_12or23, J_13or22, J_14or21, J_32, J_33, // objective function Jacobian elements
J_41, J_42, J_43, J_44, J_51, J_52, J_53, J_54, J_61, J_62, J_63, J_64; //
float SEqHatDot_1, SEqHatDot_2, SEqHatDot_3, SEqHatDot_4; // estimated direction of the gyroscope error
float w_err_x, w_err_y, w_err_z; // estimated direction of the gyroscope error (angular)
float h_x, h_y, h_z; // computed flux in the earth frame

// axulirary variables to avoid reapeated calcualtions

float halfSEq_1 = 0.5f * SEq_1;
float halfSEq_2 = 0.5f * SEq_2;
float halfSEq_3 = 0.5f * SEq_3;
float halfSEq_4 = 0.5f * SEq_4;
float twoSEq_1 = 2.0f * SEq_1;
float twoSEq_2 = 2.0f * SEq_2;
float twoSEq_3 = 2.0f * SEq_3;
float twoSEq_4 = 2.0f * SEq_4;
float twob_x = 2.0f * b_x;
float twob_z = 2.0f * b_z;
float twob_xSEq_1 = 2.0f * b_x * SEq_1;
float twob_xSEq_2 = 2.0f * b_x * SEq_2;
float twob_xSEq_3 = 2.0f * b_x * SEq_3;
float twob_xSEq_4 = 2.0f * b_x * SEq_4;
float twob_zSEq_1 = 2.0f * b_z * SEq_1;
float twob_zSEq_2 = 2.0f * b_z * SEq_2;
float twob_zSEq_3 = 2.0f * b_z * SEq_3;
float twob_zSEq_4 = 2.0f * b_z * SEq_4;
float SEq_1SEq_2;
float SEq_1SEq_3 = SEq_1 * SEq_3;
float SEq_1SEq_4;
float SEq_2SEq_3;
float SEq_2SEq_4 = SEq_2 * SEq_4;
float SEq_3SEq_4;
float twom_x = 2.0f * m_x;
float twom_y = 2.0f * m_y;
float twom_z = 2.0f * m_z;

// normalise the accelerometer measurement

norm = sqrt(a_x * a_x + a_y * a_y + a_z * a_z);
a_x /= norm;
a_y /= norm;
a_z /= norm;

// normalise the magnetometer measurement

norm = sqrt(m_x * m_x + m_y * m_y + m_z * m_z);
m_x /= norm;
m_y /= norm;
m_z /= norm;

// compute the objective function and Jacobian

f_1 = twoSEq_2 * SEq_4 - twoSEq_1 * SEq_3 - a_x;
f_2 = twoSEq_1 * SEq_2 + twoSEq_3 * SEq_4 - a_y;
f_3 = 1.0f - twoSEq_2 * SEq_2 - twoSEq_3 * SEq_3 - a_z;
f_4 = twob_x * (0.5f - SEq_3 * SEq_3 - SEq_4 * SEq_4) + twob_z * (SEq_2SEq_4 - SEq_1SEq_3) - m_x;
f_5 = twob_x * (SEq_2 * SEq_3 - SEq_1 * SEq_4) + twob_z * (SEq_1 * SEq_2 + SEq_3 * SEq_4) - m_y;
f_6 = twob_x * (SEq_1SEq_3 + SEq_2SEq_4) + twob_z * (0.5f - SEq_2 * SEq_2 - SEq_3 * SEq_3) - m_z;
J_11or24 = twoSEq_3; // J_11 negated in matrix multiplication
J_12or23 = 2.0f * SEq_4;
J_13or22 = twoSEq_1; // J_12 negated in matrix multiplication
J_14or21 = twoSEq_2;
J_32 = 2.0f * J_14or21; // negated in matrix multiplication
J_33 = 2.0f * J_11or24; // negated in matrix multiplication
J_41 = twob_zSEq_3; // negated in matrix multiplication
J_42 = twob_zSEq_4;
J_43 = 2.0f * twob_xSEq_3 + twob_zSEq_1; // negated in matrix multiplication
J_44 = 2.0f * twob_xSEq_4 - twob_zSEq_2; // negated in matrix multiplication
J_51 = twob_xSEq_4 - twob_zSEq_2; // negated in matrix multiplication
J_52 = twob_xSEq_3 + twob_zSEq_1;
J_53 = twob_xSEq_2 + twob_zSEq_4;
J_54 = twob_xSEq_1 - twob_zSEq_3; // negated in matrix multiplication
J_61 = twob_xSEq_3;
J_62 = twob_xSEq_4 - 2.0f * twob_zSEq_2;
J_63 = twob_xSEq_1 - 2.0f * twob_zSEq_3;
J_64 = twob_xSEq_2;

// compute the gradient (matrix multiplication)

SEqHatDot_1 = J_14or21 * f_2 - J_11or24 * f_1 - J_41 * f_4 - J_51 * f_5 + J_61 * f_6;
SEqHatDot_2 = J_12or23 * f_1 + J_13or22 * f_2 - J_32 * f_3 + J_42 * f_4 + J_52 * f_5 + J_62 * f_6;
SEqHatDot_3 = J_12or23 * f_2 - J_33 * f_3 - J_13or22 * f_1 - J_43 * f_4 + J_53 * f_5 + J_63 * f_6;
SEqHatDot_4 = J_14or21 * f_1 + J_11or24 * f_2 - J_44 * f_4 - J_54 * f_5 + J_64 * f_6;

// normalise the gradient to estimate direction of the gyroscope error

norm = sqrt(SEqHatDot_1 * SEqHatDot_1 + SEqHatDot_2 * SEqHatDot_2 + SEqHatDot_3 * SEqHatDot_3 + SEqHatDot_4 * SEqHatDot_4);
SEqHatDot_1 = SEqHatDot_1 / norm;
SEqHatDot_2 = SEqHatDot_2 / norm;

31
 SEqHatDot_3 = SEqHatDot_3 / norm;
SEqHatDot_4 = SEqHatDot_4 / norm;

// compute angular estimated direction of the gyroscope error

w_err_x = twoSEq_1 * SEqHatDot_2 - twoSEq_2 * SEqHatDot_1 - twoSEq_3 * SEqHatDot_4 + twoSEq_4 * SEqHatDot_3;
w_err_y = twoSEq_1 * SEqHatDot_3 + twoSEq_2 * SEqHatDot_4 - twoSEq_3 * SEqHatDot_1 - twoSEq_4 * SEqHatDot_2;
w_err_z = twoSEq_1 * SEqHatDot_4 - twoSEq_2 * SEqHatDot_3 + twoSEq_3 * SEqHatDot_2 - twoSEq_4 * SEqHatDot_1;

// compute and remove the gyroscope baises

w_bx += w_err_x * deltat * zeta;
w_by += w_err_y * deltat * zeta;
w_bz += w_err_z * deltat * zeta;
w_x -= w_bx;
w_y -= w_by;
w_z -= w_bz;

// compute the quaternion rate measured by gyroscopes

SEqDot_omega_1 = -halfSEq_2 * w_x - halfSEq_3 * w_y - halfSEq_4 * w_z;
SEqDot_omega_2 = halfSEq_1 * w_x + halfSEq_3 * w_z - halfSEq_4 * w_y;
SEqDot_omega_3 = halfSEq_1 * w_y - halfSEq_2 * w_z + halfSEq_4 * w_x;
SEqDot_omega_4 = halfSEq_1 * w_z + halfSEq_2 * w_y - halfSEq_3 * w_x;

// compute then integrate the estimated quaternion rate

SEq_1 += (SEqDot_omega_1 - (beta * SEqHatDot_1)) * deltat;
SEq_2 += (SEqDot_omega_2 - (beta * SEqHatDot_2)) * deltat;
SEq_3 += (SEqDot_omega_3 - (beta * SEqHatDot_3)) * deltat;
SEq_4 += (SEqDot_omega_4 - (beta * SEqHatDot_4)) * deltat;

// normalise quaternion
norm = sqrt(SEq_1 * SEq_1 + SEq_2 * SEq_2 + SEq_3 * SEq_3 + SEq_4 * SEq_4);
SEq_1 /= norm;
SEq_2 /= norm;
SEq_3 /= norm;
SEq_4 /= norm;

// compute flux in the earth frame

SEq_1SEq_2 = SEq_1 * SEq_2; // recompute axulirary variables
SEq_1SEq_3 = SEq_1 * SEq_3;
SEq_1SEq_4 = SEq_1 * SEq_4;
SEq_3SEq_4 = SEq_3 * SEq_4;
SEq_2SEq_3 = SEq_2 * SEq_3;
SEq_2SEq_4 = SEq_2 * SEq_4;
h_x = twom_x * (0.5f - SEq_3 * SEq_3 - SEq_4 * SEq_4) + twom_y * (SEq_2SEq_3 - SEq_1SEq_4) + twom_z * (SEq_2SEq_4 + SEq_1SEq_3);
h_y = twom_x * (SEq_2SEq_3 + SEq_1SEq_4) + twom_y * (0.5f - SEq_2 * SEq_2 - SEq_4 * SEq_4) + twom_z * (SEq_3SEq_4 - SEq_1SEq_2);
h_z = twom_x * (SEq_2SEq_4 - SEq_1SEq_3) + twom_y * (SEq_3SEq_4 + SEq_1SEq_2) + twom_z * (0.5f - SEq_2 * SEq_2 - SEq_3 * SEq_3);

// normalise the flux vector to have only components in the x and z

b_x = sqrt((h_x * h_x) + (h_y * h_y));
b_z = h_z;
}

32