See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/273383944

A robust and easy to implement method for IMU calibration without external
equipments

Conference Paper in Proceedings - IEEE International Conference on Robotics and Automation · September 2013
DOI: 10.1109/ICRA.2014.6907297

CITATIONS READS
236 19,854

3 authors:

David Tedaldi Alberto Pretto

ETH Zurich Sapienza University of Rome
2 PUBLICATIONS 340 CITATIONS 78 PUBLICATIONS 1,879 CITATIONS

SEE PROFILE SEE PROFILE

Emanuele Menegatti
University of Padova
288 PUBLICATIONS 4,720 CITATIONS

SEE PROFILE

All content following this page was uploaded by Alberto Pretto on 11 September 2017.

The user has requested enhancement of the downloaded file.

A Robust and Easy to Implement Method for IMU Calibration without
External Equipments
David Tedaldi, Alberto Pretto and Emanuele Menegatti

Abstract— Motion sensors as inertial measurement units

(IMU) are widely used in robotics, for instance in the navigation
and mapping tasks. Nowadays, many low cost micro electro
mechanical systems (MEMS) based IMU are available off the
shelf, while smartphones and similar devices are almost always
equipped with low-cost embedded IMU sensors. Nevertheless,
low cost IMUs are affected by systematic error given by
imprecise scaling factors and axes misalignments that decrease
accuracy in the position and attitudes estimation.
In this paper, we propose a robust and easy to implement
method to calibrate an IMU without any external equipment.
The procedure is based on a multi-position scheme, providing
scale and misalignments factors for both the accelerometers
and gyroscopes triads, while estimating the sensor biases. Our
method only requires the sensor to be moved by hand and
placed in a set of different, static positions (attitudes). We
describe a robust and quick calibration protocol that exploits
an effective parameterless static filter to reliably detect the static
intervals in the sensor measurements, where we assume local
stability of the gravity’s magnitude and stable temperature.
We first calibrate the accelerometers triad taking measurement
samples in the static intervals. We then exploit these results
to calibrate the gyroscopes, employing a robust numerical
integration technique.
The performances of the proposed calibration technique has
been successfully evaluated via extensive simulations and real
experiments with a commercial IMU provided with a calibra-
tion certificate as reference data.

