Improving the Precision and Speed of Euler Angles Computation from Low-Cost Rotation Sensor Data

This article compares three different algorithms used to compute Euler angles from data obtained by the angular rate sensor (e.g., MEMS gyroscope)—the algorithms based on a rotational matrix, on transforming angular velocity to time derivations of the Euler angles and on unit quaternion expressing rotation. Algorithms are compared by their computational efficiency and accuracy of Euler angles estimation. If attitude of the object is computed only from data obtained by the gyroscope, the quaternion-based algorithm seems to be most suitable (having similar accuracy as the matrix-based algorithm, but taking approx. 30% less clock cycles on the 8-bit microcomputer). Integration of the Euler angles’ time derivations has a singularity, therefore is not accurate at full range of object’s attitude. Since the error in every real gyroscope system tends to increase with time due to its offset and thermal drift, we also propose some measures based on compensation by additional sensors (a magnetic compass and accelerometer). Vector data of mentioned secondary sensors has to be transformed into the inertial frame of reference. While transformation of the vector by the matrix is slightly faster than doing the same by quaternion, the compensated sensor system utilizing a matrix-based algorithm can be approximately 10% faster than the system utilizing quaternions (depending on implementation and hardware).


Introduction
Micro-Electro-Mechanical systems (MEMS) represent the integration of mechanical elements, sensors, actuators, and electronics on a common silicon substrate through the utilization of microfabrication technology [1]. The number of MEMS used in various applications is permanently growing due to the small dimensions, light weight, lower power consumption, higher reliability, and relatively low cost which makes them commercially available. Typical MEMS-based low-cost products are accelerometers, gyroscopes, pressure sensors, microphones, digital mirror displays, micro pumps, etc. For the purpose of low-cost navigation solutions MEMS-based inertial sensors (accelerometers and gyroscopes) have been developed since orientation of an object in the three-dimensional space is key information needed for navigation, guidance and control tasks. MEMS inertial sensors may be found in variety of applications from traditional ones (navigation and positioning of various transport means and/or robots) to sensing of human body walking and movement [2][3][4][5], daily life surveillance [6] or new commercial applications available through smart phones [7]. Most studies on MEMS gyroscopes are focused on their performance, and common methods to improve the performance [8]. Unlike non-micro devices MEMS sensors experience more errors that build up over time, corrupting the precision of the measurements and eventually rendering the navigation solution useless [9,10]. Thus the first and easiest-to-measure performance criterion of a gyro is its static readout as a function of time. Accuracy is usually limited by electrical noise, systematic errors and/or mechanical thermal noise [11,12]. The static compensation of sensor inaccuracies can be enabled by proper calibration methods designed for MEMS gyroscopes and accelerometers [13]. The principle of recently developed micro-machine gyroscopes, their structures and classification can be found in [14].
Generally, gyroscopes measure rotational rate, which can be integrated to yield changes in orientation. An effective method most used to parametrize the orientation space is based on usage of so called Euler angles. Euler angles are used as a framework for formulating and solving the equations for conservation of angular momentum. This article has been written with motivation to analyze and show how precision and speed of computations of Euler angles could be improved when processing data from the MEMS gyroscope. It is organized as follows: Section 1 (Introduction) describes theoretically several methods of notation to express rotation of a body (particularly the rotation matrix, Euler angles, rotation around arbitrary axis, and quaternion). Section 2 (Experimental Section) is focused on comparison of errors occurring when algorithms utilizing described notations process data from the gyroscope. If applicable, more versions of the same algorithm are considered (focused either on accuracy or fastness of computation). At the end of the section there is discussion on how errors presented in real gyroscopes could be compensated. Section 3 (Results and Discussion) summarizes analyzed properties and gives final comparison and overview of obtained results. Finally, Section 4 gives the conclusions. The article is an extended version of the conference paper [15], elaborated and supported by the VEGA1/0453/12 grant and used with kind permission of Springer Science + Business Media. Article extensions resulted from the work under another project as stated in the Acknowledgments section.
The purpose of the inertial navigation in the 3D space is to determine six independent variables: translation of an object in three axes and its rotation in three axes, relative to the inertial frame of the reference body. In this article we describe possible ways how to express rotation (attitude) of the object and calculate it from angular velocity measured by the gyroscope.
We consider the Cartesian (orthonormal) right-hand coordinate system oriented by convention NED, i.e., North-East-Down. Moving object axes' orientations are x → forward, y → right and z → down ( Figure 1). Two reference frames are used: • Frame of reference joined with Earth (considered to be approximately inertial), marked S. All variables measured with respect to Earth will be marked without a dash.
• Frame of reference joined with rotating object, marked S'. All variables measured onboard the moving object will be marked with a dash. First we will analyze four used methods of notation that allow us to express rotation of a body. Differences among those individual approaches can be seen in data redundancy and consumption of computer time during processing of raw data from the gyroscope and during conversion from one notation to another (which has direct impact on algorithm efficiency).

