Madgwick PHD Thesis
Madgwick PHD Thesis
Sebastian O. H. Madgwick
Supervised by
Andrew Harrison
March 2014
Declaration
I declare that the work in this dissertation was carried out in accordance with the
requirements of the University’s Regulations and Code of Practice for Research
Degree Programmes and that it has not been submitted for any other academic
award. Except where indicated by specific reference in the text, the work is the
candidate’s own work. Work done in collaboration with, or with the assistance of,
others, is indicated as such. Any views expressed in the dissertation are those of
the author.
SIGNED: .............................................................
DATE: ..........................
i
Declaration
ii
Acknowledgments
iii
Acknowledgments
iv
Abstract
v
Abstract
vi
Publications arising from this
research
Journal papers
1. Sebastian O.H. Madgwick; Andrew Harrison; Paul Sharkey; Ravi
Vaidyanathan; William Harwin. “Measuring motion with kinematically
redundant accelerometer arrays: theory, simulation and implementation”.
Mechatronics, Elsevier, vol. 23, no. 5, pp. 518-529, 2013
Conference papers
1. Sebastian O.H. Madgwick, Andrew J.L. Harrison, Ravi Vaidyanathan.
“Estimation of IMU and MARG orientation using a gradient descent
algorithm”. IEEE International Conference on Rehabilitation Robotics
(ICORR), 2011
2. James Carberry, Graham Hinchly, James Buckerfield, Edward Tayler.
Thomas Burton, Sebastian Madgwick, Ravi Vaidyanathan. “Parametric
Design of an Active Ankle Foot Orthosis with Passive Compliance”. IEEE
International Conference on Computer-Based Medical Systems
3. Thomas Mitchell, Sebastian Madgwick, Imogen Heap. “Musical Interaction
with Hand Posture and Orientation: A Toolbox of Gestural Control
Mechanisms”. Proceedings of the International Conference on New Interfaces
for Musical Expression (NIME2012), 2012
4. Sebastian Madgwick, Thomas Mitchell. “x-OSC: A VersatileWireless I/O
Device For Creative/Music Applications”. Proceedings of Sound and Music
Computing Conference (SMC), 2013
5. Thomas Mitchell, Sebastian Madgwick, Simon Rankin, Geoff Hilton, Andrew
Nix, “Making the Most of Wi-Fi: Optimisations for Wireless Robust Live
Music Performance”. Proceedings of the International Conference on New
Interfaces for Musical Expression (NIME2014)
vii
Publications arising from this research
viii
Table of contents
Declaration i
Acknowledgments iii
Abstract v
Table of contents ix
List of tables xv
Glossary xxiii
1 Introduction 1
1.1 Low-cost MEMS sensors and AHRS applications . . . . . . . . . . . 1
1.1.1 Measuring motion with IMUs . . . . . . . . . . . . . . . . . 2
1.1.2 IMUs platforms . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2 Research objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.3 Contribution to knowledge . . . . . . . . . . . . . . . . . . . . . . . 5
1.4 Thesis structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2 Background 11
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2 MEMS gyroscopes, accelerometers and magnetometers . . . . . . . 11
2.2.1 Accelerometers . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.2.2 Gyroscopes . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.2.3 Magnetometers . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.2.4 Monolithic devices . . . . . . . . . . . . . . . . . . . . . . . 15
2.3 Attitude and Heading Reference System algorithms . . . . . . . . . 15
2.3.1 Complementary filter AHRS algorithms . . . . . . . . . . . . 17
2.4 Non-gyro IMUs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.5 Representation and manipulation of spatial quantities . . . . . . . . 20
2.5.1 Rotation matrices . . . . . . . . . . . . . . . . . . . . . . . . 20
2.5.1.1 Matrix transpose . . . . . . . . . . . . . . . . . . . 22
ix
Table of contents
x
Table of contents
4.8 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
xi
Table of contents
xii
Table of contents
9 Conclusions 211
References 213
xiii
Table of contents
xiv
List of tables
xv
List of tables
xvi
List of figures
xvii
List of figures
xviii
List of figures
xix
List of figures
8.1 x-IMU board alone (left) and enclosed in its plastic housing with
1000 mAh battery (right) . . . . . . . . . . . . . . . . . . . . . . . 161
8.2 x-IMU top and bottom view with labelled key hardware components 162
8.3 x-IMU GUI displaying real-time data with 2D and 3D graphics . . . 164
8.4 x-BIMU with an unattached XBee 802.15.4 module . . . . . . . . . 167
8.5 x-BIMU in plastic housing with XBee 802.15.4 module and 320
mAh battery (left) and Velcro strap designed for human motion
applications (right) . . . . . . . . . . . . . . . . . . . . . . . . . . . 168
8.6 x-BIMU top with labelled key hardware components . . . . . . . . 169
8.7 x-BIMU Terminal displaying real-time data with 2D and 3D graphics171
8.8 Example XBee-style communication modules compatible with
the x-BIMU, for: 802.15.4 (a) (b), Bluetooth (c), ZigBee (d),
Proprietary 900 MHz (e), Wi-Fi (f), Bluetooth 4.0 Low Energy (g)
and wired serial (h). . . . . . . . . . . . . . . . . . . . . . . . . . . 172
8.9 Four 802.15.4 receivers (Digi XSticks) in a USB hub provide
synchronised measurements from x-BIMUs . . . . . . . . . . . . . . 175
8.10 x-OSC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 177
8.11 x-OSC top and bottom view with labelled key hardware components 178
8.12 x-OSC settings viewed on web browser . . . . . . . . . . . . . . . . 180
8.13 Normalised distribution of measured round-trip latency . . . . . . . 183
8.14 Throughput of 1 to 15 x-OSCs sending to a single AP on one channel184
8.15 Throughput of 1 to 15 x-OSCs sending to three Access Points (APs)
on three non-overlapping channels . . . . . . . . . . . . . . . . . . . 185
8.16 SoundGrasp gloves with wrist-mounted microphone and postures
with associated control function. Image source: [170] . . . . . . . . 187
xx
List of figures
8.17 Prototype gloves system diagram (left) and being worn by musician
(right) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 188
8.18 Performance glove parts . . . . . . . . . . . . . . . . . . . . . . . . 189
8.19 Gloves being used to create and manipulate visuals projected on to
a screen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 190
8.20 Measured foot position during walking (3 steps) obtained using an
x-IMU . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 191
8.21 Mechanical model of upper limb (a), placement of x-IMUs (b) and
instrumented anthropomorphic upper limb used for evaluation (c).
Image source: [187] . . . . . . . . . . . . . . . . . . . . . . . . . . . 193
8.22 Digits wrist-worn hardware components. Image source: [188] . . . . 194
8.23 MONSUN II AUV is 60 cm long, has a hull diameter of 10 cm and
weighs <3kg. Image source: [190] . . . . . . . . . . . . . . . . . . . 194
8.24 Key internal hardware complements including an x-IMU. Image
source: [190] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 195
8.25 SMART-E AUV uses three actuated thrusters for omnidirectional
movement and a band of RGB LEDs to provide visual feedback
under water. Image source: [195] . . . . . . . . . . . . . . . . . . . 196
8.26 Internal electronics and communication networks between sensors,
actuators and controllers. Image source: [195] . . . . . . . . . . . . 196
8.27 Dimensional drawings of Tidal Turbulence Moorings deployed in
Admiralty Inlet (left) and Chacao Channel (right). Image source:
[196] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 198
8.28 Instrumented surgical tool uses x-IMU to provide orientation data.
Image source: [196] . . . . . . . . . . . . . . . . . . . . . . . . . . . 199
8.29 Virtual reality game controlled using the instrumented surgical tool.
Image source: [196] . . . . . . . . . . . . . . . . . . . . . . . . . . . 199
8.30 Instrumented equipment comprised an a soccer boot and shin
guard; each with an x-IMU mounted internally. The equipment
is shown here attached the the synchronisation trigger prior to a
data acquisition session. Image source: [198] . . . . . . . . . . . . . 200
8.31 Combined motion capture and display helmet (left) and robot head
(right). Image source: [202] . . . . . . . . . . . . . . . . . . . . . . 202
8.32 Latest generation autonomous car shown fully assembled (left) and
with shell removed (right) . . . . . . . . . . . . . . . . . . . . . . . 203
8.33 Screenshot showing view from camera mounted on car and the
plotted location of the car as estimated by the RatSLAM algorithm 204
8.34 ‘Adult size’ robot (150 cm) robot walking towards football
(Photograph courtesy of Hamidreza Kasaei) . . . . . . . . . . . . . 205
8.35 The x-IMU is being used to collect data from the Army Historic
Aircraft Flight Scout AH Mk 1 to monitor structural usage . . . . . 206
8.36 An x-BIMU mounted on a dancer’s wrist allows her to control the
real-time visualisations of the molecular simulation. Photographer:
Paul Blakemore . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 208
8.37 Changibles analysis software and resultant physical construction.
Image source: [209] . . . . . . . . . . . . . . . . . . . . . . . . . . . 208
xxi
List of figures
8.38 An individual block containing x-OSC for Wi-Fi control and 2 servos
with scissor lift mechanism. Image source: [209] . . . . . . . . . . . 209
xxii
Glossary
3D Three Dimensional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
AP Access Point . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xx
BFGS Broyden-Fletcher-Goldfarb-Shanno . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
xxiii
Glossary
I2 C Inter-Integrated Circuit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
IC Integrated Circuit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
IR Infra-Red . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
xxiv
Glossary
RF Radio Frequency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
xxv
Glossary
UK United Kingdom . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
xxvi
Chapter 1
Introduction
applications
Gyroscopes and accelerometers have traditionally been associated with the flight
control systems of an aircraft or the precision instruments within industrial
robotics. Today, they are ubiquitous; they are in our smart phones, cameras,
hard drives and toys. Advances in Microelectromechanical System (MEMS)
technology have facilitated the miniaturisation of these sensors with ever increasing
performance. In the 90s, performance of inertial MEMS was increasing by a factor
of five every year [1]. Today, performance can be achieved equally that of tactical
grade to meet the demands of the military and aerospace industry [2].
1
Chapter 1. Introduction
This project set out to demonstrate how new applications can be realised with
greater precision than previously using modern low-cost MEMS gyroscopes,
accelerometers and magnetometers. This was achieved through the development
of publicly available software and hardware platforms that have facilitated a wide
range of commercial and academic research projects exploring a diverse range of
applications.
There are several technologies that may be used for the measurement of motion.
Optical systems typically use Infra-Red (IR) light in conjunction with cameras
2
Chapter 1. Introduction
S zs
R
zr ys
d xs
yr
xr
3
Chapter 1. Introduction
these sensors provides an indicate of direction but only through their combination
in an AHRS algorithm can an absolute measurement of orientation be determined.
Until recently, commercial IMUs were limited to expensive instruments for use
in robotics and aerospace systems, and human motion capture systems. It is
only with the advent of smart phones and the proliferation of low-cost MEMS
sensors that these devices have been more widely available. Many low-cost
IMU platforms are development boards featuring sensors and a processor, for
example, the systems offered by Sparkfun, Polulu and Seeedstudio. Although
such platforms are cheap, they are often of little use as they lack the essential
sensor calibration and AHRS processing. Furthermore, these platforms do not
4
Chapter 1. Introduction
have the necessary infrastructure for application; for example, wired or wireless
communication and interfacing software. Platforms that do offer a complete
solution typically target human motion capture applications. There is an absence
of a generic IMU platforms incorporating calibrated sensors and AHRS processing
alongside a versatile communication interface.
• Design and manufacture of IMU platforms that combine the calibration and
AHRS algorithm solutions
5
Chapter 1. Introduction
AHRS algorithms
The revised AHRS algorithm presented in Chapter 7 was developed after the work
on sensor characterisation and calibration presented in Chapters 5 and 6, and
simultaneously with the IMU platforms and applications presented in Chapter 8.
The algorithm offers several improvements over the initial algorithm and other
solutions available within the subject literature.
1
http://code.google.com/p/imumargalgorithm30042010sohm/
6
Chapter 1. Introduction
This section details how this dissertation is organised along with a synopsis of each
chapter.
7
Chapter 1. Introduction
Chapter 2: Background
Some of the first work conducted in this research project was the development of
an AHRS algorithm. Although the algorithm would be significantly revised over
the proceeding years (see Chapter 7), this initial algorithm offered a novel solution
that would become widely adopted within the academic and commercial projects.
The algorithm was intended to provide an alternative solution to computationally
expensive Kalman-based solutions and explore an alternative approach to the
complementary filter solutions available. The algorithm is benchmarked against
that of a Kalman-based algorithm and demonstrated equivalent performance.
This chapter presents two schemes to measure the linear and angular kinematics of
a rigid body using a kinematically redundant array of triple-axis accelerometers.
A novel angular velocity estimation algorithm is proposed and evaluated that
can compensate for angular velocity errors using measurements of the direction
of gravity. Analysis and discussion of optimal sensor array characteristics are
provided. A damped 2 axis pendulum was used to excite all 6 Degrees of Freedom
(DOF) of the a suspended accelerometer array through determined complex motion
and is the basis of both simulation and experimental studies. The relationship
8
Chapter 1. Introduction
This chapter presents the methods and findings of the sensor characterisation
studies that form the basis of the calibration solutions presented in Chapter 6.
Several different gyroscopes, accelerometers and magnetometers from a range
of manufacturers were used throughput this project. This chapter presents the
characterisation of the two devices that would ultimately be incorporated into the
IMU platforms described in Chapter 8.
The original AHRS algorithm presented in Chapter 3 represents some of the first
work of this project. The subsequent characterisation studies presented in Chapter
5 and calibration solutions in Chapter 6 demonstrated how modern low-cost
MEMS sensors can achieve equivalent performance of such commercial IMUs. The
revised AHRS algorithm presented in this chapter was designed to operate with
9
Chapter 1. Introduction
these sensors based on the findings of Chapters 5 and 6. The revised algorithm
addresses the shortcomings of contemporary complementary filter solutions while
introducing several new features developed in conjunction with the IMU platforms
and applications presented in Chapter 8.
The chapter describes three IMU platforms; the x-IMU, the x-BIMU and the
x-OSC. Each platform fulfils a different design specification to facilitate a wide
range of applications. Two collaborative application utilising these platforms are
presented; a novel gestural music control system and a gait tracking algorithm.
This chapter also presents selected example user applications that have utilised the
platforms. The research projects chosen to appear in this chapter were selected to
demonstrate the diversity of user applications and the impact of the platforms in
peer-review publications.
Chapter 9: Conclusions
This final chapter summarises the achievements of the project, with reference to
the original research objectives. A brief description of future work is also provided.
10
Chapter 2
Background
2.1 Introduction
magnetometers
11
Chapter 2. Background
The inertial MEMS market has traditionally been dominated by the automotive
sector with almost 80% of applications being related to automotive safety features
in 2005 [15]. Today it is consumer electronics that drives the market [2] with diverse
applications ranging from hard-disk protection to complete autopilot systems
within flying toys. The Nintendo Wii games console, released in 2006, is a notable
example which attracted broad public interest in inertial MEMS and provided
researchers with a platform to innovate new applications [16, 17, 18]. The demands
of consumer electronics have focussed development on low-cost and low-power
devices. The average price of an accelerometer used in smart phones was 1.7
USD in 2006 and is predicted to be 0.2 USD in 2017 [2]. In 2012, the industries
lowest power accelerometer could provide measurements at 100 Hz while drawing
only 2 μA [19]. Smart phones and tablets typically incorporate these sensors
alongside cameras and touch screens and provide an excellent example of the new
applications that may be achieved, including: navigation, imagine stabilisation,
gestural control, augmented reality and personal health/biometrics.
12
Chapter 2. Background
2.2.1 Accelerometers
13
Chapter 2. Background
2.2.2 Gyroscopes
2.2.3 Magnetometers
14
Chapter 2. Background
Some of the most recent developments in MEMS sensors has seen the appearance
of monolithic devices that incorporate multiple sensors into a single package. What
are commonly refereed as ‘9DOF’ sensors incorporated a gyroscope, accelerometer
and magnetometer. Development is primarily driven by the low-cost demands of
consumer electronics; the budget for 9DOF sensors is less than 2.5 dollars in 2012
[2]. Several companies are also developing ‘10DOF’ devices that also incorporate a
barometer [2] though none have reached the market yet. The latest generation of
monolithic devices incorporate on-board processes to implement data processing
tasks such a calibration and AHRS data fusion.
algorithms
An AHRS algorithm is a sensor fusion algorithm that combines the data from
multiple sensors to yield a single measurement of orientation. AHRSs may
incorporate a variety of technologies, for example GPS and air speed of a UAV,
but in the context of this research project an AHRS is considered to be comprised
of a gyroscope, accelerometer and magnetometer. The accurate measurement
of orientation plays a critical role in a range of fields including: aerospace [35],
robotics [36, 37], navigation [38, 39], human motion analysis [40, 41] and machine
interaction [42].
The sensor fusion algorithm operates on the principle that each of the sensors is
able to provide some information about the orientation but not able to provide a
complete picture. The gyroscopes is the most important sensor within an AHRS.
It provides a measurement of angular velocities which may be integrated over
15
Chapter 2. Background
16
Chapter 2. Background
Although the algorithm may be derived as the sum of two filters [57, 61], most
authors express the complementary filter as a feedback loop [54, 35, 58, 55, 62,
56]. In this form, different works may be characterised by their computation of
the feedback error. For example: as the difference between angular quantities
[58, 35, 56], as the cross-product of spatial quantities [62], through rotation matrix
operations [55], or through Gauss-Newton iterations [54]. The feedback gain is
typically quantified as the crossover frequency [57, 54, 55, 35, 56]. This may
provide a convenient means of tuning the gain if the frequency characteristics of
an application are known, as in [57, 61, 56]. Most implementations will use a fixed
gain but a variable gain may enable greater accuracy if unreliable sensor data
can be detected [56]. The incorporation of integral feedback provides a means to
compensate for a gyroscope bias [62, 35]. However, this may be an undesirable
solution as it detracts from the simplicity of the complementary filter and risks
worse performance though poorly tuned second-order dynamics.
Complementary filters that use only a gyroscope and accelerometer [35, 61, 62] are
unable to provide an absolute measurement of heading. Some implementations
17
Chapter 2. Background
18
Chapter 2. Background
subject of research since the 1960s [69, 70]. For planar motion, the inverse
kinematic solution is trivial [71] and has few practical applications; for example,
the motion of a human knee [72]. Most practical applications require motion in
all six DOF to be accounted for. The so-called ‘cube configuration’ uses six single
axis accelerometers each aligned along the diagonal of a different face of a cube.
This configuration has been the focus of several studies [73, 74, 75], including
an implementation in hardware [76]. Analytical [77] and empirical [78] studies
have shown the six sensor configuration to be intrinsically unstable and limited in
application of over short time intervals. Arrays of nine accelerometers provide an
alternative, stable solution [79, 80, 81].
The arrays of six and nine accelerometer axis use specific configurations of multiple
single-axis accelerometers. Modern MEMS accelerometers are widely available
in triple-axis packages containing three mutually orthogonal accelerometer axes.
The theoretical minimum number of triple-axis accelerometers required to measure
motion in six DOF is four (incorporating a total of 12 linear accelerometers axis).
Schemes using triple-axis accelerometers [67, 82, 83, 84, 85] do not require specific
geometric configurations, only that the position of a sensor within an array is
known. As each sensor is able to provide a vector measurement of acceleration in
three dimensions, the physical orientation of the sensor is irrelevant provided that
the orientation is known. The theoretical irrelevance of each sensor’s position and
orientation within the sensor array gives triple-axis accelerometer arrays a clear
practical advantage over single-axis accelerometer arrays.
19
Chapter 2. Background
[87].
quantities
1
https://github.com/xioTechnologies/Quaternion-MATLAB-Library
20
Chapter 2. Background
ẑB ẑA
ŷB
ŷA
x̂B
x̂A
If A and B are aligned then the rotation matrix would be an identity matrix as
21
Chapter 2. Background
A
In application, BR may be a measurement of orientation provided by an AHRS
where A is a coordinate system fixed to the Earth and B represents the axes of
the IMU.
B −1 T
AR =A
BR =A
BR (2.3)
A
CR =A B
B RC R (2.4)
The rotation matrix can be obtained using Craig’s notation system by cancelling
the diagonal opposite subscript and superscript. The requirement for these
22
Chapter 2. Background
The multiplication of a rotation matrix by a vector quantity will alter the reference
coordinate system of the vector and thus achieve a rotation. For example, A v is
a vector defined in A, the vector can be equivalently described in B through the
multiplication of Equation 2.5.
A
v=A B
BR v (2.5)
2.5.2 Quaternions
23
Chapter 2. Background
ẑB ẑA
A
θ r̂
ŷB
ŷA
x̂B
x̂A
2.6 where rx , ry and rz are the elements of A r̂. The last three elements form the
imaginary parts of the quaternion.
A
B q̂ = qw qx qy qz
(2.6)
= cos 2θ −rx sin 2θ −ry sin 2θ −rz sin 2θ
A
B q̂ = 1 0 0 0 , for A = B (2.7)
The convention used throughput this dissertation positions qw as the first element,
this is consistent with Kuipers [93] and MATLAB [94]. Other implementations
may position qw as the last element; notable examples include OpenGL libraries
[95] and the Microsoft XNA Framework [96].
24
Chapter 2. Background
B A ∗
A q̂ = B q̂ = qw −qx −qy −qz (2.8)
of Equation 2.9.
A
C q̂ =B A
C q̂ ⊗ B q̂ (2.9)
The multiplication of two quaternions, q̂a and q̂b , can be determined as the
Hamilton product and defined as Equation 2.10 [93]. A quaternion product is
not commutative; that is, q̂a ⊗ q̂b 6= q̂b ⊗ q̂a .
q̂a ⊗ q̂b = qaw qax qay qaz ⊗ qbw qbx qby qbz
T
qaw qbw − qax qbx − qay qby − qaz qbz
(2.10)
qaw qbx + qax qbw + qay qbz − qaz qby
=
q q − q q + q q + q q
aw by ax bz ay bw az bx
qaw qbz + qax qby − qay qbx + qaz qbw
25
Chapter 2. Background
r11 − r22 − r33 r21 + r12 r31 + r13 r23 − r32
1
r 21 + r12 r 22 − r11 − r33 r32 + r 23 r31 − r 13
K= (2.13)
3 r +r r32 + r23 r33 − r11 − r22 r12 − r21
31 13
r23 − r32 r31 − r13 r12 − r21 r11 + r22 + r33
A
B q̂ = v3 v0 v1 v2 (2.14)
Some practical applications require that a ‘best fit’ rotation matrix is obtained
from an approximated rotation matrix derived from measurement data. This can
26
Chapter 2. Background
be achieved by first obtaining a quaternion using the method above and then
converting the quaternion into a rotation matrix as per Equation 2.12.
Although rotation matrices and quaternions provide effective tools for the
manipulation of spatial quantities, their numerical values and not immediately
meaningful. Euler angles describe an orientation as three angular quantities and
so offer a more intuitive representation. Perhaps the most common use of Euler
angles is within the aerospace sequence where the three angular quantities φ, θ
and ψ correspond to a rotations around the x, y and z axes respectively as shown
in Figure 2.3. φ, θ and ψ are commonly referred to as roll, pitch and yaw (or
heading) respectively.
θ
x φ y
Figure 2.3: The aerospace sequence of Euler angles describe orientation as the
result of angular rotations φ , θ and ψ around the x, y and z axes respectively
27
Chapter 2. Background
Figure 2.4: The Euler angle sequence describes successive rotations of around each
axis to achieve the ultimate orientation of B relative to A
28
Chapter 2. Background
2.15 to 2.17 by substituting each element of the rotation matrix with function of
the quaternion as per Equation 2.12 to yield Equations 2.18 to 2.20.
φ = atan2 2 (qy qz − qw qx ) , 2qw2 − 1 + 2qz2 (2.18)
θ = − arcsin 2 (qx qz + qw qy ) (2.19)
ψ = atan2 2 (qx qy − qw qz ) , 2qw2 − 1 + 2qx2 (2.20)
29
Chapter 2. Background
30
Chapter 3
3.1 Introduction
Some of the first work conducted in this research project was the development of an
AHRS algorithm. Although the algorithm would be significantly revised over the
proceeding years (see Chapter 7), this initial algorithm offered a novel solution
that would become widely adopted within academic and commercial projects.
The algorithm was intended to provide an alternative solution to computationally
expensive Kalman-based solutions and explore an alternative approach to the
complementary filter solutions available. The principle goals were:
31
Chapter 3. Initial Attitude and Heading Reference System algorithm
A tri-axis gyroscope will measure the angular rate about the x, y and z axes of the
sensor frame, termed ωx , ωy and ωz respectively. If these parameters (in rads−1 )
are arranged into the vector S ω defined by Equation 3.1, the quaternion derivative
describing rate of change of the Earth frame relative to the sensor frame SE q̇ can
be calculated [98] as Equation 3.2. The ⊗ operate denotes a quaternion product
32
Chapter 3. Initial Attitude and Heading Reference System algorithm
S
ω = 0 ωx ωy ωz (3.1)
S 1S S
E q̇ = E q̂ ⊗ ω (3.2)
2
The orientation of the Earth frame relative to the sensor frame at time t, E
S qω,t ,
S
can be computed by numerically integrating the quaternion derivative E q̇ω,t as
described by Equations 3.3 and 3.4, provided that initial conditions are known. In
these equations, S ωt is the angular rate measured at time t, ∆t is the sampling
S
period and E q̂est,t−1 is the previous estimate of orientation. The subscript ω
indicates that the quaternion is calculated from angular rates.
S 1S S
E q̇ ω,t = E q̂est,t−1 ⊗ ωt (3.3)
2
S
E qω,t = SE q̂est,t−1 + SE q̇ ω,t ∆t (3.4)
33
Chapter 3. Initial Attitude and Heading Reference System algorithm
orientations achieved by the rotation of the true orientation around an axis parallel
with the field. A quaternion representation requires a single solution to be found.
This may be achieved through the formulation of an optimisation problem where
an orientation of the sensor, SE q̂, is found as that which aligns a predefined reference
ˆ with the measured field in the sensor
direction of the field in the Earth frame, E d,
frame, S ŝ; thus solving 3.5 where Equation 3.6 defines the objective function.
ˆ S ŝ)
min f (SE q̂, E d, (3.5)
S q̂∈<4
E
ˆ S ŝ) = S q̂ ∗ ⊗ E dˆ ⊗ S q̂ − S ŝ
f (SE q̂, E d, (3.6)
E E
Many optimisation algorithms exist but the gradient descent algorithm is one of
the simplest to both implement and compute. Equation 3.7 describes the gradient
descent algorithm for n iterations resulting in an orientation estimation of SE q̂n+1
S
based on an ‘initial guess’ orientation E q̂0 and a variable step-size µ. Equation
3.8 computes an error direction on the solution surface defined by the objective
function, f , and its Jacobian, J .
ˆ S ŝ)
∇f T (SE q̂ k , E d,
S
E qk+1 = SE q̂ k − µ , k = 0, 1, 2...n (3.7)
∇f T (S q̂ , E d, ˆ S ŝ)
E k
ˆ S ŝ) = J T (S q̂ , E d)f
∇f (SE q̂ k , E d, ˆ (S q̂ , E d,
ˆ S ŝ) (3.8)
E k E k
Equations 3.7 and 3.8 describe the general form of the algorithm applicable to a
field predefined in any direction. However, if the reference direction of the field is
defined to only have components within 1 or 2 of the principle axes of the Earth
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
34
Chapter 3. Initial Attitude and Heading Reference System algorithm
S
E q̂ = qw qx q y qz (3.9)
E
ĝ = 0 0 0 1 (3.10)
S
â = 0 ax ay az (3.11)
2(qx qz − qw qy ) − ax
fg (SE q̂, S â) =
2(qw qx + qy qz ) − ay (3.12)
1 2 2
2( 2 − qx − qy ) − az
−2qy 2qz −2qw 2qx
Jg (SE q̂) =
2qx 2qw 2qz 2qy
(3.13)
0 −4qx −4qy 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 United Kingdom (UK)
[99]. This can be represented by Equation 3.14. Substituting E b̂ and normalised
magnetometer measurement S m̂ for E dˆ and S ŝ respectively, yields the simplified
objective function and Jacobian defined by Equations 3.16 and 3.17.
E
b̂ = 0 bx 0 bz (3.14)
35
Chapter 3. Initial Attitude and Heading Reference System algorithm
S
m̂ = 0 mx my mz (3.15)
2 2
2bx (0.5 − qy − qz ) + 2bz (qx qz − qw qy ) − mx
fb (SE q̂, E b̂, S m̂) = 2b
x x y(q q − q q
w z ) + 2b (q q
z w x + q q
y z ) − m
y (3.16)
2bx (qw qy + qx qz ) + 2bz (0.5 − qx2 − qy2 ) − mz
−2bz qy 2bz qz
Jb (SE q̂, E b̂) =
−2bx qz + 2bz qx 2bx qy + 2bz qw
2bx qy 2bx qz − 4bz qx
(3.17)
−4bx qy − 2bz qw −4bx qz + 2bz qx
2bx qx + 2bz qz −2bx qw + 2bz qy
2bx qw − 4bz qy 2bx qx
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 3.18 and 3.19. Whereas the solution surface created by the objective
functions in Equations 3.12 and 3.16 have a global minimum defined by a line,
the solution surface defined by Equation 3.18 has a minimum defined by a single
point, provided that bx 6= 0.
fg (SE q̂, S â)
fg,b (SE q̂, S â, E b̂, S m̂) = (3.18)
fb (SE q̂, E b̂, S m̂)
36
Chapter 3. Initial Attitude and Heading Reference System algorithm
JgT (SE q̂)
Jg,b (SE q̂, E b̂) = (3.19)
JbT (SE q̂, E b̂)
S ∇f T
E q ∇,t = SE q̂ est,t−1 − µt (3.20)
k∇f T k
JgT (SE q̂est,t−1 )fg (SE q̂est,t−1 , S ât )
∇f = (3.21)
Jg,b
T S
(E q̂est,t−1 , E b̂)fg,b (SE q̂est,t−1 , S ât , E b̂, S m̂t )
S
An appropriate value of µt is that which ensures the convergence rate of E q ∇,t
S
µt = α E q̇ω,t ∆t, α > 1 (3.22)
37
Chapter 3. Initial Attitude and Heading Reference System algorithm
S
In practice, E qω,t may start from incorrect initial conditions and accumulate
S
errors due to gyroscope measurement noise and E q∇,t will provide an incorrect
estimate when the accelerometer is not stationary or the magnetometer exposed
to interferences. The goal of the fusion algorithm is to provide an orientation
estimate where SE qω,t is used to filter out high frequency errors in SE q∇,t , and SE q∇,t
is used both to compensate for integral drift in SE qω,t and to provide convergence
from initial conditions.
An estimated orientation of the Earth frame relative to the sensor frame, SE qest,t ,
is obtained through the fusion of the two separate orientation calculations, SE qω,t
and SE q∇,t as described by Equation 3.23 where γt and (1 − γt ) are weights applied
to each orientation calculation.
S
E qest,t = γt SE q∇,t + (1 − γt )SE qω,t , 0 ≤ γt ≤ 1 (3.23)
µt
(1 − γt )β = γt (3.24)
∆t
β
γt = µt (3.25)
∆t
+β
38
Chapter 3. Initial Attitude and Heading Reference System algorithm
S S
The fusion process ensures the optimal fusion of E qω,t and E q∇,t assuming that
the convergence rate of SE q ∇ 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 3.22, also becomes very large and the
equations simplify. A large value of µt used in Equation 3.20 means that SE q̂est,t−1
becomes negligible and the equation can be re-written as Equation 3.26.
S ∇f
E q ∇,t ≈ −µt (3.26)
k∇f k
The definition of γt in Equation 3.25 also simplifies if the β term in the denominator
becomes negligible and the equation can be rewritten as Equation 3.27. It is
possible from Equation 3.27 to also assume that γt ≈ 0.
β∆t
γt ≈ (3.27)
µt
Substituting Equations 3.4, 3.26 and 3.27 into Equation 3.23 directly yields
Equation 3.28. It is important to note that in Equation 3.28, γt has been
substituted as both as Equation 3.26 and 0.
S β∆t ∇f S
E qest,t = −µt + (1 − 0) E q̂est,t−1 + SE q̇ ω,t ∆t (3.28)
µt k∇f k
S
Equation 3.28 can be simplified to Equation 3.29 where E q̇est,t is the estimated
orientation rate defined by Equation 3.30.
S
E qest,t = SE q̂est,t−1 + SE q̇est,t ∆t (3.29)
S ∇f
E q̇est,t = SE q̇ ω,t − β (3.30)
k∇f k
39
Chapter 3. Initial Attitude and Heading Reference System algorithm
It can be seen from Equations 3.29 and 3.30 that the algorithm calculates
S
the orientation E qest by numerically integrating the estimated rate of change
S S
of orientation E q̇ est . The algorithm computes E 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 a direction based on accelerometer and
magnetometer measurements. Figure 3.1 shows a block diagram representation
of the complete orientation estimation algorithm implementation for an IMU.
∇f
k∇f k
z−1
β
Z
1S q S
Gyroscope S ω t q̂ ⊗ S ωt .dt
kqk E q̂ est,t
2 E est,t−1
S
E q̇ est,t
z−1
40
Chapter 3. Initial Attitude and Heading Reference System algorithm
The measured direction of the Earth’s magnetic field in the Earth frame at time
t, E ĥt , can be computed as Equation 3.31. The effect of an erroneous inclination
E
of the measured direction Earth’s magnetic field, ĥt , can be corrected if the
E
algorithm’s reference direction of the Earth’s magnetic field, b̂t , is of the same
E E
inclination. This is achieved by computing b̂t as ĥt normalised to have only
components in the Earth frame x and z axes; as described by Equation 3.32.
E ∗
ĥt = 0 hx hy hz = SE q̂est,t−1 ⊗ S m̂t ⊗ SE q̂est,t−1 (3.31)
E p 2
b̂t = 0 hx + h2y 0 hz (3.32)
41
Chapter 3. Initial Attitude and Heading Reference System algorithm
S
E q̂ est,t−1 ⊗ S m̂t ⊗ SE q̂ ∗est,t−1
E
ĥt
h q i
0 h2x + h2y 0 hz
E
b̂t
S
Magnetometer 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
∇f
k∇f k
z−1
β
Z
1S q S
S
Gyroscope ω t q̂ ⊗ S ωt .dt
kqk E q̂ est,t
2 E est,t−1
S
E q̇ est,t
z−1
r
1 3
β= q̂ ⊗ 0 ω̃max ω̃max ω̃max = ω̃max (3.33)
2 4
The algorithm was tested using the Xsens MTx orientation sensor [105] containing
16-bit resolution tri-axis gyroscopes, accelerometers and magnetometers. Raw
sensor data was logged to a computer at 512 Hz and imported using the
accompanying software to provide calibrated sensor measurements which were then
processed by the proposed orientation estimation algorithm. The software also
incorporates a propriety Kalman-based orientation estimation algorithm. As both
the Kalman-based algorithm and proposed algorithm’s estimates of orientation
were computed using identical sensor data, the performance of each algorithm
could be evaluated relative to one-another, independent of sensor performance.
42
Chapter 3. Initial Attitude and Heading Reference System algorithm
The Xsens Kalman-based algorithm was used with the factory default settings.
−50
−100
42 44 46 48 50 52 54 56 58 60 62
Error
5
Kalman−based algorithm
θε (degrees)
−5
42 44 46 48 50 52 54 56 58 60 62
time (seconds)
Figure 3.3: Typical results for measured and estimated angle θ (top) and error
(bottom)
It is common [104, 105, 108, 109, 110] to quantify orientation sensor performance
as the static and dynamic Root Mean Squared (RMS) errors in the decoupled
43
Chapter 3. Initial Attitude and Heading Reference System algorithm
Table 3.1: Static and dynamic RMS error of Kalman-based algorithm and proposed
algorithm IMU and MARG implementations
The static and dynamic RMS values of φ , θ , and ψ were calculated assuming
a static state when the measured corresponding angular rate was < 5◦ /s, and a
dynamic when ≥ 5◦ /s. This threshold was chosen to be sufficiently greater than
the noise floor of the data. The results are summarised in Table 3.4 where each
value represents the mean of all eight experiments.
44
Chapter 3. Initial Attitude and Heading Reference System algorithm
1.5 1.5
β = 0.033 β = 0.041
1 1
0.5 0.5
0 0.2 0.4 0 0.2 0.4
β β
Figure 3.4: The effect of the adjustable parameter, β, on the performance of the
proposed algorithm IMU (left) and MARG (right) implementations
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 3.5: The effect of sampling rate on the performance of the proposed
algorithm IMU (left) and MARG (right) implementations
45
Chapter 3. Initial Attitude and Heading Reference System algorithm
simulate sampling rates between 1 Hz and 512 Hz. It can be seen from Figure 3.5
that the proposed algorithm achieves similar levels of performance at 50 Hz as at
512 Hz. Both algorithm implementations are able to achieve a static error < 2◦
and dynamic error < 7◦ while sampling at 10 Hz. This level of accuracy may be
sufficient for human motion applications though the sampling rate will limit the
bandwidth of the motion that may be measured.
3.5 Conclusions
46
Chapter 3. Initial Attitude and Heading Reference System algorithm
The algorithm presented in this chapter was published in 2011 [5] with the source
code initially made available on Google Code1 . The paper and its associated
technical report have since totalled over 100 citations [6]. Google Code indicates
over 10,000 downloads [7]. However, as the files have since been published in
multiple locations the total number of downloads is assumed to be higher.
1
http://code.google.com/p/imumargalgorithm30042010sohm/
47
Chapter 3. Initial Attitude and Heading Reference System algorithm
48
Chapter 4
Kinematically redundant
non-gyro IMUs
4.1 Introduction
This chapter presents two schemes of measuring the linear and angular kinematics
of a rigid body using a kinematically redundant array of triple-axis accelerometers
with potential applications in biomechanics. A novel angular velocity estimation
algorithm is proposed and evaluated that can compensate for angular velocity
errors using measurements of the direction of gravity. Analysis and discussion
of optimal sensor array characteristics are provided. A damped 2 axis pendulum
was used to excite all 6 DOF of the a suspended accelerometer array through
determined complex motion and is the basis of both simulation and experimental
studies. The relationship between accuracy and sensor redundancy is investigated
for arrays of up to 100 triple axis (300 accelerometer axes) accelerometers in
simulation and 10 equivalent sensors (30 accelerometer axes) in the laboratory
test rig. The work also reports on the sensor calibration techniques and hardware
implementation.
49
Chapter 4. Kinematically redundant non-gyro IMUs
accelerometers
ẑ i , αx
ẑ i′
ẑ B , z̈
ŷ i′
B
di
ŷ i , αy
r q x̂i , αx x̂ ′
i
p ŷ B , ÿ
x̂B , ẍ
Figure 4.1: Schematic describing the i’th accelerometer relative to the ridged body
kinematic frame
frame then it can be considered aligned to the arbitrary frame i0 defined by the
mutually orthogonal unit vectors x̂i0 , ŷi0 and ẑi0 (shown in Figure 4.1), providing
0
the measurement vector i α. The orientation of frame i0 relative to frame i is
described by the rotation matrix ii0 R so that the aligned acceleration measurement
i
α may be obtained using Equation 4.1.
i 0
α = ii0 R i α (4.1)
50
Chapter 4. Kinematically redundant non-gyro IMUs
The linear acceleration measured by the i’th sensor, i α, may be defined by the
kinematics of the body origin using Equation 4.2 where B ω and B a describe the
angular velocity and linear acceleration of the rigid body origin respectively. This
represents the sum of linear, tangential and centripetal accelerations. There is no
Coriolis term as the accelerometer position is fixed relative to the body.
i
α = B ω̇ × B di + B ω × (B ω × B di ) + B a (4.2)
T
i
α = αx αy αz (4.3)
T
B
ω= p q r (4.4)
T
B
di = dx dy dz (4.5)
T
B
a = ẍ ÿ z̈ (4.6)
0 -dx -dx 0 dz -dy dy 0 dz
Gi = I
3 -d y 0 -d y -d z 0 d x dx d z 0
(4.8)
-dz -dz 0 dy -dx 0 0 dy dx
51
Chapter 4. Kinematically redundant non-gyro IMUs
and
T
s= ẍ ÿ z̈ p2 q 2 r2 ṗ q̇ ṙ pq qr rp (4.9)
The inverse kinematic solution can be found directly from Equation 4.7 through
an inversion of Gi . If the number of sensors fixed to the body is four
then accelerometer measurements may be arranged as the vector αC and the
corresponding square matrix C constructed as described by Equations 4.10 and
4.11. C may be inverted to provide the inverse kinematic solution shown as
Equation 4.12.
T
αC = 1
αT 2
αT 3
αT 4
αT (4.10)
T
C = G1 T
G2 T
G3 T
G4 T (4.11)
s = C -1 αC (4.12)
n n!
m= = (4.13)
4 24(n − 4)!
52
Chapter 4. Kinematically redundant non-gyro IMUs
1 X
m
s= sj (4.14)
m j=0
If the positions of any sensors are equal or if all sensors exist on a plane (for a
non-planar system) then Gi is singular and cannot be inverted. For n > 4, a
unique value of C -1 must be computed for each combination. This method may
result in a considerable computational load as an array of n sensors requires m
12×12 matrix inversions and multiplications.
T
P = G1 T G2 T G3 T ... Gn T (4.16)
s = P + αP (4.17)
If the positions of any sensors are equal or if all sensors exist on a plane (for a
non-planar system) then the rows of P are not linearly independent and Equation
53
Chapter 4. Kinematically redundant non-gyro IMUs
4.18 may not be used. The computational requirement of this method considerably
is less than that of the combinatorial method as only one n×12 matrix inversion
and multiplication is required for n sensors.
The state vector s does not directly yield individual angular velocity terms. This
is a problem common to all accelerometer-only inertial measurement methods and
many solutions have proposed. Cardou et al. provide a discussion and analysis of
the different approaches to this problem and present a method more robust than
many existing solutions [112].
54
Chapter 4. Kinematically redundant non-gyro IMUs
The available angular velocity terms may be used to construct an objective function
f (ω̃, s) as shown in Equation 4.19 where ω̃ is the estimated angular velocity. If
f (ω̃, s) = 0 then ω̃ = ±ω. As each angular velocity term is only available as
the product of itself or another, the sign of each quantity is lost and cannot be
recovered directly.
2 2
p̃ − p
q̃ 2 − q 2
2
r̃ − r2
f (ω̃, s) =
(4.19)
p̃q̃ − pq
q̃r̃ − qr
r̃p̃ − rp
T
ω̃ = p̃ q̃ r̃ (4.20)
55
Chapter 4. Kinematically redundant non-gyro IMUs
ηp
η = JfT (ω̃)f (ω̃, s) (4.21)
q
ηr
2 2
p̃ − p
T q̃ 2 − q 2
2p̃ 0 0 q̃ 0 r̃ 2
r̃ − r2
=
0 2q̃ 0 p̃ r̃ 0
(4.22)
p̃q̃ − pq
0 0 2r̃ 0 q̃ p̃
q̃r̃ − qr
r̃p̃ − rp
2 2
2p̃(p̃ − p ) + q̃(p̃q̃ − pq) + r̃(r̃p̃ − rp)
=
2q̃(q̃ 2
− q 2
) + p̃(p̃q̃ − pq) + r̃(q̃r̃ − qr)
(4.23)
2r̃(r̃2 − r2 ) + q̃(q̃r̃ − qr) + p̃(r̃p̃ − rp)
The direction of this vector may be combined with a magnitude computed as the
difference between the magnitudes of the estimated angular velocities and actual
angular velocities. The error in the angular velocity is therefore estimated, ẽ, as
the dot product shown in Equation 4.24. Due to the loss of sign in the available
angular velocity terms the error will be computed relative to the closest value of
±ω.
p
sign(η )
p |p̃| − p2
ẽ = sign(η ) · |q̃| − pq 2 (4.24)
q
√
sign(ηr ) |r̃| − r2
The estimated angular velocity is computed from the angular acceleration and
estimated error in the angular velocity using Equation 4.25. This is represented
as the block diagram shown as Figure 4.2.
Z
ω̃ = (ω̇ − K ẽ).dt (4.25)
56
Chapter 4. Kinematically redundant non-gyro IMUs
R
ω̇ ω̃
K
2
p
q2
2
r
Calculate ẽ
pq
qr
rp
As the estimated error may only be computed relative to the closest value of ±ω,
ω̃ will converge to the closest value of either ±ω. A correct sign assertion may be
assured if the dynamics in ω̇ are of a sufficiently large magnitude relative to the
convergence rate of ω̃ governed by K, at the zero-crossing point of ω. A sufficient
magnitude of ω̇ relative to the convergence rate of ω̃ would mean that integral
drift would not dominate the ω̃ dynamics and the value of ω̃ would not drift across
zero and hence cause an incorrect sign.
57
Chapter 4. Kinematically redundant non-gyro IMUs
drift compensation
In many applications it can be assumed that the mean linear acceleration measured
over an extended period of time represents gravity. The direction of gravity
measured relative to the accelerometer array may therefore be used to estimate
the attitude of the array relative to the Earth’s surface. As will be shown in
Section 4.6.3, sensor measurement and alignment errors result in a bias error in
the estimated angular velocity at low velocities. Knowledge of the accelerometer
array orientation may be used to correct for steady state errors in the estimation
angular velocity and so compensate for integral drift.
Z
B B
Eq = E q̇ (4.27)
B 1B
E q̇ = E q̂ ⊗ 0 (B ω − δ)T (4.28)
2
Z
δ = KP + KI (4.29)
58
Chapter 4. Kinematically redundant non-gyro IMUs
2qx qz − qw qy
= Ba ×
2qw qx + qy qz
(4.30)
2 2 2 2
qw − qx − qy + qz
The error in the estimated angular velocity is represented by the integral feedback
term. The compensated estimated angular velocity is therefore obtained as
Equation 4.31.
Z
0
ω̃ = ω̃ − Ki (4.31)
A 2 axis pendulum was the chosen subject of simulation and experimental studies
as it approximates a number of human and robotic manipulation joints and the
oscillatory motion is analogous to gait and excites all 6 DOF of the suspended
body. The damped motion results in both low and high magnitude rates and the
complex motion of the 2 axis joint results in rotational kinematics in all three axes
of the suspended body. The pendulum is represented by the schematic shown as
Figure 4.3 where θ1 and θ2 represent the angle in each axis of the joint. The 2
axis pendulum equations of motion have been derived as Equations 4.32 and 4.33
where c1 and c2 are the damping of each axis of the joint.
g sin θ1 c1 θ̇1
θ̈1 = 2θ̇1 θ̇2 tan θ2 − − (4.32)
l cos θ2 m cos2 θ2
g c2 θ̇2
θ̈2 = −θ̇12 cos θ2 sin θ2 − cos θ1 sin θ2 − (4.33)
l m
59
Chapter 4. Kinematically redundant non-gyro IMUs
{E}
{A} z zE
A
g
ϴ1 yA
yE
xE xA ϴ2
l
m
θ̇1 cos θ2
B
ω=
θ̇ 2
(4.35)
θ̇1 sin θ2
θ̈1 cos θ2 − θ̇1 θ̇2 sin θ2
B
ω̇ =
θ̈ 2
(4.36)
θ̈1 sin θ2 + θ̇1 θ̇2 cos θ2
cos θ2 0 sin θ2
E
B R = sin θ1 sin θ2 cos θ1 - sin θ1 cos θ2 (4.37)
- cos θ1 sin θ2 sin θ1 cos θ1 cos θ2
60
Chapter 4. Kinematically redundant non-gyro IMUs
Simulation studies were used for the initial testing of the proposed methods
and to investigate the relationship between system accuracy and redundancy.
Simulations used the damped 2 axis pendulum model with an array of n triple-axis
accelerometers positioned evenly around a 360◦ helix concentric with the length of
the pendulum and with a radius of 0.05 m. Figure 4.4 shows the pendulum and
triple-axis accelerometers as drawn by the simulation software. For convenience,
all sensors were aligned with the pendulum body coordinate frame; arbitrary
orientations were used and calibrated for in the experimental studies discussed
in Section 4.6.
Earth frame
Pendulum frame and arm
Pendulum axes
Accelerometer
Accelerometer x axes
0.05
Accelerometer y axes
0 Accelerometer z axes
−0.05
Z (m)
−0.25
−0.2
0
−0.2 −0.1 0.2
0 0.1 0.2 0.3
X (m)
Y (m)
61
Chapter 4. Kinematically redundant non-gyro IMUs
The integral drift inherent to the angular velocity estimation algorithm would
introduce a bias in the estimated velocity at low rates. This random steady
state error meant that the apparent relationship between the angular velocity
error and redundancy was non-monotonic. To compensate for this, the mean
result of 5 simulations was computed for each value of redundancy. The results
for angular acceleration and velocity are shown as Figure 4.5 and Figure 4.6.
It can be seen that the combinatorial method results in an greater error than
the pseudoinverse method at greater computational expense. Simulations of the
combinatorial method for arrays > 33 sensors were abandoned as the time taken
to process arrays of this size became impractical. The non-monotonic relationship
observed in the estimated angular velocity is due to the ‘random walk’ in this
state that comes about at low velocities and is constrained by the feedback in
the angular velocity estimation filter. This ‘random walk’ is due to the noise
62
Chapter 4. Kinematically redundant non-gyro IMUs
Oscillatory motion of the pendulum and the frequency response of the sensors
results in a minimum achievable error represented by the oscillatory error
components due to the phase-lag of the sensor outputs (due to a simulated
anti-aliasing filter). This minimum achievable error was computed by running a
simulation using an accelerometer error model comprised of only the anti-aliasing
filter, i.e. zero noise. The resultant error is the same for arrays of all sizes and
represents the minimum achievable error. The respective minimum achieve errors
are indicated in Figure 4.5 and Figure 4.6.
0.26 0.038
mean[RM S[pǫ ], RM S[qǫ ], RM S[rǫ ] (rads−1 )
mean[RM S[ṗǫ ], RM S[q̇ǫ ], RM S[ṙǫ ] (rads−2 )
0.2 0.032
0.18 0.03
0.16 0.028
0.14 0.026
0.12 0.024
0.1 0.022
5 20 40 60 80 100 5 20 40 60 80 100
Redundancy (number of tri-axis sensors) Redundancy (number of tri-axis sensors)
Figure 4.5: Angular acceleration error vs.Figure 4.6: Angular velocity error vs.
sensor redundancy sensor redundancy
The pendulum testbed allows more complex motion to be evaluated. If the initial
conditions use a non-zero velocity vector and dissimilar damping coefficients, all
cross terms in Equation 4.7 can be elicited at a range of frequencies. The spectrum
of any individual accelerometer axis will contain frequencies up to the natural
frequency of the pendulum and hence the results can be demonstrated to work
across this frequency range. The upper graphs in Figures 4.7 and 4.8 demonstrate
63
Chapter 4. Kinematically redundant non-gyro IMUs
ṗ
50 q̇
ṙ
rads−2
−50
−100
3 3.5 4 4.5 5 5.5 6
5 ṗǫ
q̇ǫ
ṙǫ
rads−2
−5
3 3.5 4 4.5 5 5.5 6
Time (s)
p̃′
5 q̃ ′
r̃ ′
rads−1
−5
0.8
p̃′ǫ
0.4 q̃ǫ′
r̃ǫ′
rads−1
−0.4
−0.8
2 2.5 3 3.5 4 4.5 5 5.5 6 6.5
Time (s)
Experimental studies were conducted to verify the above simulation results. The
studies were designed to be as similar to the simulation studies as possible and
64
Chapter 4. Kinematically redundant non-gyro IMUs
4.6.1 Hardware
1
https://github.com/xioTechnologies/DAQ32
65
Chapter 4. Kinematically redundant non-gyro IMUs
the accelerometer package to minimise size and weight. The probe’s 0.4 m
lead consisted of wires twisted together so that the power supply wires would
provide a level of electromagnetic shielding to the accelerometer signals. Fine wire
connections were sunk in a silicon sealant to provide mechanical durability. The
complete probe (shown in Figure 4.9) occupied a volume of a 5 mm cube and
weighed less than 1 g (excluding lead).
66
Chapter 4. Kinematically redundant non-gyro IMUs
Figure 4.10: 2 axis pendulum with attached accelerometer probes and annotated
potentiometer axes θ1 and θ2
4.6.2 Calibration
67
Chapter 4. Kinematically redundant non-gyro IMUs
individual axis bias and sensitivity respectively. The skew matrix S accounts
for non-orthogonality and coupling between the accelerometer axes. The gain smn
defines the proportion of axis n affecting axis m.
α = SKu − b
1 sxy sxz kx 0 0 ux bx (4.38)
= s
yx 1 s 0 k 0 u − b
yz y y y
szx szy 1 0 0 kz uz bz
g = kSKu − bk (4.39)
X
min (g − kSKuj − bk)2 (4.40)
S,K,b∈<
j
68
Chapter 4. Kinematically redundant non-gyro IMUs
pitch and roll axes through their full range in 10◦ steps so that every combination
of pitch and roll angle was achieved. At each orientation, the gimbal was allowed
10 seconds to come to rest and then the mean of each accelerometer output was
taken for a 15 second period.
A value of g was obtained as 9.812 ms−2 for the geographical location of the
studies [113]. Equation 4.40 was solved using MATLAB’s Optimisation toolbox’s
fminunc function and the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method.
Initial values for gains and biases were chosen as the typical values detailed in
the AXDL335 datasheet. Calibration was verified by comparing the measured
magnitude of gravity with the known magnitude at each orientation within the
calibration dataset. Figure 4.11 shows typical calibrated measurements of one
accelerometer for the full calibration dataset.
Calibrated accelerometer output
10
αx
αy
αz
5
ms−2
−5
−10
Calibrated accelerometer measured magnitude of gravity
9.85 q
α2x + αy2 + αz2
ms−2
9.8 g
9.75
100 200 300 400 500 600 700
Dataset sample
4.6.2.2 Alignment
69
Chapter 4. Kinematically redundant non-gyro IMUs
0
accelerometer measurement, i α, is related to the misaligned measurement i α by
the rotation matrix ii0 R as described by Equation 4.1. The alignment may be
equivalently represented by the quaternion ii0 q. If the array is stationary, Equation
4.41 will be true for sensors at all orientations. Therefore ii0 q may be found as the
0
solution to 4.42 where i αj and B aj represent the measured acceleration of the i’th
accelerometer and the acceleration of the body at the j’th orientation respectively.
As the acceleration of the body is unknown, one fiducial sensor must be chosen
B
to provide a measurement equal to a. An alternative kinematic minimisation
approach for alignment calibration was also implemented but performance was
found to be no better than simply choosing a convenient sensor to be the fiducial.
0
B
a = ii0 q̂ ∗ ⊗ i α ⊗ ii0 q̂ (4.41)
Xh i2
i ∗ i0 i B
min
i q∈<
i0 q̂ ⊗ αj ⊗ i0 q̂ − aj (4.42)
i0 j
Equation 4.42 was solved using the same calibration dataset and method described
in Section 4.6.2.1. The fiducial sensor was chosen as that judged to be best
aligned to the pendulum on visual inspection. This gravitational based attitude
calibration scheme is similar to previously proposed techniques; however, the use of
a quaternion representation avoids problems of non-orthogonality associated with
a rotation matrix based approach [114].
4.6.3 Results
The pendulum was left to swing with motion in both axes remaining in phase
until rest while data from each of the 10 triple-axis accelerometers and 2 joint
potentiometers was recorded. The ‘true’ suspended body kinematics were then
estimated for the first 30 seconds of the dataset using the potentiometer data and
70
Chapter 4. Kinematically redundant non-gyro IMUs
Figure 4.12 and Figure 4.13 show the pendulum kinematic states as estimated
using the pseudoinverse method and data from all 10 sensors. The angular
velocity estimation algorithm used a gain of K = 2 and the orientation estimation
algorithm used gains of KP = 5 and KI = 1.0. Each figure consists of 2 plots: the
upper plot represents the measured or estimated state and the lower plot represents
the corresponding error. Small errors in alignment calibration and the constant
acceleration of gravity result in an angular velocity bias error and low velocities
(see Figure 4.13). The orientation estimation algorithm described in Section 4.3
compensates for this error in p and q but errors r cannot be compensated for as
this axis of rotation eventually aligns with the direction of gravity.
Angular acceleration
40
ṗ
20 q̇
ṙ
rads−2
−20
−40
0 5 10 15 20 25 30
Angular acceleration error
5 ṗǫ
q̇ǫ
ṙǫ
rads −2
−5
0 5 10 15 20 25 30
Time (s)
To investigate the relationship between error and redundancy, the suspended body
kinematics were computed using the pseudoinverse method for each combination
of sensors for the available array of ten sensors, thus results from four to ten
sensors were possible. Since the number of combinations for a given number of
sensors is defined by 10
n
, 210 results were possible for 4 sensors, rising to 252 for
5 sensors and then falling to give only the single result when the fully redundant
array of 10 sensors was computed. Performance is quantified as the mean of the
71
Chapter 4. Kinematically redundant non-gyro IMUs
−5
−10
0 5 10 15 20 25 30
Compensated estimate of angular velocity error
0.8
p̃′ǫ
0.4 q̃ǫ′
r̃ǫ′
rads −1
−0.4
−0.8
0 5 10 15 20 25 30
Time (s)
RMS error in the x, y and z DOF for linear acceleration, angular acceleration and
angular velocity. These results are shown as box plots in Figure 4.14 The results
compare well with the simulation studies given in Section 4.5. The anomalous
rise apparent in the third graph of Figure 4.14 when all 10 sensors are used is an
artefact of the imperfect alignment of sensors causing the information of one sensor
to dominate the calculation. This would suggest a practical system could make a
confidence estimate of individual sensors and down grade sensors that might be
close to saturation, become mechanically misaligned (or knocked), or operating at
a low signal to noise ratio.
3.5
0.9 1.6
3
0.8 1.4
0.6 1
2
0.5 0.8
0.3 0.4
1
0.2 0.2
0.1 0.5 0
4 5 6 7 8 9 10 4 5 6 7 8 9 10 4 5 6 7 8 9 10
Redundancy (number of tri-axis sensors) Redundancy (number of tri-axis sensors) Redundancy (number of tri-axis sensors)
Figure 4.14: Relationship between error and sensor redundancy for linear
acceleration (left), angular acceleration (middle), angular velocity (right).
Whiskers indicate 1.5× interquartile range.
72
Chapter 4. Kinematically redundant non-gyro IMUs
4.7 Discussion
Due to sensor alignment errors and the constant acceleration of gravity, the
estimate of angular velocity r was subject to a bias error at low velocities
(< 0.5 rads−1 ). Observations of the direction of gravity enabled the orientation
estimation algorithm to compensate for such errors in p and q; however,
compensation is not possible when the axis of r is parallel with the gravity
vector. Such compensation may be achieved if the algorithm is extended to
incorporate information from a magnetometer [59]. An accelerometer array cannot
be used to determine directly the direction of an angular velocity. The angular
velocity estimation algorithm presented can potentially assert an incorrect sign
when angular velocities and accelerations are both of low magnitude. The velocity
estimation method appeared to be robust in dealing with this problem although a
sign switch is possible.
Calibration has been limited to being off-line using a specific calibration dataset
73
Chapter 4. Kinematically redundant non-gyro IMUs
4.8 Conclusion
This work analyses the benefits of using large arrays of accelerometers. Results
include both simulation studies of up to 100 tri-axial sensors and empirical results
of an array of up to 10 MEMS accelerometers. A comparison is given between
simulation and empirical data based on a 2 axis pendulum simulating a knee.
Two methods are reported for the reconstruction of an acceleration state vector
74
Chapter 4. Kinematically redundant non-gyro IMUs
to represent the limb. Both methods allow the pre-computation of the necessary
matrix inverse. The method based on computing a pseudoinverse (using singular
value decomposition) is shown to achieve better accuracy than that based on
combinatorial averaging of the sensors, and at less computational expense.
75
Chapter 4. Kinematically redundant non-gyro IMUs
76
Chapter 5
Characterisation of low-cost
MEMS sensors
5.1 Introduction
The objective of this project is to demonstrate how new applications can be realised
using modern low-cost MEMS gyroscopes, accelerometers and magnetometers. A
key part of achieving this is the development of calibration solutions so that each
sensor is able to provide precise measurements with a determined accuracy. For
example, it is shown in Chapter 6 that calibration of the accelerometer bias and
sensitivity may reduce orientation measurement errors from 8◦ to <0.5◦ .
Low-cost MEMS are typically targeted towards consumer electronics that do not
demand a high level of accuracy. A competitive market means that many modern
devices are supported by the manufacture for only a few years before being
superseded by the next generation. As a result, detailed information about sensor
characteristics and performance is often unavailable. It was therefore necessary for
this project to determine significant sensor characteristics as a first step towards
the development of calibration solutions. This chapter presents the methods and
77
Chapter 5. Characterisation of low-cost MEMS sensors
findings of the sensor characterisations studies that form the basis of the calibration
solutions presented in the next chapter.
The investigations presented in this the chapter focus on the sensor performance
over an operating temperature range or at a fixed temperature. For gyroscopes
and accelerometers, this requires precise controlled motion while simultaneously
regulating the device temperature. Ideally, this would be achieved using an
industrial three-axis thermal chamber. For example, the Ideal Aerosmith 2003HP
shown in Figure 5.1 is able to achieve an accuracy in angular velocities of 0.0001%
and orientations of ±1 arcsecond (∼0.00027◦ ) for temperatures between -65◦C to
+85◦C [115]. Temperature is regulated using electric heating and Liquid Nitrogen
(LN2 ) or Carbon Dioxide (CO2 ) cooling.
Figure 5.1: Ideal Aerosmith 2003HP three-axis thermal chamber of achieve precise
velocities and orientations for temperatures between -65◦C to +85◦C. Image source:
[115]
78
Chapter 5. Characterisation of low-cost MEMS sensors
Figure 5.2: MPU-6050 (top IC) and HMC5883L (bottom IC) mounted on the
x-IMU PCB. The larger IC is 4×4 mm in size.
79
Chapter 5. Characterisation of low-cost MEMS sensors
Standards such as IEEE 1431-2004 [119] and 1293-1998 [120] offer an extensive
list of gyroscope and accelerometer characteristics that may be quantified. The
investigations presented in this chapter are limited to the characteristics indicated
as being most significant in the the device datasheets [121, 122] with a particular
focus on the gyroscope bias. The scope of the investigations is summarised by the
following.
• Sensitivity and bias variation between parts for gyroscope, accelerometer and
magnetometer.
The following sections present the methods and findings of each investigation in
the order above.
The MPU-6050 datasheet specifies the gyroscope bias “variation over temperature”
as ±20◦/s over -40◦C to 80◦C [121]. Further information about characteristic, such
as linearity or repeatability, is not specified in the datasheet or any other available
documentation. The absence of such information about this potentially very large
measurement error exemplifies the necessity for these characterisation studies. The
gyroscope bias thermal response was characterised for eight devices by logging
80
Chapter 5. Characterisation of low-cost MEMS sensors
the gyroscope output when stationary within a thermal chamber to expose the
gyroscope to a controlled temperature range.
Figure 5.3: Thermal chamber internal components including the prototype x-IMU
being tested
During operation, the fan continuously blows on to the thermal mass to circulate
air within the chamber and transfer heat produced by the x-IMU to the thermal
mass. The fan heater blows directly on to the thermal mass and away from the
81
Chapter 5. Characterisation of low-cost MEMS sensors
x-IMU to prevent rapid thermal dynamics. Separate variable power supplies were
used to power the fan and heater and a thermalcouple meter used to monitor
internal temperature. Software was written to log sensor from the x-IMU and
display temperature measurements in real-time to guide manual control of the
heater. Figure 5.4 shows the sealed thermal chamber with separate variable power
supplies and thermocouple meter.
Figure 5.4: Sealed thermal chamber with thermocouple meter and separate
variable power supplies for fan and heater
Eight x-IMUs were tested individually. For each test, the thermal chamber and
x-IMU were initially left in an industrial freezer at -20◦C for >12 hours. The
chamber was then sealed within the freezer before being removed to conduct the
experiment. The heat produced by the fans and x-IMU caused the temperature to
rise of ∼5◦C over the first 5 minutes. The heater power supply was then gradually
increased to maintain a fixed ascent of ∼0.5◦C per minute over 3 hours. The
temperature ascent for each of the eight devices (named ‘A’ to ‘H’) is shown in
Figure 5.5.
82
Chapter 5. Characterisation of low-cost MEMS sensors
A B C D E F G H
80
70
60
50
Temperature (◦ C) 40
30
20
10
−10
−20
0 15 30 45 60 75 90 105 120 135 150 165 180 195
Time (minutes)
Figure 5.5: Measured gyroscope temperature during experiments for each of the
eight devices tested
The measured gyroscope bias for each of the eight devices over a -20◦C to 70◦C
temperature range is shown in Figure 5.6. The three plots indicate the bias for
the x, y and z axes termed bωx , bωy and bωz . The gyroscope output originally in
units of lsb was converted to degrees per second using the typical sensitivity of
16.4 lsb/◦/s specified in the datasheet.
A B C D E F G H
3 3 3
2 2 2
1 1 1
0 0 0
bωx (◦ /s)
bωy (◦ /s)
bωz (◦ /s)
−1 −1 −1
−2 −2 −2
−3 −3 −3
−4 −4 −4
−25 −15 −5 5 15 25 35 45 55 65 75 −25 −15 −5 5 15 25 35 45 55 65 75 −25 −15 −5 5 15 25 35 45 55 65 75
Temperature ( C) ◦
Temperature ( C) ◦
Temperature (◦ C)
Figure 5.6: Gyroscope bias thermal response of eight devices over a -20◦C to 70◦C
temperature range
The bias variation over temperature is within the ±20◦/s range specified in the
83
Chapter 5. Characterisation of low-cost MEMS sensors
datasheet though the relationship was revealed to be both non-linear and unique to
each device. A bias calibration solution that accounts for this variation would need
to incorporate a non-linear model evaluated uniquely for each device. Although
some devices may be approximated by a linear model over a limited temperature
range, this would not be sufficient for all; notably, bωy of device ‘E’ approximates
a parabola centred at 25◦C.
A calibration solution may model the thermal response of the bias as an nth
order polynomial function of temperature. To determine the necessary order, a
polynomial fit was performed for each device for first to fifth order polynomials
and the resultant maximum error of all devices summarised in Table 5.1. Each
polynomial was computed using the MATLAB polyfit function with data points
computed as the mean of each consecutive 2◦C bin.
Table 5.1: Gyroscope bias error assuming an nth order temperature calibration
model
Table 5.1 suggests at least a fourth order polynomial is necessary to account for
the non-linear response over the investigated temperature range. At least five data
points would be required to approximate such a polynomial. In practice, using the
minimum number of data points will result in a less accurate approximation; the
above accuracy was achieved using thousands of data points obtained over a three
hour temperature ascent. Evaluating calibration parameters in this way for a large
number devices would represent a significant practical challenge.
84
Chapter 5. Characterisation of low-cost MEMS sensors
The gyroscope bias thermal response presented in Section 5.2.2 was obtained for a
slow, monotonic temperature ascent. These experiments were extended to observe
the bias behaviour when exposed to a range of thermal dynamics that could be
achieved with the chamber, including: slow cooling, rapid cooling, rapid heating
and thermal cycling. These investigations exposed a hysteresis characteristic of
the gyroscope that was not described within any of the available manufacturer’s
documentation. Figure 5.7 shows the measured temperature during one such
experiment that exemplifies the observations. The annotated markers ‘a’, ‘b’,
‘c’ and ‘d’ indicate key events during the experiment.
C
80
(b)
70
(c)
60
50
Temperature (◦ C)
40
30 (d)
20
10
0
(a)
−10
−20
0 15 30 45 60 75 90 105 120 135 150 165 180 195 210 225
Time (minutes)
The first ∼180 minutes of the experiment, from point ‘a’ to point ‘b’, show the
temperature ascent to ∼70◦C presented in Section 5.2.2. The heater was then
switched off and the temperature slowly descended for ∼15 minutes until point
‘c’. The thermal chamber lid was then removed and the temperature descended
rapidly until stabilising at ∼25◦C at point ‘d’. The bias response is shown as
5.8. The annotated markers correspond to the events previously described within
Figure 5.7. These plots indicate thermal hysteresis in all three axes.
Although this thermal hysteresis is a known characteristic [123, 27], there is little
85
Chapter 5. Characterisation of low-cost MEMS sensors
2.2 0.1 3
(a) (a)
0.05 2.9 (b)
2
2.8
0
1.8 2.7 (c)
−0.05
(d) 2.6
1.6 −0.1
(b) 2.5
1.4 −0.15 (c) 2.4
−0.2
bωx (◦ /s)
bωy (◦ /s)
bωz (◦ /s)
2.3
1.2
−0.25 2.2
1
−0.3 2.1
0.8 2 (d)
−0.35
1.9
0.6 −0.4 (d)
1.8
−0.45
0.4 (c) 1.7
−0.5
1.6
0.2
(a)
(b) −0.55 1.5
0 −0.6 1.4
−25 −15 −5 5 15 25 35 45 55 65 75 −25 −15 −5 5 15 25 35 45 55 65 75 −25 −15 −5 5 15 25 35 45 55 65 75
◦ ◦
Temperature ( C) Temperature ( C) Temperature (◦ C)
Figure 5.8: Thermal response of the gyroscope bias indicating thermal hysteresis
86
Chapter 5. Characterisation of low-cost MEMS sensors
50 50 50
40 40 40
30 30 30
∆bωx (0.001◦ /s)
20 20 20
10 10 10
0 0 0
Figure 5.9: Gyroscope bias deviation for six different orientations relative to
gravity. Box plots show the distribution for 27 devices (whiskers indicate 1.5×
interquartile range).
The deviation of the gyroscope bias for each axis, ∆bωx , ∆bωy and ∆bωz , was
calculated as the difference between the gyroscope bias at a given orientation
from the mean for all orientations. The box plots in Figure 5.9 shows the
distribution of ∆bωx , ∆bωy and ∆bωz for all six orientations of the gyroscope
relative to gravity. The data was scaled to units of degrees per second using
87
Chapter 5. Characterisation of low-cost MEMS sensors
the typical sensitivity specified in the datasheet. Although the specific nature of
the acceleration sensitivity is unknown, the box plots were expected to indicate a
correlation between the direction of gravity and an offset error of with an order
of magnitude of ∼0.1◦/s. Instead, the box plots indicate a mean deviation of
approximately zero and a variance representative of noise. These results suggest
that, contrary to the manufacturer’s specification, the gyroscope bias does not
have a significant sensitivity to static acceleration.
The MPU-6050 datasheet does not provide any information about the gyroscope
bias time-variant random walk and so an investigation was conducted to quantify
this characteristic. Four x-IMUs were sealed in the thermal chamber, each set up
to log the gyroscope output and thermometer measurements to the on-board SD
card at a rate of 512 samples per second. The thermal chamber was intended only
to isolate the internal temperature from the ambient environment; the internal
fan and heater were not used. Data was logged over a 48 hour period and then
cropped to select a period with minimal temperature variations.
0.1
∆T (◦ C)
−0.1
0 1 2 3 4 5 6 7 8 9 10
X Y Z
150
100
∆bω (◦ /hour)
50
−50
−100
0 1 2 3 4 5 6 7 8 9 10
Time (hours)
88
Chapter 5. Characterisation of low-cost MEMS sensors
Figure 5.10 shows the measured temperature variation, ∆T , and gyroscope bias
variation, ∆bω , for one device over a 10 hour period with a temperature variation of
±0.025◦C. The measurement were average for each 30 second window to attenuate
the high-frequency noise. Bias measurements were scaled to units of degrees per
hour using the typical sensitivity specified in the datasheet. The results shown
above are representative of the observations made for all devices and indicate a
bias variation that does not correlate to temperature. This random walk resulted
in a peak error of 190◦ /hour (∼0.053◦/s) for the four devices over the 10 hour
period.
Plotting the gyroscope bias variation over time clearly demonstrates random walk
but does not provide an effective means for quantifying this characteristic. Allan
variance is the accepted standard for evaluating gyroscope bias stability [125].
Stockwell provides a concise introduction to the implementation of the method in
a Crossbow white paper [126].
The Allan variance, AVAR, is calculated by splitting the dataset into multiple
consecutive windows of time and evaluating the variance in the gyroscope bias
averaged for each window of time. The AVAR is computed for a range of different
window sizes, termed τ . The maximum value of τ must not split the dataset into
more than nine windows to maintain statistical significance [125]. Equation 5.1
describes the calculation where AVR2 (τ ) is the squared Allan variance for period
τ , n is the number of windows created by the given value of τ , and y(τ )i is the
averaged gyroscope bias calculated for the it h window.
1 X 2
AVR2 (τ ) = y(τ )i+1 − y(τ )i (5.1)
2(n − 1) i
89
Chapter 5. Characterisation of low-cost MEMS sensors
Figure 5.11 shows the AVAR for each axis calculated using values of τ between 0.01
and 1000 seconds. The linear descent created on a logarithmic plot for small values
of τ indicates a reduction in variance due to the averaging of mean-zero random
noise. The AVAR reaches a minimum at which point the length of τ is sufficiently
long to incorporate the changes due to the random walk and the variance beings
to increase again. The stability of the gyroscope bias is quantified as the minimum
Allan variance and corresponding value of τ [126]. In this case the gyroscope bias
stability is quantified as ∼4◦ /hour for τ = 100 seconds.
X Y Z
3
10
2
10
AVAR (◦ /hour)
1
10
0
10
−2 −1 0 1 2 3
10 10 10 10 10 10
Averaging time, τ (seconds)
Figure 5.11: Allan variance plot of the gyroscope bias random walk
90
Chapter 5. Characterisation of low-cost MEMS sensors
most of the specified quantities are ambiguous. The gyroscope bias investigations
in Section 5.2 have already revealed the potential significance of such ambiguities.
This provides a motivation to investigate the thermal response of the bias and
sensitivity for all three types of sensor.
The sensitivity and bias and may be evaluated for each sensor by exposing each
axis to both a positive and negative reference. For example, this could be a ±1
g reference for each accelerometer axis. The axis bias may then be calculated as
Equation 5.2 where b is the axis bias in lsb, u+r is the axis output when exposed
to the positive reference and u−r is the axis output when exposed to the negative
reference.
u+r + u−r
b= (5.2)
2
The axis sensitivity may be calculated as Equation 5.3 where s is the axis bias
in lsb per units of the reference signal, u+r is the axis output when exposed to
the positive reference and u−r is the axis output when exposed to the negative
reference, and r is the magnitude of the reference signal.
|u+r | + |u−r |
s= (5.3)
2r
91
Chapter 5. Characterisation of low-cost MEMS sensors
The three-axis thermal chamber was created from a polystyrene box and a cube
structure designed to slot within the internal volume to mount four battery
powered prototype x-IMUs in alignment with the outer dimensions of the chamber.
A 1 kg copper plate with aluminium fins and battery powered fan were mounted
internally to blow air directly on to the x-IMUs. The chamber can be orientated at
six orientations on a flat surface to expose each accelerometer axis to +1 g and -1 g;
and at six orientations relative to the rotational axis of a rate table to expose each
gyroscope axis to a positive and negative reference of angular velocity. A Technics
SL-1200MK2 turntable was used to provide a repeatable reference of 200◦/s (331/3
Rotations Per Minute (RPM)). Details on the evaluation of this turntable are
discussed in Chapter 6. Figure 5.12 shows the components that make up the
three-axis thermal chamber.
Each x-IMUs was set up to log sensor data to on-board SD cards at a rate of
256 samples per second while simultaneously streaming the data via Bluetooth
to be viewed in real-time. The HMC5883L was mounted only a few millimetres
from the MPU-6050 and so assumed to be at the same temperature indicated
by the MPU-6050 internal thermometer. The thermal chamber was stored in
a -30◦C industrial freezer for >12 hours prior to the investigation. The x-IMUs
and fan were switched on and thermal chamber sealed before being removed
from the freezer. The thermal chamber was then continuously cycled through
six static orientations relative to gravity and six orientations on the turntable
92
Chapter 5. Characterisation of low-cost MEMS sensors
while the internal temperature increased with the heat generated by the internal
electronics from approximately -25◦C to 40◦C over 2.5 hours. Figure 5.13 shows
this temperature ascent for each of the four devices: ‘A’, ‘B’, ‘C’ and ‘D’. The
necessity for the chamber to be wireless prevented the use of active heating which
both limited the maximum temperature temperature and prevented control of the
rate of ascent. Nevertheless, the temperature range achieved was sufficient to
investigate the thermal response of each of the sensors.
A MATLAB script was written to process the logged data and extract the periods
at which the each gyroscope axis was detected as being at the reference speeds
of -200◦/s and +200◦/s. The average gyroscope output was computed for discrete
93
Chapter 5. Characterisation of low-cost MEMS sensors
A B C D
40
35
30
25
20
Temperature (◦ C)
15
10
5
0
−5
−10
−15
−20
−25
0 15 30 45 60 75 90 105 120 135 150
Time (minutes)
Figure 5.13: Measured temperature ascent of sensors within the three-axis thermal
chamber over 2.5 hours
temperature intervals of 2◦C to yield a value of u+r and u−r for each axis at each
temperature interval. The corresponding sensitivities, sωx , sωy and sωz ; and biases,
bωx , bωy and bωz , were calculated using Equations 5.3 and 5.2 to yield the results
shown in Figure 5.14. The bias was scaled to units of degrees per second using the
typical sensitivity specified in the datasheet.
A B C D
bωy (◦ /s)
bωz (◦ /s)
1.5 1
−2
1 0
−3
0.5
−4 −1
0
−0.5 −5 −2
−20 −10 0 10 20 30 40 −20 −10 0 10 20 30 40 −20 −10 0 10 20 30 40
Temperature ( C) ◦
Temperature ( C)◦
Temperature (◦ C)
Figure 5.14: Thermal response of gyroscope sensitivity and bias from -15◦C to
40◦C
The gyroscope biases bωx , bωy and bωz for the four devices has already been
investigated in greater detail in Section 5.2. The thermal responses shown in Figure
94
Chapter 5. Characterisation of low-cost MEMS sensors
5.14 match those previously demonstrated for devices ‘A’, ‘B’, ‘C’ and ‘D’. The
gyroscope sensitivity temperature coefficient was calculated for the four devices
as the ratio of the best-fit gradient to the mean sensitivity over the temperature
range. The results are summarised in Table 5.3.
As with the gyroscope analysis, accelerometer data was extracted to yield averaged
measurements of u+r and u−r for each axis at each temperature interval. The
corresponding sensitivities sax , say and saz , and biases bax , bay and baz , were
calculated using Equations 5.3 and 5.2 to yield the results shown in Figure 5.15.
The bias was scaled to units of mg using the typical sensitivity specified in the
datasheet.
The accelerometer sensitivity temperature coefficient was calculated for the four
devices as the ratio of the best-fit gradient to the mean magnitude of the measured
sensitivity over the temperature range. The results are summarised in Table 5.4.
95
Chapter 5. Characterisation of low-cost MEMS sensors
A B C D
say (lsb/g)
saz (lsb/g)
2060 2060 2070
bay (mg)
baz (mg)
14 8 50
12 6 0
10 4 −50
8 2 −100
6 0 −150
−20 −10 0 10 20 30 40 −20 −10 0 10 20 30 40 −20 −10 0 10 20 30 40
Temperature (◦ C) Temperature (◦ C) Temperature (◦ C)
Figure 5.15: Thermal response of accelerometer sensitivity and bias from -15◦C to
40◦C
96
Chapter 5. Characterisation of low-cost MEMS sensors
The magnetometer sensitivity temperature coefficient was calculated for the four
devices as the ratio of the best-fit gradient to the mean magnitude of the measured
97
Chapter 5. Characterisation of low-cost MEMS sensors
A B C D
smy (lsb/Ga)
smz (lsb/Ga)
500 475 500
bmy (lsb)
bmz (lsb)
8.5 9.5
8 9 −0.5
7.5 8.5
−1
7 8
−1.5
6.5 7.5
6 7 −2
−20 −10 0 10 20 30 40 −20 −10 0 10 20 30 40 −20 −10 0 10 20 30 40
Temperature (◦ C) Temperature (◦ C) Temperature (◦ C)
Figure 5.16: Thermal response of magnetometer sensitivity and bias from -15◦C
to 40◦C
98
Chapter 5. Characterisation of low-cost MEMS sensors
bmx , bmy and bmz result in the increased variance of the bias temperature coefficient
relative to the sensitivity temperature coefficients. In the case of bmz , calculations
approach divisions by zero and so cannot be expected to provide meaningful results.
The absence of a bias variation means that the magnetometer thermal response
can be accounted for by scaling all measurements by approximately -0.3%/◦C,
though the variation in above results indicate that the specific coefficient would
need to be evaluated for each device. An alternative solution would be to use
the magnetometer internal field generators to calculate the sensitivity at discrete
time or temperature intervals; this technique is suggested in the magnetometer
datasheet.
99
Chapter 5. Characterisation of low-cost MEMS sensors
Figure 5.17 shows the normalised distribution of the gyroscope axis sensitivities
sωx , sωy and sωz and biases bωx , bωy and bωz measured for the 150 MPU-6050 parts.
The biases have been scaled using the typical sensitivity specified in the datasheet.
Normalised distribution
0 0 0
15.8 16 16.2 16.4 16.6 16.8 15.8 16 16.2 16.4 16.6 16.8 15.8 16 16.2 16.4 16.6 16.8
◦ ◦
sωx (lsb/ /s) sωy (lsb/ /s) sωz (lsb/◦ /s)
0.35 0.35 0.35
Normalised distribution
Normalised distribution
Normalised distribution
The datasheet specifies that the gyroscope sensitivity of 16.4 lsb/◦/s ±3% at 25◦C;
corresponding to a range of 15.9 lsb/◦/s to 16.9 lsb/◦/s. The histograms confirm
a mean of ∼16.4 lsb/◦/s with no outliers exceeding these tolerances. The bias is
specified as being within the range ±20◦/s at 25◦C. The lower histograms of Figure
5.17 reveal that the true distribution is more limited. Whereas the sensitivity
100
Chapter 5. Characterisation of low-cost MEMS sensors
distributions are centred on the specified nominal value, the mean bias for each
axis is not centred on zero. For example, the mean value of bωx is approximately
-3◦/s. Although this mean suggests a more accurate nominal value, this may be
specific to the batch.
Figure 5.18 shows the normalised distribution of the accelerometer axis sensitivities
sax , say and saz and biases bax , bay and baz measured for the 150 MPU-6050 parts.
The biases have been scaled using the typical sensitivity specified in the datasheet.
Normalised distribution
Normalised distribution
Normalised distribution
Normalised distribution
101
Chapter 5. Characterisation of low-cost MEMS sensors
5.6.3 Magnetometer
Normalised distribution
Normalised distribution
Normalised distribution
Normalised distribution
Normalised distribution
Normalised distribution
102
Chapter 5. Characterisation of low-cost MEMS sensors
5.7 Conclusions
103
Chapter 5. Characterisation of low-cost MEMS sensors
104
Chapter 6
6.1 Introduction
105
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
The calibration solutions presented in this chapter account for the following
characteristics.
• Sensor misalignment
Notable characteristics not accounted for by the calibration models include the
gyroscope bias acceleration sensitivity and the thermal response of the sensors.
Previous investigations indicated could not demonstrate a correlation between the
gyroscope bias and static acceleration. Thermal response investigations reveal a
significant sensitivity in the magnetometer that may be accounted for using the
manufactures recommended technique; and of the gyroscope bias which can be
accounted for by the AHRS algorithm.
The next section summaries the commercial IMUs chosen as benchmarks. Section
6.4 describes the calibration solution for the gyroscope. Section 6.5 describes
the calibration solution for the accelerometer along side an empirical equation
to contrast uncalibrated an calibrated AHRS performance. Section 6.6 describes
three different magnetometer calibration solutions and demonstrates the AHRS
performance of each solution for both ‘typical’ and ‘extreme’ magnetic distortion
scenarios. Section 6.7 presents an optimisation that can be made to simplify
the sensor alignment arithmetic, and significantly reduce the computation load
in application. The chapter concludes with a summary of proposed calibration
solutions.
106
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
Four commercial IMUs were chosen to benchmark the performance of the proposed
calibration solutions. Each product appeared on the market within the last 5 years
and is priced less than 3000 USD. The chosen devices, shown Figure 6.1, were
the MicroStrain 3DM-GX3-25, Xsens MTw, VectorNav VN-100 and CH Robotics
UM6. MicroStrain and Xsens are recognised as long-established and leading
manufacturers whereas VectorNav and CH Robotics represent new manufacturers
whom have taken advantage of recent advances in MEMS technology to offer
lower-end devices at a significantly lower price.
Table 6.1 summarises the calibrated performance of each IMU as specified in the
manufacturer’s documentation [108, 128, 109, 52, 129]. Prices and release dates
were either obtained from the manufacturer’s website or confirmed by contacting
the manufacturer directly. The MPU-6050 [121] and HMC-5883L [122] have
also been included in the table for comparison, these devices function only as
an uncalibrated gyroscope, accelerometer and magnetometer and do not have an
associated AHRS performance.
There is a clear correlation between the AHRS performance and price with the
more expensive devices providing greater accuracy. Both MicroStrain and Xsens
specify a dynamic error. However, the meaning of this quantity is ambiguous and it
cannot be used for the purposes of benchmarking; Xsens specify the dynamic error
along side the disclaimer: “May depend on type of motion” [128]. The proposed
107
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
MPU-6050 /
3DM-GX3-25 MTw VN-100 UM6
HMC5883L
Price $2550 $2990 $500 $199 <$10
Year of release (circa) 2009 2011 2008 2009 2010
AHRS performance
Static error (pitch/roll) ±0.5◦ <0.5◦ 0.5◦ <2◦ N/A
◦ ◦
Static error (heading) ± 0.5 1 2◦ < ◦
5 N/A
Dynamic error (pitch/roll) ±2◦ 2◦ - - N/A
Dynamic error (heading) ±2◦ 2◦ - - N/A
Operating temperature -40 to 70◦C 0 to 55◦C 25◦C - N/A
Gyroscope
Range (◦/s) 300 1200 2000 2000 2000
Linearity (%)
√ 0.03 0.1 0.1 0.2 0.2
Noise (◦/s/ Hz) 0.03 0.05 0.005 0.03 0.005
Accelerometer
Range (g) 5 16 16 2 16
Linearity (%)
√ 0.1 0.2 0.5 - 0.5
Noise (μg/ Hz) 80 0.0003 400 - 400
Magnetometer
Range (Gauss) 2.5 1.5 2.5 2.5 8.1
Linearity (%) √ 0.4 0.2 0.1 - 0.1
Noise (μGauss/ Hz) 100 150 140 - -
calibration solutions were therefore benchmarked against the static errors only.
The specified heading accuracy of each device also represents a potential ambiguity.
It can be shown (see Appendix C) that a calibration offset error corresponding to
a 0.5◦ heading error at the equator would result in a 1.5◦ heading error in the
UK due to a varying inclination of the Earth’s magnetic field. This should be
considered when benchmarking the heading performance.
108
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
Figure 6.2: Calibration cube with nine prototype x-BIMUs held in alignment with
the cube axes
The calibration cube was laser cut from 3 mm Medium-Density Fibreboard (MDF)
and assembled through press-fit joints to avoid metal components that may
interfere with magnetic calibration. The structure is comprised of three horizontal
levels separated by vertical spacers. The middle and top levels feature grids of
mounting holes and cut-outs for securing the sensors and additional equipment
such as batteries. The bottom level is an empty frame to provide rigidity to the
structure. Sensors mounted on the middle level are held in alignment by a laser cut
MDF plate with cut-outs matching the dimensions of the x-BIMU. Each corner of
the cube functions as a socket to ensure alignment when the cube is mounted on the
109
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
xs , ra xc
r
γ
zs
zc
ys
yc
Figure 6.3: A misalignment between the calibration cube and sensor x axes of
angle γ prevents the sensor axis from being exposed to the full magnitude of the
reference signal r.
The relationship between the reference signal, r, and the attenuated reference
signal, ra , is defined by Equation 6.1 where γ is the angular misalignment between
the sensor axis and the direction of the reference signal. ra can be approximated
as a truncated Taylor series.
γ2
ra = r cos γ ≈ r 1 − (6.1)
2
The cosine relationship means that the error in the reference signal has greater
tolerance to misalignment errors of a small magnitude. The attenuation of the
1
https://github.com/xioTechnologies/Calibration-Cube
110
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
reference signal for misalignments between -10◦ and +10◦ is shown in Figure 6.4.
A sensor calibrated against an attenuated reference would incorporate a sensitivity
error equal to this attenuation if the misalignment is not accounted for within the
calibration model.
Reference signal attention (%)
0
−0.2
−0.4
−0.6
−0.8
−1
−1.2
−1.4
−1.6
−10 −8 −6 −4 −2 0 2 4 6 8 10
Alignment error, γ (degrees)
Figure 6.4: Relationship between calibration cube alignment error and attenuation
of reference signal
The calibration cube protractor is used in conjunction with the calibration cube
to evaluate the static pitch/roll and heading error achieved by an accelerometer
and magnetometer. The protractor is not required for the calibration. Figure 6.5
shows the calibration cube mounted on the protractor.
The protractor is comprised of a 290 mm diameter cog with 72 teeth, and a base
plate incorporating a socket for the cog. The cog can be mounted in the socket
at 72 different angular positions corresponding to a 360◦ rotation in 5◦ steps. The
cog has four sockets that mate with any four corners of the calibration cube. The
protractor may therefore be used to rotate the cube around all three axes.
111
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
Figure 6.5: Calibration cube mounted on the protractor (left) and close up of cog
teeth and engraved angular increments (right)
The manufacture process was essential in achieving the accuracy of the protractor.
The cog and socket were laser cut from a single sheet of 5 mm acrylic plastic and
the cog flipped within the socket so that the natural chamfers created by the laser
cutter formed parallel surfaces. Theoretically, the protractor accuracy would be
determined by the backlash between the cog and the socket. In practise, the finite
accuracy of the laser cutter results in an elliptical cog and socket that achieve either
a flush fit or a ‘press fit’ as the cog is rotated. This potentially improves accuracy
by distributing geometric errors between the 72 teeth. An accelerometer was
demonstrated as achieving 0.1◦ static pitch/roll accuracy through a 360◦ rotation
using the protractor in Section 6.5.4; it can therefore be inferred that the protractor
accuracy is <0.1◦ .
The gyroscope calibration model is described by Equations 6.2 and 6.3 where ω
is the calibrated gyroscope measurement in degrees per second, C
Ω R is a rotation
matrix describing the alignment of the gyroscope relative to the calibrated frame,
uω is the uncalibrated gyroscope output in lsb, Sω is the gyroscope sensitivity
in lsb per degrees per second, and bω is gyroscope bias in lsb. Sω is a diagonal
112
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
matrix as the calibration model does not account for cross-axis coupling or axis
non-orthogonality.
−1
ω=C
Ω RSω (uω − bω ) (6.2)
−1
ωx rω11 rω12 rω13 sωx 0 0 uωx bωx
ω = r
y ω21 rω22 rω23 0 sωy 0 uωy − bωy (6.3)
ωz rω31 rω32 rω33 0 0 sωz uωz bωz
1. Calculate the bias, bω , as the the mean gyroscope output, uω , while the
gyroscope is stationary.
3. Calculate each axis sensitivity, sωx , sωy and sωz , as the average magnitude
of the axis output when exposed the +ω and −ω reference.
C
4. Calculate the misalignment, Ω R, using the previously determined values
bω and Sω and three orthogonal quantities of uω within the dataset that
represent each axis begin exposed to a +ω reference.
113
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
The gyroscope bias parameters bωx , bωy and bωz are calculated as the mean
uncalibrated gyroscope output over a period of several seconds while the device is
stationary. The gyroscope bias is known to vary over time and with temperature
though this is not accounted for in the calibration model. It is assumed that an
AHRS algorithm will provide a means of updating bω during operation.
The dataset is collected by orientating the calibration cube on the turntable so that
so that each axis is aligned with the vertical rotational axis and pointing either up
or down; exposing each axis to a reference of +ω and −ω. This is represented by
the six orientations shown in Figure 6.6. The mean uncalibrated gyroscope output,
114
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
uω , is taken over a period of several seconds for each orientation to yield a dataset
of six vector quantities. Each pair of measurements representing +ω and −ω for a
given axis should be made contiguously with the calibration cube rotated around
only a single axis relative to the turntable between the two orientations. This
ensures that unexpected time-variant or temperature deviations are minimised
and that a misalignment between the plater and rotational axis results in an equal
attenuation of the +ω and −ω reference.
Figure 6.6: The gyroscope calibration dataset is represented by the six orientations
of the calibration cube on the turntable to expose each axis to a +ω and −ω
reference
.
The gyroscope axis sensitivities sωx , sωy and sωz are calculated as the average
magnitude of the axis output when exposed to the +ω and −ω reference. This is
described by Equation 6.4 where sω is the axis sensitivity in lsb per degrees per
second, u+ω and u−ω are the axis output in lsb when exposed to the +ω and −ω
reference respectively, and ω is the magnitude of the reference angular velocity
115
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
known to be 200◦/s.
|u+ω | + |u−ω
sω = (6.4)
2ω
using the values of the Sω and bω previously determined in Sections 6.4.1 and 6.4.3
and the three orthogonal vector quantities of uω within the calibration dataset that
C
represent each axis being exposed to a positive reference of angular velocity. ΩR is
initially approximated as C
G R̃ as described by Equation 6.5 where uωx , uωy and uωz
are the gyroscope outputs when the positive reference is applied to the x, y and
z axes respectively. The three columns of C
G R̃ are constructed as each measured
C
G R̃ can only be assumed to be an approximation of a rotation matrix because
C
measurement errors may prevent it from being orthogonal. ΩR is obtained from
C
G R̃ by computing the best-fit rotation matrix as described in Chapter 2.
The accelerometer calibration model is described by Equations 6.6 and 6.7 where
C
a is the calibrated accelerometer measurement in g, AR is a rotation matrix
describing the alignment of the accelerometer relative to the calibrated frame,
ua is the uncalibrated accelerometer output in lsb, Sa is the accelerometer
sensitivity in lsb per g, and ba is accelerometer bias in lsb. Sa is a diagonal
116
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
matrix as the calibration model does not account for cross-axis coupling or axis
non-orthogonality.
−1
a=C
A RSa (ua − ba ) (6.6)
−1
ax ra11 ra12 ra13 sax 0 0 uax bax
a = r
y a21 ra22 ra23 0 say 0 uay − bay (6.7)
az ra31 ra32 ra33 0 0 saz uaz baz
2. Calculate each axis bias, bax , bay and baz , as the average of the axis output
when exposed the +1 g and −1 g reference.
3. Calculate each axis sensitivity, sax , say and saz , as the average magnitude of
the axis output when exposed the +1 g and −1 g reference.
C
4. Calculate the misalignment, A R, using the previously determined values
Sa and ba and the three orthogonal quantities of ua within the dataset
representing each axis begin exposed to a +1 g.
The following sections describe the accelerometer calibration process in detail and
present a practical evaluation of the calibrated accelerometer performance. The
calibrated performance in contrasted with the uncalibrated performance and that
of the commercially available IMUs.
117
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
118
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
u+g + u−g
ba = (6.8)
2
The sensitivity of each axis is calculated as the average magnitude of the axis
output when exposed to the +1 g and −1 g reference. This is described by Equation
6.9 where sa is the axis sensitivity in lsb per g; and u+g and u+g are the axis output
in lsb when exposed to the +1 g and −1 g reference respectively.
|u+g | + |u−g |
sa = (6.9)
2g
calculated using the values of the Sa and ba previously determined in Section 6.5.2
and the three orthogonal vector quantities of ua within the calibration dataset that
represent each axis begin exposed to a +1 g reference. C
A R is initially approximated
119
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
T
−1
kSa (uax − ba )k
C −1
A R̃ = kSa (uay − ba )k
(6.10)
−1
kSa (uaz − ba )k
C
A R̃ can only be assumed to be an approximation of a rotation matrix because
C
measurement errors may prevent it from being orthogonal. AR is obtained from
C
A R̃ by computing the best-fit rotation matrix as described in Chapter 2.
120
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
The calibration cube protractor was used to rotate the cube 360◦ around each
axis in 5◦ steps. At each step, the mean accelerometer output was obtained
for approximately 2000 samples over 30 seconds to yield both uncalibrated and
calibrated measurements for each accelerometer. The corresponding values of θx ,
θy or θz for each step would then be computed and compared with the angle
indicated by the protractor to yield the error in the angular measurement of
inclination provided by an accelerometer. Figure 6.8 shows the calibration cube
protractor being used to collect the accelerometer evaluation dataset.
A misalignment between the zero position of the protractor and the direction of
gravity introduced an offset to each accelerometer’s measurement of inclination.
This offset was removed by subtracting the mean angular error of all nine
121
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
accelerometers from each. Figure 6.9 shows the performance of the uncalibrated
accelerometers. The top three plots indicate the error in each accelerometer’s
measurement of inclination for a 360◦ rotation around each of accelerometer axes.
The distribution of these errors are shown for each accelerometer as bottom three
box plots. Rotations around the z axis can be seen to result in the lowest
error. This is expected because rotations around the x axis or y axis incorporate
measurements from the accelerometer z axis which is indicated as providing a
greater initial offset error in the MPU-6050 data sheet [121]. This is due to the
physical structure of the MEM sensor. The worst case inclination measurement
errors are ±6.5◦ , 6.9◦ and ±2.6◦ for rotations around the x, y and z axes respectively.
These errors are three times greater than that of the lowest performing commercial
IMU and demonstrate the need for calibration.
A B C D E F G H I
10 10 10
8 8 8
θx error (degrees)
θy error (degrees)
θz error (degrees)
6 6 6
4 4 4
2 2 2
0 0 0
−2 −2 −2
−4 −4 −4
−6 −6 −6
−8 −8 −8
−10 −10 −10
0 90 180 270 360 0 90 180 270 360 0 90 180 270 360
Inclination (degrees) Inclination (degrees) Inclination (degrees)
10 10 10
8 8 8
θx error (degrees)
θx error (degrees)
θx error (degrees)
6 6 6
4 4 4
2 2 2
0 0 0
−2 −2 −2
−4 −4 −4
−6 −6 −6
−8 −8 −8
−10 −10 −10
A B C D E F G H I A B C D E F G H I A B C D E F G H I
Accelerometer Accelerometer Accelerometer
Figure 6.10 shows the performance of the calibrated accelerometers. The worst case
inclination measurement errors are ±0.35◦ , ±0.44◦ and ±0.15◦ for rotations around
the x, y and z axes respectively. This performance matches that of the highest
performing commercial IMU; the 3DM-GX3-25 which is indicated as achieving an
122
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
accuracy of ±0.5◦ .
A B C D E F G H I
θy error (degrees)
θz error (degrees)
0.3 0.3 0.3
0.2 0.2 0.2
0.1 0.1 0.1
0 0 0
−0.1 −0.1 −0.1
−0.2 −0.2 −0.2
−0.3 −0.3 −0.3
−0.4 −0.4 −0.4
−0.5 −0.5 −0.5
0 90 180 270 360 0 90 180 270 360 0 90 180 270 360
Inclination (degrees) Inclination (degrees) Inclination (degrees)
0.5 0.5 0.5
0.4 0.4 0.4
θx error (degrees)
θx error (degrees)
θx error (degrees)
0.3 0.3 0.3
0.2 0.2 0.2
0.1 0.1 0.1
0 0 0
−0.1 −0.1 −0.1
−0.2 −0.2 −0.2
−0.3 −0.3 −0.3
−0.4 −0.4 −0.4
−0.5 −0.5 −0.5
A B C D E F G H I A B C D E F G H I A B C D E F G H I
Accelerometer Accelerometer Accelerometer
m = Su − h (6.15)
mx sxx sxy sxz ux hx
m = s
y yx syy syz uy − hy (6.16)
mz szx szy szz uz hz
123
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
hard-iron biases, h, are proposed in the following sections. Each solution represents
a different level of complexity with associated practical implications.
Figure 6.11: Magnetometer in the ‘typical’ (left) and ‘worst case’ (right) magnetic
distortion arrangements. The magnetometer chip is indicated by the orange circle.
124
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
for 2D rotations around each calibrated axis would form three circles on the
surface of the sphere, each concentric with the axis of rotations. However, soft-iron
distortions will distort the sphere into an ellipsoid and hard-iron biases will shift
the centre away from the origin. This visual method of evaluation is demonstrated
in some magnetometer manufacturer application notes [133].
Y X Y X
125
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
T
a m
vx vy vz = × (6.17)
kak kmk
The calibration cube protractor was used to rotate the cube 360◦ around
each axis in 5◦ steps and the mean accelerometer and magnetometer outputs
obtained for approximately 2000 samples over 30 seconds at each step. The
corresponding values of θx , θy or θz for each step were then be computed and
compared with the angle indicated by the protractor to yield the error in the
angular measurement of heading provided by an magnetometer. A misalignment
126
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
between the protractor and magnetic North introduced an offset error to each
magnetometer’s measurement of heading. This offset was removed by subtracting
the mean error of the nine measurements of heading for a given axis of rotation.
Figure 6.13 shows the calibration cube protractor being used to collect this
magnetometer evaluation data.
The top three plots of 6.14 show the error in each uncalibrated magnetometer’s
measurement of heading for a 360◦ rotation around each of axes for the ‘typical’
magnetic distortion scenario. The distribution of these errors are shown for each
magnetometer as bottom three box plots. The distribution of the errors between
±180 seen in Figure 6.14 demonstrates that magnetometer calibration is essential
for the calculation of heading. The error in uncalibrated measurements of heading
obtained for the ‘worst case’ dataset were all uniformly distributed between ±180.
127
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
A B C D E F G H I
θy error (degrees)
θz error (degrees)
90 90 90
45 45 45
0 0 0
−45 −45 −45
−90 −90 −90
−135 −135 −135
−180 −180 −180
0 90 180 270 360 0 90 180 270 360 0 90 180 270 360
Heading (degrees) Heading (degrees) Heading (degrees)
180 180 180
135 135 135
θx error (degrees)
θx error (degrees)
θx error (degrees)
90 90 90
45 45 45
0 0 0
−45 −45 −45
−90 −90 −90
−135 −135 −135
−180 −180 −180
A B C D E F G H I A B C D E F G H I A B C D E F G H I
Magnetometer Magnetometer Magnetometer
Magnetometer offset calibration accounts for only the hard-iron distortions. The
soft-iron matrix, S is assumed to be an identity so that the magnetometer
calibration model described by Equation 6.16 simplifies to Equation 6.21. The
calibrated magnetometer measurement and hard-iron bias are both represented in
units of lsb for this model.
mx ux hx
m = u − h (6.21)
y y y
mz uz hz
128
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
is the magnitude of the magnetic field in lsb. This can be expressed in vector form
as shown by Equation 6.23.
hx
h2
x
hy
1 −2mz 2 2 2 2
−2mx 1 −2my 2 = v − mx − my − mz (6.23)
hy
h
z
h2z
Mh = v (6.25)
The hard-iron biases can therefore the calculated as the solution to Equation 6.25
for a specified value of v and numerical dataset of n measurements. The value of
h may be obtained using the pseudoinverse of M as described by Equation 6.31
for a calibration dataset of dataset of n > 6. Elements hx , hy and hz of vector h
129
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
h = M +v
−1 (6.26)
T T
= M M M v
Y X Y X
130
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
A B C D E F G H I
10 10 10
8 8 8
θx error (degrees)
θy error (degrees)
θz error (degrees)
6 6 6
4 4 4
2 2 2
0 0 0
−2 −2 −2
−4 −4 −4
−6 −6 −6
−8 −8 −8
−10 −10 −10
0 90 180 270 360 0 90 180 270 360 0 90 180 270 360
Heading (degrees) Heading (degrees) Heading (degrees)
10 10 10
8 8 8
θx error (degrees)
θx error (degrees)
θx error (degrees)
6 6 6
4 4 4
2 2 2
0 0 0
−2 −2 −2
−4 −4 −4
−6 −6 −6
−8 −8 −8
−10 −10 −10
A B C D E F G H I A B C D E F G H I A B C D E F G H I
Magnetometer Magnetometer Magnetometer
Figure 6.16: Heading measurement error provided by the nine magnetometers for
the ‘typical’ magnetic distortion scenario after offset calibration.
Figure 6.17 summarises the heading measurement error provided by the nine
magnetometers for the ‘worst case’ magnetic distortion scenario after offset
calibration. The simplified offset calibration model is clearly insufficient for this
scenario. The distortion of a sphere into the ellipsoid seen in Figure 6.15 causes
the x axis to fall outside of the circle of data points describing θx ; consequently it
is not possible to resolve a heading and θx is distributed between ±180◦ .
Ellipsoid calibration requires the fitting an ellipsoid to the data points provided by
the uncalibrated magnetometer when rotated through many different orientation
in a homogeneous magnetic field. The soft-iron matrix, S, and hard-iron vector,
h, may then be obtained as the geometric solution that equated the ellipsoid with
a sphere. Ellipsoid fitting is the basis for many existing magnetometer calibration
solutions [134]. The magnetometer calibration solution presented here builds on
131
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
A B C D E F G H I
θy error (degrees)
θz error (degrees)
90 90 90
45 45 45
0 0 0
−45 −45 −45
−90 −90 −90
−135 −135 −135
−180 −180 −180
0 90 180 270 360 0 90 180 270 360 0 90 180 270 360
Heading (degrees) Heading (degrees) Heading (degrees)
180 180 180
135 135 135
θx error (degrees)
θx error (degrees)
θx error (degrees)
90 90 90
45 45 45
0 0 0
−45 −45 −45
−90 −90 −90
−135 −135 −135
−180 −180 −180
A B C D E F G H I A B C D E F G H I A B C D E F G H I
Magnetometer Magnetometer Magnetometer
Figure 6.17: Heading measurement error provided by the nine magnetometers for
the ‘worst case’ magnetic distortion scenario after offset calibration.
an ellipsoid fitting algorithm shared by Yury Petrov on the MATLAB file exchange
[135] that evaluates the ellipsoid origin, radii and principle axes for a given dataset.
The calibration parameters S and h are derived from these parameters. The
ellipsoid fitting algorithms describes an ellipsoid using with the polynomial shown
as Equations 6.27 and 6.28.
132
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
v1
v
2
v3
v4
x2 y 2 z 2 2xy 2xz 2yz 2x 2y 2z
v5 = 1 (6.28)
v6
v7
v
8
v9
A
B
C
2 2 2
x1 y1 z1 2x1 y1 2x1 z1 2y1 z1 2x1 2y1 2z1
1
D
x2 y 2 z 2 2x2 y2 2x2 z2 2y2 z2 2x2 2y2 2z2 1
2 2 2
. .. .. .. .. .. .. .. .. E = . (6.29)
.. .
. . . . . . . ..
F
x2n yn2 zn2 2xn yn 2xn zn 2yn zn 2xn 2yn 2zn
1
G
H
I
Mh = v (6.30)
The ellipsoid can therefore the calculated as the solution to Equation 6.30 for a
numerical dataset of n measurements. The value of h may be obtained using the
133
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
h = M +v
−1 (6.31)
T T
= M M M v
The radii of the ellipsoid can be determined from the Eigenvalues of Rc as described
by Equation 6.35. Tc describes the T having been translated to the origin. λ1 to
λ2 are the first three Eigenvalue of Rc .
T
t11 t12 t13 t14 1 0 0 1 0 0
t21 t22 t23 t24 1 0 c 0 1 0 c
0
Tc = = T (6.34)
t t32 t33 t34 0 1 0 0 1
31 0
t41 t42 t43 t44 0 0 0 1 0 0 0 1
t11 t12 t13
1
Rc =
t21 t22 t23 t44 (6.35)
t31 t32 t33
134
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
R = λ1 λ2 λ3 (6.36)
The above ellipsoid fitting algorithm describes the ellipsoid as an origin, c, radii,
r and a rotation R. The soft-iron matrix, S, and hard-iron vector, h, can be
obtained using Equations 6.37 and 6.38. r is the intensity of the Earth’s field in
Gauss within the calibration dataset, known to be approximately 0.5 Gauss in the
UK [30].
r
rx 0 0
T
S = R
0
r
rx
0
R (6.37)
r
0 0 rx
h = Sc (6.38)
Figure 6.19 summarises the heading measurement error provided by the nine
magnetometers for the ‘typical’ magnetic distortion scenario after ellipsoid
135
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
Z
Z
Y X Y X
calibration. The error for rotations around each axis is approximately ±5◦ ,
matching the performance of the worst performing commercial IMU, the UM6.
A B C D E F G H I
10 10 10
8 8 8
θx error (degrees)
θy error (degrees)
θz error (degrees)
6 6 6
4 4 4
2 2 2
0 0 0
−2 −2 −2
−4 −4 −4
−6 −6 −6
−8 −8 −8
−10 −10 −10
0 90 180 270 360 0 90 180 270 360 0 90 180 270 360
Heading (degrees) Heading (degrees) Heading (degrees)
10 10 10
8 8 8
θx error (degrees)
θx error (degrees)
θx error (degrees)
6 6 6
4 4 4
2 2 2
0 0 0
−2 −2 −2
−4 −4 −4
−6 −6 −6
−8 −8 −8
−10 −10 −10
A B C D E F G H I A B C D E F G H I A B C D E F G H I
Magnetometer Magnetometer Magnetometer
Figure 6.19: Heading measurement error provided by the nine magnetometers for
the ‘typical’ magnetic distortion scenario after ellipsoid calibration.
The heading measurement errors for the ‘worst case’ magnetic distortion scenario
are summarised in Figure 6.20. The misalignment around the y axis observed in
Figure 6.18 indicates that the variance of errors in θx and θz will be greater than
136
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
A B C D E F G H I
40 40 40
30 30 30
θx error (degrees)
θy error (degrees)
θz error (degrees)
20 20 20
10 10 10
0 0 0
−10 −10 −10
−20 −20 −20
−30 −30 −30
−40 −40 −40
0 90 180 270 360 0 90 180 270 360 0 90 180 270 360
Heading (degrees) Heading (degrees) Heading (degrees)
40 40 40
30 30 30
θx error (degrees)
θx error (degrees)
θx error (degrees)
20 20 20
10 10 10
0 0 0
−10 −10 −10
−20 −20 −20
−30 −30 −30
−40 −40 −40
A B C D E F G H I A B C D E F G H I A B C D E F G H I
Magnetometer Magnetometer Magnetometer
Figure 6.20: Heading measurement error provided by the nine magnetometers for
the ‘worst case’ magnetic distortion scenario after ellipsoid calibration.
The proposed alignment calibration method requires that the magnetometer has
already been calibrated using the ellipsoid method to yield an initial soft-iron
matrix, Ŝ, and hard-iron vector, ĥ. The calibrated datasets describing rotations
137
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
around the calibration cube x, y and z axis form the red, green and blue circles seen
Figure 6.18 are used to calibrated for misalignment. The data points of each circle
lie on a plane orthogonal to the axis of rotation. A rotation matrix describing the
calibration cube relative to the ellipsoid-calibrated magnetometer can therefore be
found as the solutions to these planes.
Equations 6.39 and 6.40 describe a plane where α, β, γ and δ are fixed
coefficients and mx , bmy and cmz represents an ellipsoid-calibrated magnetometer
measurement that sits on the plane.
α
β
mx my mz 1 =0 (6.40)
γ
δ
mx1 my1 mz1
α 1
−δ
mx2 my2 mz2 1
β
. .. .. −δ = .. (6.41)
.. . .
.
γ
−δ
mxn myn mzn 1
T
Mv = 1 1 . . . 1 (6.42)
The axis orthogonal to the plane, v, can be found as the solution to Equation 6.43
for a dataset of n measurements where n > 3 where M + is the pseudoinverse of
138
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
M.
T
+
v=M 1 1 ... 1
T (6.43)
T
−1 T
= M M M 1 1 ... 1
The vectors orthogonal to each planes are found as vx , vy and vz . These vectors
can be normalised and combined to approximate a rotation matrix as shown by
Equation 6.44. A true rotation matrix, R, can be obtained from R̃ by computing
the best-fit rotation matrix as described in chapter 2.
T
vy
R̃ = vx
kvx k kvy k
vz
kvz k
(6.44)
The aligned soft-iron matrix, S, and hard-iron vector, h, can be obtained from
the initial ellipsoid-calibrated parameters, Ŝ and ĥ by correcting each using the
calculated misalignment, R, as described by Equations 6.45 and 6.46 .
S = RT Ŝ (6.45)
h = RT ĥ (6.46)
139
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
Z
Z
Y X Y X
Figure 6.22 summarises the heading measurement error provided by the nine
magnetometers for the ‘typical’ magnetic distortion scenario after alignment
calibration. The error remains appropriately equal for each axis of rotation but is
the peak-to-peak error through a 360◦ is now half that of the ellipsoid-calibrated
performance. The performance is now approximates that of the 2◦ accuracy
achieved by the VN-100 commercial IMU. Geographic location is known to affect
the heading accuracy achieved by a magnetometer; for example, a ±0.5◦ at the
equator would result in a ±1.5◦ heading error in the UK (see Appendix C). This
suggests that the approximate ±2◦ heading accuracy demonstrated here may be
closer to that of the high-end MTw and 3DM-GX3-25 if the test were to be repeated
at the equator. However, it was not possible to verify this empirically, nor was
it possible to confirm the test conditions of the heading accuracy specified by the
commercial IMU.
The heading measurement errors for the ‘worst case’ magnetic distortion scenario
are summarised in Figure 6.23. The error for rotations around each axis are now
appropriately equal though the errors remain greater than those for the ‘typical’
magnetic distortions. Although some devices achieve the 2◦ accuracy achieved by
140
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
A B C D E F G H I
5 5 5
θx error (degrees) 4 4 4
θy error (degrees)
θz error (degrees)
3 3 3
2 2 2
1 1 1
0 0 0
−1 −1 −1
−2 −2 −2
−3 −3 −3
−4 −4 −4
−5 −5 −5
0 90 180 270 360 0 90 180 270 360 0 90 180 270 360
Heading (degrees) Heading (degrees) Heading (degrees)
5 5 5
4 4 4
θx error (degrees)
θx error (degrees)
θx error (degrees)
3 3 3
2 2 2
1 1 1
0 0 0
−1 −1 −1
−2 −2 −2
−3 −3 −3
−4 −4 −4
−5 −5 −5
A B C D E F G H I A B C D E F G H I A B C D E F G H I
Magnetometer Magnetometer Magnetometer
Figure 6.22: Heading measurement error provided by the nine magnetometers for
the ‘typical’ magnetic distortion scenario after alignment calibration.
the VN-100, the performance is generally closer to that of the UM6 which may
be sufficient for many applications. The contrast between the ‘typical’ and ‘worst
case’ performance demonstrates the need to minimise magnetic distortions within
the physical design of a device.
A B C D E F G H I
10 10 10
8 8 8
θx error (degrees)
θy error (degrees)
θz error (degrees)
6 6 6
4 4 4
2 2 2
0 0 0
−2 −2 −2
−4 −4 −4
−6 −6 −6
−8 −8 −8
−10 −10 −10
0 90 180 270 360 0 90 180 270 360 0 90 180 270 360
Heading (degrees) Heading (degrees) Heading (degrees)
10 10 10
8 8 8
θx error (degrees)
θx error (degrees)
θx error (degrees)
6 6 6
4 4 4
2 2 2
0 0 0
−2 −2 −2
−4 −4 −4
−6 −6 −6
−8 −8 −8
−10 −10 −10
A B C D E F G H I A B C D E F G H I A B C D E F G H I
Magnetometer Magnetometer Magnetometer
Figure 6.23: Heading measurement error provided by the nine magnetometers for
the ‘worst case’ magnetic distortion scenario after alignment calibration.
141
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
Heading error
Calibration solution
‘Typical’ ‘Worst case’
Hard-iron ±5.5◦ ±55.5◦
Ellipsoid ±3.3◦ ±21.1◦
Ellipsoid + alignment ±1.9◦ ±3.9◦
Table 6.2: Heading error for each magnetometer calibration solution achieved for
‘typical’ and ‘worst case’ magnetic distortion scenarios
Where as hard-iron and ellipsoid calibration require only the arbitrary rotation
dataset, alignment calibration requires the additional three single-axis rotation
datasets. This increases the complexity the calibration process but enables
significantly improved performance. In the case of ‘typical’ distortions, the ±1.9◦
accuracy achieved matches that of the VN-100. The ±3.9◦ accuracy achieved for
‘worst case’ distortions is less accurate but still sufficient for many applications.
None of the calibration solutions were able to match the ±0.5◦ heading accuracy
of the 3DM-GX3-25. Although it may be assumed that the greatest accuracy
142
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
achieved, ±1.9◦ , is limited by the geographical location of the tests, there was not
an opportunity to demonstrate this empirically or to clarify the test conditions of
the 3DM-GX3-25.
In this simplified calibration model, the common calibrated frame becomes that
of the physical accelerometer and the magnetometer calibration parameters are
adjusted accordingly. The gyroscope is assumed to be physically aligned with the
accelerometer so that the calibration models for the three sensors may be redefined
as Equations 6.47 to 6.49.
m = S̃um − h̃ (6.49)
143
Chapter 6. Calibration of low-cost MEMS sensors for AHRS applications
The adjusted soft-iron and hard-iron parameters S̃ and h̃ in Equation 6.49 are
calculated as Equations 6.50 and 6.51 where R is the rotation matrix describing
the calibrated frame relative to the physical accelerometer, previously determined
in Section 6.5.3.
T
S̃ = C
AR S (6.50)
T
h̃ = C
AR h (6.51)
6.8 Conclusion
This chapter has presented a complete calibration solution for IMU incorporating a
gyroscope, accelerometer and magnetometer. Practical performance was evaluated
for the low-cost (<10 USD) MPU-6050 and HM-5883L sensors and benchmarked
against four commercial IMUs price between 199 USD and 3000 USD. The results
indicate a static pitch and roll accuracy matching that to the highest performing
commercial IMU and a static heading accuracy matching that of a mid-range
commercial IMU.
144
Chapter 7
7.1 Introduction
145
Chapter 7. Revised AHRS algorithm
Although the calculation of ‘zero-g’ and global accelerations are not required for
the computation of attitude or heading, these accelerations can only be obtained
in conjunction with an AHRS algorithm and are often the ultimate quantities
required by an application.
146
Chapter 7. Revised AHRS algorithm
Z
I I
Eq = E q̇.dt (7.1)
I 1I T
E q̇ = E q̂ ⊗ 0 (ω − Ke) (7.2)
2
2qx qz − 2qw qy
a
ea = â × 2q
y z q + 2q q
w x , where â = (7.3)
kak
2 2
2qw − 1 + 2qz
147
Chapter 7. Revised AHRS algorithm
matrix obtained from IE q̂. The cross product of â and m̂ yields a vector that is
orthogonal to both gravity as assumed by the accelerometer, and the measured
direction of the magnetic field thus providing a vector measurement on a plane
parallel with the Earth’s surface that is insensitive to the orientation of the IMU
and the inclination of the magnetic field, as in [137]. This is a feature of Digital
Magnetic Compasses (DMCs) sometimes refereed to as tilt compensation.
2 2
−2qw + 1 − 2qx
m
em = â × m̂ ×
−2q q
x y + 2q q
w z , where m̂ = (7.4)
kmk
−2qx qz − 2qw qy
The initial value of IE q̂, assumed to be an identity, will likely be incorrect upon
initialisation of the algorithm. The low value of K and corresponding slow
convergence of IE q̂ desired for normal operation will cause this error to be sustained
for a period of several seconds. This initialisation time can be reduced by ramping
148
Chapter 7. Revised AHRS algorithm
down K from a high value to the value intended during normal operation. This
process is described by Equation 7.6 where tinit is the initialisation period in
seconds, Kinit is the initial large value of K and Knormal is the value intended
during normal operation where Knormal > Kinit .
Knormal + tinit −t
tinit
(Kinit − Knormal ) if t < tinit
K= (7.6)
K else
normal
The sensor conditioning processes described in this section are performed prior to
the main complementary filter described in Section 7.2.
ω 0 = ω − ωbias (7.7)
149
Chapter 7. Revised AHRS algorithm
Z
ωbias = 2πfc pω.dt (7.8)
p is determined as Equation 7.9 where the function fb (ω, ωmin ) that computes the
time that the magnitude of each element of ω has been below ωmin . ωmin is the
minimum angular velocity below which the bias estimator can be enabled. tb is
the minimum stationary period after which the bias estimator will be enabled.
1 if fb (ω, ωmin ) > tb
p= (7.9)
0 else
Magnetic measurements that deviate from the expected intensity of the Earth’s
magnetic field are rejected to prevent magnetic distortions from corrupting the
algorithm output. This is achieved by substituting m0 in place of m in Equations
7.4 and 7.5. m0 is defined by Equation 7.10 where mmin and mmax are the specified
minimum and maximum magnitudes of the Earth’s magnetic field in units of
the magnetometer. The selection of m0 as zero works in conjunction with the
conditional selection of e in Equation 7.5 to omit the magnetometer from the
algorithm update.
m if mmin < kmk < mmax
0
m = T (7.10)
0 0 0 else
150
Chapter 7. Revised AHRS algorithm
period of time. The process is described by 7.11 where the function f (a, gmin , gmax )
calculates the period of time that the condition of −gd < (kak − 1) < gd has been
false. T
0 0 0 if f (a, gd ) > ta
a0 = (7.11)
a else
2qx qz − 2qw qy
azero =a−
2qy qz + 2qw qx (7.12)
2qw2 − 1 + 2qz2
The global acceleration, aglobal , is the acceleration of the IMU in the Earth
coordinate system as defined by the ENU convention. It is calculated as Equation
7.13.
0 aTglobal =
I
E q̂ ⊗ 0 aTzero ⊗ IE q̂ ∗ (7.13)
Figure 7.1 shows a block diagram of the complete AHRS algorithm and indicates
the connection between the separate processes described in Sections 7.2 to 7.4.
151
Chapter 7. Revised AHRS algorithm
Gyroscope bias Sensor conditionaing Complementary filter estimator ‘Zero-g’/global acceleration calculation
ωbias
Gyroscope
Gyroscope bias compensation
ω0 + ω 0 − Ke Calculate IE q̇
I q̇
E
R I q
E I q
Orientation
ω E I q̂
(Equation 7.7) (Equation 7.2) kI q
E k E
−
Ke
Initialisation behaviour
K (Equation 7.6)
Magnetometer m0 e
152
Error calculation
(Equations 7.3 to 7.5)
Accelerometer
Linear acceleration rejection
a (Equation 7.11) a0
Global acceleration
Zero-g acceleration calculation Global acceleration calculation
aglobal
(Equation 7.12) (Equation 7.13)
Zero-g acceleration
azero
7.6 Discussion
The algorithm computes the error in ω 0 from the accelerometer and magnetometer
through a cross product operation, the reason for this may not be immediately
obvious. In the case of the accelerometer, the error is determined as the difference
between the direction of gravity measured by the accelerometer, â, and the
direction of gravity assumed by IE q̂. The cross product yields ea as a vector with a
magnitude equal to the sine of the angular error in IE q̂ and a direction orthogonal
to both â and the direction of gravity assumed by IE q̂. ea therefore describes the
necessary ‘adjustment’ to IE q̂ applicable to the dimensions of ω 0 . In the case of
the magnetometer, em represents an equivalent adjustment though its direction
is always aligned with gravity to prevent the magnetometer from influencing the
pitch and roll components of orientation.
The sine relationship defining kek suggests a potential flaw in the algorithm.
Angular errors exceeding 90◦ would yield an increasingly small feedback term. An
error of 180◦ would result in no feedback and the algorithm would never converge.
In practice, it is may assumed that angular errors in IE q̂ will always be <90◦ after
initialisation and an appropriate selection of Kinit will avoid the hazard associated
with an initial error of 180◦ .
153
Chapter 7. Revised AHRS algorithm
The gyroscope bias is estimated each time the IMU is detected as being stationary.
Some manufactures suggest similar mechanisms where measurements below a
threshold are zeroed [139] or steady-state signals are sampled as the bias [140].
Such solutions would risk corrupting small amplitude signals and result in
undesirable step changes of the gyroscope bias. The bias compensation proposed in
this chapter is not subject to these shortcomings. Some applications may involve no
154
Chapter 7. Revised AHRS algorithm
natural stationary periods. In such scenarios, bias compensation may still effective
if an initial stationary period can be imposed to allow the system to ‘warm up’
prior to use.
The gyroscope bias compensation parameters include the filter corner frequency,
fc , the minimum angular velocity, ωmin , below which the bias estimator can
enabled and the minimum stationary period, tb , after which the bias estimator
will be enabled. A corner frequency of fc = 0.05 Hz achieves convergence (within
2% of final value) within 13 seconds yet is relatively low to reduce the risk of
attenuating low amplitude motion that may falls below ωmin . A value of ωmin may
be determined as the expected variation in the gyroscope bias during operation.
For example, the thermal response and random walk investigations presented in
Chapter 5 indicate that a threshold of ωmin = 4◦/s would be sufficient. The value
should include a margin to account unexpected noise within the application; for
example, mechanical vibrations. A period threshold of tb = 2 seconds would be
sufficient to prevent the bias compensation of activating during zero-crossing in
most applications.
The magnetic distortion rejection parameters, mmin and mmax define the minimum
and maximum valid intensity for the ambient magnetic field. Default values of
mmin = 0.22 Gauss and mmax = 0.67 Gauss may be assumed as these represent
the natural limits of the Earth’s magnetic field [30]. More robust performance
can be achieved if a reduced range can be specified representative of a specific
operating environment.
155
Chapter 7. Revised AHRS algorithm
The algorithm provides IE q̂ describing the orientation of the Earth relative to the
IMU. Many applications will require a measurement of the IMU relative to the
I ∗
Earth, E
I q̂, which is obtained as the conjugate E q̂ . An equivalent Euler angle or
The algorithm also provides the estimated gyroscope bias, ωbias . This
may be useful to processes external to the AHRS algorithm that require a
bias-compensated gyroscope measurement. The provision of ωbias also enables
the gyroscope bias to be saved to non-volatile memory so that the ‘tracked’ value
can be maintained even when the system is shut down.
156
Chapter 7. Revised AHRS algorithm
7.7 Conclusion
The revised AHRS algorithm presented in this chapter builds on the findings of
the sensor characterisation and calibration work presented in Chapters 5 and 6.
It was developed simultaneously with the IMU platforms and applications and
applications in the next chapter and incorporates features that address practical
requirements, verified through these empirical works.
157
Chapter 7. Revised AHRS algorithm
158
Chapter 8
8.1 Introduction
This chapter describes the development of IMU platforms and the broad
applications they have facilitated. These platforms take advantage of modern
MEMS devices, the calibration solutions presented in Chapter 6, and the AHRS
algorithm presented in chapter 7. The following three sections (8.2 to 8.3)
summarise the development of three IMU platforms developed through this
research project; the: x-IMU, x-BIMU and x-OSC. Each platform fulfils a different
design specification and together have facilitated a wide range of applications and
research projects. Section 8.5 summarises projects that have utilised the platforms
for specific applications in collaboration with other academic researchers. Section
8.6 section presents a selection of user applications that have utilised the IMU
platforms within commercial or academic research projects. The chapter concludes
with a description of the future developments scheduled for for the IMU platforms
and applications.
159
Chapter 8. IMU platforms and applications
The x-IMU was developed in 2010 to provide a versatile IMU and data acquisition
platform for use in a wide range of potential applications. In the four years since,
the hardware and firmware has continued to evolve to take advantage of modern
MEMS sensors and introduce new functionality based on user feedback. The
original design was intended to meet the following specification.
160
Chapter 8. IMU platforms and applications
Figure 8.1: x-IMU board alone (left) and enclosed in its plastic housing with 1000
mAh battery (right)
parameters saved to on-board memory. During normal operation, the sensors are
continuously sampled and the on-board AHRS algorithm updated at 512 Hz. An
extensive set of internal settings can be managed using the x-IMU software to
configure the device behaviour and the send rate of individual packet types up to
512 Hz.
The x-IMU board dimensions are 33×42 mm and it weights 12 g. Figure 8.2 shows
an annotated top and bottom view with key hardware components are described
below.
(b) Crystal oscillator: 32.768 kHz crystal provides a time base for the RTC
sample timer.
(c) USB bridge: Future Technology Devices International (FTDI) FT232R USB
bridge. See Section 8.2.2 for more information.
161
Chapter 8. IMU platforms and applications
(o) (e)
(n)
(f)
(m)
(g)
(l)
(h)
(k)
(j)
(i)
Figure 8.2: x-IMU top and bottom view with labelled key hardware components
(f) Bluetooth LED: Blue Light-Emitting Diode (LED) that indicates the
connection status of Bluetooth.
(g) Battery connector: Socket for 3.7 V Lithium-ion cell. Battery voltage is
measured on-board.
(h) Power switch: Slide switch for switching the device on and off or selecting
162
Chapter 8. IMU platforms and applications
(i) SD card socket: Micro SD card socket supports FAT16 and FAT32 formats
for capacities up to 32 GB.
(j) Bluetooth: Microchip RN-41 Bluetooth module. See Section 8.2.2 for more
information.
(k) Auxiliary port: Eight channel configurable auxiliary port also includes power
I/O pins. See Section 8.2.4 for more information.
(l) Status LED: Green LED indicates the status of the x-IMU.
(n) SD card LED: Amber LED that indicates the connection status of SD card.
The x-IMU Graphical User Interface (GUI) is a Windows application that provides
an interface for configuring internal settings and displaying real-time data. Figure
8.3 shows a screen shot of the GUI displaying internal sensor measurements in 2D
plots and the AHRS data as a 3D representation. The software also provides tools
for synchronisation of real-time clock, firmware updates, inertial and magnetic
calibration tools and data logging to Comma-Separated Values (CSV) files for
exporting data to software such as MATLAB and Microsoft Excel. Figure 8.3
shows the x-IMU GUI with gyroscope and accelerometer measurements being
displayed as 2D plots, and AHRS data being displayed as a 3D representation.
163
Chapter 8. IMU platforms and applications
Figure 8.3: x-IMU GUI displaying real-time data with 2D and 3D graphics
The x-IMU enables wireless communication via USB, Bluetooth or serial through
the the axillary port. USB communication is enabled through a FTDI FT232R
chip to maximise compatibility with FTDI drivers are available for Windows, OS
X, Linux [143] and Android [144]. The RN-41 Bluetooth module is a class 1 device
with a range of 100 m [145]. It uses the Serial Port Profile (SPP) Bluetooth profile
eliminates the need for host drivers for compatibly with a wide range of platforms
1
https://github.com/xioTechnologies/x-IMU-GUI
2
https://github.com/xioTechnologies/x-IMU-Arduino-Example
3
https://github.com/xioTechnologies/x-IMU-Android-Example
164
Chapter 8. IMU platforms and applications
such laptops and smart phones. The auxiliary port can be configured as a serial
interface for communication with embedded systems. Settings support hardware
flow control and standard baud rates.
The auxiliary port is a 12 pin socket that provides access to power and eight
general purpose analogue/digital I/O pins. The auxiliary port can be configured
in several different modes including: analogue inputs, digital I/O, Pulse-Width
Modulation (PWM), serial and external sleep/wake mode. Analogue and digital
input modes allows data from external sensors to be sent in real-time or logged
to the SD card. PWM and digital output modes allow real-time control of wide
range of device, for example, LEDs or vibration motors for haptic feedback. The
165
Chapter 8. IMU platforms and applications
real-time applications
The x-BIMU was developed in 2012 to provide a more versatile wireless alternative
to the x-IMU. The x-IMU uses Bluetooth which imposes several limitations. The
Bluetooth standard limits the maximum number of active devices to seven [147]
which may be insufficient for many applications. Throughput and latency become
unreliable when using multiple devices with a single Bluetooth receiver which may
comprise real-time applications. Bluetooth modules also consume more power that
some other wireless technologies and so reduce the potential battery life.
• Low-power operation
166
Chapter 8. IMU platforms and applications
Size and weight are often important for wireless applications and so the x-BIMU
was designed for a compact form factor. When enclosed within its plastic housing
with a 320 mA battery and wireless module, the x-BIMU has a total size of 30 ×
38 × 22 mm and weight of 22 g. Figure 8.5 shows the x-BIMU in plastic housing
alone (left) and combined with a Velcro body strap (right) designed to secure the
plastic housing for human motion applications.
167
Chapter 8. IMU platforms and applications
Figure 8.5: x-BIMU in plastic housing with XBee 802.15.4 module and 320 mAh
battery (left) and Velcro strap designed for human motion applications (right)
168
Chapter 8. IMU platforms and applications
(b) (c)
(a) (d)
(e)
(j)
(f)
(i)
(g)
(h)
(a) RGB status LED: The colour indicates the selected wireless channel of
the XBee 802.15.4 module of the battery charging status if the x-BIMU
is connected to a charger.
(b) Battery connector: Socket for 3.7 V Lithium-ion cell. Battery voltage is
measured on-board and can be transmitted via the connected communication
module to enable to user to monitor the battery level.
(c) USB connector: A socket charging the connected battery, it does not provide
USB communication.
(d) Power button: Used to turn the device on and off. The button also provides
a fail-safe mechanism of recovery if the device is configured with incorrect
169
Chapter 8. IMU platforms and applications
settings.
(j) XBee socket: Connection for any XBee-style module. Provides interface to
serial communication as well as control signals such as sleep and reset.
4
https://github.com/xioTechnologies/x-BIMU-Terminal
5
https://github.com/xioTechnologies/x-BIMU-Arduino-Example
6
https://github.com/xioTechnologies/x-BIMU-Android-Example
170
Chapter 8. IMU platforms and applications
Figure 8.7: x-BIMU Terminal displaying real-time data with 2D and 3D graphics
The x-BIMU is optimised for XBee 802.15.4 modules but compatible with any
XBee-style communication module. Upon start up, the x-BIMU will automatically
detect if an XBee 802.15.4 module is present and configure it to operate in API
mode. This gives the x-BIMU explicit control over packetisation at the 802.15.4
Medium Access Control (MAC) layer [151] to enable low-latency and low-power
performance. If an XBee 802.15.4 module is not detected, communication reverts
to a generic 115200 baud serial interface compatible with any XBee-style module
with a matching configuration.
(a) Digi XBee 802.15.4 (1 mW): Standard 802.15.4 module with a maximum
transmit power for 1 mW for a range of 100 m [151]. The module shown
incorporates a low-profile PCB antenna.
171
Chapter 8. IMU platforms and applications
(b) Digi XBee 802.15.4 (60 mW): High-power variant of the 802.15.4 module;
operates identically but has a maximum transmit power for 60 mW for a
range of 1 km [151]. The module shown has a wire antennae to provide
better omnidirectional performance.
(c) Microchip RN41XV Bluetooth: Class 1 Bluetooth 2.1 module for 100 m
range [152] Supports SPP to avoid the need for host drivers.
(d) Digi XBee ZB: Supports licensed ZigBee for complex network configurations
such as mesh networks, and compatibility with third-party ZigBee systems.
The module shown is a 2 mW variant for 120 m range, 63 mW variants are
also available for 3.2 km range [153].
(e) Digi XBee-PRO 900: Uses proprietary 900 MHz communication protocol
with range of 15.5 km [154]. The module requires an external antennae.
(f) Digi XBee Wi-Fi: 802.11n for high-throughput User Datagram Protocol
(UDP) or Transmission Control Protocol (TCP) transmission over ad hoc
or existing infrastructure Wi-Fi networks [155].
(g) Seeed Studio BLEbee: Bluetooth 4.0 Low Energy module for direct
172
Chapter 8. IMU platforms and applications
(h) x-BIMU Serial Breakout: Open-source8 serial adapter for direct wired
interface to microelectronics or USB and RS-232 platforms when combined
with a corresponding adapter.
The 802.15.4 standard is specially designed for low-power devices to spend large
amount of time asleep [157]. The XBee 802.15.4 idle current is 50 mA but also
includes a <10 μA sleep mode with a typical wake up time of 10 ms [151]. The
x-BIMU has explicit control of packetisation of the 802.15.4 MAC layer and so is
able to cycle this sleep mode between packet transmission to achieve significant
power savings for low data-rate applications.
7
https://github.com/michaelkroll/BLEbee
8
https://github.com/xioTechnologies/x-BIMU-Serial-Breakout
173
Chapter 8. IMU platforms and applications
In some applications it may be desirable for the wireless IMU to power down
during periods of inactivity to extend the battery life. The x-BIMU supports this
through a sleep timer and motion trigger wake up. The sleep timer will trigger
the x-BIMU to enter sleep mode if no motion is detected for a specified period
of time. The motion trigger wake up will cause the x-IMU to immediately wake
up again once motion is detected. An important feature of this mechanism is the
low-latency wake up as any missed motion data could compromise measurements.
To minimise power consumption, motion is detected using the MPU-6050 built-in
Digital Motion Processor (DMP) [121]. This allows the all other electronics
to completely power down until a hardware interrupt pulse is generated by the
MPU-6050.
Table 8.1: Measured current consumption of x-BIMU and XBee 802.15.4 module
174
Chapter 8. IMU platforms and applications
Figure 8.9: Four 802.15.4 receivers (Digi XSticks) in a USB hub provide
synchronised measurements from x-BIMUs
9
https://github.com/xioTechnologies/x-BIMU-Logger
175
Chapter 8. IMU platforms and applications
I/O interface
The x-OSC was developed in 2012 and intended to address the need for a combined
high-performance wireless I/O interface and IMU. The x-OSC was introduced
in a 2013 paper [9] which provides a survey of existing wireless I/O devices and
discussion on their shortcomings. Recent years have seen the appearance of several
commercial wireless IMU systems, for example the Xsens MTw [128] and APDM
opal [158]. However, these systems utilise low-bandwidth wireless technologies
and so are typically limited to an update rate of ∼120 Hz for multiple IMUs.
Furthermore, none of these systems offer an I/O interface for external sensors and
output devices.
176
Chapter 8. IMU platforms and applications
The x-OSC, shown in Figure 8.10, is a wireless I/O board that provides just
about any software with access to 32 high-performance analogue/digital channels
and on-board IMU measurements via Open Sound Control (OSC) messages over
Wi-Fi. There is no user programmable firmware and no software or drivers to
install making x-OSC immediately compatible with any Wi-Fi-enabled platform.
177
Chapter 8. IMU platforms and applications
(b)
(a)
(l)
(k) (j) (i)
Figure 8.11: x-OSC top and bottom view with labelled key hardware components
(b) Ping button: Used to trigger network ‘pings’ to facilitate device discover.
Also provides a fail-safe mechanism of recovery if the device is configured
with incorrect settings.
(c) 16 Input channels: See Section 8.4.2 for more information. Includes
regulated 3.3 V output to power user electronics.
178
Chapter 8. IMU platforms and applications
(e) RGB LED: Indicates the networking mode and connection status. Can also
be controlled by user to provide visual feedback within an application.
(g) 16 Input channels: See Section 8.4.2 for more information. Includes
unregulated battery voltage output to power user electronics.
(j) 74LVC125A Quad buffers : Output buffers enable each of the 16 outputs to
sink/source up to 50 mA [166].
(k) Battery connector: Socket for 3.7 V Lithium-ion cell. Battery voltage is
measured on-board and available to the user as an OSC message.
(l) Crystal oscillator: 16 MHz, ±20 ppm crystal drives system clock [167].
An embedded web server enables all internal settings to be configured using a web
browser, see Figure 8.12. Settings may be viewed and modified during run-time
without interrupting the acOSC messages. This eliminates the need for any
device-specific software and allows the x-OSC to be configured using any laptop,
smart phone or tablet.
179
Chapter 8. IMU platforms and applications
180
Chapter 8. IMU platforms and applications
In addition to modes described above, the first four inputs and outputs can be
configured to serial mode with each transmit and receive pair utilising a dedicated
hardware Universal Asynchronous Receiver/Transmitter (UART) module. Each
serial channel supports baud rates in the range 9600 to 1 M baud and incorporates
a 2 kB buffer to ensure high throughput without loss of data. Received serial
data is framed before being sent as OSC-blob messages. Framing boundaries are
determined by a user defined buffer size, timeout and optional framing character.
181
Chapter 8. IMU platforms and applications
The round-trip latency was evaluated by wiring an x-OSC digital input and output
together and then measuring the time between an output toggle message being
sent by the host computer and the input change message being received from
x-OSC. The digital input and digital output OSC messages are 100 and 32 bytes
long respectively, each incorporating an address pattern and either 16 or 1 int32
argument/s. A software application was written to send a message to toggle the
digital output every 50 ms and WireShark was used to log the time of each
packet sent and received. This method of evaluating latency is different from
that previously proposed [9] and has the advantage of eliminating the software
application from the loop, which may impose an unknown latency specific to
182
Chapter 8. IMU platforms and applications
the software, Operating System (OS) or processing load of the computer. The
infrastructure network was created by connecting the computer and x-OSC to a
single router/AP.
Ideal (latency test-packets only) Under load (with additional 200 packets per second)
0.6
Normalised distribution
0.5
0.4
0.3
0.2
0.1
0
0 2 4 6 8 10 12 14 16 18
Round-trip latency (ms)
Figure 8.13 shows the round-trip latency distribution of over 13000 samples
achieved for either ‘ideal’ conditions, where communication is limited to only
the digital input/output messages; and ‘under load’ conditions, where x-OSC
was configured to simultaneously send analogue input data at 200 messages per
second. Table 8.2 summarises the results. It may be assumed that the latency
for communication in either direction is approximately half that of the observed
round-trip latency, for example, <3 ms for ‘ideal’ conditions.
8.4.4.2 Throughput
183
Chapter 8. IMU platforms and applications
OSC message including 16 float32 arguments. WireShark was used to log the
time of arrival of each packet and packets per second was calculated as the
number of packets arriving from each x-OSC within each one second window.
Each experiment starts with a single running x-OSC. At one minute intervals, an
additional x-OSC is activated for a period of 15 minutes, to yield a recording of
throughput for 1 to 15 x-OSCs. Tests were conducted with the 15 x-OSCs sharing
a single channel and evenly distributed between three non-overlapping channels to
investigate the benefit.
x−OSC 1 x−OSC 2 x−OSC 3 x−OSC 4... Sum
4000
3600
3200
Packets per second
2800
2400
2000
1600
1200
800
400
0
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
Minutes
Figure 8.15 shows the throughput for 1 to 15 x-OSCs connected to three APs, each
operating on a separate non-overlapping channel. Distribution between multiple
channels can be seen to produce an increased net throughput of ∼4800 packets per
second. The first group of five x-OSCs were configured on channel 1, the next on
184
Chapter 8. IMU platforms and applications
4800
4000
2400
1600
800
0
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
Minutes
channel 11 and the final five on channel 6. This specific order demonstrates that the
channel 1 and 11 groups are able to operate simultaneously without interference.
After 10 minutes, the inclusion of the the final group (channel 6) increases the net
throughput proportionally but with significantly increased variance. Crucially, it
can be seen that the the throughput of an x-OSC in the channel 11 group falls
below that of the rest once the channel 6 group appears. This demonstrates the
potential for channel 6 to interfere with channels 1 or 11 [168].
Although the above results show that the x-OSC is able to achieve both
high-throughput and low latency, the large channel access overhead involved in
the transmission of 802.11 packets does produce a trade-off between throughput
and latency [169]. For example, the x-OSC firmware was modified to demonstrate
a maximum throughput of ∼400 packets per second with a 1 byte payload (0.4
kB/s), and ∼295 packets with a 1472 byte payload (∼434 kB/s). This corresponds
to a 1000 times increase in throughput at the cost of a ∼25% reduction in send rate.
This indicates the potential for the x-OSC to achieve significantly higher sample
rates by allowing the user to sacrifice latency in favour of bundling multiple OSC
185
Chapter 8. IMU platforms and applications
The development of the IMU platforms has coincided with exploration of numerous
applications over the past four years. However, the largest projects were achieved
through collaboration with other people. This section present two selected
projects.
“The gloves” are a pair of data gloves that allow the wearer interact with software
though the natural motion and gesticulation of their hands. Although the gloves
have potential for a wide variety of applications, development has focussed on
live music performance where the gloves provide the wearer with an unlimited
number of virtual musical instruments and access to all the tools of a professional
production desk; all played and controlled live, through precise motion and
gestures.
The project started in 2010 as a partnership between a musician Imogen Heap, and
a computer scientist Dr Tom Mitchell of the University of the West of England.
Together they created ‘SoundGrasp’ [170]; a single glove with a wrist-mounted
mic that would allow the wearer to record and play back audio samples. The
glove was a 5DT 14 Ultra glove which incorporated 14 fibre optic bend sensors to
provide measurements finger flexion and abduction via USB. An ANN was used
to identify predefined postures that would either trigger events or switch modes.
The continuous motion of flexion also provided a means to control variable audio
effects including reverb and filtering. Figure 8.16 shows the glove with the postures
and their associated control audio function.
186
Chapter 8. IMU platforms and applications
Figure 8.16: SoundGrasp gloves with wrist-mounted microphone and postures with
associated control function. Image source: [170]
The fingertip-less gloves did not obstruct the wearer from playing most musical
instruments and the wrist-mounted microphone could be used to sample both voice
and acoustic instruments. SoundGrasp was an exciting proof of concept but is was
clear that much more could be achieved with the adoption of IMU technology. The
author of this dissertation joined the partnership to help achieve this.
The prototype gloves built upon the SoundGrasp system by using a pair of gloves
and mounting an x-IMU to each wrist. Each x-IMU used a USB connection along
side that of the 5DT gloves. The x-IMU auxiliary ports were configured in PWM
mode to drive an Red, Green, Blue (RGB) LED on each hand. A dedicated
headset microphone was also added. Figure 8.17 shows the prototype gloves system
diagram (left) and being worn by musician (right).
187
Chapter 8. IMU platforms and applications
Figure 8.17: Prototype gloves system diagram (left) and being worn by musician
(right)
The performance gloves represented a significant advance to the system. The 5DT
gloves were abandoned in order to achieve a wireless system with custom designed
components to meet technical requirements. The complete performance gloves
system is shown in Figure 8.18 with labelled key components.
• Arm bands: Contained an x-IMU with a wired connection to the RGB LED
and vibration motor on the hand.
188
Chapter 8. IMU platforms and applications
Arm bands
Gloves
• Harness and back hub: Mounts wireless radio units for the two
wrist-mounted microphones, headset microphone and audio earpiece. Also
mounts a x-IMU an custom designed hub that routes wired serial
communication between the five x-IMUs and the Bluetooth radio within
the head piece.
• Head piece: Mounts a Bluetooth radio that provides a single wireless link
between the five x-IMUs and the host software.
The performance gloves software utilised the same control mechanisms as outlines
for the prototype gloves but was extended to provide a direct interface to existing
professional music production software. The use of five x-IMUs distributed
between the arms and torsos provided an opportunity to implement more detailed
measurement of the body’s motion. A key use was to use the torso orientation
189
Chapter 8. IMU platforms and applications
measurement as a ’zero reference’ for each glove so that orientation controls would
be relative to the wearer and not the Earth.
The performance gloves have provided a unique tool for the composition and
performance of music. They have also have potential for a use in a range of
applications not yet fully explored. For example, Figure 8.19 shows the gloves
being used to create and manipulate visuals projected on to a screen.
Figure 8.19: Gloves being used to create and manipulate visuals projected on to
a screen
190
Chapter 8. IMU platforms and applications
This project explored the potential application of the x-IMU to human gait studies.
The autonomous data logging and motion triggered sleep/wake functionality make
the x-IMU an idea tool ambulatory monitoring with the potential to address
concerns outlined by Zhoua.
foot y axis
foot z axis
0.2
0
3 2.8 2.6 2.4 2.2 2 1.8 1.6 1.4 1.2 0.2
1 0.8 0
0.6 0.4 0.2
0.2 0
Figure 8.20: Measured foot position during walking (3 steps) obtained using an
x-IMU
191
Chapter 8. IMU platforms and applications
The x-IMU has been on public sale since 2010 and the x-BIMU and x-OSC since
2013. More than 500 of these platforms have been sold to customers world-wide;
including many commercial and academic research institutions. This section
presents 15 selected example user applications that have utilised the platforms
developed. The research projects chosen to appear in this section were selected
to demonstrate the diversity of user applications and impact of the platforms in
peer-reviewed publications.
Multi IMUs mounted on an upper limb can be used to determine the limb joint
kinematics and provide a valuable tool for rehabilitation of stroke patients. Daniel
Galinski and Bruno Dehez used four x-IMUs to explore this in work conducted at
Center for Research in Mechatronics, Université catholique de Louvain, Belgium.
In a 2012 paper [187] they present an evaluation of initialisation procedures for
this application of IMUs. The evaluation was achieved through a comparison of
the measurement provided by four x-IMU fixed to an anthropomorphic upper limb
and its built in absolute encoders. Figure 8.21 shows the corresponding mechanical
model of upper limb, placement of x-IMUs and instrumented anthropomorphic
192
Chapter 8. IMU platforms and applications
upper limb. Their investigation validated the tested initialisation procedures and
suggested a potential advantage in one procedure for physically weak patients.
Figure 8.21: Mechanical model of upper limb (a), placement of x-IMUs (b) and
instrumented anthropomorphic upper limb used for evaluation (c). Image source:
[187]
Digits is a wrist-worn sensor that recovers the full 3D pose of the user’s hand
without requiring any external sensing infrastructure or covering the hand itself.
It is a Microsoft Research project in collaboration with the Culture Lab, Newcastle
University; and FORTH, University of Crete. The hardware, shown in Figure is
comprised on an x-IMU, an IR laser line projector, diffused IR LEDs, an and
IR camera. The x-IMU was connected via USB and used to measure the wrist
orientation while the IR LEDs line projector and LEDs project on to the hand to
be seen by the camera. A wide range of signal processing techniques are used with
a kinematic model of the hand to determine a full 3D posture from the data. The
techniques and an evaluation of performance are discussed in detail in a 2012 paper
[188] published by the team. A range of applications have been demonstrated on
the Microsoft Research website10 with digits providing an eyes-free interaction for
10
http://research.microsoft.com/en-us/projects/digits/
193
Chapter 8. IMU platforms and applications
underwater swarms
194
Chapter 8. IMU platforms and applications
The MONSUN project started in 2009 and has been a source of numerous
publications covering the hardware and software design process [191] including
novel components such as an acoustic modem for underwater communication
[192]; simulation studies exploring practical applications such as environmental
monitoring [190]; and significant towards fault-tolerant and energy-efficient swarms
of underwater robots [193]. The most recent publication presented a preview of
the forthcoming third generation MONSUN in 2013 [194] which may incorporate
the x-BIMU.
195
Chapter 8. IMU platforms and applications
on-board electronics. The x-IMU is used along side the Tritech Micron Scanning
Sonar for localisation and navigation of the craft. SMART-E won the innovation
award at the European SAUC-E competition 2012.
Figure 8.25: SMART-E AUV uses three actuated thrusters for omnidirectional
movement and a band of RGB LEDs to provide visual feedback under water.
Image source: [195]
196
Chapter 8. IMU platforms and applications
The team tested two versions of the TTM, each comprised of a steel float tethered
to an anchor with several instruments mounted along the tether as shown in Figure
8.27. The mounted instruments include: ADVs, current meters, x-IMUs, and an
Acoustic Wave And Current Profiler. The x-IMU logged data to its on-board SD
card to provide orientation data necessary for motion correction. The TTMs were
deployed at different locations through 2012 and 2013. Results showed although
raw data saw sufficient to evaluate turbulence intensities and reasonable turbulence
spectra, Motion correction would more detailed aspects of the turbulence.
training system
197
Chapter 8. IMU platforms and applications
and place them on small cylinders within a virtual environment. Figure 8.29 shows
a screen shot of this software. The authors suggest that future development could
improve the system through the creation of virtual human organs.
198
Chapter 8. IMU platforms and applications
Figure 8.28: Instrumented surgical tool uses x-IMU to provide orientation data.
Image source: [196]
Figure 8.29: Virtual reality game controlled using the instrumented surgical tool.
Image source: [196]
199
Chapter 8. IMU platforms and applications
Figure 8.30: Instrumented equipment comprised an a soccer boot and shin guard;
each with an x-IMU mounted internally. The equipment is shown here attached
the the synchronisation trigger prior to a data acquisition session. Image source:
[198]
‘low-cost’ IMUs
Standalone ‘low-cost’ IMUs, such as the x-IMU with its on-board data logging
to an SD card, could facilitate large-scale studies into the orthopaedic deficits
(lameness) of horses. In a 2013 paper [199] submitted to the journal of Computer
Methods in Biomechanics and Biomedical Engineering , Brighton et al. prevent
their investigation to the accuracy and limits of agreement of the x-IMU compared
with an established IMU-based gait analysis system, the Xsens MTx [48]. The
200
Chapter 8. IMU platforms and applications
study involved six horses, both IMUs were mounted on two anatomical landmarks,
the sacrum and sternum. The vertical velocity and position of each location were
estimated by combining the orientation and accelerometer measurements and the
symmetry of the horse gait evaluated. The x-IMU was found to provide sufficient
precision when assessing symmetry but was not able to provide reliable data for
an alternative analytical method that as tested. The authors suggest that further
investigations into specific calibration and processing algorithms could address
these shortcomings.
robot head
In a recent publication [202], the research team presented a system that allows
a human operator to control Flobi by wearing a combined motion capture and
display helmet. The system allows the the direct and live transfer of human
facial expressions, gaze and head movements while providing the operator with
a real-time display of the scene as perceived by the robot’s vision sensors.
The operators helmet, shown in Figure 8.31, incorporates a video projector,
face camera, x-IMU and earphone. The x-IMU is used to provide real-time
measurement of the head orientation to control Flobi’s neck actuators. The system
201
Chapter 8. IMU platforms and applications
also records human motion datasets for analysis or playback. Other applications
of Flobi have included its role as a fitness instructor [203] and an investigation into
the effects of loneliness on psychological anthropomorphism [204].
Figure 8.31: Combined motion capture and display helmet (left) and robot head
(right). Image source: [202]
202
Chapter 8. IMU platforms and applications
for another; and ‘intimacy’, where physical contact between performers is felt by
others and also represented through visualisations to an audience. The authors
conclude that the paper aims to invite further discussion on if such technology can
enhance the experience for both performers and audience.
Figure 8.32: Latest generation autonomous car shown fully assembled (left) and
with shell removed (right)
The x-BIMU is used with a serial breakout board to connect to the main processor
via serial. The angular rate data provided by the x-BIMU is used along side
images provided by the camera to estimate the location of the vehicle on a map
using a RatSLAM algorithm. Simultaneous Localisation And Mapping (SLAM)
encompasses all the processes required for autonomous navigation of a robot in a
203
Chapter 8. IMU platforms and applications
Figure 8.33: Screenshot showing view from camera mounted on car and the plotted
location of the car as estimated by the RatSLAM algorithm
compete in RoboCup
204
Chapter 8. IMU platforms and applications
the PERSIA platform is described in a paper currently under review for the Elsiver
journal of Mechatronics [206]. The robot, shown in Figure 8.34, is actuated by 20
servo motors and senses the environment though a camera and the x-IMU. The
main processor implements an ANN to maintain balance and gait control. An
x-IMU installed on the hip is used to provide feedback to facilitate the training of
the ANN. Other software mechanisms for path planing and decision making are
also described in the paper. In the recent RobotCup IranOpen 2013 competition,
the team were awarded first place for the ‘adult size’ robot (150 cm) category.
Figure 8.34: ‘Adult size’ robot (150 cm) robot walking towards football
(Photograph courtesy of Hamidreza Kasaei)
205
Chapter 8. IMU platforms and applications
military helicopters
Monitoring the in-service usage of military aircraft and comparing this usage with
design assumptions is essential for assuring structural integrity throughout the life
of a fleet. Traditional approach to these validation programmes can be expensive
and so there is a demand for alternative low-cost solutions. A team at the Defence
Science and Technology Laboratory, a part of the UK Ministry of Defence, has
been evaluating the potential of the x-IMU to provide a cost-effective solution for
monitoring the usage of historic and small-fleet rotary-wing aircraft. The x-IMU
is used as a self-contained data logger to record gyroscope and accelerometer data
to the on-aboard SD card. The sleep timer and motion trigger wake up allow the
x-IMU to remain in situ for extended periods of time for minimal maintenance.
The x-IMU has been used to capture flight data from the Army Historic Aircraft
Flight Scout AH Mk 1 aircraft shown in Figure 8.35. Currently, over 20 flying
hours of data have been captured and is being analysed and compared with the
design usage assumptions for the aircraft.
Figure 8.35: The x-IMU is being used to collect data from the Army Historic
Aircraft Flight Scout AH Mk 1 to monitor structural usage
206
Chapter 8. IMU platforms and applications
207
Chapter 8. IMU platforms and applications
Prototype hardware was developed using x-OSC, which provided the team with
a versatile wireless interface board with the potential to explore a wide range of
208
Chapter 8. IMU platforms and applications
sensor and actuator concepts. The current block incorporate two servo motors
to actuate the one face of the block about two degrees of freedom. Each block
connects to a router to achieve a single Wi-Fi connection to the computer running
the software... x-OSC/Wi-Fi offers the potential for large number of nodes.
The prototype block design, shown in Figure 8.38 incorporates two servo motors
to acute the one face of the block about two degrees of freedom. The paper
demonstrates these six blocks being assembled together to form pulsing heart.
Figure 8.38: An individual block containing x-OSC for Wi-Fi control and 2 servos
with scissor lift mechanism. Image source: [209]
Future work will see the IMU platforms continue to evolve based on user feedback
and to take advantage of increasingly high-performance MEMS technology.
Specific tasks have already be outlined for 2014. The modular communication
of the x-BIMU as proved to be a versatile design paradigm but the device is
limited by the processor; A new hardware revision will utilise the latest 32-bit
processors to achieve equivalent low-power performance with superior same rates
and signal processing capabilities. Planned x-OSC firmware developments will
take full advantage of Wi-Fi bandwidth to enable sample rates in the order kHz
with the potential to facilitate new applications.
The gloves project presented in Section 8.5.1 is an ongoing collaboration and will
see the development of next generation hardware and software in 2014 using the
x-OSC.
209
Chapter 8. IMU platforms and applications
210
Chapter 9
Conclusions
This project set out to demonstrate how new applications can be realised using
modern low-cost MEMS gyroscopes, accelerometers and magnetometers. This was
achieved through the development of publicly available software and hardware
platforms that have facilitated a wide range of commercial and academic research
projects exploring a diverse range of applications. In achieving this, the project
has met the following objectives:
• Design and manufacture of IMU platforms that combine the calibration and
AHRS algorithm solutions
In achieving its objectives, this project has made several contributions. Chapters
3, 4 and 7 present sensor fusion algorithms applicable to conventional IMUs
211
Chapter 9. Conclusions
and non-gyro IMUs. The conventional AHRS algorithms address the specific
requirements of modern low-cost MEMS sensors and have been proven in the IMU
platforms and applications presented in Chapter 8. The impact of these algorithms
is also demonstrated by the numerous citations and code downloads. The non-gyro
IMU was not demonstrated in application but represents an interesting alternative.
The novel work presented in Chapter 4 demonstrated the potential practical
benefits of such kinematically redundant sensor arrays.
The AHRS algorithm and sensor calibration works were brought together in
the development of the Chapter IMU platforms presented in Chapter 8. Each
of these platforms each addressed a specific design specification and together
facilitated a wide range of new applications; demonstrated by the numerous
scientific publications that resulted from collaborative projects and user projects.
This project achieved its original objectives but MEMS technology is continuing
to advance and there are undoubtedly many more new and exciting applications
that can be facilitated by the work conducted within this project. Future work
will see the ongoing development of the IMU platforms to incorporate improved
MEMS sensors and new features to meet the needs of an expanding user base.
Work will also continue on the AHRS algorithm and sensor calibration solutions
to meet the specific needs of improved MEMS sensors and new applications.
212
References
213
References
214
References
215
References
2013.
[34] M. Tanenhaus, D. Carhoun, T. Geis, E. Wan, and A. Holland, “Miniature
imu/ins with optimally fused low drift mems gyro and accelerometers
for applications in gps-denied environments,” in Position Location and
Navigation Symposium (PLANS), 2012 IEEE/ION, pp. 259–264, April 2012.
[35] S. K. Hong, “Fuzzy logic based closed-loop strapdown attitude system for
unmanned aerial vehicle (uav),” Sensors and Actuators A: Physical, vol. 107,
no. 2, pp. 109 – 118, 2003.
[36] B. Barshan and H. Durrant-Whyte, “Inertial navigation systems for
mobile robots,” Robotics and Automation, IEEE Transactions on, vol. 11,
pp. 328–342, Jun 1995.
[37] 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, vol. 1, pp. 317–322, May
11–15, 2002.
[38] D. H. Titterton and J. L. Weston, Strapdown inertial navigation technology.
The Institution of Electrical Engineers, 2004.
[39] S. Beauregard, “Omnidirectional pedestrian navigation for first responders,”
in Proc. 4th Workshop on Positioning, Navigation and Communication
WPNC ’07, pp. 33–36, Mar. 22–22, 2007.
[40] H. J. Luinge and P. H. Veltink, “Inclination measurement of human
movement using a 3-d accelerometer with autocalibration,” vol. 12,
pp. 112–121, Mar. 2004.
[41] H. Zhou and H. Hu, “Human motion tracking for rehabilitation–a survey,”
Biomedical Signal Processing and Control, vol. 3, no. 1, pp. 1 – 18, 2008.
[42] 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, pp. 98–102, May 22–24, 2006.
[43] S. Moafipoor, D. A. Grejner-Brzezinska, and C. Toth, “Adaptive
calibration of a magnetometer compass for a personal navigation system,”
in International Global Navigation Satellite Systems Society (IGNSS
Symposium 2007), (The University of New South Wales, Sydney, Australia),
December 2007.
[44] M. Cordoba, “Attitude and heading refernce system i-ahrs for the efigenia
autonomous unmanned aerial vehicles uav based on mems sensor and a
neural network strategy for attitude estimation,” in Control Automation,
2007. MED ’07. Mediterranean Conference on, pp. 1–8, June 2007.
[45] P. Martin and E. Salaun, “Invariant observers for attitude and heading
estimation from low-cost inertial and magnetic sensors,” in Decision and
216
References
Control, 2007 46th IEEE Conference on, pp. 1039–1045, Dec 2007.
[46] R. E. Kalman, “A new approach to linear filtering and prediction problems,”
Journal of Basic Engineering, vol. 82, pp. 35–45, 1960.
[47] R. G. Brown and P. Y. C. Hwang, Introduction to Random Signals and
Applied Kalman Filtering. John Wiley & Sons, Inc., third edition ed., 1997.
[48] Xsens Technologies B.V., MTi and MTx User Manual and Technical
Documentation, mt0100p, revision O ed., October 2010.
[49] LORD Corporation, MicroStrain Sensing Systems, 459 Hurricane Lane,
Suite 102, Williston, VT 05495 USA, 3DM-GX3-45 - Miniature GPS-Aided
Inertial Navigation System (GPS/INS), version 8400-0035 revision 003 ed.,
2014.
[50] InterSense, LLC, 700 Technology Park Drive, Suite 102, Billerica, MA 01821,
USA, InertiaCube4 - Precision inertial orientation sensor, 072-00136-0k12
revision 1.0 ed.
[51] VectorNav Technologies, LLC, VN-100 User Manual, revision 1.2.10 ed.,
2009.
[52] CH Robotics LLC, UM6 Ultra-Miniature Orientation Sensor Datasheet,
revision 2.4 ed., October 2013.
[53] YEI Technology, 630 Second Street, Portsmouth, Ohio 45662, 3-Space Sensor
Embedded Ultra-Miniature Attitude & Heading Reference System User’s
Manual, August 2013.
[54] R. B. McGhee, E. R. Bachmann, X. Yun, and M. J. Zyda, “Real-time
tracking and display of human limb segment motions using sourceless sensors
and a quaternion-based filtering algorithm - part i: Theory,” tech. rep.,
NAVAL POSTGRADUATE SCHOOL, Monterey, California, 2000.
[55] R. Mahony, T. Hamel, and J.-M. Pflimlin, “Complementary filter design
on the special orthogonal group so(3),” in Decision and Control, 2005 and
2005 European Control Conference. CDC-ECC ’05. 44th IEEE Conference
on, pp. 1477–1484, Dec 2005.
[56] T. S. Yoo, S. K. Hong, H. M. Yoon, and S. Park, “Gain-scheduled
complementary filter design for a mems based attitude and heading reference
system,” Sensors, vol. 11, no. 4, pp. 3816–3830, 2011.
[57] A.-J. Baerveldt and R. Klang, “A low-cost and low-weight attitude
estimation system for an autonomous helicopter,” in Intelligent Engineering
Systems, 1997. INES ’97. Proceedings., 1997 IEEE International Conference
on, pp. 391–395, Sep 1997.
[58] J. Roberts, P. Corke, and G. Buskey, “Low-cost flight control system
for a small autonomous helicopter,” in Robotics and Automation, 2003.
Proceedings. ICRA ’03. IEEE International Conference on, vol. 1,
pp. 546–551 vol.1, Sept 2003.
217
References
218
References
219
References
220
References
221
References
222
References
interferometric fiber optic gyros,” IEEE Std 952-1997, pp. i–, 1998.
[126] D. W. Stockwell, “Bias stability measurement: Allan variance,” tech. rep.,
Crossbow Technology, Inc., 2003.
[127] I. Skog and P. Händel, “Calibration of a mems inertial measurement unit,”
in in Proc. XVII IMEKO WORLD CONGRESS, (Rio de Janeiro, 2006.
[128] Xsens Technologies B.V., MTw User Manual, revision F ed., January 2013.
[129] “UM6 orientation sensor product webpage.” http://www.chrobotics.com/
shop/orientation-sensor-um6, 2014. [Accessed: February 25, 2014].
[130] Ideal Aerosmith, Inc., East Grand Forks, MN USA, Grand Forks, ND USA,
1270VS Series Single-Axis Rate Table Datasheet, rev: q ed.
[131] Panasonic Company Division of Matsushita Electric Corporation of America,
One Panasonic Way, Secaucus, New Jersey 07094, SL-1200MK2 Turntable
System Service Manual (M), (MC).
[132] Thurlby Thandar Instruments Ltd., Glebe Road, Huntingdon,
Cambridgeshire, PE29 7DR, England (United Kingdom), TF960 6GHz
Universal Counter Instruction Manual, 48581-1440 issue 2 ed.
[133] “Application note AN3192: Using LSM303DLH for a tilt compensated
electronic compass,” Tech. Rep. Doc ID 17353 Rev 1, STMicroelectronics,
August 2010.
[134] J. Vasconcelos, G. Elkaim, C. Silvestre, P. Oliveira, and B. Cardeira,
“A geometric approach to strapdown magnetometer calibration in sensor
frame,” Aerospace and Electronic Systems, IEEE Transactions on, vol. 47,
pp. 1293–1306, April 2011.
[135] Y. Petrov, “Matlab central - ellipsoid fit.” http://www.mathworks.
co.uk/matlabcentral/fileexchange/24693-ellipsoid-fit, July 2009.
[Accessed: March 26, 2014].
[136] InvenSense, Inc., MPU-6000 and MPU-6050 Register Map and Descriptions
Revision 4.2, revision: 4.2 ed., 2013.
[137] Y. Xiaoping, E. Bachmann, and R. McGhee, “A simplified quaternion-based
algorithm for orientation estimation from earth gravity and magnetic field
measurements,” Instrumentation and Measurement, IEEE Transactions on,
vol. 57, pp. 638 –650, mar. 2008.
[138] N. S. Nise, Control Systems Engineering. John Wiley & Sons, Inc., fourth
edition ed., 2004.
[139] STMicroelectronics, TA0343 Technical article - Everything about
STMicroelectronics’ 3-axis digital MEMS gyroscopes, doc id 022032
rev 1 ed., July 2011.
[140] D. Vargha and M. Maia, “An overview of motion processing solutions for
consumer products,” tech. rep., InvenSense, I n c., 1197 Bor regas Avenue,
223
References
Sunnyvale, CA 94089.
[141] Microchip Technology Inc., dsPIC33FJ32GP302/304,
dsPIC33FJ64GPX02/X04, and dsPIC33FJ128GPX02/X04 Data Sheet,
DS70292G ed., 2012.
[142] Maxim Integrated Products, MAX1811 - USB-Powered Li+ Charger, rev
2 ed., 2003.
[143] Future Technology Devices International Limited, FT232R USB UART IC
Datasheet, version 2.10 ed., 2012.
[144] Future Technology Devices International Ltd., TN 132 Adding FTDI Devices
VCP Driver Support to Android, version 1.1 ed., 2011.
[145] Roving Networks, Inc., RN41/RN41N Class 1 Bluetooth Module, version
3.42r ed., November 2013.
[146] x-io Technologies Limited, x-IMU User Manual, version 5.2 ed., November
2013.
[147] bluetooth.org, BLUETOOTH SPECIFICATION Version 2.1 + EDR, July
2007.
[148] Microchip Technology Inc., MCP73831/2 Datasheet, DS20001984F ed., 2013.
[149] Micrel Inc., MIC5319 Datasheet, M9999-071812 ed., July 2012.
[150] Microchip Technology Inc., PIC24FJ64GA104 Family Data Sheet,
DS39951C ed., 2010.
[151] Digi International Inc., XBee XBee-PRO RF Modules: Product Manual
v1.xEx - 802.15.4 Protocol, 90000982 M ed., July 2013.
[152] Roving Networks, Inc., RN41XV & RN42XV Bluetooth Module, version
1.0 ed., December 2013.
[153] Digi International Inc., XBee XBee-PRO ZB RF Modules, 90000976 P ed.,
April 2013.
[154] Digi International Inc., XBee-PRO 900HPXBee-PRO XSC RF Modules,
90002173 H ed., September 2013.
[155] Digi International Inc., XBee Wi-Fi RF Module, 90002180 G ed., September
2013.
[156] B. Ivey, AN1416 - Low-Power Design Guide. Microchip Technology Inc.,
ds01416a ed., 2011.
[157] “IEEE standard for local and metropolitan area networks–part 15.4:
Low-rate wireless personal area networks (lr-wpans),” IEEE Std
802.15.4-2011 (Revision of IEEE Std 802.15.4-2006), pp. 1–314, 2011.
[158] APDM Inc., 2828 Southwest Corbett Avenue, Portland, OR 97201, opal -
Technical Specifications, 2012.
224
References
225
References
the knee during walking after stroke,” Physiotherapy, vol. 92, no. 3, pp. 159
– 165, 2006.
[174] A. Vega-Gonzàlez and M. H. Granat, “Continuous monitoring of upper-limb
activity in a free-living environment,” Archives of Physical Medicine and
Rehabilitation, vol. 86, no. 3, pp. 541 – 548, 2005.
[175] K. Aminian, B. Najafi, C. Büla, P. F. Leyvraz, and P. Robert,
“Spatio-temporal parameters of gait measured by an ambulatory system
using miniature gyroscopes,” Journal of Biomechanics, vol. 35, no. 5, pp. 689
– 699, 2002.
[176] F. B. V. M. M. S. T. E. Uswatte G, Miltner WH, “Objective measurement
of functional upper-extremity movement using accelerometer recordings
transformed with a threshold filter,” Stroke, pp. 662–667, 2000.
[177] G. Uswatte, W. L. Foo, H. Olmstead, K. Lopez, A. Holand, and L. B.
Simms, “Ambulatory monitoring of arm movement using accelerometry: An
objective measure of upper-extremity rehabilitation in persons with chronic
stroke,” Archives of Physical Medicine and Rehabilitation, vol. 86, no. 7,
pp. 1498 – 1501, 2005.
[178] E. Bernmark and C. Wiktorin, “A triaxial accelerometer for measuring arm
movements,” Applied Ergonomics, vol. 33, no. 6, pp. 541 – 547, 2002.
[179] A. V. Rowlands, P. W. M. Thomas, R. G. Eston, and R. Topping, “Validation
of the RT3 triaxial accelerometer for the assessment of physical activity,”
Medicine and Science in Sports and Exercise, vol. 36, pp. 518–524, March
2004.
[180] S. Y. Kumahara H, Tanaka H, “Daily physical activity assessment: What is
the importance of upper limb movements vs whole body movements,” Int J
Obes, pp. 1105–1110, 2004.
[181] P. Veltink and H. Franken, “Detection of knee unlock during stance by
accelerometry,” Rehabilitation Engineering, IEEE Transactions on, vol. 4,
pp. 395 –402, Dec. 1996.
[182] X. Yun, E. R. Bachmann, H. Moore, and J. Calusdian, “Self-contained
position tracking of human movement using small inertial/magnetic sensor
modules,” in ICRA, pp. 2526–2533, 2007.
[183] E. Foxlin, “Pedestrian tracking with shoe-mounted inertial sensors,” IEEE
Comput. Graph. Appl., vol. 25, no. 6, pp. 38–46, 2005.
[184] H. M. Schepers, H. Koopman, and P. H. Veltink, “Ambulatory assessment of
ankle and foot dynamics,” IEEE Transactions on Biomedical Engineering,
vol. 54, no. 5, pp. 895–902, 2007.
[185] R. Stirling, K. Fyfe, and G. Lachapelle, “Evaluation of a new method
of heading estimation for pedestrian dead reckoning using shoe mounted
sensors,” The Journal of Navigation, vol. 58, no. 01, pp. 31–45, 2005.
226
References
227
References
2013.
[197] W.-H. Chen, B.-Y. Huang, S.-R. Chen, S.-F. Yang, C.-C. Lee, D.-S. Yu,
, and Y.-H. Lin, “FPGA and virtual reality based minimally invasive
surgery training system,” in Proceedings of 2013 FPGA Workshop an Design
Contest, (Southeast UNiversity, Nanjing, China), November 2013.
[198] J. S. Akins, DEVELOPMENT AND EVALUATION OF INSTRUMENTED
SOCCER EQUIPMENT TO COLLECT ANKLE JOINT KINEMATICS IN
THE FIELD. PhD thesis, Swanson School of Engineering, University of
Pittsburgh, 2013.
[199] C. Brighton, E. Olsen, and T. Pfau, “Is a standalone inertial measurement
unit accurate and precise enough for quantification of movement symmetry
in the horse?,” Computer Methods in Biomechanics and Biomedical
Engineering, vol. 0, pp. 1–6, 2013.
[200] F. Hegel, F. Eyssel, and B. Wrede, “The social robot ‘flobi’: Key concepts
of industrial design,” in RO-MAN, 2010 IEEE, pp. 107–112, 2010.
[201] I. Lutkebohle, F. Hegel, S. Schulz, M. Hackel, B. Wrede, S. Wachsmuth,
and G. Sagerer, “The bielefeld anthropomorphic robot head “Flobi”,” in
Robotics and Automation (ICRA), 2010 IEEE International Conference on,
pp. 3384–3391, 2010.
[202] S. Schulz, F. Lier, I. Lutkebohle, and S. Wachsmuth, “Robot reality - a
motion capture system that makes robots become human and vice versa,”
in Robotics and Automation (ICRA), 2013 IEEE International Conference
on, pp. 2126–2133, 2013.
[203] L. Sussenbach, K. Pitsch, I. Berger, N. Riether, and F. Kummert, ““Can you
answer questions, Flobi?”: Interactionally defining a robot’s competence as
a fitness instructor,” in RO-MAN, 2012 IEEE, pp. 1121–1128, 2012.
[204] F. Eyssel and N. Reich, “Loneliness makes the heart grow fonder (of
robots) - on the effects of loneliness on psychological anthropomorphism,”
in Human-Robot Interaction (HRI), 2013 8th ACM/IEEE International
Conference on, pp. 121–122, 2013.
[205] T. Michailidis, D. Polydorou, and J. Bullock, “A multimodal integration of
sensory feedback modalities for dance performers,” in Proceedings of the 1st
Fascinate Conference, (Falmouth University, Cornwall, UK), 28-30 August
2013.
[206] S. Kasaei, S. Kasaei, and S. Kasaei, “Persia: A fully autonomous humanoid
robot based on hybrid ann controller and vision system,” Mechatronics, 2013.
[under review].
[207] Glowacki, D. R., Tew, P., Mitchell, T., Kriefman, L, Hyde, J., Malcolm,
L. J., Price, J., McIntosh-Smith, and S., “Sculpting molecular dynamics in
real-time using human energy fields,” Molecular Aesthetics, 2013.
228
References
[208] Glowacki, D. R., Tew, P., Mitchell, T., Price, J., McIntosh-Smith, and
S., “danceroom spectroscopy: Interactive quantum molecular dynamics
accelerated on gpu architectures using opencl,” in UK Many Core
Development Conference 2012 (UKMAC ’12), (Bristol), 2012.
[209] A. Roudaut, R. Reed, T. Hao, and S. Subramanian, “Changibles: Analyzing
and designing shape changing constructive assembly,” in SIGCHI Conference
on Human Factors in Computing Systems (CHI ’14), (ACM, New York, NY,
USA), 2014. under review.
[210] P. Cardou and J. Angeles, “Singularity analysis of accelerometer strapdowns
for the estimation of the acceleration field of a planar rigid-body motion,”
in 12th IFToMM World Congress, 2007.
[211] T. R. Williams and K. R. Fyfe, “Planar accelerometer configurations,”
Journal of Applied Mechanics, vol. 71, no. 1, pp. 10–14, 2004.
[212] B. Zappa, G. Legnani, A. J. van den Bogert, and R. Adamini, “On
the number and placement of accelerometers for angular velocity and
acceleration determination,” Journal of Dynamic Systems, Measurement,
and Control, vol. 123, no. 3, pp. 552–554, 2001.
[213] G. Baselli, G. Legnani, P. Franco, F. Brognoli, A. Marras, F. Quaranta, and
B. Zappa, “Assessment of inertial and gravitational inputs to the vestibular
system,” Journal of Biomechanics, vol. 34, no. 6, pp. 821 – 826, 2001.
[214] P. Spevak and P. Forstner, “Msp430 32-khz crystal oscillators,” Tech. Rep.
SLAA322B, Texas Instruments Incorporated, August 2006.
[215] Micro Crystal Switzerland, Micro Crystal AG, Mühlestrasse 14, CH-2540
Grenchen, Switzerland, Micro Crystal Switzerland - CC7V-T1A Tuning Fork
Crystal 32.768 kHz, version 5 ed., 2009.
[216] K. Kudapali, AN1155 - Run-Time Calibration of Watch Crystals, 2008.
[217] S. Bible, AN826 - Crystal Oscillator Basics and Crystal Selection for
rfPICTM and PICmicro Devices. Microchip Technology Inc., ds00826a ed.,
2002.
229
References
230
Appendix A
T
i
α = αx αy (A.1)
1 0 -dx -dy
Gi = (A.2)
0 1 -dy dx
T
s= ẍ ÿ ω 2
ω̇ (A.3)
231
Appendix A. Non-gyro IMU planar simplification
For the simplified planar system, the objective function is simply the difference
between the squared estimated angular velocity term and the squared angular
velocity term yielded by the state vector s. The error originally defined in Equation
4.24 in the estimated angular velocity may be redefined as the scalar quantity of
Equation A.4.
√
ẽ = sign(ω̃(ω̃ 2 − ω 2 )) |ω̃| − ω2 (A.4)
For ideal accelerometers, all valid array configurations perform equally. In practice,
errors due to noise and signal saturation in accelerometer measurements will mean
that some array configurations perform better than others. We wish to evaluate
the effect of individual array characteristics to understand what determines an
optimal array. An accelerometer array may be characterised by the following:
Qin et al. [87] use inspection of matrix singular values to determine the suitability
of a sensor configuration and determined that an optimal design is represented by
an isotropic matrix; that is, a matrix with a condition number of one where all
singular values are identical and nonzero. This is achieved by a sensor constellation
represented by the vertices of a Platonic solid. The analysis of other authors exist
for planar [210, 211] and 3D [212] [213] accelerometer arrays.
232
Appendix A. Non-gyro IMU planar simplification
To analyse the effect of the kinematic origin location and array volume we will
consider a simplified planar array of 2 dual-axis accelerometers represented by the
measurements α10 to α40 , each incorporating an error (δ1 to δ4 ). The kinematic origin
of the array is equidistant from each sensor and may be displaced by a distance p
along the origin y axis (remaining equidistant from each sensor). The volume of
the array is defined by the distance of separation, l. The system is described by
the schematic shown in Figure A.1. Equation A.5 describes the forward kinematic
solution.
ÿ
ω
α4 α2
ẍ
p
α3 α1
1 1
2l 2l
α0 = Hs + δ (A.5)
0
α1 1 0 − 2l p ẍ δ1
α0 0 1 p l
ÿ
2 2 δ2
= + (A.6)
α0 1 l
p
3 0 2 ω 2 δ3
α40 0 1 p −2l
ω̇ δ4
The inverse kinematic solution may be arranged to provide the measured kinematic
states of the body origin in terms of the true kinematic states and measurement
0
errors as shown by Equations A.7 and A.8. In these equations, ẍ0 , ÿ 0 , ω 2 and ω̇ 0
233
Appendix A. Non-gyro IMU planar simplification
s0 = H −1 α0
= H −1 (Hs + δ) (A.7)
= s + H −1 δ
0
ẍ ẍ
1
2
− pl 1
2
p
l δ1
p
ÿ 0 ÿ 1
− pl 1
δ
l 2 2 2
= + (A.8)
ω 2 0 ω 2 − 1 0 1
0
l l δ3
ω̇ 0 ω̇ 0 + 1l 0 − 1l δ4
For this simplified equidistant configuration it can be seen from Equation A.8 how
the distance from the origin (p) and sensor separation (l) may be manipulated to
minimise the error in the measured kinematics states. The error in the measured
0
angular states (ω 2 and ω̇ 0 ) is independent of the origin displacement and inversely
proportional to the sensor separation. Error components in the measured linear
states (ẍ0 , ÿ 0 ) are proportional to the origin displacement and inversely proportional
to the sensor separation. This method extends to assessing redundant planar
sensor arrays, and leads to a weighted average (linear accelerations) or a weighted
sum of difference (angular acceleration). Thus it can be demonstrated that an
optimal array would maximise the sensor separation. The kinematics origin should
be located at the point where the expected mean acceleration is minimised. Where
more specific information is not available this could be considered as the centre of
mass of the limb. Other authors have arbitrarily chosen the centre of volume as
the kinematic origin of a sensor array [86, 87], so assuming a homogeneous body
the methodology we describe gives a justification for their decision.
Analysis of the equidistant planar structure supports Qin’s assertion [87] that a
reasonable sensor array would be to distribute the sensors on the vertices of a
Platonic solid. This only applies if the body is assumed to be rotating around the
234
Appendix A. Non-gyro IMU planar simplification
centre of the sensor constellation. Where more sensors than vertices are used in
the IMU, a consideration of the singular values of the reconstruction matrix is an
appropriate method to optimise sensor placement.
235
Appendix A. Non-gyro IMU planar simplification
236
Appendix B
B.1 Introduction
These investigations focusses specially on the 32.678 kHz crystal used by the
x-IMU. The is the clock source for both the RTC and sample clock. Crystal
oscillators typically have an error of around ±20 ppm though this will increase with
changes in temperature. Temperature Compensated Crystal Oscillator (TCXO)
and Oven-Controlled Crystal Oscillator (OCXO) provide far greater accuracy over
a temperature range but are large, expensive and require allot of power. IMU
237
Appendix B. Crystal oscillator thermal response
Equation B.1 [214] describes the crystal error model where e is the period error in
ppm, e0 is the error at the turnover frequency, c is the temperature coefficient, T
is the temperature of the crystal and T0 is the turnover temperature. Table B.1
summarises the values of each of these parameters are stated in the datasheet of
the crystal used by the x-IMU [215].
e = e0 + c(T − T0 )2 (B.1)
A typical error of ±20 ppm at 25◦C can only be achieved if the design of the
crystal oscillator circuit accounts for electrical characteristics of the PCB layout.
For example, a poorly matched load capacitance can result in an error of almost
400 ppm [216]. Component values were initially selected according to manufacturer
application notes [217] but were found to provide errors of the order of 150 ppm
at 25◦C Instead, component values were selected as those found to provide a error
centred around 0 ppm 25◦C for the eight tested devices.
238
Appendix B. Crystal oscillator thermal response
The thermal response of eight crystals were characterised using the thermal chapter
described in Chapter 5. Each device was exposed to a monotonic ascent from -20◦C
to 70◦C over 3 hours. A TTi TF930 frequency counter was used to measure the
crystal frequency with an accuracy of ±0.2 ppm [132]. Software was written to log
the device temperature along side the crystal error. Figure B.1 shows measured
period error for each device over the temperature range.
A B C D E F G H
100
90
80
Crystal period error (ppm)
70
60
50
40
30
20
10
−10
−20 −15 −10 −5 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70
Temperature (◦ C)
A second-order polynomial was fit was performed for each device to determine the
model parameters. The results are summarised in Table B.2.
Table B.2: Crystal error model parameters for eight devices over temperature
range
239
Appendix B. Crystal oscillator thermal response
The simplest calibration solution is to use the generic error model defined by
Equation B.1 and the typical parameters specified in Table B.1. The resultant
error of the eight devices is shown as Figure B.2.
A B C D E F G H
30
25
Crystal period error (ppm)
20
15
10
−5
−10
−20 −15 −10 −5 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70
Temperature (◦ C)
Figure B.2: Crystal oscillator error using the generic calibration model and typical
parameters specified in the datasheet
This simple calibration solution is able to reduce the crystal error from 90 ppm to
25 ppm. This solution can be implemented in software and does not require any
specific values t be evaluated through testing.
240
Appendix B. Crystal oscillator thermal response
A significant reduction in error can be achieved if the generic error model is used
in conjunction with the specific specific parameters evaluated for each device and
summarised in Table B.2. The resultant error of the eight devices is shown as
Figure B.3. This calibration solution reduces the errors to within ppm.
A B C D E F G H
1.5
1
Crystal period error (ppm)
0.5
−0.5
−1
−1.5
−2
−2.5
−20 −15 −10 −5 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70
Temperature (◦ C)
Figure B.3: Crystal oscillator error using the generic calibration model and specific
parameters evaluated empirically
The curve of the error shown in Figure B.3 indicates that the second-order generic
error model is too low an order. This may due to an inaccuracy in the model form
or may be due to a non-linearity in the thermometer. The source of the error is
irreverent provide that a calibration solution accounts for the characteristics of
both the thermometer and crystal. Figure B.4 shows the crystal oscillator error
using a third-order calibration model.
This third-order calibration model is able to achieve an error of 0.2 ppm. Initial
large error is thought to be because of an initial temperature difference between
the thermometer device and the crystal oscillator.
241
Appendix B. Crystal oscillator thermal response
A B C D E F G H
0.5
0.4
0.3
Crystal period error (ppm)
0.2
0.1
−0.1
−0.2
−0.3
−0.4
−0.5
−20 −15 −10 −5 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70
Temperature (◦ C)
B.5 Conclusion
242
Appendix C
The Earth’s magnetic field is not constant around the world and it is of interest
to consider how this will affect the accuracy of a heading measurement obtained
from a magnetometer-based AHRS. The field varies in declination (deviation from
geodetic north), inclination (deviation from the horizontal plane) and in intensity.
Given that a magnetometer-based AHRS is intended to measure heading relative
to magnetic north and not geodetic north, declination variations are irreverent.
Maus et al. provide a comprehensive explanation and quantification of the Earth’s
magnetic field in The US/UK World Magnetic Model for 2010-2015 [30]. It is
reported that the intensity can be expected to between 0.22 Gauss and 0.67 Gauss
around the world. The varying inclination is presented as the Mercator projection
overlay shown as Figure C.1. The contours indicate that non-polar regions can
expect a maximum inclination of ∼80◦ , with the 0◦ inclination contour remaining
with ±15◦ latitude of the equator. The inclination in the UK is indicted as being
between 66◦ and 70◦ .
243
Appendix C. Effect of geographic location on AHRS heading accuracy
80 80
60°N 60°N
60
45°N 45°N
60 60
30°N 30°N
40 40 40
20 20
15°N 15°N
20
0
0
0 -20
0° -20 0°
-20 -40
-40
15°S 15°S
-40 -60
-60
30°S 30°S
-60
45°S 45°S
-80
-60
60°S 60°S
kj
-80
-60
70°S 70°S
180° 135°W 90°W 45°W 0° 45°E 90°E 135°E 180°
Figure C.1: Contour interval: 2 degrees, red contours positive (down); blue
negative (up); green zero line. Image source: [30]
x y
θ = atan2 , (C.2)
||m|| ||m||
244
Appendix C. Effect of geographic location on AHRS heading accuracy
will be subject to measurement errors from a variety of sources. For the purposes of
this investigation we will consider an error of e in the magnetometer magnetometer
x axis to yield the measurement m̃ as described by Equation C.3.
x̃ e
m̃ = = m + (C.3)
ỹ 0
Using Equations C.1, C.3 and C.4, it is possible to determine θerror for selected
values of b and γ. The relationship of Equation C.1 suggests that θerror will have
a greater sensitivity to variations in γ and so we will first consider a numerical
example for only variations in γ. Fixed values of b = 0.5 Gauss and e = 0.0044
Gauss were chosen as this ratio was found to provided a maximum heading error
±0.5◦ for an inclination angle of 0◦ ; this is heading accuracy specified by high-end
commercial IMUs such as the MicroStrain 3DM-GX3-25 [108]. Figure C.2 shows
the corresponding heading error achieved for an inclination of 0◦ (on the equator),
30◦ , 70◦ (in the UK) and 80◦ .
245
Appendix C. Effect of geographic location on AHRS heading accuracy
2
θerror (degrees)
−1
−2
−3
0 60 120 180 240 300 360
θ (degrees)
Figure C.2: Heading error achieved for an inclination of 0◦ (on the equator), 30◦ ,
70◦ (in the UK) and 80◦ and fixed values of b = 0.5 Gauss and e = 0.0044 Gauss
Location b γ θerror
Phuket, Thailand 0.41 Gauss 0◦ ±0.6◦
Mina, Peru 0.26 Gauss 0◦ ±1.0◦
Edinburgh, UK 0.5 Gauss 70◦ ±1.5◦
Trua, Russia 0.6 Gauss 80◦ ±2.5◦
Table C.1: Extreme values of b and γ selected from [30] to demonstrate the least
disparity in θerror for differing geographic locations
The results in Table C.1 confirm that the variation in heading performance for
differing geographic locations may be less than that assumed by an inclination-only
model but that the specific geographic location will still have a significant effect on
the heading performance. A heading error performance observed at the equator
may correspond to an error almost times greater in the UK and equally a heading
accuracy observed in the UK may be expected to improve by a factor of three
when operating on the equator.
246