I. I NTRODUCTION
IMUs (Inertial Measurement Units) are very popular
sensors in robotics: among others, they have been exploited Fig. 1. Calibration protocol: (top) Some examples of the Xsens MTi IMU
disposed in different attitudes as required by our method; (down) Diagram
for inertial-only navigation [1], attitude estimation [2], and of our calibration protocol
visual-inertial navigation [3], [4], also using a smartphone
device [5]. IMUs used in robotics are usually based on
MEMS (micro electro mechanical systems) technology.
sensitivities, and non zero biases. The IMU calibration
They are composed by a set of tri-axial clusters: an
refers to the process of identifying these quantities.
accelerometers, a gyros and often a magnetometer cluster.
In an ideal IMU, the tri-axial clusters should share the same
Many commercial IMU in the cost range form 1000 $
3D orthogonal sensitivity axes that span a three dimensional
to 2000 $, such as the Xsens MTi [6] exploited in the
space, while the scale factor should convert the digital
experiments (Sec. V), are factory calibrated 1 . Each sensor
quantity measured by each sensor into the real physical
is sold with its own calibration parameter set stored into the
quantity (e.g, accelerations and gyro rates). Unfortunately,
firmware or inside a non-volatile memory, providing accurate
low cost MEMS based IMU are usually affected by non
measurements off the shelf. Unfortunately, the overhead
accurate scaling, sensor axis misalignments, cross-axis
cost for the factory calibration is predominant: usually the
This research has been partially supported by Consorzio Ethics with sensor hardware (sensors, chips, embodiment, . . . ) is likely
the grant ”Rehabilitation Robotics”, by Univ. of Padova with the grant to be only a fraction of the final device cost. Actually, the
“TIDY-UP: Enhanced Visual Exploration for Robot Navigation and Ob-
ject Recognition”, and by the European Commission under FP7-600890-
factory calibration is usually performed using standard but
ROVINA. The authors are with the Department of Information Engineering, effective methods, where the device outputs are compared
University of Padova, Italy. Email: [email protected], with known references: this process requires time for each
[email protected]. Pretto is also with the Department of Computer,
Control, and Management Engineering “Antonio Ruberti“, Sapienza Univer-
sity of Rome, Italy. Email: [email protected] 1 Often they are also compensated over temperature
sensor and a high cost equipment. On the other hand, low- drifts in a period estimated using the Allan variance
cost IMUs (20-100 $) and the IMU sensors that equip current • We introduce a simple but effective static detector
smartphones are usually poorly calibrated, resulting in mea- that exploits the sensor noise magnitude, a fixed-time
surements coupled with not negligible systematic errors. For sampling window and a cutting threshold, automatically
instance, state-of-the-art visual-inertial navigation systems estimated inside the optimization framework
such the one presented in [5], that exploits a smartphone as • We employ the Runge-Kutta numerical integration
experimental platform, while performing so well in forward, method to improve the accuracy of the gyroscope cali-
almost regular, motion 2 , shows lower performances in more bration.
“exciting” motions, i.e. in motions that quickly change linear We have extensively tested our system using synthetic data
acceleration and rotational axes. affected by variable biases, misalignments, scale factor er-
In this paper, we propose an effective and easy to imple- rors, and noise. In all of the cases, we have obtained
ment calibration scheme, that only requires to collect IMU stable and accurate results. Moreover, we have performed
data with the simple procedure described in the flow chart the calibration of a commercial, factory calibrated Xsens
reported in Fig. 1. After an initial initialization period with MTi IMU, using its raw, uncalibrated data as input. Our
no motion, the operator should move the IMU in different calibration results are comparable to the factory parameters
positions, in order to generate a set of distinct, temporarily reported by the device’s calibration certificate.
stable, rotations. The collected dataset is used to calibrate the
scale and misalignments factors for both the accelerometers A. Related Works
and gyroscopes triads, while estimating the sensor biases. Traditionally the calibration of an IMU has been done by
As other calibration technique, we neglect the effect of the using special mechanical platforms such as a robotic ma-
cross-axis sensitivities, since for minor misalignments and nipulator, moving the IMU with known rotational velocities
minor cross-axis sensitivities errors it is usually difficult to in a set of precisely controlled orientations [9], [10], [11].
distinguish between them. At each orientation, the output of the accelerometers are
Our procedure exploits the basic idea of the multi-position compared with the precomputed gravity vector while during
method, firstly presented in [7] for accelerometers calibra- the rotations the output of gyroscopes are compared with the
tion: in a static position, the norms of the measured acceler- precomputed rotational velocity. However, the mechanical
ations is equal to the magnitudes of the gravity plus a multi- platforms used for calibration are usually very expensive,
source error factor (i.e., it includes biases, misalignment, resulting in a calibration cost that often exceeds the cost of
noise,...). All these quantities can be estimated via minimiza- the IMU’s hardware.
tion over a set of static attitudes. After the calibration of the In [12] a calibration procedure that exploits a marker-
accelerometer triad, we can use the gravity vector positions based optical tracking system has been presented, while in
measured by the accelerometers as a reference to calibrate the [13], the GPS readings are used to calibrate initial biases
gyroscope triad. Integrating the angular velocities between and misalignments. Clearly, the accuracy of these methods
two consecutive static positions, we can estimate the gravity strongly depends on the accuracy of the employed kinematic
positions in the new orientation. The gyroscopes calibration reference (i.e., the motion capture system or the GPS). The
is finally obtained minimizing the errors between these multi-position method was firstly introduced by Lotters et al.
estimates and the gravity references given by the calibrated [7]: they proposed to calibrate the biases and the scale factor
accelerometers. of the accelerometers using the fact that the magnitude of
In this procedure the gyroscopes calibration accuracy the static acceleration must equal to the gravity’s magnitude.
strongly depends on the accuracy of the accelerometers cali- This technique has been extended in [14] and [15] to include
bration, being used as a reference. Moreover, signal noise and the accelerometer axis misalignment. The error model they
biases should negatively affect both the calibration accuracy proposed for the gyroscopes is similar to the one used for
and the reliability of the algorithms used to detect the actual the accelerometers, but the calibration procedure in this
static intervals used in the calibration. Finally, a consistent case requires a single axis turntable to provide a strong
numerical integration process is essential to mitigate the rotation rate signal, providing high calibration accuracy. Un-
effect of the signal discretization, usually sampled at 100 fortunately, these approaches not only require a mechanical
Hz. In our approach, we face these problems introducing equipment, but the two triads are independently calibrated,
the following modifications to the standard multi-position and the misalignment between them can’t be detected. In
method: [8] and [16] authors presented two calibration schemes that
• The proposed calibration protocol exploits a larger do not require any external mechanical equipment. Similarly
number of static states with reduced periods, in order to to our approach, in the first work the authors calibrate
increase the cardinality of the dataset while preserving the accelerometers exploiting the high local stability of the
the assumption of local stability of the sensors biases gravity vector’s magnitude, and then gyroscopes calibration
• As proposed in [8], we characterize the gyroscope bias is obtained comparing the gravity vector sensed by the
calibrated accelerometer with the gravity vector obtained by
2 Actually, during an almost regular motion miscalibration errors may integrating the angular velocities. In the second work the
easily be assimilated by the biases included in the system state authors also exploit the local stability of the magnetic field.
Hwangbo et al. [17] recently proposed a self-calibration
technique based on an iterative matrix factorization. They
use gravity as accelerometers reference, and a camera as
gyroscopes reference.
II. S ENSOR E RROR M ODEL
For an ideal IMU, the 3 axes of the accelerometers
triad and the 3 axes of the gyroscopes triad define a
single, shared, orthogonal 3D frame. Each accelerometer
senses the acceleration along of one distinct axis, while
each gyroscope measures the angular velocity around the
same axis. Unfortunately in real IMUs, due to assembly
inaccuracy, the two triads form two distinct (i.e., misaligned),
non-orthogonal, frames. Also the single sensors are not
perfect: typically the scaling factors used to convert the
digital outputs of the sensors in real physical quantity Fig. 2. Non-orthogonal sensor (accelerometers or gyroscopes) axes (xS ,
are different for different instances of the same sensors, y S , z S ), and body frame axes (xB , y B , z B ).
while only a default, nominal scaling factor is provided
by manufacturers. Moreover, the output signals are almost
always affected by non zero, variable biases. where we have changed letter β, referring to the general
case, with the letter α, referring to the accelerometer case,
As introduced above, both the accelerometers frame (AF) while aO and aS denote the specific acceleration in AOF
and the gyroscopes frame (GF) are usually non-orthogonal. and AF, respectively3 .
We can define two associated orthogonal, ideal frames (AOF
and GOF, respectively) in the following way: As mentioned before, gyroscope and accelerometer mea-
• The x-axis of the AOF and the one of the AF coincide surements should refer to the same reference frame, in our
• The y-axis of the AOF lies in the plan spanned by the case AOF. Then, using Eq. 1, for the gyroscopes, we have:
x and y axes of the AF.  
1 −γyz γzy
For the gyroscopes case, it is sufficient to substitute the AF ω O = Tg ω S , Tg =  γxz 1 −γzx  (3)
and AOF acronyms with GF and GOF, respectively. Finally, −γxy γyx 1
we define a body frame (BF), which is an orthogonal frame
that represents, for example, the coordinate frame of the where ω O and ω S denote the specific angular velocities in
IMU’s chassis. The body frame usually differs from the AOF and in GF , respectively.
AF and GF frames by small angles but, in general, there Both the accelerometers and the gyroscopes are affected
is no direct relation between them. For small angles, a by biases and scale errors. Two scaling matrix are introduced
measurement sS in a non-orthogonal frame (AF or GF) can  a   g 
sx 0 0 sx 0 0
be transformed in the orthogonal body frame as (for details
Ka =  0 say 0  , Kg =  0 sgy 0  (4)
of the derivation, see [18]): a
  0 0 sz 0 0 sgz
