Introduction

This article is a continuation of my IMU Guide, covering additional orientation kinematics topics. I will go through some theory first and then I will present a practical example with code build around an Arduino and a 6DOF IMU sensor (acc_gyro_6dof). The scope of this experiment is to create an algorithm for fusing gyroscope and accelerometer data in order to create an estimation of the device orientation in space. Such an algorithm was already presented in part 3 of my “IMU Guide” and a practical Arduino experiment with code was presented in the “Using a 5DOF IMU” article and was nicknamed “Simplified Kalman Filter”, providing a simple alternative to the well known Kalman Filter algorithm. In this article we’ll use another approach utilizing the DCM (Direction Cosine Matrix). For the reader that is unfamiliar with MEMS sensors it is recommended to read Part 1 and 2 of the IMU Guide article. Also for following the experiments presented in this text it is recommended to acquire an Arduino board and an acc_gyro_6dof sensor.

Prerequisites

No really advanced math is necessary. Find a good book on matrix operations, that’s all you might need above school math course. If you would like to refresh your knowledge below are some quick articles:
Cartesian Coordinate System – http://en.wikipedia.org/wiki/Cartesian_coordinate_system
Rotation – http://en.wikipedia.org/wiki/Rotation_%28mathematics%29
Vector scalar product – http://en.wikipedia.org/wiki/Dot_product
Vector cross product – http://en.wikipedia.org/wiki/Cross_product
Matrix Multiplication – http://en.wikipedia.org/wiki/Matrix_multiplication
Block Matrix – http://en.wikipedia.org/wiki/Block_matrix
Transpose Matrix – http://en.wikipedia.org/wiki/Transpose
Triple Product – http://en.wikipedia.org/wiki/Triple_product

Notations

Vectors are marked in bold text – so for example “v” is a vector and “v” is a scalar (if you can’t distinguish the two there’s problem with the text formatting wherever you’re reading this).

Part 1. The DCM Matrix

Generally speaking orientation kinematics deals with calculating the relative orientation of a body relative to a global coordinate system. It is useful to attach a coordinate system to our body frame and call it Oxyz, and another one to our global frame and call it OXYZ. Both the global and the body frames have the same fixed origin O (see Fig. 1). Let’s also define i, j, k to be unity vectors co-directional with the body frame’s x, y, and z axes – in other words they are versors of Oxyz and let I, J, K be the versors of global frame OXYZ.

Figure 1

Thus, by definition, expressed in terms of global coordinates vectors I, J, K can be written as:

I^G = {1,0,0} ^T, J^G={0,1,0} ^T , K^G = {0,0,1} ^T

Note: we use {…} ^T notation to denote a column vector, in other words a column vector is a translated row vector. The orientation of vectors (row/column) will become relevant once we start multiplying them by a matrix later on in this text.

And similarly, in terms of body coordinates vectors i, j, k can be written as:

i^B = {1,0,0} ^T, j^B={0,1,0} ^T , k^B = {0,0,1} ^T

Now let’s see if we can write vectors i, j, k in terms of global coordinates. Let’s take vector i as an example and write its global coordinates:

i^G = {i_x^G , i_y^G , i_z^G} ^T

Again, by example let’s analyze the X coordinate i_x^G, it’s calculated as the length of projection of the i vector onto the global X axis.

i_x^G = |i| cos(X,i) = cos(I,i)

Where |i| is the norm (length) of the i unity vector and cos(I,i) is the cosine of the angle formed by the vectors I and i. Using the fact that |I| = 1 and |i| = 1 (they are unit vectors by definition). We can write:

i_x^G = cos(I,i) = |I||i| cos(I,i) = I.i

Where I.i. is the scalar (dot) product of vectors I and i. For the purpose of calculating scalar product I.i it doesn’t matter in which coordinate system these vectors are measured as long as they are both expressed in the same system, since a rotation does not modify the angle between vectors so: I.i = I^B.i^B = I^G.i^G = cos(I^B.i^B) = cos(I^G.i^G) , so for simplicity we’ll skip the superscript in scalar products I.i , J.j , K.k and in cosines cos(I,i), cos(J,j), cos(K,k).