Euler Angles
Euler angles are expressing rotation of the object as a sequence of three rotations around objects' local coordinate axes. This way of rotation expression is most interpretative and has zero data redundancy because only three real numbers are needed. Different sequence of axis rotation produces different resultant rotation; therefore Euler angles are defined according to chosen sequence (convention). In aviation the most used convention is z-y-x convention (sometimes called Yaw-Pitch-Roll convention or 3-2-1, see Figure 2): 1. Rotate the object around its z-axis by angle Yaw (marked γ); 2. Rotate the object around its new y1-axis by angle Pitch (marked β); 3. Rotate the object around its new x2-axis by angle Roll (marked α).
Rotation order of z-y-x convention can be expressed by the following operator: Inverse rotation is given by the reversed rotation order by inverted angles: Main disadvantage of representing object's rotation by Euler angles is a lack of the simple algorithm for vector transformation. This can be realized by transferring Euler angles to the rotation matrix by Equation (10) and following application of Equation (4). Trivial chaining (adding) of two rotations represented by Euler angles is not possible.

The Rotational Matrix
The rotation matrix defines change of coordinates of the object in the coordinate system S during rotational movement. It is a typical representation of object's attitude (very often used, e.g., in computer graphics). It is clear that this form has the greatest data redundancy due to needs of saving nine real numbers: Transformation of coordinates from the system S to the system S' can be done by multiplication of the position column vector r by the rotation matrix: Result of rotation R1 followed by R2 is given by matrix multiplication: Inverse rotation is given by the transposed matrix: While the original vector r has the same length as the resultant vector r', the rotational matrix has to be orthogonal with its determinant equal to 1. The matrix is orthogonal when all its row or column vectors are perpendicular to each other. The following algorithm can be used to normalize the matrix to be pure rotational [16]: 1. Calculate deviations eik from orthogonality of the matrix columns: Conversion from 3-2-1 Euler angles to the rotational matrix is given by the following formula [17]: where: Conversion from the rotational matrix to 3-2-1 Euler angles can be done by the following algorithm: The function atan2(y, x) is a four quadrant inverse tangent function, i.e., arctangent function extended to the output angle interval from -π to π. Inputs x and y are coordinates of any point in 2D plane, output is an oriented angle between x-axis and the vector [x, y]. Function is supported by many programming languages by standard (e.g., C-language), having two arguments. The purpose of using two arguments instead of one is to gather information on the signs of the inputs in order to return the appropriate quadrant of the computed angle, which is not possible for the single-argument arctangent function.

Rotation around Arbitrary Axis
According to the Euler theorem it is possible to replace every rotation representation by simple rotation around angle θ around the arbitrary axis given by the unit vector n = n' = [nx, ny, nz] (length of the axis vector is |n| = 1). Note that the axis vector has the same coordinates in the inertial system S and the body-fixed system S'.
Transformation of the vector r from the system S to S' is expressed by the Rodriguez rotation formula: (13) Inverse rotation is expressed by the identical axis n and opposite angle −θ. Chaining of two rotations around non-parallel axes of rotation is impossible to implement trivially, transformation to another type of expression is needed.

Quaternion
Quaternion (invented by sir William Rowan Hamilton in 1843) is a modification of rotation around arbitrary axis expression utilizing algebra of complex numbers expanded to three imaginary dimensions with the complex units i, j, k, for which it is valid: Based on the expanded Euler's formula, the rotation for quaternion around the axis ] , , [ While the axis n is a unit 3D vector, quaternion must follow unit constraint to be pure rotational: Normalization of quaternion is done by the similar way like normalization of any vector. An approximate formula (like matrix normalization) can be used only if normalization is performed after each update of the quaternion: The advantage of quaternions is quick computing of chaining of rotation q1 followed by q2 utilizing Hamilton's product: (18) There are two basic variants of vector transformation utilizing quaternion. The first one takes the transformed vector as a quaternion 0 Concerning speed it is better to use the following formula: where q is a vector part of quaternion: Conversion from 3-2-1 Euler angles to unit quaternion is given by the following formula:

Experimental Section
In this section we compare errors of Euler angle estimation caused by algorithms processing gyroscopic data and being based on different rotation notations. These errors increase during run-time and depend on sampling frequency. The gyroscope firmly joined with the moving object S' is measuring angular velocity as a tri-component vector . These data are sampled with given sample frequency fsample = 1/ΔT. The sensor system has to process data sample by sample in real-time ( Figure 3). As mentioned above, outputs of the algorithm are Euler angles α, β, γ, the system should also provide utility of the transformation of the vector from the S to S' coordinate system.
In order to eliminate influence of the sensor itself a model of the ideal digital-output gyroscope with the following properties was used for algorithm testing: • Gyroscope output in each axis is a signed integer with 16-bit precision (like in many of available low-cost gyroscopes). Full-scale range of the output angular rate is ±500°/s.
• No noise is present at gyroscope output; also sampling frequency is absolutely precise (we want to examine errors of data processing algorithms, not precision of data itself). Therefore, data simulation was used instead of real experiment. In order to obtain comparable results, simulated movement of the object has to be exactly the same for all experiments. Therefore the pre-defined non-random movement has to be simulated. As a test input for algorithms we used a model of precession motion with perpendicular precession axis (see Figure 4). Such rotational movement is easy to define and also it is possible to analytically compute object's attitude (Euler angles) at any time.
Angular velocity of primary rotation and precession was chosen A = 1 rad·s −1 . Simulated angular velocity of the object (measured in its frame of reference) is then given by following: Simulation time corresponds to 20 turns (tend = 40π/A ≈ 2 min). Euler angles during simulated movement are shown in Figure 5. Initial rotation is {α0 = 0°, β0 = 60°, γ0 = 0°}. Euler angles (3-2-1 convention) during defined movement are given by following: ( ) Note that the difference between two angles has to be computed as angular difference (e.g., difference between 180° and −180° is zero) and the maximal shown error is 180°. Figure 5. Euler angles during one turn of the simulated movement.

The Algorithm Based on Updating of the Rotational Matrix
The first version of the algorithm for processing of measured angular velocity is utilizing a matrix as a primary expression of rotation. The principle of this method is shown in Figure 6. The update matrix defines rotation of the object between 2 recent samples of the angular velocity vector ω' (samples ωn-1 and ωn) with time span ΔT. It is possible to create the update matrix from angular velocity by two ways-precise and fast. Fast version uses linear approximation of sine and cosine functions which significantly reduces computational demands; precise version uses non-linear goniometric functions. There is a possibility of using Taylor series of higher order as an approximation of sine and cosine functions.

Precise Version
We can assume that between 2 samples there is constant angular velocity, so its direction defines rotation axis and magnitude multiplied by sample period ΔT defines the angle of rotation: The corresponding update matrix is: Because of the linearity of equations (there is no need for calculation of trigonometric functions or normalization of the axis vector) this is the fastest of all mentioned methods (it is about 3-times faster than precise version, depending on the used hardware). However, the main disadvantage is low accuracy, which constrains this algorithm for systems with high sampling frequency. Figure 7 compares the fast and precise versions by their relative errors with respect to sampling frequency. Expression of rotation based on rotational matrices does not contain any singularities; therefore it is working with constant precision for every tilt. The advantage is also the quick algorithm of vector transformation. In order to maintain rotation matrix orthogonality, normalization is strongly recommended if fast version of the matrix-based algorithm is used. Shown results are computed after normalization in each step.

The Algorithm Based on the Integration of the Euler Angle Rates
Using this algorithm it is possible to avoid intermediate expression of rotation (e.g., by the matrix) and following need for conversion to Euler angles. The principle is shown in Figure 8. This version uses relation between angular velocity ω' measured in the coordinate system S' and time derivations of Euler angles (Euler angle rates): , we get resulting Euler angles. There are two algorithms of numerical integration used in real-time processing: Step integration: Although trapezoidal integration is usually more precise than simple step integration, according to Figure 9 step integration is in case of Euler angle rates little more precise. This is caused by non-linearity of transformation Equation (34). The algorithm is precise enough only at high sampling frequency. The main disadvantage of this algorithm is singularity of expression Equation (34) in case of cosβ = 0 called gimbal-lock, which is representing the state, when x-axis is pointing downwards or upwards (β = 90° or β = −90° respectively). In surroundings of this singularity numerical error is rising. In case that position reaches this singularity, information about two DoF is lost (see Figure 10).  This error can be avoided by early conversion to another Euler convention which reaches singularity in other points (for example conversion to 1-2-1, 1-3-1, 2-3-1, 3-1-2, 3-1-3 or 3-2-3 Euler angle convention [18]). After calculation of Euler angles in substitute convention, they are transformed back to the primary convention. Accuracy is then achieved in the whole angle range. This is computation demanding non-linear operation [18].