1 −βyz βzy
We introduce also two bias vector
sB = TsS , T =  βxz 1 −βzx  (1)  a  g
−βxy βyx 1 bx bx
ba = bay  , bg = bgy  (5)
where sB and sS denote the specific force (acceleration),
baz bgz
or equivalently the rotational velocity, in the body frame
coordinates and accelerometers (or gyroscopes) coordinates, The complete sensor error model is
respectively. Here βij is the rotation of the i-th accelerometer
aO = Ta Ka (aS + ba + ν a ) (6)
or gyroscope axis around the j-th BF axis, see Fig. 2.
On the other hand, the two orthogonal frames BF and for the accelerometers, and
AOF (and, equivalently, BF and GOF) are relate by a pure
rotation. ω O = Tg Kg (ω S + bg + ν g ) (7)
In the presented calibration method, we assume that the for the gyroscope, where ν g and ν g are the accelerometer
body frame BF coincides with the accelerometers orthogonal measurement noise and the gyroscope measurement noise,
frame AOF: in such case, the angles βxz , βxy , βyx become respectively.
zero, so in the accelerometers case Eq. 1 becomes:
  3 To relate the obtained calibration with a different body frame (e.g.
1 −αyz αzy
O a S a BF’), it is sufficient to estimate the rotation matrix that relate AOF to BF’,
a =T a , T = 0 1 −αzx  (2) for instance using the accelerometers outputs in three different orthogonal
0 0 1 orientations.
 III. BASIC C ALIBRATION F RAMEWORK calibrated accelerometer readings in the k-th static interval,
In order to calibrate the accelerometers triad, we need to and ug,k is the acceleration versor computed using Eq. 11
estimate the following unknown parameter vector: (i.e., integrating the angular velocities between the k − 1-th
and the k-th static intervals). We obtain θ gyro minimizing
θ acc = αyz , αzy , αzx , sax , say , saz , bax , bay , baz
 