Similarly we can show that:

i_y^G = J.i , i_z^G=K.i , so now we can write vector i in terms of global coordinate system as:

i^G= { I.i, J.i, K.i}^T

Furthermore, similarly it can be shown that j^G= { I.j, J.j, K.j} ^T , k^G= { I.k, J.k, K.k} ^T.

We now have a complete set of global coordinates for our body’s versors i, j, k and we can organize these values in a convenient matrix form:

(Eq. 1.1)

This matrix is called Direction Cosine Matrix for now obvious reasons – it consists of cosines of angles of all possible combinations of body and global versors.

The task of expressing the global frame versors I^G, J^G, K^G in body frame coordinates is symmetrical in nature and can be achieved by simply swapping the notations I, J, K with i, j, k, the results being:

I^B= { I.i, I.j, I.k}^T , J^B= { J.i, J.j, J.k}^T , K^B= { K.i, K.j, K.k}^T

and organized in a matrix form:

(Eq. 1.2)

It is now easy to notice that DCM^B = (DCM^G)^T or DCM^G = (DCM^B)^T , in other words the two matrices are translates of each other, we’ll use this important property later on.

Also notice that DCM^B. DCM^G = (DCM^G)^T .DCM^G = DCM^B. (DCM^B)^T = I₃ , where I₃ is the 3×3 identity matrix. In other words the DCM matrices are orthogonal.

This can be proven by simply expanding the matrix multiplication in block matrix form:

(Eq. 1.3)

To prove this we use such properties as for example: i^GT. i^G = | i^G|| i^G|cos(0) = 1 and i^GT. j^G = 0 because (i and j are orthogonal) and so forth.

The DCM matrix (also often called the rotation matrix) has a great importance in orientation kinematics since it defines the rotation of one frame relative to another. It can also be used to determine the global coordinates of an arbitrary vector if we know its coordinates in the body frame (and vice versa).

Let’s consider such a vector with body coordinates:

r^B= { r_x^B, r_y^B, r_z^B} ^T and let’s try to determine its coordinates in the global frame, by using a known rotation matrix DCM^G.

We start by doing following notation:

r^G = { r_x^G , r_y^G , r_z^G } ^T.

Now let’s tackle the first coordinate r_x^G:

r_x^G = | r^G| cos(I^G,r^G) , because r_x^G is the projection of r^G onto X axis that is co-directional with I^G.

Next let’s note that by definition a rotation is such a transformation that does not change the scale of a vector and does not change the angle between two vectors that are subject to the same rotation, so if we express some vectors in a different rotated coordinate system the norm and angle between vectors will not change:

| r^G| = | r^B| , | I^G| = | I^B| = 1 and cos(I^G,r^G) = cos(I^B,r^B), so we can use this property to write

r_x^G = | r^G| cos(I^G,r^G) = | I^B || r^B| cos(I^B,r^B) = I^B. r^B = I^B. { r_x^B, r_y^B, r_z^B} ^T , by using one the two definition of the scalar product.

Now recall that I^B= { I.i, I.j, I.k}^T and by using the other definition of scalar product:

r_x^G = I^B. r^B = { I.i, I.j, I.k}^T . { r_x^B, r_y^B, r_z^B} ^T = r_x^B I.i + r_y^B I.j + r_z^B I.k

In same fashion it can be shown that:

r_y^G = r_x^B J.i + r_y^B J.j + r_z^B J.k
r_z^G = r_x^B K.i + r_y^B K.j + r_z^B K.k

Finally let’s write this in a more compact matrix form:

(Eq. 1.4)

 

Thus the DCM matrix can be used to covert an arbitrary vector r^B expressed in one coordinate system B, to a rotated coordinate system G.

We can use similar logic to prove the reverse process:

(Eq. 1.5)

Or we can arrive at the same conclusion by multiplying both parts in (Eq. 1.4) by DCM^B which equals to DCM^GT, and using the property that DCM^GT.DCM^G = I₃ , see (Eq. 1.3):

DCM^B r^G = DCM^B DCM^G r^B = DCM^GT DCM^G r^B = I₃ r^B = r^B

 