The Algorithm Based on Quaternion
The third possibility is to utilize primary expression of rotation using quaternion. The principle is expressed by Figure 11. Similarly as in the case of the rotational matrix, two variants of calculation are possible.

Precise Version
It is an analogy of the precise matrix-based algorithm. The form of update quaternion is following: Then it is valid: 1 . − = n update n q q q (38) Figure 11. Principle of the quaternion-based algorithm.

Fast Version
Neglecting higher order members, using approximations: We obtain update quaternion in the form: By integration of quaternion derivation by time we get resulting rotation quaternion. Figure 12 compares precision of fast and precise versions of the algorithm. Like the fast matrix-based algorithm also the fast quaternion-based algorithm requires normalization of quaternion after each step. Normalization of rotation quaternion is described by Equation (17). Presented results are obtained with normalization. Figure 12. Errors of the quaternion-based algorithms during simulated movement.

Compensation of MEMS Gyroscope Data Using a MEMS Accelerometer and Magnetic Compass
Results given above are valid in an ideal case when gyroscope data are absolutely precise. Real MEMS gyroscope readings are noisy and sensitive to vibrations. The greatest impact on precision of Euler angles estimation has offset of the gyroscope. Due to variance of parameters of an electro-mechanical system with temperature the offset is also temperature dependent. The aim is to use secondary sensor (accelerometer, magnetic compass) to compensate increasing (offset-caused) error of the gyroscope-only system.
The accelerometer is sensing its acceleration (3D vector) relative to inertial frame of reference. In gravitational field the accelerometer is sensing gravity as acceleration upwards. Reading of the accelerometer is (see Figures 13 and 14): where a' is own acceleration of the object expressed in the coordinate system S', g' is a vector of gravitational acceleration (depending on locality near Earth) transformed to the coordinate system of the object S' based on data concerning object rotation and anoise is the noise caused by: • Vibrations of this object • Thermal noise of the sensor • Quantization noise of the A/D converter Figure 13. Accelerometer and magnetic compass readings at non-zero pitch β and yaw γ. Acceleration aacc is measured by the on-board accelerometer as a sum of the gravity acceleration g and object's acceleration a. Earth's magnetic field induction B has inclination θ, declination δ and its horizontal complement points to magnetic North. According to Figure 14 for roll and pitch angle we obtain: If we assume that noise anoise has zero mean value and lower limiting frequency fmin, then the noise can be effectively suppressed by the low pass filter.
Since we cannot determine the rotation around vertical z'-axis (yaw γ) from accelerometer data, it is necessary to add a magnetic compass to the sensor system. For ensuring proper function of the system for all rotations of the object, the magnetic sensor has to determine magnetic induction B' of the Earth's magnetic field in all three axes (compass output is the vector ] , , [ ). For yaw rotation calculated from readings of the magnetic sensor it is valid: where δm is magnetic declination (offset between magnetic and geographic north direction, depending on actual position on Earth), Bx1 and By1 are components of measured magnetic induction after transformation to the coordinate system S1 (inverted x-and following y-rotation, see definition of Euler angles) according to the formula: In terms of avoiding preparation of partial inverse rotation, it is more convenient to determine the difference between yaw γgyro calculated from gyroscope data and yaw from the magnetic compass γmag as: where Bx a By are components of measured magnetic induction after transformation to the coordinate system S (inverted x-, y-and z-rotation), which are: For fusion of Euler angles measured by the gyroscope as a primary sensor and accelerometer and magnetic compass as secondary sensors we can use the algorithm as shown in Figure 15. Gain K << 1 expresses relative weight of the accelerometer with respect to the gyroscope (if K = 0, the accelerometer does not affect output Euler angles). Delay block and gain forms the first order discrete low pass filter in the accelerometer signal path with cutoff frequency: The fusion schema does not filter out any noise from gyroscope reading; it suppresses increasing error of estimated Euler angles in long term caused mainly by offset.
While the schematics in Figure 15 contains the reverse conversion block from Euler angles to the rotation matrix, normalization of the matrix is no longer needed.