(8) Eq. 13 with LM.
We define the following function: IV. C ALIBRATION P ROCEDURE
aO = h(aS , θ acc ) = Ta Ka (aS + ba ) (9) As introduced in Sec. I, the proposed calibration frame-
work requires to collect a dataset with the stream of the
Here we can neglect the measurements noise thanks to raw accelerometers and gyroscopes readings, taken while the
the fact that in our calibration procedure we apply signal operator moves the IMU in different static positions, in order
averaging in each static interval. to generate a set of distinct, temporarily stable, rotations.
As in the conventional multi-position scheme, we move the A simple diagram of our calibration protocol is reported in
IMU in a set of M distinct, temporarily stable, rotations. Fig. 1. As mentioned in Sec. III, to mitigate the noise effect
We can extract M acceleration vectors aSk (measured in the in the minimization of Eq. 10 and Eq. 13, we need to average
non-orthogonal AF), averaging the accelerometers readings the signals over a suitable time interval. This imposes a lower
in a temporal window inside each static interval. The cost bound in the length of the static interval (twait in Fig. 1).
function we use to estimate accelerometers’ parameters is: A initialization period (Tinit in Fig. 1) with no motion is
M
X essential as well: this will be exploited to characterize the
acc
L(θ )= (||g||2 − ||h(aSk , θ acc )||2 )2 (10) gyroscopes biases (Sec. IV-C) and the static detector operator
k=1 (Sec. IV-A).
where ||g|| is the actual magnitude of the local gravity A. Static Detector
vector that can easily recovered from specific public tables
(e.g., knowing latitude, longitude and altitude of the location
where we are performing the calibration). In order to
minimize Eq. 10, we employ the Levenberg-Marquardt
(LM) algorithm.

In order to calibrate the gyroscope triad, we can assume

the system as bias-free simply averaging the static gyroscope
signals over a suitable initial period of no motion. This
is justified by the following discussion about the Allan
variance (see Sec. IV-C). Moreover, since we need to use the
accelerometers as known references, we use the calibration
parameters θ acc computed above, correcting the accelerom- Fig. 3. An example of the static detector applied to the accelerometers
eters readings with Eq. 9. data: the static detector is represented by the black square wave, its high
We define the operator Ψ, that takes as input a sequence of n level classify an interval as static.
gyroscopes readings ω Si and an initial gravity versor ua,k−1
(i.e., a unit vector representing the gravity direction) given The accuracy of the calibration strongly depends on the
by the calibrated accelerometers, and return the final gravity reliability in the classification between static and motion
versor ug,k , computed using the gyroscopes measurements intervals: to calibrate the accelerometers we use static inter-
between the k − 1-th and the k-th static intervals: vals, while for gyroscopes calibration we also include the
motion intervals between two consecutive static intervals.
ug,k = Ψ ωiS , ua,k−1
 
(11) In our experience, band-pass filter based operators, like the
Ψ can be any integration algorithm that computes the final quasi-static detector used in [8], perform poorly with real
orientation through integrating the input angular velocities. datasets: detected static intervals frequently includes some
The unknown parameter vector we need to estimate to small portion of motion. Moreover, they require a fine tuning,
calibrate the gyroscope is: since they depend on three parameters.
We propose instead to use a variance based static detector
θ gyro = γyz , γzy , γxz , γzx , γxy , γyx , sgx , sgy , sgz
 
(12) operator, that exploits the lower bound in the lengths of the
In this case, we can define the cost function as: static interval introduced above. We base our detector on
the accelerometer signals: given a time interval of length
M
X twait seconds (see Fig. 1), for each accelerometers sample
L(θ gyro ) = ||ua,k − ug,k ||2 (13) (atx , aty , atz ) at time t, we compute the variance magnitude,
k=2
i.e. the magnitude of the variance, as:
where M is the number of static intervals, ua,k is the accel- q
eration versor measured averaging in a temporal window the ς(t) = [vartw (atx )]2 + [vartw (aty )]2 + [vartw (atz )]2 (14)
where vartw (at ) is an operator that compute the variance of C. Allan Variance
a general signal at in a time interval of length tw seconds
centered in t. We classify between static and motion intervals
simply checking if the square of ς(t) is lower or greater then
a threshold. As a threshold, we consider an integer multiple
of the square of the variance magnitude ςinit , computed over
all the initialization period Tinit . In all the experiments, we
use tw = 2 sec, while Tinit is estimated using the Allan
variance (see Sec. IV-C). It is important to note that our
static detector does not require any parameter tuning: the
integer multiplier used in the classification is automatically
estimated by our calibration algorithm (see Sec. IV-D). Fig. 3
reports an example of how our static filter works on real data:
in this case the estimated integer multiplier is 6.
Fig. 4. Allan Variance computed for the Xsens MTi gyroscopes triad.
B. Runge-Kutta Integration
As reported in Eq. 11, in the gyroscopes calibration we We characterize the random gyroscope bias drifts using
need to perform a discrete time angular velocity integration: the Allan variance ([20], [8]), which measures the variance
a robust and stable numerical integration method is desir- of the difference between consecutive interval averages. The
able since it can improve the calibration accuracy. Given Allan variance σa2 is defined as:
a common instruments rate of 100 Hz (like the Xsens
MTi IMU used in the experiments) and since we represent 1
rotations using quaternion arithmetic, with this setup a proper σa2 = (x(t̃, k) − x(t̃, k − 1))2 =
2
integration algorithm choice [19] is the Runge-Kutta 4th K
1 X
order normalized method (RK4n). In our experience (experi- = (x(t̃, k) − x(t̃, k − 1))2 (22)
ments not reported for space constraints) RK4n outperforms 2K
k=1
the standard linear integration procedure providing higher
where x(t, k) is the k-th interval average which spans t̃
accuracy results. Let Eq. 15 be the differential equation
seconds, and K is the number of interval which the total
describing the quaternion kinematics:
considered time is segmented in. We compute the Allan
1 variance for each gyroscope axis, with t0 ≤ t̃ ≤ tn . We fix t0
f(q, t) = q̇ = Ω(ω(t))q, (15)
2 = 1s, tn = 225s. The time in which the Allan variances of the
where Ω(ω) is the operator which turns the considered tri- three axis converge to a small value represents a good choice
dimensional angular velocity into the real skew symmetric for initialization period Tinit (Fig. 1). In this initialization
matrix representation, that is: period, we compute the average of the static gyroscope
 