Part 2. Angular Velocity

So far we have a way to characterize the orientation of one frame relative to another rotated frame, it is the DCM matrix and it allows us to easily convert the global and body coordinates back and forth using (Eq. 1.4) and (Eq. 1.5). In this section we’ll analyze the rotation as a function of time that will help us establish the rules of updating the DCM matrix based on a characteristic called angular velocity. Let’s consider an arbitrary rotating vector r and define it’s coordinates at time t to be r(t). Now let’s consider a small time interval dt and make the following notations: r = r (t) , r’= r (t+dt) and dr = r’ – r:

Figure 2

Let’s say that during a very small time interval dt → 0 the vector r has rotated about an axis co-directional with a unity vector u by an angle dθ and ended up in the position r’. Since u is our axis of rotation it is perpendicular to the plane in which the rotation took place (the plane formed by r and r’) so u is orthogonal to both r and r’. There are two unity vectors that are orthogonal to the plane formed by r and r’, they are shown on the picture as u and u’ since we’re still defining things we’ll choose the one that is co-directional with the cross product r x r’, following the rule of right-handed coordinate system. Thus because u is a unity vector |u| = 1 and is co-directional with r x r’ we can deduct it as follows:

u = (r x r’) / |r x r’| = (r x r’) / (|r|| r’|sin(dθ)) = (r x r’) / (|r|² sin(dθ)) (Eq. 2.1)

Since a rotation does not alter the length of a vector we used the property that| r’| = |r|.

The linear velocity of the vector r can be defined as the vector:

v = dr / dt = ( r’ – r) / dt (Eq. 2.2)

Please note that since our dt approaches 0 so does dθ → 0, hence the angle between vectors r and dr (let’s call it α) can be found from the isosceles triangle contoured by r , r’ and dr:

α = (π – dθ) / 2 and because dθ → 0 , then α → π/2

What this tells us is that r is perpendicular to dr when dt → 0 and hence r is perpendicular to v since v and dr are co-directional from (Eq. 2.2):

v ⊥ r (Eq. 2.21)

We are now ready to define the angular velocity vector. Ideally such a vector should define the rate of change of the angle θ and the axis of the rotation, so we define it as follows:

w = (dθ/dt ) u (Eq. 2.3)

Indeed the norm of the w is |w| = dθ/dt and the direction of w coincides with the axis of rotation u. Let’s expand (Eq. 2.3) and try to establish a relationship with the linear velocity v:

Using (Eq. 2.3) and (Eq. 2.1):

w = (dθ/dt ) u = (dθ/dt ) (r x r’) / (|r|² sin(dθ))

Now note that when dt → 0, so does dθ → 0 and hence for small dθ, sin(dθ) ≈ dθ , we end up with:

w = (r x r’) / (|r|² dt) (Eq. 2.4)

Now because r’ = r + dr , dr/dt = v , r x r = 0 and using the distributive property of cross product over addition:

w = (r x (r + dr)) / (|r|² dt) = (r x r + r x dr)) / (|r|² dt) = r x (dr/dt) / |r|²

And finally:

w = r x v / |r|²­ (Eq. 2.5)

This equation establishes a way to calculate angular velocity from a known linear velocity v.

We can easily prove the reverse equation that lets us deduct linear velocity from angular velocity:

v = w x r ­ (Eq. 2.6)

This can be proven simply by expanding w from (Eq. 2.5) and using vector triple product rule (a x b) x c = (a.c)b – (b.c)a. Also we’ll use the fact that v and r are perpendicular (Eq. 2.21) and thus v.r = 0

w x r = (r x v / |r|²­) x r = (r x v) x r / |r|²­ = ((r.r) v + (v.r)r) / |r|²­ = ( |r|²­ v + 0) |r|² = v

So we just proved that (Eq. 2.6) is true. Just to check (Eq. 2.6) intuitively – from Figure 2 indeed v has the direction of w x r using the right hand rule and indeed v ⊥ r and v ⊥ w because it is in the same plane with r and r’.

 

Part 3. Gyroscopes and angular velocity vector

