This is the ac cepted version of the   manuscript :   D. Nemec, A. Janota, M. Hruboš and V. Šimák, "Intelligent Real - Time MEMS Sensor Fusion and Calibration,"  in IEEE Sensors Journal, vol. 16, no.   19, pp. 7150 - 7160, Oct.1, 2016,   doi: 10.1109/JSEN.2016.2597292 ,   URL:   http://ieeexplore.ieee.org/document/7529156    Abstract — This   paper   discusses   an   innovative   adaptive  heterogeneous fusion algorithm based on estimation of the mean  square error of all va riables used in real time processing. The  algorithm   is   designed   for   a   fusion   between   derivative   and  absolute sensors and is explained by the fusion of the 3 - axial  gyroscope, 3 - axial accelerometer and 3 - axial magnetometer into  attitude and heading estimatio n. Our algorithm has similar error  performance   in   the   steady   state   but   much   faster   dynamic  response   compared   to   the   fixed - gain   fusion   algorithm.   In  comparison   with   the   extended   Kalman   filter   the   proposed  algorit hm   converges   faster   and   takes less computatio nal time . On  the other hand, Kalman filter has smaller mean square output  error in a steady state but becomes unstable if the estimated state  changes   too   rapidly.   Additionally, the noisy fusion deviation can  be used in the process of calibration. The paper   proposes and  explains   a   real - time   calibration   method   based   on   machine  learning working in the online mode during run - time. This allows  compensation   of   sensor   thermal   drift   right   in   the   sensor’s  working   environment   without   need   of   re - calibration   in   the  lab oratory .  Index Terms — calibration, inertial navigation,   m ean square  error methods , sensor fusion.  I.   I NTRODUCTION  NERTIAL   SENSORS   manufactured   by   the   MEMS   (Micro  Electro - Mechanical Systems) technology are the core of  modern low - cost AHRSs (Attitude and He ading Reference  Systems).   The   purpose   of   these   systems   is   to   determine  rotation of the measured object with respect to the horizontal  plane and northern direction   which is a crucial task in mobile  robotics,   aviation, automated car navigation and many other s.  These   sensor   systems   use   nonlinear   discrete   numerical  integration of the measured angular velocity   which   is a typical  © 20 16   IEEE. Personal use of this material is permitted. Permission from  IEEE must be obtained for all other uses, in any current or future media,  including reprinting/republishing this material for advertising or promotional  purpo ses, creating new collective works, for resale or redistribution to servers  or lists, or reuse of any copyrighted compone nt of this work in other works.  The paper was submitted for review in February   22 , 2016.   This work was  supported by the European Regional Development Fund and the Ministry of  Education of the Slovak Republic, within the project ITMS 26220220089  “New   methods   of   mea surement   of   physical   dynamic   parameter s   and  interactions of motor vehicles, traffic flow and road”.  Authors are   with Dep t.   of Control and Information Systems, Faculty of  Electrical   Engineering,   University   of   Žilina,   Univerzitná   8215/1,   Žilina  01026, Slovakia.   e - mail:  { dusan.nemec ;   ales.janota;   marian.hrubos; vojtech.simak} @fel.uniza.sk  example   of   velocity   measurement   when   the   sensor   is  measuring   time   derivative   of   desired   variable.   Main  disadvantage of this method is g reat sensitivity of the output  quality to the precision of the sensor measurements (especially  sensor bias will cause increasing drift of the integrated result).  First step of elimination of the integrated error is calibration of  the sensor. A standard way   of calibration measures raw output  of the sensor as a response to stimulus with known amplitude.  Relation between raw and real sensor outputs is formed into  the transfer function (calibration curve) and its parameters are  obtained from the measurements du ring calibration in offline  mode. It is possible to calibrate by:  -- One point (zero order transfer function   -   bias only) .  -- Two points (first order transfer function   -   bias and  gain),  -- Multiple points (calibration curve is a polyline or higher  order cu rve) .  In order to eliminate influence of the sensor random noise  each calibration point has to be computed as an average of  multiple   measurements   in   the   same   conditions   [1] .   This  requires special laboratory equipment which provides accurate  and steady simu lation of different sensor stimuli. Zhang et al.  proposed a method of estimation of the calibration constants  for   the   3 - axial inertial sensor (gyroscope, accelerometer)   [2] .  Gyroscope bias is determined directly in   a   steady state and  accelerometer bias i s computed after multiple steps when the  acceleration sensor is oriented vertically along each of its axes  one by one or the sensor has to be exposed to precisely known  stimuli   [3][4] . All these methods are working in the offline  mode. Wang and Hao propose d a method utilizing an artificial  neural network combined with the Kalman filter for estimation  of nonlinear calibration parameters   [ 5 ] . For online calibration  it is necessary to detect steady state of the object ,   e.g. by lower  vibrations   [ 6 ] .  When MEMS s ensors are used, their calibration parameters  tend to drift with temperature   [ 5 ]   [ 7 ] . Transfer function is  therefore two - dimensional   –   one input   corresponds to   raw  sensor data and the second input is sensor temperature. Most  of   commercially   available   integ rated   MEMS   sensors  incorporate   a   temperature   sensor   which   allows   usage   of  advanced temperature compensation techniques .  In order to compensate integrated error during run - time it is  necessary to use a secondary absolute sensor and provide a  data fusion. Th e secondary sensor   may   be much noisier and  have slower response but its error has to be kept inside fixed  bounds. A modification of the Kalman filter can be used as a  Intelligent real - time MEMS sensor fusion and  calibration  Dušan Nemec ,   Aleš   Janota ,   Marián   Hruboš ,   and   Vojtech Šimák  I
This is the ac cepted version of the   manuscript :   D. Nemec, A. Janota, M. Hruboš and V. Šimák, "Intelligent Real - Time MEMS Sensor Fusion and Calibration,"  in IEEE Sensors Journal, vol. 16, no.   19, pp. 7150 - 7160, Oct.1, 2016,   doi: 10.1109/JSEN.2016.2597292 ,   URL:   http://ieeexplore.ieee.org/document/7529156  core of the sensor fusion algorithm   [3][ 8 ][ 9 ] [23] , however it  might   be   difficult   to   estim ate   parameters   of   the   filter  (covariance   matrix,   state   model)   for   a   standalone   sensor  system because the Kalman filter parameters depend on the  measured system. Another sensor fusion utilizes Bayesian  networks and   the   stochastic approach   [ 1 0 ][ 1 1 ][ 1 2 ] . We h ave  proposed   a   heterogeneous   sensor   fusion   method   for   one  differential sensor   and one absolute sensor which requires  only minimum count of parameters independently from the  measured system.   The performance of our algorithm will be  compared with the perform ance of the extended Kalman filter  (EKF) used in direct form described in [23].  Our   real - time   calibration   method   utilizes error   estimate  obtained as a side output from the sensor fusion algorithm.  This approach eliminates the need of steady state detectio n  and offline calibration. Since it can be run ning   all the time  when the sensor is   in use   our method should compensate long -  term drifts continuously.  II.   H ETEROGENEOUS   F USION   A LGORITHM   C ONSIDERING  Q UALITY  The method will be explained on the example of the fu sion  of   the   3 - axial   gyroscope   (velocity   sensor),   3 - axial  accelerometer   (absolute   attitude   sensor)   and   3 - axial  magnetometer   (absolute   heading   sensor).   Sensor   axes   are  orientated   according   to   the   NED   convention   ( x - North   or  forward,   y - East or right,   z - Down),   Euler angles are computed  in the ZYX convention ( α   –   Roll,   β -   Pitch,   γ -   Yaw ).   Attitude  of object is   then expressed by roll and pitch angles; heading is  expressed by yaw angle.  In order to express the quality of estimation we will use   the  mean square error (MSE). In general the error model of the  attitude   estimation is nonlinear   [ 1 3 ] . MSE of the directly  measured data is considered constant and depends on the used  sensor; MSE of a computed variable   y   =   f ( x 1 ,   x 2 , ..., x N )   is  approximated by: ). ( MSE ) ( MSE   1 2  k N k   k  x x f y                   (1)  A.   Estimating   Euler angles from gyroscope rea dings  In the   inertial   navigation the   object’s   Euler   angles   are  primary computed from   the   angular velocity   ω   measured by  the gyroscope. There are several methods of angular velocity  integration into Euler angles; we usually use the matrix - based  algorithm. Rotation is expressed as the 3D transformation  matrix   R   updated by the infinitesimal update matrix computed  from each sample of   the   angular velocity   [ 14 ] :                        1 1 1 where ,   update update u  t t t t t t  x y x z y z         R R R R  (2)  and   ω x ,   ω y ,   ω z   are the Cartesian components of the angular  velocity vector and   Δ t   is a sampling period of the gyroscope.  The error of the gyroscope can be considered the same for  each axis and constant: ) ( MSE ) ( MSE   gyro   const E t   i i          (3)  MSEs of the updated uncompensated rotational matrix   R u  are: ), ( MSE ) ( MSE ) ( MSE ) ( MSE  , 2 2 , 3 2 gyro 2 , 3 gyro 2 , 2 , 1 , 1 u  k z k y k k k k  R R E R E R R R             (4.1) ), ( MSE ) ( MSE ) ( MSE ) ( MSE  , 3 2 , 1 2 gyro 2 , 1 gyro 2 , 3 , 2 , 2 u  k x k z k k k k  R R E R E R R R             (4.2) ), ( MSE ) ( MSE ) ( MSE ) ( MSE  , 1 2 , 2 2 gyro 2 , 2 gyro 2 , 1 , 3 , 3 u  k y k x k k k k  R R E R E R R R             (4.3)  where index   k   = 1  3 . As can be seen, errors of all matrix  elements are increasing with new samples.   Uncompensated  Euler angles are then   [15 ] :     3 , 3 u 3 , 2 u gyro   , atan2   R R    ,   (5.1)     3 , 1 u gyro   arcsin   R     ,   (5.2)     1 , 1 u 2 , 1 u gyro   , atan2   R R    .   (5.3)  Corresponding MSEs of the Euler angles are:     2 2 u3,3 2 u2,3 3 , 3 u 2 3 , 2 u 3 , 2 u 2 u3,3 gyro  ) ( MSE ) ( MSE ) ( MSE  R R R R R R       ,   (6.1) 2 u1,3 3 , 1 u gyro   1 ) ( MSE ) ( MSE   R R      ,   (6.2)     2 2 u1,2 2 u1,1 1 , 1 u 2 2 , 1 u 2 , 1 u 2 u1,1 gyro  ) ( MSE ) ( MSE ) ( MSE  R R R R R R       .   (6.3)  These formulas are undefined at   the   gimbal lock ( cos     = 0  and   R u 1,3 =    1 ). If such condition occurs it is impossible to  determine both roll and yaw (one has to be chosen) and MSE  estimation is very imprec ise.  B.   Estimating roll and pitch from accelerometer readings  Secondary ,   the attitude of the object can be obtained from  acceleration readings by formulas   [ 14 ] : ) , ( atan2 acc   z y   a a        ,   (7) ) , ( atan2   2 2 acc   z y x   a a a       ,   (8)  where   a x ,   a y ,   a z   are   the   components   of   the   acceleration  measured by   the   accelerometer bound with a moving object.  MSE of attitude estimation depends on the dynamics of the  system (the vector   a   measured by the accelerometer is a sum  of   the   gravitational   acceleration   g   and   the   object’s   own  accelerat ion including vibrations which might be useful in  different applications   [ 1 6 ] ). If the object is steady, variance of  the vector   a   is smaller and attitude estimation is more precise:     2 2 2 2 2 acc  ) ( MSE ) ( MSE ) MSE(  z y y z z y  a a a a a a       ,   (9)                   2 2 2 2 2 2 2 2 2 2 2 2 acc  ) ( MSE ) ( MSE ) ( MSE ) MSE(  z y x z y z z y y x x z y  a a a a a a a a a a a a a            (10)  In order to decrea se error of estimation the acceleration can  be   averaged   from   multiple   samples   (oversampled
This is the ac cepted version of the   manuscript :   D. Nemec, A. Janota, M. Hruboš and V. Šimák, "Intelligent Real - Time MEMS Sensor Fusion and Calibration,"  in IEEE Sensors Journal, vol. 16, no.   19, pp. 7150 - 7160, Oct.1, 2016,   doi: 10.1109/JSEN.2016.2597292 ,   URL:   http://ieeexplore.ieee.org/document/7529156  measurement).   Then   MSEs   of   average   acceleration  components are:      . ) ( MSE ] [ 1 ] [ 1 ) ( MSE ] [ 1 ) ( MSE  2 1 1 2 1 2  N a k a N k a N N a a k a N a  i N k   i N k   i i N k   i i i                       (11)  The second expression allows processing of accelerometer  samples   by   batch es   with   size   N   with   lower   memory  requirements.   Accelerometer’s   own   errors   MSE( a i )   are  negligible at higher   N   with respect to the errors caused by  vibrations. However, the simulations have shown that using  batches is causing relatively large step changes in   the resultant  estimated   variable.   Therefore   it   is   better   to   use   online  approximation of average and MSE; then the fusion can be  performed in each step   [ 1 7 ] :      N n a n s N n s   i i i  2  ] [ ] 1 [ 1 ] [        ,   (12.1)      N n a n a N n a   i i i  ] [ ] 1 [ 1 ] [        ,   (12.2)  where ] [ n s i   is a mean   square of the   i - th acceleration component  in the   n - th step. Formulas (12) are first order low - pass IIR  filters with cut - off frequency:                    1 2 1 1 arccos 2  sample cut   N N f f    .   (13)  MSE of the estimated average ] [ n a i   is then: N n a n a n s n a   i i i i  ]) [ ( MSE ] [ ] [ ]) [ ( MSE   2       .   (14)  I n order to compute   yaw   from the magnetic induction vector  B , it has to be rotated from objects’ local coordinates to the  global horizontal plane by following:                                       z y x z y x  B B B B B B               cos cos cos sin sin sin cos 0 sin cos sin sin cos  .   (15)  Yaw   is t hen computed by the formula: ) , ( atan2 mag   x y   B B         .   (16)  The vertical component of the magnetic induction z B    is not  used; therefore the third row of the matrix in (15) can be  omitted in the algorithm. Since the transformation (15) is  using estimated roll and pitch, quality of the   yaw   estimati on  depends on the quality of the attitude estimation.                          ) ( MSE cos cos cos sin sin ) ( MSE sin sin sin cos ) ( MSE sin cos ) ( MSE sin sin ) ( MSE cos ) ( MSE  2 2 2 2 2                   z y x z y z y x x  B B B B B B B B B               (17.1)                  ). ( MSE sin cos ) ( MSE sin ) ( MSE cos ) ( MSE  2 2 2        y z z y y  B B B B B         (17.2)  Then MSE of the yaw estimated by   the   magnetometer is:     2 2 2 2 2 mag  ) ( MSE ) ( MSE ) ( MSE  x y y x x y  B B B B B B               .   (18)  Error of the magnetometer is also the same for each axis  an d can be considered constant: mag ) ( MSE   E B i     .   (19)  C.   Heterogeneous fusion algorithm  The scheme   shown   in Fig. 1describes the algorithm that can  be   used   for   fusion   of   the   Euler   angles   computed   from  gyroscope, accelerometer and magnetometer readings.  Th e fusion algorithm is based on incremental compensation  of difference between incrementally integrated Euler angles   gyro ,   gyro ,   gyro   and absolute but noisy Euler angles   acc ,   acc ,   mag . The adaptive gain block (see Fig. 1) optimizes the impact  of eac h data source in order to increase resultant precision. A  value with lower MSE has a higher effect to the result   [1 8 ] . If  the gain block has unit gain, the values obtained by gyroscope  integration are not taken into account and result is equal to  Euler ang les obtained from the accelerometer and gyroscope.  Resultant compensated Euler angles are equal to:                                                                  gyro mag gyro acc gyro acc gyro gyro gyro  1 1 1                         K K K K K K d d d  ,   (20)  where   K   ,   K   ,   K     are variable fusion gain coefficients adjusted  according to the estimation erro rs and their values vary from 0  to 1. The vector   [ d   ,   d   ,   d   ]   represents fusion deviations:                                      gyro mag gyro acc gyro acc                K K K d d d  (21)  In case when fusion deviation   d i   is not available due to  different sampling frequencies, it is simply considered zero.  Usually   the   magne tometer   has   much   lower   sampling  frequency than gyroscope or accelerometer. When   the   new  sample   from   magnetometer   is   not   available,   the   fusion  Fig. 1.   The fusion algorithm scheme.
This is the ac cepted version of the   manuscript :   D. Nemec, A. Janota, M. Hruboš and V. Šimák, "Intelligent Real - Time MEMS Sensor Fusion and Calibration,"  in IEEE Sensors Journal, vol. 16, no.   19, pp. 7150 - 7160, Oct.1, 2016,   doi: 10.1109/JSEN.2016.2597292 ,   URL:   http://ieeexplore.ieee.org/document/7529156  deviation   d γ   is null.  Fusion   deviations are important inputs for the automatic  calibration algorithm described in   the next chapter. Their  errors are:     ) ( MSE ) ( MSE ) ( MSE   gyro acc 2             K d  ,   (22.1)     ) ( MSE ) ( MSE ) ( MSE   gyro acc 2             K d  ,   (22.2)     ) ( MSE ) ( MSE ) ( MSE   gyro mag 2             K d  .   (22.3)  The fusion gain coefficients have to be adjusted in order to  minimize the output error which is equal to: ) ( MSE ) 1 ( ) ( MSE ) ( MSE   gyro 2 acc 2                 K K  .   (23)  Note that this formula is also valid for pitch and   yaw   angles  with   corresponding   coefficients.   Output   error   is   minimal  when: 0 ) ( MSE 0 ) ( MSE  2 2                K K  .   (24)  Solving these conditions we obtain the optimal gain: ) ( MSE ) ( MSE ) ( MSE  acc gyro gyro           K  .   (25)  D.   Proof of Pr ecision Improvement  We can substitute the optimal gain in the formula (23): ) ( MSE ) ( MSE ) ( MSE ) ( MSE ) ( MSE  acc gyro acc gyro             .   (26)  By dividing by   MSE(   gyro )   we   can   get: 1 ) ( MSE ) ( MSE ) ( MSE ) ( MSE ) ( MSE  acc gyro acc gyro           .   (27)  Since MSE is always higher than zero   MSE(   )   is always  smaller   than   MSE(   gyro ) .   Sin ce   the   formula   (26)   is  symmetrical, the same is valid for   MSE(   acc ) . Therefore the  theoretical   fusion output error is always smaller than errors of  any single estimation.   Note that   the   output error directly  depends on estimation of the MSEs during   executio n of the  algorithm.  E.   Error of the   R esultant   R otational   M atrix  Fig. 1 contains the conversion block from compensated  Euler angles to the rotational matrix   [ 15 ] :                                                c c c s s s c s s c s c c s c c s s s s c c s s s s c c c  R  ,   (28)  where   s     = sin   ,   c     = cos   , etc.   Note that   the   matrix   R   =   R u  when   no   magnetometer   and   accelerometer   samples   are  available due to the different sampling frequency.  When   a   new   sample   is   processed   and   errors   of   the  compensated Euler angles are computed according to (26), it  is possible to compute MSE of each element in the   rotational  matrix which is used in next sample processing (used in (4)):             ) ( MSE sin cos ) ( MSE cos sin ) ( MSE   2 2 1 , 1             R  ,   (29.1)             ) ( MSE cos cos ) ( MSE sin sin ) ( MSE   2 2 2 , 1            R  ,   (29.2)       ) ( MSE cos ) ( MSE   2 3 , 1     R  ,   (29.3)               ), ( MSE sin sin sin cos cos ) ( MSE cos cos sin ) ( MSE cos sin cos sin sin ) ( MSE  2 2 2 1 , 2                         R  (29.4)               ), ( MSE cos sin sin sin cos ) ( MSE sin cos sin ) ( MSE sin sin cos cos sin ) ( MSE  2 2 2 2 , 2                         R  (29.5)             ) ( MSE sin sin ) ( MSE cos cos ) ( MSE   2 2 3 , 2            R  ,   (29.6)               ), ( MSE sin sin cos cos sin ) ( MSE cos cos cos ) ( MSE cos sin sin sin cos ) ( MSE  2 2 2 1 , 3                         R  (29.7)               ), ( MSE cos sin cos sin sin ) ( MSE sin cos cos ) ( MSE sin sin sin cos cos ) ( MSE  2 2 2 2 , 3                         R  (29.8)             ). ( MSE sin cos ) ( MSE cos sin ) ( MSE   2 2 3 , 3            R  (29 .9)  By comparing expressions in brackets with elements of the  matrix   R   in (28) it is possible to simplify   the   formulas (29):       ), ( MSE ) ( MSE cos sin ) ( MSE   2 2 , 1 2 1 , 1        R R       (30.1)       ), ( MSE ) ( MSE sin sin ) ( MSE   2 1 , 1 2 2 , 1        R R      (30.2)       ), ( MSE cos ) ( MSE   2 3 , 1     R  (30.3)      ), ( MSE ) ( MSE cos cos sin ) ( MSE ) ( MSE  2 2 , 2 2 2 1 , 3 1 , 2         R R R       (30.4)      ), ( MSE ) ( MSE sin cos sin ) ( MSE ) ( MSE  2 1 , 2 2 2 2 , 3 2 , 2         R R R       (30.5)       ), ( MSE sin sin ) ( MSE ) ( MSE   2 2 3 , 3 3 , 2            R R  (30.6)      ), ( MSE ) ( MSE cos cos cos ) ( MSE ) ( MSE  2 2 , 3 2 2 1 , 2 1 , 3         R R R       (30.7)      ) ( MSE ) ( MSE sin cos cos ) ( MSE ) ( MSE  2 1 , 3 2 2 2 , 2 2 , 3         R R R       (30.8)       ) ( MSE sin cos ) ( MSE ) ( MSE   2 2 3 , 2 3 , 3          R R  .   (30.9)  Initial MSE of Euler angles s hould be initialized to a large  value representing low quality, e.g.:     2 0 2 0 0   2 ) ( ) ( ) (              MSE MSE MSE  .   (31)  In order to avoid extreme values of MSE (e.g. around  gimbal lock singularity, see (6)), it is convenient to limit MSE  of the rotational matrix to the interva l ( - 1, 1). Maximal values  for Euler angles’ MSEs can be set according to (31).  III.   T HE   R EAL - T IME   C ALIBRATION   A LGORITHM  In this chapter we will discuss how fusion deviation ( a   side  product of the fusion of two or more sensors) can be used for  sensor calibration . It is convenient when error estimate of  deviation is also available. The algorithm will be explained on  the   example   of   the   gyroscope,   accelerometer   and  magnetometer fusion proposed in the previous chapter.
This is the ac cepted version of the   manuscript :   D. Nemec, A. Janota, M. Hruboš and V. Šimák, "Intelligent Real - Time MEMS Sensor Fusion and Calibration,"  in IEEE Sensors Journal, vol. 16, no.   19, pp. 7150 - 7160, Oct.1, 2016,   doi: 10.1109/JSEN.2016.2597292 ,   URL:   http://ieeexplore.ieee.org/document/7529156  A.   Calibration mode l  A universal relation between the raw value   w   measured by  the MEMS sensor and the actual value   ω   of the measured  variable is following: ) , (   T w C    ,   (32)  w here   T   is sensor’s temperature and   C ( w , T )   is a calibration  transfer function   [ 19 ] . Note   that errors caused by the sensor’s  hysteresis are neglected   [ 20 ] .  Our   algorithm   assumes   first   order   calibration   transfer  function at any given temperature, therefore there are two one -  dimensional functions: gain   G ( T )   and bias   B ( T )   for which it is  valid: ) ( ) (   T B w T G        (33)  Since the   C ( w , T )   function is continuous, partial functions  G ( T )   and   B ( T )   are also continuous and they can be expressed  by polynomials:        N k k k T g T G   0  ) (  ,   (34)        N k k k T b T B   0  ) (  .   (35)  where   g k   and   b k   are constants which   change only by the long -  term drift (e.g. by aging). In order to suppress sensitivity of  the higher order coefficients, full scale of the sensor raw data  and temperature should be normalized to the interval (0,1) or  ( - 1,1). The proposed learning algorithm r equires floating - point  number implementation   [ 21 ] . In order to maintain precision in  full   range   and   avoid   ignoring   small   incremental   steps  (especially by 32 - bit floating point IEEE754 format), the gain  function should be shifted by one. Bias es   of the MEMS  sensors drift faster than gain ,   therefore it is convenient to  make the bias function independent from the gain function.  The modified relation is:         ) ( 1 ) (   T B w T G         .   (36)  If the calibration formula (36) is used, the polynomial  coefficients   g k   and   b k   c an be initialized to zero.  In inertial navigation it is required to measure vectors of  linear acceleration and angular velocity in all three axes (the  3 - axial MEMS vibrational gyroscope and accelerometer are  usually used for this purpose   [ 4 ] ).   Each axis h as its own gain  and the offset calibration functions (independently from each  other) according to (33). However, if the vector variable is  measured,   it   is   necessary   to   take   sensor   orientation   into  account. The most general case occurs when each axis is  mea sured by a single physical sensor. The relation between  the vector of measured raw data   w   and the calibrated vector   ω  is following (in the orthonormal coordinate system):     ) ( ) (   T T   B w G A ω        ,   (37)  where   G ( T )   is   a   diagonal   matrix   of   the   shifted   gain  functions and   B ( T )   is a vector of the bias functions:                1 ) ( 0 0 0 1 ) ( 0 0 0 1 ) ( ) (  T G T G T G T  z y x  G  ,   (38.1             ) ( ) ( ) ( ) (  T B T B T B T  z y x  B  .   (38.2)  Each   of   G i ( T )   and   B i ( T )   functions   are   polynomials  according   to   (34)   and   (35)   respectively.   For   example,  c oefficients of the polynomial   G x ( T )   are marked   g x , 0   g x , 1   ... ,  coefficients of the polynomial   B x ( T )   are marked   b x , 0   b x , 1   ... etc.  The   matrix   A   is   a   constan t   alignment   matrix   which  considers misalignment between sensor’s axes and the object’s  axes   [22] . Since the gain matrix   G ( T )   normalizes each axis  separately, the alignment matrix has to be column - normalized:                       2 3 , 2 2 3 , 1 2 , 3 1 , 3 3 , 2 2 2 , 3 2 2 , 1 1 , 2 3 , 1 2 , 1 2 1 , 3 2 1 , 2  1 1 1  a a a a a a a a a a a a  A  .   (39)  B.   Backpropagation of   the fusion error  Parameters   a i,k   g i,k   b i,k   have to be obtained by learning.   Since  the true value of the measured variable is unavailable du ring  runtime and   the   data processing algorithm is nonlinear and  recursive,   we cannot minimize the error between raw r eadings  of the   sensor and the actual value as the static (offline)  calibration does.   However ,   it is possible to minimize the  fusion deviation (see (21)) in long term.   Fusion deviation is  therefore used as an output error for learning.   Implementation  of   the   deterministic   least - squares   method   for   non - linear  recursive   system   would   require   remembering   all   previous  samples and would be very complicated. Such approach would  be difficult to implement into low - cost hardware.   Considering  mentioned drawbacks,   we have   chosen to use a stochastic  learning algorithm .  In order to decrease memory requirements and decrease the  computational time for one iteration of the learning algorithm,  we have decided to avoid recurrent   learning.   To be able to do  that it is necessary to   back - propagate the f usion deviation  vector   [ d   ,   d   ,   d   ] T   through the data processing algorithm   to  the angular velocity error vector   e   =   [ e x ,   e y ,   e z ] T . The data  processing algorithm   disregarding sensor fusion   can   also   be  modelled   by   the   following   continuou s non - linear differential  equation   [15] :  Fig.   2 .   The   real - time calibration scheme in the context of the sensor data  processing and fusion.
This is the ac cepted version of the   manuscript :   D. Nemec, A. Janota, M. Hruboš and V. Šimák, "Intelligent Real - Time MEMS Sensor Fusion and Calibration,"  in IEEE Sensors Journal, vol. 16, no.   19, pp. 7150 - 7160, Oct.1, 2016,   doi: 10.1109/JSEN.2016.2597292 ,   URL:   http://ieeexplore.ieee.org/document/7529156                                                  z y x  dt d dt d dt d                     cos cos cos sin 0 sin cos 0 cos sin cos cos sin sin 1  ,   (4 0 )  The 3x3 matrix in (40) represents sensitivity of the output  (Euler   angles)   to   the   calibrated   gyroscope   readings.   The  i nverted matrix can be used to transform errors of the Euler  angles   ( fusion   deviations)   to   the   errors   of   the   calibrated  angular velocity :                                                 d d d t e e e  z y x  cos cos sin 0 cos sin cos 0 sin 0 1 1  cal  ,   (4 1 )  w here   Δ t cal   is a period of the calibration procedure (can be  larger   than   the   sampling   period   of the gyroscope). The relation  (41) is valid only if the Euler   angles   α, β, γ   are estimated  precisely, which is achieved by compensation   of the errors by  sensor fusion. In order to avoid invalid adjustments of the  calibration   parameters   cau sed by initial low quality of the  estimated Euler angles, the calibration proce dure is applied  only   if   the   fusion   deviation   is   below   some   predefined  threshold.   This approach also eliminates the need for recurrent  learnin g and the learning algorithm optimizes the error vector  of the gyroscope (instead of direct optimization of the fus ion  deviation).  C.   Adaptive Scaling of the Fusion Deviation  Since the smaller calibration error will cause the smaller  parameter change we can scale the fusion deviations according  to   their   MSEs   (22).   The   scaling   function   applied   before  backpropagation and l earning   can be e.g.:                0 , ) ( MSE 1 max   2 max e d d d  ,   (4 2 )  where   e max   is the   maximal   RMS   of the   fusion deviation  acceptable for usage in calibration. This avoids degradation of  calibration by initial low - quality data.  D.   The Learning Algorithm  It is possible to   choose   from many incremental learning  algorithms .   W e have chosen   the   Adam algorithm   [24] , which  is an extension of the well - known gradient descent algorithm.  Its advantage over   the standard   gradient descent algorithm is  its build - in 1 st   order infinite response f ilter for gradient and  adaptive   learning rate .   As a result, the calibration parameters  c i,k   develop smoother than by   the standard   gradient descent  method.  Goal of the learning is to minimize   the loss function   E :       2 2 real   2 1 2 1   e ω ω      E  .   ( 4 3 )  The   Adam algor ithm requires gradients of the loss function by  each parameter which is given by: k i k i   c c E  , ,            ω  e  (4 4 )  By differentiation of (37) we get:          T 2 , 2 2 , 1 2 , 1  0 1 ) ( 1 ) (                 a a T B w T G a   y y y  ω  ,   ( 45 .1)         T 3 , 3 3 , 1 3 , 1  0 1 ) ( 1 ) (                 a a T B w T G a   z z z  ω  ,   ( 45 .2)         T 1 , 1 1 , 2 1 , 2  0 1 ) ( 1 ) (                 a a T B w T G a   x x x  ω  ,   ( 45 .3)         T 3 , 3 3 , 2 3 , 2  1 0 ) ( 1 ) (                 a a T B w T G a   z z z  ω  ,   ( 45 .4)         T 1 , 1 1 , 3 1 , 3  1 0 ) ( 1 ) (                 a a T B w T G a   x x x  ω  ,   ( 45 .5)          T 2 , 2 2 , 3 2 , 3  1 0 ) ( 1 ) (                 a a T B w T G a   y y y  ω  ,   ( 45 .6)           T , 3 , 2 , 1 ,  ) (   i i i k i i k i  a a a T T B w g         ω  ,   ( 46 )           T , 3 , 2 , 1 ,  1 ) (   i i i k i k i  a a a T T G b          ω  .   ( 47 )  where   G i ( T )   is the value of the polynomial   G i   at the current  normalized   temperature   T   etc.   Speed of learni ng depends on  the setting of learning rate. Because the bias of MEMS sensors  drifts rapidly, learning rate of the coefficients   b i,k   should be  highest. Since the alignment matrix is constant, learning rate  of   a i,k   coefficients should be the smallest. This d istribution  avoids false changes in the alignment matrix according to the  short - term drift of   the   sensor bias.  IV.   E XPERIMENTAL   R ESULTS  For   evaluation   of   the   proposed   fusion   and   calibration  algorithms we have used both the simulation model of the  sensor and t he real MEMS sensor. Simulation allows us to  compare real bias and gain parameters with those obtained by  learning.   We have compared performance of the extended  Kalman filter implemented according to [23], fusion with  fixed gain (also known as   the   compleme ntary filter) and our  proposed   fusion   algorithm   with   adaptive   gain   based   on  estimation of MSE .  A.   Sensor Fusion Simulation  The first   simulation   analyzed   a steady - state at   non zero  attitude     = 30°,     =   - 45°,     = 60° . The simulated sensors’ raw  readings include d noise according to Tab le   1 ; all sensors have  used   the   output   sample   rate   512 Hz .   MSE   of   the   sensor  readings were considered to be equal to the square of RMS  and time - invariant.  TABLE   I  S IMULATED   N OISE   P ARAMETERS  Sensor   RMS   Bias  Gyroscope  Accelerometer  Magnetometer  0.5 °/s  1.0 m/s 2  1 0 % of full range  20 °/ s  none  none
This is the ac cepted version of the   manuscript :   D. Nemec, A. Janota, M. Hruboš and V. Šimák, "Intelligent Real - Time MEMS Sensor Fusion and Calibration,"  in IEEE Sensors Journal, vol. 16, no.   19, pp. 7150 - 7160, Oct.1, 2016,   doi: 10.1109/JSEN.2016.2597292 ,   URL:   http://ieeexplore.ieee.org/document/7529156  The   value   of   the   fixed   fusion   gain   has   been   cho sen  according to the steady value of the adaptive gain (approx.  0.05 according to the Fig.   4 ).   As can be seen in Fig. 3   the  fusion algorithm with adaptive gain has shorter convergence  time than the other methods.   Faster convergence is caused by  the peak of   fusion gain (see Fig. 4) which is a result of initial  low quality of the output   (see (31))   and system   increased  fusion gain in order to compensate rapidly the initial errors by  accelerometer.   Results of fusion   using   fixed gain converge to  the results of t he fusion with adaptive gain after approx. 150  milliseconds.   Output of the Kalman filter has much smoother  response   and lower error   in   a   steady state .  RMS errors of   the estimated   Euler   angle s   in a   steady attitude  (100   seconds   of   experiment)   are   compared   in   Table   2 .  Execution time of the algorithm depends on implementation  and hardware, therefore the   values shown in Table   2   can be  used only for relative comparison between discussed methods.  All   algo rithms   were   implemented   in   the   MATLAB  environment; it is possible to decrease the execution time by  using   a   compiled programming   language like C and   the   code  optimization.  In order to verify dynamic parameters of the sensor fu sion  algorithm we have simulated harmonic rotation around axis   x  (see Fig.   5 ).   Due to the definition of Euler angles in Z - Y - X  convention, the rotation around x - axis will affect only roll  angle.   Noise parameters of the sensors are the same as those  used in   the   previous simulation.   In order to visualize the  difference between compared algorithms, Fig. 5 displays only  the   beginning   of   the   experiment.   Fig.   6   illustrates   one  important advantage of our proposed fusion algorithm over  extended   Kalman   filter   –   stabi lity   in   highly   dynamic  conditions. As can be seen the Kalman filter   becomes instable  after 6 periods of harmonic banking   at frequency 1   Hz , but our  proposed algorithm maintains its stability   and precision . If we  decrease the frequency of the banking (rotat ion)   below 0.5 Hz ,  the instability of the extended Kalman filter disappears.  TABLE   II  C OMPARISON OF   FUSION   ALGORITHMS   IN A STEADY ATTITUD E  Fusion  algorithm  Roll  RMS  [deg]  Pitch  RMS  [deg]  Yaw  RMS  [deg]  Execution time  [ms/sample]  Fixed gain  Adaptive gain  Kalman filter  1. 17  1. 09  0. 58  1. 06  0.93  0. 53  2.3 1  1.56  1.06  0.21  0 .2 3  0. 29  Fig.   6 . Instability of the extended Kalman filter during periodic banking  after longer period.  Fig. 5.   The roll angle estima tion during simulated periodic banking from  - 60° to +60° .  Fig. 4.   The adaptive gain of the roll estimation fusion in a steady state  after reset .  Fig.   3 .   The estimated roll angle in a steady state from noisy data with the  fusion after reset.
This is the ac cepted version of the   manuscript :   D. Nemec, A. Janota, M. Hruboš and V. Šimák, "Intelligent Real - Time MEMS Sensor Fusion and Calibration,"  in IEEE Sensors Journal, vol. 16, no.   19, pp. 7150 - 7160, Oct.1, 2016,   doi: 10.1109/JSEN.2016.2597292 ,   URL:   http://ieeexplore.ieee.org/document/7529156  B.   Sensor Fusion Experiment  Experiments   with   the   real   MEMS   sensor   utilized   a  combined sensor board containing the 3 - axial gy roscope and  accelerometer   IMU - 3000   and   the   3 - axial   magnetometer  LSM303DLH).   The   noise   parameters   of   the   used   sensor  module are   shown   in Tab le   3 .  All noise parameters are estimated in a steady state. The bias  of the acceleromete r and magnetometer is considered to be  zero (due to the previous static calibration).  The sensor board was banked by a servomotor from zero  roll to approx.   -   54 degrees (see Fig.   7 ). Filter parameter   N  used in (12) and (14) was set to   N   = 5   which correspon ds to  the cut - off frequency   f cut   ≈   18 Hz at the sampling rate 512 Hz.  The absolute errors of the constant gain and the adaptive gain  fusion are compared in Fig.   8 .   MSE of the Kalman filter is  slightly higher (0,5%)   than   MSE of our proposed algorithm.  Due   to   the   vibrations   caused   by   the   servomotor   during  movement   the   adaptive   fusion   gain   is   lower   while   the  servomotor is running (see Fig.   9 ). Note that the gain is  stabilized in a steady state in the value   K   = 0.1 which was also  used   in the constant gain fusion in above comparison.  C.   Experimental   V alidation of the   C alibration   A lgorithm  The proposed real - time calibration algorithm was used to  estimate the bias and gain matrices of our gyroscope sensor  module.   Initial calibration parameters   g i,k   and   b i,k   were null.  Since   all   sensors   are   placed   on   one   board   the   sensor  misalignment was neglected; therefore the   A   matrix according  to (39) was not needed to be adjusted. The thermal drifts of  calibration   parameters   we re   not   considered   during   the  experiment because temperature of the sensor module was  stable. If the general temperature - depend e nt version of the  calibration   algorithm   is   used   the   higher - order   thermal  coefficients   g i,k >0   and   b i,k >0   has to be adjusted slowly   because  their real values do not change rapidly during lifetime of the  sensor. Mentioned higher - order calibration parameters allow  faster adaptation to different temperature which can   be   useful  in   some   applications   (e.g.   indoor - outdoor   transition   of   a  mob ile robot).  Fig.   8.   Error of the estimated roll during   the   experiment .  TABLE   III  R EAL   N OISE   P ARAMETERS  Sensor   Full scale   RMS   [ x , y , z ]  Gyroscope    500 °/s   [0.43, 0.41, 0.48]   °/ s  Accelerometer    8 g ≈ 78.5 m/s 2   [0.65, 0.53, 0.51] m / s 2  Magnetometer   4 Gauss   [0.12, 0.19, 0.09] Gauss  Fig. 10.   Raw measured angular velocity around the x - axis during the  experiment .  Fig. 9.   The adaptive fusion gain during   the   experiment .  Fig. 7.   Roll movement during   the   experiment compared with the values  estimated from the accelerometer only and from the adaptive fusion .
This is the ac cepted version of the   manuscript :   D. Nemec, A. Janota, M. Hruboš and V. Šimák, "Intelligent Real - Time MEMS Sensor Fusion and Calibration,"  in IEEE Sensors Journal, vol. 16, no.   19, pp. 7150 - 7160, Oct.1, 2016,   doi: 10.1109/JSEN.2016.2597292 ,   URL:   http://ieeexplore.ieee.org/document/7529156  As   an   experimental   input   we   have   used   a   previously  mentioned combined MEMS sensor module. The module was  moved randomly with no pre - defined pattern (see Fig.   10   for  an illustration) therefore the precise attitude a nd   yaw   are  unknown.   The   raw   (non - calibrated)   accelerometer   and  magnetometer readings as a secondary input have been used in  order to demonstrate robustness of the learning algorithm. The  learning parameters were following:  The   given learning rate is close to the upper limit for usage  in real systems; since the whole fusion algorithm is iterative  (closed loop) higher learning rates could cause instability. The  lower learning rate will result in smoother but also slower bias  deve lopment.   Resultant bias development during automated  learning is shown in Fig.   11 .  Comparison with the actual bias values measured after the  experiment by static calibration is in Tab le   5 . As can be   seen,  the real - time calibration method converges to the real values.  The used gyroscope module has internally compensated  gain therefore the gain deviance is very small (note scale of  the   y - axis in Fig.   12 ).  V.   D ISCUSSION  As can   be seen in Fig.   3   gyroscope bias is effectively  suppressed by the sensor fusion. The fusion with the constant  fusion gain causes smoother output but higher estimation  error.   The adaptive gain allows very fast reactions and start -  up.   According to   Fig.   4   the   adaptive algorithm reflects initial  low quality of roll estimation; therefore the gain (weight of the  attitude estimated by the accelerometer) is initially high and  then rapidly decreases with the roll estimation error.   If we  compare   the   proposed   algorith m   with   the   widely - used  extended Kalman filter, our algorithm   converg es   faster   and it  is stable even during very dynamic changes. On the other  hand,   the   extended Kalman filter has smoother response   and  lower RMS .   Another advantage of our algorithm comes fr om  its lower execution time.  Experiments with   the   real MEMS sensors approved results  obtained by simulations. The MSE estimation algorithm is  able to detect rapid changes in movement including vibrations  (see Fig.   9 ) which adaptively decreases the influenc e of the  accelerometer   (absolute   but   noisy   sensor)   to   the   result.  Estimation of MSE during algorithm execution also provides  additional valuable information about the quality of the result.  If additional information about the object’s state is available,   it  is possible to change input MSE of the accelerometer to reflect  known systematic errors.  Second   experiment   series   evaluated   the   automated  calibration algorithm. According to Fig.   11   the bias learning  works also during movement and slowly converges to th e  precise bias values. Disadvantage of the intelligent calibration  in comparison with laboratory calibration is its lower precision  which can be improved by decreasing of the learning rate.   The  l ower learning rate will however require longer learning time.  The   learning rate   of the bias   B i   has to be much   higher   (100 -  times)   than the learning rate   of the gain   G i , otherwise the  calibration   stability   might   be   corrupted   and   the   overall  calibration parameters would diverge.  VI.   C ONCLUSION  In this paper we have propose d the improved sensor fusion  algorithm. As an explanatory example we have used the fusion  between   the   3 - axial gyroscope (measuring angular velocity),  3 - axial accelerometer (measuring acceleration of the local  system   including   gravity)   and   3 - axial   magnetome ter  (measuring Earth’s magnetic field induction) into estimation  of attitude and   yaw   (AHRS system). The algorithm is based on  Fig. 12. Real - time learning of the gyroscope ’s   gain .  TABL EV  C OMPARISON   B ETWEEN   L EARNED   B IAS AND   V ALUES   E STIMATED BY  S TATIC   C ALIBRATION  Gyroscope  axis  Learned  bias  [ rad/s ]  Bias obtained by  static calibration  [rad/s]  Error of th e learned  bias  [% of the full scale]  x   - 2.0·10 - 2   - 2.52·10 - 2   0.05   %  y   - 0.9·10 - 2   - 1.19·10 - 2   0.03 %  z   1.5·10 - 2   1.26·10 - 2   0.03 %  TABLE   IV  L EARNING   P ARAMETERS  Learning Parameter   Symbol   Valu e  Bias learning rate   λ bias   10 - 6  Gain learning rate   λ gain   10 - 8  1 st   momentum filter [24]   β 1   0.999  2 nd   momentum filter [24]   β 2   0.9999  Maximal fusion MSE   e max   5 °/s  Fig.   1 1 .   Gyroscope bias real - time learning.
This is the ac cepted version of the   manuscript :   D. Nemec, A. Janota, M. Hruboš and V. Šimák, "Intelligent Real - Time MEMS Sensor Fusion and Calibration,"  in IEEE Sensors Journal, vol. 16, no.   19, pp. 7150 - 7160, Oct.1, 2016,   doi: 10.1109/JSEN.2016.2597292 ,   URL:   http://ieeexplore.ieee.org/document/7529156  the estimation of the mean square error during run - time. Then  the theoretical optimal fusion gain is computed. According to  Fig. 4   and Fig. 5 the adaptive gain has   comparable results with  the carefully chosen fixed gain fusion, but has much better  dynamic characteristics (at least 5 - times shorter rise time).  Additionally ,   the parameters of adaptive systems are easier to  measure direc tly (e.g. noise parameters of the used sensors are  usually available from their manufacturer) comparing with the  difficulty of proper fixed fusion gain selection. However,  output MSE should be considered only as a qualitative factor  since the formula used   for estimation of the mean square error  is only the first order approximate. Main disadvantage of the  adaptive fusion algorithm is the higher CPU load (approx. 2 -  times more CPU time needed for MSE estimation compared  with the fixed - gain based fusion algori thm).  Second part of the paper proposes the innovative run - time  calibration method based on processing of fusion deviation  data.   The system was designed to utilize any incremental  stochastic optimization (learning) method; in this article we  have deployed   learning   method called Adam, which is an  extension of   the   gradient descent   method   [24].   Because the  learning algorithm is using small learning rate ,   it is resistant to  the occasional high - power noise contained in fusion data. This  feature is supported by e valuation of the mean square error of  the   fusion.   In   the   discussed   case   of   inertial   attitude  measurement our calibration algorithm is most suitable if the  measured movement is not continuous because in a steady  state the quality of   estimated Euler angles   r ises with time. The  algorithm automatically recognizes such a state and measured  data ha ve   greater impact on the calibration parameters. Since  the fusion is heterogeneous (absolute sensor data are merged  with derivative sensor data) the bias of the absolut e sensor  does not affect the calibration procedure.  Although the fusion method proposed by this manuscript is  derived for this special case, authors believe it can be used in  many other applications as an alternative to the extended  Kalman   filter.   The   prop osed   fusion   method   is   especially  suitable ,   if the sensor readings are processed by non - linear  functions   in   a   recursive   way.   The   proposed   fusion   and  calibration methods can be adjusted to other sensor fusion  scenarios (e.g. combination of the GNSS system   –   absolute  velocity   and   a   position   sensor   and   accelerometer   –   a  differential velocity sensor; an impulse volume sensor and a  flow sensor and many others).  A CKNOWLEDGMENT  This work has been supported by the European Regional  Development Fund and the Ministry   of Education of the  Slovak Republic, within the project ITMS 26220220089 “New  methods of measurement of physical dynamic parameter s   and  interactions of motor vehicles, traffic flow and road . ”  R EFERENCES  [1]   W.   Yuanxin ,   and   S.   Wei ,   “On   Calibration   of   Three - Axis  Magnetometer , ”   Sensors Journal, IEEE , vol .   15,   no.   11, pp. 6424 -  6431 .  [2]   Z.   Zhang ,   and   G.   Yang,   "Calibration   of Miniature   Inertial   and  Magnetic   Sensor   Units   for   Robust   Attitude   Estimation,"   IEEE  Transactions on   Instrumentation and Measurement ,   vol.6 3, no.3,  pp.711 - 718 .  [3]   Q.   Gan ,   and   C.J.   Harris, "Comparison of two measurement fusion  methods   for   Kalman - filter - based   multisensor   data   fusion,"   IEEE  Transactions on   Aerospace and Electronic Systems,   vol.37, no.1,  pp.273 - 279 .  [4]   H.   Qasem , O.   Gorgis ,   and   L.   Reindl, "Design and calibration of an  inertial   sensor   system   for   precise   vehicle   navigation,"   in   Proc.  WPNC ,   2008 , pp.229 - 231 .  [5]   L.   Wang ,   Y.   Hao ,   and F.   Wang, "Calibration of low cost MEMS  inertial Measurement Unit for an FPGA - based navigation system, "   in  ICIA , 2011 ,   pp.181 - 186 .  [6]   L.   You ,   J.   Georgy,   N.   Xiaoji ,   L. Qingli ,   and   N.   El - Sheimy,  “Autonomous Calibration of MEMS Gyros in Consumer Portable  Devices , ”   Sensors Journal ,   IEEE , vol .   15,   no.   7, pp. 4062 - 4072 .  [7]   S.A.   Zotov,   B.R.   Simon,   A.A.   Trusov ,   a nd   A.M.   Shkel, “High  Quality Factor Resonant MEMS Accelerometer   w ith Continuous  Thermal Compensation , ”   Sensors Journal, IEEE , vol .   15,   no.   9, pp.  5045 - 5052 .  [8]   Z.   Yingwei, “Cubature + Extended Hybrid Kalman Filtering Method  and   i ts   Application   in   PPP/IMU   Tightly   Coupled   Navigation  Systems , ”   Sensors Journal, IEEE , vol .   15,   no.   12, pp. 6973 - 6985 .  [9]   A.   Assa ,   and F.   Janabi - Sharifi, “A Kalman Filter - Based Framework  for Enhanced Sensor Fusion , ”   Sensors Journal, IEEE , vol .   15,   no.   6,  pp. 3281 - 3292 .  [10]   G.   Mirzaei ,   M. M.   Jamali,   J.   Ross,   P.V.   Gorsevski ,   and   V.P.  Bingman, “Data Fusion of Acoustics, Infrared, and Marine Radar for  Avian Study , ”   Sensors Journal, IEEE , vol .   15,   no.   11, pp. 6625 - 6632 .  [11]   P.   Chongyan ,   and   S.   Shuli ,   “Fusion   Predictors   for   Multisensor  Stochasti c   Uncertain   Systems   w ith   Missing   Measurements   and  Unknown Measurement Disturbances , ”   Sensors Journal, IEEE , vol .  15,   no.   8, pp. 4346 - 4354 .  [12]   L .   Chang ,   K .   Li ,   and   B .   Hu ,   “Huber’s M - Estimation - Based Process  Uncertainty Robust Filter for Integrated INS/GPS , ”   Se nsors Journal,  IEEE , vol .   15,   no.   6, pp. 3367 - 3374 .  [13]   M .   Zhong ,   D .   Guo ,   and   J .   Guo ,   “PMI - Based   Nonlinear   H ∞  Estimation   of   Unknown   Sensor   Error   for   INS/GPS   Integrated  System , ”   Sensors Journal, IEEE , vol .   15,   no.   5, pp. 2785 - 2794 .  [14]   A.   Janota,   V.   Šimák,   D.   Nemec,   andJ.   Hrbček ,   „ Improving precision  and speed of Euler angles computing from low cost sensor data , “  Sen sors ,   vol.   15 , no.   3,   pp.   7016 - 7039 .  [15]   J.   Diebel .   (2006).   Representing   Attitude:   Euler   Angles,   Unit  Quaternions, and Rotation Vectors .   Stanford University . [Online].  Available:  http://www.astro.rug.nl/software/kapteyn/_downloads/attitude.pdf  [16]   R.   Hostettler ,   W.   Birk,   M.L.   Nordenvaad, “Joint Vehicle Trajectory  and Model Parameter Estimation Using Road Side Sensors , ”   Sensors  Journal, IEEE , vol .   15,   no.   9, pp. 5075 - 5086 .  [17]   V.   Kreinovich ,   et al. , “Towards Combining Probabilistic and Interval  Uncertainty in Engine ering Calculations: Algorithms for Computing  Statistics   under   Interval   Uncertainty,   and   Their   Computational  Complexity , ”   Reliable Computing ,   v ol .   12,   no.   6, pp. 471 - 501 .  [18]   J.   Wei , L.   Yong , and C.   Rizos,   “Optimal Data Fusion Algorithm for  Navigation   Using   Tri ple   Integration   of   PPP - GNSS,   INS,   and  Terrestrial Ranging System , ”   Sensors Journal, IEEE , vol .   15,   no.   10,  pp. 5634 - 5644 .  [19]   L.   Meyer,   A.   Buhmann,   D.C.   Meisel,   and   J.G.   Korvink,  “Relationship Between Zero - Rate Output and the MEMS Element in
This is the ac cepted version of the   manuscript :   D. Nemec, A. Janota, M. Hruboš and V. Šimák, "Intelligent Real - Time MEMS Sensor Fusion and Calibration,"  in IEEE Sensors Journal, vol. 16, no.   19, pp. 7150 - 7160, Oct.1, 2016,   doi: 10.1109/JSEN.2016.2597292 ,   URL:   http://ieeexplore.ieee.org/document/7529156  a Closed - Loop Sys tem , ”   Sensors Journal, IEEE , vol .   15,   no.   12, pp.  7200 - 7207 .  [20]   S.   Luczak,   “Experimental   Studies   of   Hysteresis   in   MEMS  Accelerometers: A Commentary , ”   Sensors Journal, IEEE , vol .   15,  no.   6, pp. 3492 - 3499 .  [21]   M.   Hoffman,   P.   Bauer,   B.   Hemrnelman,   and A.   Hasan, "Ha rdware  synthesis of artificial neural networks using field programmable gate  arrays and fixed - point numbers,"   in   Region 5   IEEE   Conference,  2006 ,   pp.324 - 328 .  [22]   F.   Jiancheng ,   and   L.   Zhanchao, “In - Flight Alignment of POS Based  on State - Transition Matrix , ”   Senso rs Journal, IEEE , vol .   15,   no.   6,  pp. 3258 - 3264 .  [23]   R. Munguia and A. Grau, “A Practical Method for Implementing an  Attitude and Heading Reference System,”   International Journal of  Advanced Robotic Systems ,   vol. 11, no. 62, pp. 1 – 12, 2014 .  [24]   D.   P.   Kingma,   J.   L .   Ba,   “Adam:   A   Method   for   Stochastic  Optimization,”   International   Conference   on   Learning  Representations ,   2015,   arXiv:1412.6980  Dušan Nemec   was born in   Ž ilina, Slovakia in 1991.  He received his MSc. degree in Automation from  the University of Ž ilina in 2015.   He is currently  pursuing the Ph.D. degree in   automation at the Dep .  of Control and Information Systems of the Faculty  of Electrical Engineering at the   University of   Ž ilina.  His research area is focused on mobile robotics and  sensor systems, especially inertial navigation and  flying mobile robots. In 2015 he became a laureate  of the award: “Student personality of the Slovakia of  the   academic   year   2014/20 15   in   the   area   of  Electrical Engineering and Industrial Technologies”.  Aleš Janota   was born in 1963 and received his MSc.  degree from Technical University of Transport and  Communications,   Ž ilina, Czechoslovakia in 1981.  He   post - graduated   from   Universit y   of   Zilina  (UNIZA),   Slovakia,   in   1998   in   the   field   of  Telecommunications.   From 2003 to   2009 he   acted as  an   Assistant Professor for Information and Safety -  related systems. Since 2010 he has worked as a  professor in Control Engineering at the Dep t.   of  Contr ol   and   Information   Systems ,   UNIZA .   His  research   interests   include   sensors,   artificial  intelligence and   intelligent transportation systems.   From   2010   to   2014 he was  a national delegate in the Domain Committee for Transport and Urban  Development of the COST   program.  Mari á n   Hruboš   was   born   in   1987   in   Martin,  Slovakia. Currently   works   as a   postdoc   researcher at  the Dep t.   of Control and Information Systems of the  Faculty of Electrical Engineering at the University of  Ž ilina   (UNIZA) .   His   research   is   focused   on  industrial   automation,   electrical   engineering   and  programming,   particularly   on   data   fusion   from  multiple   sensors   for   creat ionof   3D macro - textured  model s   of the real space sand their   use   in both virtual  and   real   applications.   During   his   Bc.   and   MSc.  s tudies he was with   the Institute of Competitiveness and Innovations ,   UNIZA .  H e   was   awarded   a   Certificate   o f   Merit   from   the   Association   of  Electrotechnical Industry   of the Slovak Republic   as a designer of the year  2014.  Vojtech Šimák   was b orn in   Ž ilina, Czechoslovakia,  in 1980. He received the M.S. degree in information  and safety - related systems (2004) and the Ph.D.  degree   in   Control   Engineering   (2008)   from   the  Faculty   of   Electrical   Engineering,   University   of  Ž ilina   (UNIZA) , Slovakia.   Duri ng his   Ph.D. study  he stayed for   5 months at the Helsinki University of  Technology in the Control Engineering Laboratory  (2007). From 2007 to 2009, he was a researcher  with the Dep t.   of Industrial Engineering, Faculty of  Mechanical Engineering,   UNIZA . Sinc e 2009 he   has   worked with the Dep t.   of  Control and Information Systems at the Faculty of Electrical Engineering,  UNIZA ,   Slovakia.   His   research   interests   include   industrial   and   mobile  robotics and positioning systems, particularly navigation, inertial syste ms and  computer vision.