0 −ωx −ωy −ωz signals to correctly determine the gyroscopes biases used in
ωx 0 ωz −ωy  the calibration. In the case of the Xsens MTi IMU, a good
Ω(ω) =  ωy −ωz
. (16) value for Tinit is 50 seconds (see Fig. 4).
0 ωx 
ωz ωy −ωx 0
D. Complete Procedure
The RK4n integration algorithm is:
To avoid unobservability in the calibration parameters
1 estimation, a minimum of nine different attitudes [15] has to
qk+1 = qk + ∆t (k1 + 2k2 + 2k3 + k4 ) , (17)
6 be collected (e.g., Fig. 1). In our experience, a higher number
ki = f(q(i) , tk + ci ∆t), (18) N of distinct attitudes are required to get better calibration
(i)
q = qk , f or i = 1, (19) results, while keeping reduced the duration of each static
i−1
X interval in order to preserve the assumption of temporal
q(i) = qk + ∆t aij kj , f or i > 1. (20) stability of the gyroscopes biases. With 36 ≤ N ≤ 50 and
j=1 1 sec ≤ twait ≤ 4 sec, we obtain a good trade-off between
where all the coefficients needed, ci and aij , are calibration accuracy, biases stability, and noise reduction.
The duration of the initialization period Tinit is given by
c1 = 0, c2 = 12 , c3 = 12 , c4 = 1,
the Allan variance analysis (see Sec. IV-C). The calibration
1
a21 = 2, a31 = 0, a41 = 0, protocol is summarized in Fig. 1, while in Algorithm 1 the
1 pseudo-code of the calibration algorithm is reported.
a32 = 2, a42 = 0, a43 = 1.
Finally, for each step, we also need to normalize the (k + 1)- V. E XPERIMENTS
th quaternion:
We have tested our method with both synthetic and
qk+1
qk+1 → . (21) real data. In the simulations, we can compare the results
||qk+1 || with a perfect ground truth, i.e., a noise-free, undistorted
Algorithm 1 IMU Calibration
Require: Tinit , twait ; aS and ω S (accelerometer’s and −1
gyroscope’s dataset collected according to Fig. 1). ω Ssynth = (Tg Kg ) ωO
synth − b
g
(24)
for the angular velocities. aSsynth and aO synth are the
bg ← average gyroscope signals over Tinit ; synthetic acceleration in the non-orthogonal sensor frame
ω Sbiasf ree ← ω S - bg ; and in the associated orthogonal frame, respectively. ω Ssynth
Minf ← empty matrix; and ω O synth are the synthetic angular velocities in the non-
ςinit ← Eq. 14, with tw = Tinit ; orthogonal sensor frame and in the associated orthogonal
for i = 1 : k frame, respectively. Eq. 23 and Eq. 24 are obtained from
2
threshold ← i ∗ ςinit ; models proposed in Eq. 6 and in Eq. 7. The metrics we use
s intervals ← static detector computed to evaluate the quality of the results are:
using twait and threshold;
1) Comparing the estimated values to the real ones;
[Residual, P aramsacc ] ← optimize Eq. 10 using
2) Comparing the average difference between the perfect,
s intervals and aS , averaging with twait ;
noise-free and undistorted signal with the noisy signal
Minf (i) ← [Residual, P aramsacc , threshold,
before and after calibration;
s intervals];
3) For the accelerometers only: we consider the magni-
end
tude and the angular error between the sensed acceler-
indexopt ← index of the minimum residual in Minf ;
ation and the applied one during each static intervals.
P aramsacc ← from Minf using indexopt ;
Since the magnitude of the gravity vector is assumed
s intervalsopt ← from Minf using indexopt ;
to be the only quantity known, the angular error
aO ← calibrate aS using P aramsacc ;
here is calculated for the worst case, where the full
P aramsgyro ← optimize Eq. 13 using s intervalsopt ,
error appears on a single accelerometer axis which is
ω Sbiasf ree and aO , averaging with twait .
perfectly horizontal, i.e. perpendicular to the gravity
acc
vector. An error of g · sin(θdiv ) will result in the pitch
acc
or roll angle being measured as θdiv radiants instead
signal, while the calibration matrices are known. With the of zero;
real dataset, we compare the estimated calibration param- 4) For the gyroscopes only: we consider the magnitude
eters with the calibration parameters reported in the IMU and the angular error between the acceleration sensed
datasheet. by the calibrated accelerometer and the acceleration
computed integrating the angular velocities given by
A. Synthetic Data the gyroscope.
We first generate a set of ideal, noise-free signals. The We have generated 1200 different signals, combining 40
accelerometers readings are generated starting from a three- different distortion parameter sets with a set of 30, randomly
dimensional signal based on three different-pulsation sinu- generated, ideal signals. Due to reduced space, we report
soids randomly modulated. At the beginning we add 5000 average and worst case results for one distortion parameter
zero samples (the initialization period) and every time the set, applied to the whole set of ideal signals. Results for the
three signals are simultaneously zero we introduce 400 zero other distortion parameter sets are very similar4 .
samples (the static intervals). The three dimensional gravity
vector projected onto the three axis has been added as well. TABLE I
For the angular velocities sensed by the gyroscope, the idea First metric: ACCELEROMETERS PARAMETERS
is to consider a tri-dimensional angular velocity vector ω, Real Mean RMS Mean Error RMS
which describes the perceived rotation of the aforementioned Value value x10−3 x10−3 x10−3
gravity vector, and then project ω onto the three axis of the αyz 0.0049 0.0049 0.0481 0.0398 0.0275
gyroscope. In this way we correlate the measurement of the αzy -0.0055 -0.0055 0.0401 0.0334 0.0214
αzx 0.0079 0.0079 0.0296 0.0248 0.0190
two different clusters of sensors. For each motion interval sax 0.9908 0.9908 0.0327 0.0265 0.0191
different zenith and azimuth velocities are randomly chosen, sa
y 1.0068 1.0068 0.0304 0.0258 0.0199
while for the rest of the time these velocities are considered sa
z 1.0066 1.0066 0.0215 0.0178 0.0151
a
bx 0.0793 0.0793 0.1369 0.1163 0.0819
to be equal to zero. The sampling frequency of the whole a
by -0.0024 -0.0024 0.2138 0.1760 0.1178
synthetic data has been fixed to 100 Hz. ba
z 0.0636 0.0636 0.1332 0.0953 0.0919
For each ideal signal, we add a white gaussian noise and
finally we distort the data with random generated distortion From Tab. I and Tab. II, we report a negligible average
parameters, i.e.: error, with an order of magnitude always equal or less
than 10−4 . From Tab. IV and Tab. V, we report significant
−1
aSsynth = (Ta Ka ) aO
synth − b
a
(23) 4 All synthetic results, the software for real calibration and a real
dataset are available following the link: robotics.dei.unipd.it/
for the accelerations, and: ˜pretto/software/tpm_icra2014.zip
 TABLE II