A 3-axis MEMS gyroscope is a device that senses rotation about 3 axes attached to the device itself (body frame). If we adopt the device’s coordinate system (body’s frame), and analyze some vectors attached to the earth (global frame), for example vector K pointing to the zenith or vector I pointing North – then it would appear to an observer inside the device that these vector rotate about the device center. Let w_x , w_y , w_z be the outputs of a gyroscope expressed in rad/s – the measured rotation about axes x, y , z respectively. Converting from the raw output of the gyroscope to physical values is discussed for example here: http://www.starlino.com/imu_guide.html . If we query the gyroscope at regular, small time intervals dt, then what gyroscope output tells us is that during this time interval the earth rotated about gyroscope’s x axis by an angle of dθ_x = w_xdt, about y axis by an angle of dθ_y = w_ydt and about z axis by an angle of dθ_z = w_zdt. These rotations can be characterized by the angular velocity vectors: w_x = w_x i = {w_x , 0 , 0 }^T , w_y = w_y j = { 0 , w_y , 0 }^T , w_z = w_z k = { 0 , 0, w_z }^T , where i,j,k are versors of the local coordinate frame (they are co-directional with body’s axes x,y,z respectively). Each of these three rotations will cause a linear displacement which can be expressed by using (Eq. 2.6):

dr₁ = dt v₁ = dt (w_x x r) ; dr₂ = dt v₂ = dt (w_y x r) ; dr₃ = dt v₃ = dt (w_z x r) .

The combined effect of these three displacements will be:

dr = dr₁ + dr₂ + dr₃ = dt (w_x x r + w_y x r + w_z x r) = dt (w_x + w_y + w_z) x r (cross product is distributive over addition)

Thus the equivalent linear velocity resulting from these 3 transformations can be expressed as:

v = dr/dt = (w_x + w_y + w_z) x r = w x r , where we introduce w = w_x + w_y + w_z = {w_x , w_y , w_z }

Which looks exactly like (Eq. 2.6) and suggests that the combination of three small rotations about axes x,y,z characterized by angular rotation vectors w_x , w_y , w_z is equivalent to one small rotation characterized by angular rotation vector w = w_x + w_y + w_z = {w_x , w_y , w_z }. Please note that we’re stressing out that these are small rotations, since in general when you combine large rotations the order in which rotations are performed become important and you cannot simply sum them up. Our main assumption that let us go from a linear displacement to a rotation by using (Eq. 2.6) was that dt is really small, and thus the rotations dθ and linear displacement dr are small as well. In practice this means that the larger the dt interval between gyro queries the larger will be our accumulated error, we’ll deal with this error later on. Now, since w_x , w_y , w_z are the output of the gyroscope, then we arrive at the conclusion that in fact a 3 axis gyroscope measures the instantaneous angular velocity of the world rotating about the device’s center.

Part 4. DCM complimentary filter algorithm using 6DOF or 9DOF IMU sensors

In the context of this text a 6DOF device is an IMU device consisting of a 3 axis gyroscope and a 3 axis accelerometer. A 9DOF device is an IMU device of a 3 axis gyroscope, a 3 axis accelerometer and a 3 axis magnetometer. Let’s attach a global right-handed coordinate system to the Earth’s frame such that the I versor points North, K versor points to the Zenith and thus, with these two versors fixed, the J versor will be constrained to point West.

Figure 3

Also let’s consider the body coordinate system to be attached to our IMU device (acc_gyro used as an example),

Figure 4

We already established the fact that gyroscopes can measure the angular velocity vector. Let’s see how accelerometer and magnetometer measurements will fall into our model.