Results and Discussion
Effect of sensor fusion is more significant after longer time (especially at low angular velocities). Figure 16 shows effect of using fusion of gyroscope, accelerometer and magnetic sensor readings. Simulated rotation was slowed down 100-times (A = 0.01, compare with Equation (24)). Fusion gain was K = 0.01, noise in secondary sensor data has SNR = 0dB. The precise quaternion-based algorithm at sampling frequency 1 kHz was used. Due to gyroscope offset the estimation error continuously increases with time. The low pass filter within the data fusion algorithm suppresses noise in secondary sensor data and roll angle obtained by fusion slightly oscillates around actual roll.
As can be seen in Figure 17, sensor fusion with weak bound of secondary absolute but noisy sensor can effectively suppress error of estimation caused by sensor offsets. Fusion gain K has to be set according to offset variance of the gyroscope (the more precise the gyroscope is the lesser fusion gain can be obtained).  The second great aspect of the algorithm is its computational time. Two types of reference hardware were used: • 8-bit low-cost microprocessor (Atmel ATmega1284P running at 20 MHz); • 32-bit microprocessor with FPU and DSP support (Atmel UC3C1512C running at 48 MHz). Table 1 compares computational time of algorithms in terms of the CPU cycles of 8-bit low-cost microprocessor. The mentioned cycle counts are average values from 1000 random inputs, using the mathematical library optimized for AVR 8-bit microcontrollers. Algorithms are using software-implemented single precision floating point arithmetic (according to IEEE 754) due to the fact that AVR microcontrollers do not contain the floating point unit (FPU). Using highly optimized implementation of the matrix-based algorithm including fusion of the gyroscope with accelerometer and magnetic compass allows algorithm sampling rate up to approximately 200 Hz (running on AVR 8-bit core @ 20 MHz). Table 2 shows the same algorithms running on the 32-bit microprocessor. Utilization of the 32-bit microcontroller with FPU significantly reduces the count of needed clock cycles (in case of adding and multiplication of real numbers approx. 30 times depending on the used processor). While the representation of numbers is the same for all architectures (32-bit floating point number), accuracy of the algorithm does not depend directly on the used microcontroller. However, decreasing time needed for one cycle of the algorithm allows higher maximal sample rate (up to maximal sample rate of gyroscope itself). Increasing sample rate will improve accuracy significantly. For example, increasing sampling rate from 200 Hz to 1 kHz will decrease error caused by the algorithm by approx. 50% (see Figures 5, 7 and 10). Table 1. Comparison of methods in terms of 8-bit AVR processor clock cycles.
If the microcontroller with hardware support of floating-point calculations is used, linear (fast) versions of algorithms are much faster than the precise non-linear algorithms. Results given in Tables 1  and 2 strongly depend on implementation of the discussed algorithms (execution speed can be improved by using optimized mathematical libraries for hardware supporting floating-point calculations). Number of the clock ticks is shown mainly for simple comparison purposes.
The Euler angle rates integration can be faster than the remaining two algorithms but it has significantly worse accuracy at the same sampling frequency and also has intrinsic singularity. Therefore the choice should be between the matrix-and quaternion-based algorithms. If the sensor system should be able to quickly transform many vectors between inertial and local frame of reference the matrix-based algorithm can be a better choice (2-3 times faster vector transformation than by quaternion).

Conclusions
By comparing relative errors of each mentioned algorithm we can see that the worst algorithm is direct Euler angles integration due to its singularity. Precise version of the quaternion-based algorithm is slightly faster than the precise matrix-based algorithm. Fast version of the quaternion-based algorithm at lower sampling frequency is also more accurate than the matrix-based algorithm (see Table 3). Difference in accuracy between fast and precise versions of the same algorithm decreases with sampling frequency (see Figures 7 and 12). The choice of the proper algorithm depends on: • Available computational power (CPU) and maximal sampling frequency of the sensors (which reflects in overall cost of the sensor system and its accuracy). At lower sampling frequency the fast quaternion-based algorithm is more precise than the fast matrix-based algorithm.
• Precision requirements (in order to achieve long-term stability the compensated system with sensor fusion has to be used).
• Amount of vectors transformed from the non-rotated coordinate system to rotated coordinates and vice versa (transformation performed by the matrix is faster).