reductions in the absolute errors. Finally, from Tab. VI we
First metric: G YROSCOPES PARAMETERS
report that the divergence’s magnitude is reduced by a factor
Real Mean RMS Mean Error RMS of 11.9 and the angular error by a factor of 12.7, while
value value x10−3 x10−3 x10−3 in Tab. VII the magnitude of the considered divergence is
γyz 0.0112 0.0110 0.8547 0.6392 0.5920
γzy -0.0211 -0.0210 0.4419 0.3468 0.2669 reduced by a factor of 21.4 and the angular error by a factor
γxz 0.0040 0.0039 1.0630 0.9080 0.5266 of 21.9. In Fig. 5, we show some examples of the calibrated
γzx -0.0010 -0.0011 0.4102 0.3386 0.2302 signal compared with the ideal one.
γxy 0.0270 0.0270 0.8154 0.6375 0.4944
γyx 0.0151 0.0155 0.7250 0.7315 0.3958
sgx 0.8786 0.8785 0.4121 0.3366 0.2299
sgy 0.9703 0.9704 0.4059 0.3353 0.2237 B. Real Data
sgz 1.0460 1.0460 0.4216 0.3410 0.2397
As introduced, we have exploited an Xsens MTi IMU
TABLE III (Fig. 1) as experimental platform. The device datasheet
First metric: W ORST C ASE - PARAMETERS provides the factory calibrated misalignment matrices that
align the accelerometers (AF) and gyroscopes (GF) frames
Real value Est value Real value Est Value
to the body frame BF, while we estimate the matrices that
αyz 0.0049 0.0049 γyz 0.0112 0.0099
αzy -0.0055 -0.0055 γzy -0.0211 -0.0207 align that AF and GF to AOF. In order to compare our
αzx 0.0079 0.0079 γxz 0.0040 0.0030 results with the results of the factory calibration, we need to
sax 0.9908 0.9908 γzx -0.0010 -0.0011 know the matrix Rb that relates AOF to BF. Given Rb , we
sa
y 1.0068 1.0068 γxy 0.0270 0.0252
sa
z 1.0066 1.0066 γyx 0.0151 0.0166 can express our calibration vectors in BF. Given the factory
a
bx 0.0793 0.0792 sgx 0.8786 0.8790 calibration, Rb can be easily estimated using data contained
a
by -0.0024 -0.0026 sgy 0.9703 0.9701 in misalignment matrix (see discussion about Eq. 1) and
ba 0.0636 0.0636 sgz 1.0460 1.0460
z
following a simple method proposed in [14]5 .
TABLE IV We acquired a dataset as described in Fig. 1, with an initial
Second metric: A BSOLUTE ERRORS ALONG THE AXIS - static period of about 50 seconds, followed by a set of 37
ACCELEROMETERS rotations separated by static intervals of 2-4 seconds. We
uses as initial guess for the optimization the ideal values for
x-axis y-axis z-axis the accelerometer, that are (see Eq. 8):
m/s2 m/s2 m/s2
Uncalib 0.0842 0.0564 0.0635
Calib 0.0055 0.0056 0.0056  
1 0 0 0 1 0 0 0 1 . (25)
TABLE V
While for the gyroscope we used (see Eq. 12):
Second metric: A BSOLUTE ERRORS ALONG THE AXIS - G YROSCOPES
1 1 1
 (2n − 1)