Accelerometers are devices that can sense gravitation. Gravitation vector is pointing towards the center of the earth and is opposite to the vector pointing to Zenith K^B. If the 3 axis accelerometer output is A = {A_x , Ay , A_z } and we assume that there are no external accelerations or we have corrected them then we can estimate that K^B = –A. (See this IMU Guide for more clarifications http://www.starlino.com/imu_guide.html).

Magnetometers are devices that are really similar to accelerometers, except that instead of gravitation they can sense the Earth’s magnetic North. Just like accelerometers they are not perfect and often need corrections and initial calibration. If the corrected 3-axis magnetometer output is M = {M_x , My , M_z }, then according to our model I^B is pointing North , thus I^B = M.

Knowing I^B and K^B allows us calculate J^B = K^B x I^B.

Thus an accelerometer and a magnetometer alone can give us the DCM matrix , expressed either as DCM^B or DCM^G

DCM^G = DCM^BT = [I^B, J^B, K^B]^T

The DCM matrix can be used to convert any vector from body’s(devices) coordinate system to the global coordinate system. Thus for example if we know that the nose of the plane has some fixed coordinates expressed in body’s coordinate system as r^B = {1,0,0}, the we can find where the device is heading in other words the coordinates of the nose in global coordinate systems using (Eq. 1.4):

r^G = DCM^G r^B

So far you’re asking yourself if an accelerometer and a magnetometer gives us the DCM matrix at any point in time, why do we need the gyroscope ? The gyroscope is actually a more precise device than the accelerometer and magnetomer are , it is used to “fine-tune” the DCM matrix returned by the accelerometer and magnetometer.

Gyroscopes have no sense of absolute orientation of the device , i.e. they don’t know where north is and where zenith is (things that we can find out using the accelerometer and magnetometer), instead if we know the orientation of the device at time t, expressed as a DCM matrix DCM(t) , we can find a more precise orientation DCM(t+dt) using the gyroscope , then the one estimated directly from the accelerometer and magnetometer direct readings which are subject to a lot of noise in form of external (non-gravitational) inertial forces (i.e. acceleration) or magnetically forces that are not caused by the earth’s magnetic field.

These facts call for an algorithm that would combine the readings from all three devices (accelerometer, magnetometer and gyroscope) in order to create our best guess or estimate regarding the device orientation in space (or space’s orientation in device’s coordinate systems), the two orientations are related since they are simply expressed using two DCM matrices that are transpose of one another (DCM^G = DCM^BT ).

We’ll now go ahead and introduce such an algorithm.

We’ll work with the DCM matrix that consists of the versors of the global (earth’s) coordinate system aligned on each row:

If we read the rows of DCM^G we get the vectors I^B, J^B, K^B. We’ll work mostly with vectors K^B (that can be directly estimated by accelerometer) and vector I^B (that can be directly estimated by the magnetometer). The vector J^B is simply calculated as J^B = K^B x I^B , since it’s orthogonal to the other two vectors (remember versors are unity vectors with same direction as coordinate axes).

Let’s say we know the zenith vector expressed in body frame coordinates at time t₀ and we note it as K^B_0. Also let’s say we measured our gyro output and we have determined that our angular velocity is w = {w_x , w_y , w_z }. Using our gyro we want to know the position of our zenith vector after a small period of time dt has passed we’ll note it as K^B_1G . And we find it using (Eq. 2.6):

K^B_1G ≈ K^B₀ + dt v = K^B₀ + dt (w_g x K^B₀) = K^B₀ + ( dθ_g x K^B₀)

Where we noted dθ_g = dt w_g. Because w_g is angular velocity as measured by the gyroscope. We’ll call dθ_g angular displacement. In other words it tells us by what small angle (given for all 3 axis in form of a vector) has the orientation of a vector K^B changed during this small period of time dt.

Obviously, another way to estimate K^B is by making another reading from accelerometer so we can get a reading that we note as K^B_1A .

In practice the values K^B_1G will be different from from K^B_1A. One was estimated using our gyroscope and the other was estimated using our accelerometer.

Now it turns out we can go the reverse way and estimate the angular velocity w_a or angular displacement dθ_{a ­}= dt w_a , from the new accelerometer reading K^B_1A­ , we’ll use (Eq. 2.5):

w­_a = K^B₀ x v_a / | K^B₀|²­

Now v_a = (K^B_1A­ – K^B₀) / dt , and is basically the linear velocity of the vector K^B₀. And | K^B₀|²­­ = 1 , since K^B₀ is a unity vector. So we can calculate:

dθ_{a ­}= dt w_a = K^B₀ x (K^B_1A­ – K^B₀)

The idea of calculating a new estimate K^B_{1 ­} that combines both K^B_1A and K^B_1G is to first estimate dθ as a weighted average of dθ_a and dθ_g :

dθ = (s_a dθ_a + s_g dθ_g) / (s_a + s_g­), we’ll discuss about the weights later on , but shortly they are determined and tuned experimentally in order to achieve a desired response rate and noise rejection.

Note for advanced readers: To improve algorithm we could divide dϴ_g into two components: one orthogonal to K^B₀ and one parallel to it. We’d then only apply averaging formula on the component orthogonal to K^B₀ , the idea being that the accelerometer cannot sense rotation  about K^B₀ (azimuth to zenith) axis, thus the equivalent parallel component  of dϴ_a will always be zero.

And then K^B_{1 ­} is calculated similar to how we calculated K^B_1G:

K^B₁ ≈ K^B₀ + ( dθ x K^B₀)

Why we went all the way to calculate dθ and did not apply the weighted average formula directly to K^B_1A and K^B_1G ? Because dθ can be used to calculate the other elements of our DCM matrix in the same way:

I^B₁ ≈ I^B₀ + ( dθ x I^B₀)

J^B₁ ≈ J^B₀ + ( dθ x J^B₀)

The idea is that all three versors I^B, J^B, K^B are attached to each other and will follow the same angular displacement dθ during our small interval dt. So in a nutshell this is the algorithm that allows us to calculate the DCM₁ matrix at time t₁ from our previous estimated DCM₀ matrix at time t­­₀. It is applied recursively at regular small time intervals dt and gives us an updated DCM matrix at any point in time. The matrix will not drift too much because it is fixed to the absolute position dictated by the accelerometer and will not be too noisy from external accelerations because we also use the gyroscope data to update it.

So far we didn’t mention a word about our magnetometer. One reasons being that it is not available on all IMU units (6DOF) and we can go away without using it, but our resulting orientation will then have a drifting heading (i.e. it will not show if we’re heading north, south, west or east), or we can introduce a virtual magnetometer that is always pointing North, to introduce stability in our model. This situation is demonstrated in the accompanying source code that used a 6DOF IMU.

Now we’ll show how to integrate magnetometer readings into our algorithm. As it turns out it is really simple since magnetometer is really similar to accelerometer (they even use similar calibration algorithms), the only difference being that instead of estimating the Zenith vector K^B vector it estimates the vector pointing North I^B. Following the same logic as we did for our accelerometer we can determine the angular displacement according to the updated magnetometer reading as being:

dθ_{m ­}= dt w_m = I^B₀ x (I^B_1M­ – I^B₀)

Now let’s incorporate it into our weighted average:

dθ = (s_a dθ_a + s_g dθ_g + s_m dθ_m) / (s_a + s_g +_­ s_m)

From here we go the same path to calculate the updated DCM_1­

I^B₁ ≈ I^B₀ + ( dθ x I^B₀) , K^B₁ ≈ K^B₀ + ( dθ x K^B₀) and J^B₁ ≈ J^B₀ + ( dθ x J^B₀),

In practice we’ll calculate J^B₁ = K^B₁ x I^B_1, after correcting K^B₁ and I^B₁ to be perpendicular unity vectors again , note that all our logic is approximated and dependent on dt being small, the larger the dt the larger the error we’ll accumulate.

So if vectors I^B₀, J^B₀, K^B₀ form a valid DCM matrix , in other words they are orthogonal to each other and are unity vectors, then we can’t say the same about I^B₁, J^B₁, K^B₁ , the formulas used for calculating them does not guarantee the orthogonality or length of the vector to be preserved , however we will not get a big error if dt is small, all we need to do is to correct them after each iteration.

First let’s see how we can ensure that two vectors are orthogonal again. Let’s consider two unity vectors a and b that are “almost orthogonal” in other words the angle between these two vectors is close to 90°, but not exactly 90°. We’re looking to find a vector b’ that is orthogonal to a and that is in the same plane formed by the vectors a and b. Such a vector is easy to find as shown in Figure 5. First we find vector c = a x b that by the rules of cross product is orthogonal to both a and b and thus is perpendicular to the plane formed by a and b. Next the vector b’ = c x a is calculated as the cross product of c and a. From the definition of cross product b’ is orthogonal to a and because it is also orthogonal to c – it end up in the plane orthogonal to c , which is the plane formed by a and b. Thus b’ is the corrected vector we’re seeking that is orthogonal to a and belongs to the plane formed by a and b.

Figure 5

We can extend the equation using the triple product rule and the fact that a.a = |a| = 1:

b’ = c x a = (a x b) x a = –a (a.b) + b(a.a) = b – a (a.b) = b + d , where d = – a (a.b) (Scenario 1, a is fixed b is corrected)

You can reflect a little bit on the results … So we obtain corrected vector b’ from vector b by adding a “correction” vector d = – a (a.b). Notice that d is parallel to a. Its direction is dependent upon the angle between a and b, for example in Figure 5 a.b = cos (a,b) > 0 , because angle between a and b is less than 90°thus d has opposite direction from a and a magnitutde of cos(a,b) = sin(b,b’).

In the scenario above we considered that vector a is fixed and we found a corrected vector b’ that is orthogonal to a. We can consider the symmetric problem – we fix b and find the corrected vector a’:

a’ = a – b (b.a) = a – b (a.b) = a + e, where e = – b (a.b) (Scenario 2, b is fixed a is corrected)

Finally in the third scenario we want both vectors to move towards their corrected state, we consider them both “equally wrong”, so intuitively we apply half correction to both vectors from scenario 1 and 2:

a’ = a – b (a.b) / 2 (Scenario 3, both a and b are corrected)
b’ = b – a (a.b) / 2

Figure 6

This is an relatively easy formula to calculate on a microprocessor since we can pre-compute Err = (a.b)/2 and then use it to correct both vectors:

a’ = a – Err * b
b’ = b – Err * a

Please note that we’re not proving that a’ and b’ are orthogonal in Scenario 3, but we presented the intuitive reasoning why the angle between a’ and b’ will get closer to 90°if we apply the above corrective transformations.

Now going back to our updated DCM matrix that consists of three vectors I^B₁, J^B₁, we apply the following corrective actions before reintroducing the DCM matrix into the next loop:

Err = ( I^B₁ . J^B₁ ) / 2

I^B₁^’ = I^B₁ – Err * J^B₁
J^B₁^’ = J^B₁ – Err * I^B₁
I^B₁^’’ = Normalize[I^B₁^’]
J^B₁^’’ = Normalize[J^B₁^’]
K^B₁^’’ = I^B₁^’’ x J^B₁^’’

Where Normalize[a] = a / |a| , is the formula calculating the unit vector co-directional with a.

So finally our corrected DCM₁ matrix can be recomposed from vectors I^B₁^’’, J^B₁^’’, K^B₁^’’ that have been ortho-normalized (each vector constitutes a row of the updated and corrected DCM matrix).

We repeat the loop to find DCM_{2 ,} DCM_{3 ,} or in general DCM _n , at any time interval n.

References

1. Theory of Applied Robotics: Kinematics, Dynamics, and Control (Reza N. Jazar)

2. Linear Algebra and Its Applications (David C. Lay)

3. Fundamentals of Matrix Computations (David S. Watkins)

4. Direction Cosine Matrix IMU: Theory (W Premerlani)

Additional Notes

For the implementation of the algorithm for now see my quadcopter project in particular releases 6/7 have a nice Processing program for visual display of the DCM matrix and a model plane. The entire code is on SVN repository:

https://code.google.com/p/picquadcontroller/source/browse/trunk/

The code is in imu.h file:

https://code.google.com/p/picquadcontroller/source/browse/trunk/imu.h

A PDF Version of this article is available here

DCM Tutorial – An Introduction to Orientation Kinematics by Starlino (PDF, Rev 0.1 Draft)

Please mention and link to the source when using information in this article:

http://www.starlino.com/dcm_tutorial.html

Starlino Electronics  // Spring , 2011