x-axis y-axis z-axis
r 0 0 0 r 0 0 0 r , r= (26)
(rad/s) (rad/s) (rad/s) 2y
Uncalib 0.1043 0.1097 0.0345
Calib 0.0035 0.0039 0.0042 where n is the numbers of bit of the A/D converter, and the
gyroscope full-scale from datasheet is [-y, +y] rad/s.
TABLE VI
Third metric: ACCELEROMETERS DIVERGENCE ERROR The achieved calibration is comparable to the factory cali-
Average Max obser- Worst case Worst case
bration parameters given by the datasheet, see tables starting
error ved error average error max error from VIII. For each table, we report the factory calibration
m/s2 (rad) m/s2 (rad) m/s2 (rad) m/s2 (rad) parameters (on the left), and our calibration parameters (on
Uncalib 0.0665 0.2133 0.0623 0.2098 the right). It is important to point out that in our results we
( 0.0114) ( 0.0226) ( 0.0115) ( 0.0240)
Calib 0.0056 0.0299 0.0056 0.0298 are implicitly including an error that can’t be attributed to
( 0.0009) ( 0.0035) ( 0.0009) ( 0.0038) our calibration method. This is the propagated error caused
by the IMU’s datasheet rounded values (see table) when we
TABLE VII compute Rp . In spite of this problem, the results we obtain
Fourth metric: G YROSCOPES DIVERGENCE ERROR are very close to the parameters reported in the datasheet4 .
Average Max obser- Worst case Worst case
TABLE VIII
error ved error average error max error
m/s2 (rad) m/s2 (rad) m/s2 (rad) m/s2 (rad) S CALING - ACCELEROMETER
Uncalib 4.7125 9.2930 5.2859 8.5822
( 0.5494) ( 0.5494) ( 0.6029) ( 0.6029) 415 0.00 0.00 414.41 0 0
Calib 0.2208 0.4469 0.5102 0.8597 0.00 413 0.00 0 412.05 0
( 0.0256) ( 0.0256) ( 0.0569) ( 0.0569) 0.00 0.00 415 0 0 414.61

5 Due to reduced space, we can’t report here other details.

 (a) (b)

Fig. 5. Calibration Improvement: (a) accelerometers, (b) gyroscopes (the small box is zoomed in the larger box).

TABLE IX TABLE XI
S CALING - G YROSCOPE M ISALIGNMENT - G YROSCOPE

4778 0.00 0.00 4778.0 0 0 1.00 -0.01 -0.02 0.9998 -0.0149 -0.0218
0.00 4758 0.00 0 4764.8 0 0.00 1.00 0.04 0.0003 1.0007 0.0433
0.00 0.00 4766 0 0 4772.6 -0.01 0.01 1.00 -0.0048 0.0121 1.0004

TABLE X TABLE XII

M ISALIGNMENT - ACCELEROMETER O FFSET - ACCELEROMETER

1.00 0.00 -0.01 1.0000 -0.0066 -0.0110 33123 33276 32360 32768 32466 32485
0.01 1.00 0.01 0.0102 1.0001 0.0114
0.02 0.01 1.00 0.0201 0.0098 0.9998

[7] J. Lotters, J. Schipper, P. Veltink, W. Olthuis, and P. Bergveld,

“Procedure for in-use calibration of triaxial accelerometers in medical
VI. C ONCLUSIONS AND F UTURE W ORKS applications,” Sensors and Actuators A: Physical, vol. 68, no. 1-3, pp.
In this paper, we have presented an effective and semi 221–228, 1998.
[8] W. Fong, S. Ong, and A. Nee, “Methods for in-field user calibration of
automatic IMU calibration method. The calibration protocol an inertial measurement unit without external equipment,” Measure-
is quick, reliable and it does not need any external equipment. ment Science and Technology, vol. 19, pp. 1–11, 2008.
The proposed algorithms are easy to implement and they [9] R. Rogers, Applied Mathematics in Integrated Navigation Systems,
ser. AIAA education series. American Institute of Aeronautics and
do not require any parameters tuning. Results achieved with Astronautics, 2003.
both synthetic and real data show the real possibility to [10] A. B. Chatfield, Fundamentals of high accuracy inertial navigation,
improve at no cost the performances of low cost, poor ser. Progress in astronautics and aeronautics. Reston, VA. American
Institute of Aeronautics and Astronautics, Inc., 1997.
calibrated, IMUs, such the ones that equip the current smart- [11] J. J. Hall and R. L. Williams II, “Inertial measurement unit calibration
phones. Actually, we are currently testing our approach on platform,” Journal of Robotic Systems, vol. 17, no. 11, pp. 623–632,
smartphone-grade IMU sensors. We also plan to integrate the 1998.
[12] A. Kim and M. Golnaraghi, “Initial calibration of an inertial mea-
proposed procedure in a calibration framework for camera- surement unit using an optical position tracking system,” in Position
IMU systems, that includes IMU, camera-IMU frame and Location and Navigation Symposium, 2004. PLANS 2004, 2004, pp.
temporal calibrations. 96–101.
[13] E. M. Nebot and H. F. Durrant-Whyte, “Initial calibration and align-
R EFERENCES ment of low-cost inertial navigation units for land vehicle applica-
tions.” Journal of Robotic Systems, vol. 16, no. 2, pp. 81–92, 1999.
[1] V. Kubelka and M. Reinstein, “Complementary filtering approach to [14] I. Skog and P. Hndel, “Calibration of a mems inertial measurement
orientation estimation using inertial sensors only.” in Proc. of. IEEE unit,” in Proc. of XVII IMEKO WORLD CONGRESS, 2006.
International Conference on Robotics and Automation (ICRA), 2012, [15] Z. F. Syed, P. Aggarwal, C. Goodall, X. Niu, and N. El-Sheimy, “A
pp. 599–605. new multi-position calibration method for mems inertial navigation
[2] T. Hamel and R. E. Mahony, “Attitude estimation on so[3] based systems,” Measurement Science and Technology, vol. 18, pp. 1897–
on direct inertial measurements.” in Proc. of. IEEE International 1907, 2007.
Conference on Robotics and Automation (ICRA), 2006, pp. 2170– [16] C. M. Cheuk, T. K. Lau, K. W. Lin, and Y. Liu, “Automatic calibration
2175. for inertial measurement unit,” in in Proc. of International Conference
[3] K. Tsotsos, A. Pretto, and S. Soatto, “Visual-inertial ego-motion esti- on Control, Automation, Robotics and Vision, 2012.
mation for humanoid platforms,” in Proc. of: IEEE-RAS International [17] M. Hwangbo, J.-S. Kim, and T. Kanade, “Imu self-calibration using
Conference on Humanoid Robots, 2012, pp. 704–711. factorization,” IEEE Transactions on Robotics, vol. 29, no. 2, pp. 493–
[4] K. Konolige and M. Agrawal, “Frameslam: From bundle adjustment 507, 2013.
to real-time visual mapping,” IEEE Transactions on Robotics, vol. 24, [18] C. Jekeli, Inertial Navigation Systems with Geodetic Applications. De
no. 5, pp. 1066–1077, Oct. 2008. Gruyter, 2001.
[5] M. Li, B. Kim, and A. I. Mourikis, “Real-time motion estimation on [19] M. Andrle and J. Crassidis, “Geometric integration of quaternions,” in
a cellphone using inertial sensing and a rolling-shutter camera,” in in Proc. of AIAA/AAS Astrodynamics Specialist Conference, 2012.
Proceedings of the IEEE International Conference on Robotics and [20] A. M. Sabatini, “A wavelet-based bootstrap method applied to inertial
Automation, Karlsruhe, Germany, May 2013, pp. 4697–4704. sensor stochastic error modelling using the allan variance,” Measure-
[6] XSens, http://www.xsens.com. ment Science and Technology, vol. 17, no. 11, pp. 2980–2988, 2006.

View publication stats