Update: This post is part 1 of a series on sensor calibration. The next part, for gyroscopes and accelerometers, can be found here.
Fitting an ellipsoid
Suppose we have a $n \times 3$ matrix containing readings from a MEMs magenetometer:
$M=\begin{bmatrix}x_1 & y_1 & z_1\\x_2 & y_2 & z_2\\\vdots & \vdots & \vdots\\x_n & y_n & z_n\end{bmatrix}$
$x_i$,$y_i$ and $z_i$ denote the $x$,$y$ and $z$ readings of the magnetometer at time index $i$.
A scatter plot of $M$ would create a rough shape of an ellipsoid.
The general equation form an ellipsoid is:
$Ax^2 + By^2 + Cz^2 + 2Dxy + 2Exz + 2Fyz + 2Gx + 2Hy + 2Iz = 1 \tag{1.0}$
Therefore if the magnetometer outputs all lay on the surface of an ellipsoid, then (1.0) would be satisfied for each of it's outputs for some $A,B,...,I$. The magnetometer readings would also form the following system of equations:
$\begin{align*} Ax_1^2 &+ By_1^2 + Cz_1^2 + 2Dx_1y_1 + 2Ex_1z_1 + 2Fy_1z_1 + 2Gx_1 + 2Hy_1 + 2Iz_1=1\\ Ax_2^2 &+ By_2^2 + Cz_2^2 + 2Dx_2y_2 + 2Ex_2z_2 + 2Fy_2z_2 + 2Gx_2 + 2Hy_2 + 2Iz_2=1\\ \vdots\\ Ax_n^2 &+ By_n^2 + Cz_n^2 + 2Dx_ny_n + 2Ex_nz_n + 2Fy_nz_n + 2Gx_n + 2Hy_n + 2Iz_n=1\\ \end{align*}$
Which can be re-written in matrix form: \[ \begin{bmatrix} x_1^2 & y_1^2 & z_1^2 & 2x_1y_1 & 2x_1z_1 & 2y_1z_1 & 2x_1 & 2y_1 & 2z_1\\ x_2^2 & y_2^2 & z_2^2 & 2x_2y_2 & 2x_2z_2 & 2y_2z_2 & 2x_2 & 2y_2 & 2z_2\\ \vdots & \vdots & \vdots & \vdots & \vdots &\vdots & \vdots & \vdots & \vdots\\ x_n^2 & y_n^2 & z_n^2 & 2x_ny_n & 2x_nz_n & 2y_nz_n & 2x_n & 2y_n & 2z_n\\ \end{bmatrix} \textbf{P} = \begin{bmatrix} 1\\ 1\\ \vdots\\ 1 \end{bmatrix} \tag{1.1} \] where $\textbf{P} = [A, B, C, D, E, F, G, H, I]^T$
Hence finding an ellipsoid whose surface contains the magnetometers output rest on solving the 9 parameters of $\textbf{P}$ that satisfy the system of equations. MEMs magnetometers are not perfect however. So instead we find a $\textbf{P}$ which "best fits" the set of magnetometer readings.
How could we pick desirable parameters for $\textbf{P}$? By letting $D$ be the leftmost matrix in (1.1) and $b$ be an $n \times 1$ vector of ones, a suitable set of parameters $\textbf{P}$ is one that satisfies:
\[\underset{\textbf{P} \in \mathbb{R}^9 }{\min} \|D\textbf{P}-b\|_2^2 \tag{1.2} \]
The above simply states that $\textbf{P}$ "best fits" an ellipsoid if $\textbf{P}$ minimizes the distance from the surface of the ellipsoid and 1. In other words the system of equations is as close to 1 as possible for each $i$.
As it happens, finding $\textbf{P}$ amounts to solving the following equation for $\textbf{P}$:
$D^TD\textbf{P} = D^Tb \tag{1.3}$
Interested readers can find a proof to why that is here.
Hard Iron Estimations
Lets begin by considering how the $x$,$y$ and $z$ variables are grouped in (1.0):
$\underbrace{Ax^2 + By^2 + Cz^2}_{\text{Squared} } + \underbrace{2Dxy + 2Exz + 2Fyz}_{\text{Combination} } + \underbrace{2Gx + 2Hy + 2Iz}_{\text{Linear} }\underbrace{-1}_\text{Constant} = 0$
An ellipsoid may be arbitrarily rotated or have it's center translated. For example, compare the form of the above equation to the general form of an unrotated ellipsoid centered at the origin :
$\frac{x^2}{p_x^2} + \frac{y^2}{p_y^2} + \frac{z^2}{p_z^2} -1 = 0 \tag{2.0}$
$p_x$,$p_y$ and $p_z$ each denote lengths of the principle axes. Though (2.0) contains squared and constant terms, it doesn’t contain any combination or linear terms. This suggests that the combination and linear terms have a relationship with the rotation and translation of an ellipsoid.
To illustrate this relationship, let $E_f$ be the ellipsoid which has been fitted to the magnetometer data using Section 1. Suppose $E_0$ is an unrotated ellipsoid that is centred at the origin and has principle axis of equal length to those of $E_f$. Let:
\[ R = \begin{bmatrix} a &b &c\\ d &e &f\\ g &h &i\\ \end{bmatrix} \] Where $R$ is a rotation matrix which rotates $E_0$ to the orientation of $E_f$. That is, $R$ aligns the principle axis of $E_0$ to those of $E_f$. To see how a rotation affects the form of (2.0) let $x\prime, y\prime$ and $z\prime$ denote the respective results of rotating points $x,y$ and $z$ using $R$:
\[ \begin{bmatrix} a &b &c\\ d &e &f\\ g &h &i\\ \end{bmatrix} \begin{bmatrix} x\\ y\\ z\\ \end{bmatrix} = \begin{bmatrix} ax + by + cz\\ dx + ey + fz\\ gx + hy + iz\\ \end{bmatrix} = \begin{bmatrix} x\prime \\ y\prime\\ z\prime\\ \end{bmatrix} \tag{2.1} \] Substituting $x\prime, y\prime$ and $z\prime$ into the $x$, $y$ and $z$ variables of (2.0), allows us to track how the rotation affects the general form. After expanding and grouping:
\begin{align*} (\frac{a^2}{p_x^2} &+ \frac{d^2}{p_y^2} + \frac{c^2}{p_z^2})x^2 +(\frac{b^2}{p_x^2} + \frac{e^2}{p_y^2} + \frac{h^2}{p_z^2})y^2 +(\frac{c^2}{p_x^2} + \frac{f^2}{p_y^2} + \frac{i^2}{p_z^2})z^2 + \\ 2(\frac{ab}{p_x^2} &+ \frac{de}{p_y^2} + \frac{gh}{p_z^2})xy +2(\frac{ac}{p_x^2} + \frac{df}{p_y^2} + \frac{gi}{p_z^2})xz +2(\frac{bc}{p_x^2} + \frac{ef}{p_y^2} + \frac{hi}{p_z^2})yz = 1 \end{align*}
Don't worry about the actual equation, instead note the form that is obtained. $a,b,...,i$ are all constants, as are the principle axes $p_x,p_y$ and $p_z$. Hence the above equation can be relabeled as:
$A_rx^2 + B_ry^2 + C_rz^2 + 2D_rxy + 2E_rxz + 2F_ryz = 1 \tag{2.2}$
It's still unknown if the coefficients in (2.2) are equal to the squared and combination terms in (1.0). This is why the subscript of $r$ , meaning "rotation only", is used. Regardless, this does show that the combination terms are related to how the ellipsoid is rotated as those are the only new terms. To see how translation affects the terms, substitute $(x-x_0),(y-y_0)$ and $(z-z_0)$ for the respective $x,y$ and $z$ variables in (2.2). Expanding and grouping like terms gives:
$ A_rx^2 + B_ry^2 + C_rz^2 + 2D_rxy + 2E_rxz + 2F_ryz +\\ 2(-A_rx_0 -D_rz_0 - E_ry_0)x + 2(-B_ry_0 -E_rx_0 - F_rz_0)y + 2(-C_rz_0 -D_rx_0 - F_ry_0)z +\\ (A_rx_0^2 + B_ry_0^2 + C_rz_0^2 + 2D_rx_0z_0 + 2E_rx_0y_0 + 2F_ry_0z_0) = 1 $
Since $x_0,y_0$ and $z_0$ are all constant the above can be relabelled:
$A_rx^2 + B_ry^2 + C_rz^2 + 2D_rxy + 2E_rxz + 2F_ryz + 2G_tx + 2H_ty + 2I_tz + l - 1 = 0 \tag{2.3}$
Again it is uncertain that $G_t, H_t$ and $I_t$ are equal to $G,H$ and $I$ respectively, so they are given the subscript $t$ to denote "translation." The above expression is nearly in the same form as (1.0) however there is one exception. (1.0) did not have a constant term in it other than $-1$, yet (2.3) has $l - 1$. To arriving at (1.0)'s form we must divide (2.3) by $-w$, where $w = l - 1$. This results in:
\begin{alignat*}{1} \frac{A_r}{-w}x^2 &+ \frac{B_r}{-w}y^2 &+ \frac{C_r}{-w}z^2 &+ 2\frac{D_r}{-w}xy &+ 2\frac{E_r}{-w}xz &+ 2\frac{F_r}{-w}yz &+ 2\frac{G_t}{-w}x &+ 2\frac{H_t}{-w}y &+ 2\frac{I_t}{-w}z &- 1 &= 0\\ Ax^2 &+ By^2 &+ Cz^2 &+ 2Dxy &+ 2Exz &+ 2Fyz &+ 2Gx &+ 2Hy &+ 2Iz &- 1 &= 0 \end{alignat*}
When each coefficient: $A_r,B_r,...,I_t$, are divided by $-w$, they become equal to the respective coefficients: $A,B,C,..., I$ in (1.0).
This allows us to derive an equation to find the center of the ellipsoid. Recall that (2.4) gave an equation for an ellipsoid, $E_0$, centred at the origin with equal principle axes to that of the fitted ellipsoid, $E_f$. By letting $v = [x, y, z]$, (2.2) can be put into matrix form.
\[ 0 = vUv^T - 1 \text{ , } U = \begin{bmatrix} A_r &D_r &E_r\\ D_r &B_r &F_r\\ E_r &F_r &C_r \end{bmatrix} \tag {2.4} \]
(1.0) can be put into matrix form as well:
$vSv^T +2v^Tb - 1$, where \[ S = \begin{bmatrix} A &D &E\\ D &B &F\\ E &F &C \end{bmatrix} , b = \begin{bmatrix} G\\ H\\ I\\ \end{bmatrix} \]
Note $A ... I$ are given from the results of section 1. Letting $C_0 = [x_0,y_0,z_0]$ be the vector containing the $x,y$ and $z$ values for the center of $E_f$, $E_0$ can be translated to the center of $E_f$:
$(v-C_0)U(v-C_0)^T - 1 = vUv^T - 2C_0Uv^T + C_0UC_0^T - 1 = 0 \tag{2.5}$
(2.5) is just (2.3) in matrix form, in fact $C_0UC_0^T - 1 = w$ and it was determined that dividing each of the coefficients in (2.3) by $-w$ will give the coefficients in (1.0) i.e:
$S = \frac{1}{-w}U$
Therefore dividing (2.5) by $-w$ gives:
\begin{align*} vSv^T - 2C_0Sv^T -1 &= 0 = vSv^T +2v^Tb - 1 &\Rightarrow\\ - 2C_0Sv^T &= 2v^Tb &\Rightarrow\\ C_0 &= -S^{-1}b \tag{2.6} \end{align*}
Hard iron offsets are defined as how far the ellipsoid has been translated from the origin. This is precisely what $C_0$ contains, so we can use (2.6) to derive $C_0$ after $S$ and $b$ are given via section 1.
Soft Iron Estimations
After the center of the ellipsoid is known the remanding soft iron errors can be compensated for. First the fitted ellipsoid, $E_f$, will be translated to the origin because this will make correcting soft iron distortions easier. In Section 2, (2.3) was obtained by substituting $(x-x_0),(y-y_0)$ and $(z-z_0)$ into the $x,y$ and $z$ variables of (2.2). This had the effect of translating an ellipsoid, $E_0$, from the origin to the center of $E_f$. Therefore to translate $E_f$ to the origin substitute $(x+x_0),(y+y_0)$ and $(z+z_0)$ for the $x,y$ and $z$ variables of (1.0). After expanding and grouping liked terms the following is derived:$Ax^2 + By^2 + Cz^2 + 2Dxy + 2Exz + 2Fyz + \\ 2(G + Ax_0 + Dy_0 + Ez_0)x + 2(H + Dx_0 + By_0 + Fz_0)y + 2(I + Ex_0 + Fy_0 + Gz_0)z + w^{\prime} = 0$
where $w^{\prime}$ is defined as:
$w^{\prime} = (Ax_0^2 + By_0^2 + Cz_0^2 + 2Dx_0y_0 + 2Ex_0z_0 + 2Fy_0z_0 + 2Gx_0 + 2Hy_0 + 2Iz_0 -1)$
The above equation describes an ellipsoid centered at the origin. Therefore it's linear terms are zero:
$Ax^2 + By^2 + Cz^2 + 2Dxy + 2Exz + 2Fyz + w^{\prime} = 0 \tag{3.0}$
To match (3.0)'s form to that of (2.2), divide (3.0) by $-w^{\prime}$. The fitted ellipsoid $E_f$'s matrix form now becomes:
\begin{align*} -\frac{1}{w^{\prime}}\textbf{x}^{T}S\textbf{x} = 1 = \textbf{x}^{T}U\textbf{x} \tag{3.1}\\ \text{$S$ is as defined in section 2 and } \textbf{x} = \begin{bmatrix} x\\ y\\ z \end{bmatrix} \end{align*}
Another way to think about $w$ and $w^{\prime}$ is that they are coefficients which translate the fitted ellipsoid $E_f$. $w$ translates the fitted ellipsoid from the origin to some center as seen in the relation:
$S = -\frac{1}{w}U$
Where as $w^{\prime}$ translates the fitted ellipsoid from it's current center to the origin as seen in the relation:
$-\frac{1}{w^{\prime}}S = U$
However $w^{\prime}$ and $w$ are distinct.
Now that $E_f$ has been translated to the origin soft iron errors can be corrected. First we'll need to characterize how soft iron errors act on a magnetometer. For simplicity, lets examine how the soft iron effects are modeled in 2D (the same principles will apply in 3D).
If a 2-axis magnetometer was not affected by soft iron or hard iron distortions, then it's output would lay on the border of a unit circle centered around the origin as shown in Figure 1. This is also assuming their is no random noise associated with the device.
 |
Figure 1. Noiseless mag. output absent of any soft iron or hard iron distortion
|
It is assumed that soft iron distortion are an linear transformation of that unit circle [1]. For the sake of the example, suppose the soft iron errors are given by the linear transformation:
\[ A = \begin{bmatrix} 5 & 1.5\\ 1.5 & 3 \end{bmatrix} \]
If $\textbf{x} = [x,y]^T$ represent a 2D vector, then $A$ transforms $\textbf{x}$ under the following linear equation:
\begin{align*} A\textbf{x} &= b \hfill \Rightarrow\\ \textbf{x}&=A^{-1}b \tag{3.2} \end{align*}
The matrix form of the equation for a unit circle is:
$\textbf{x}^T\textbf{x} = 1$
Substitute (3.2) into the above equation:
\begin{align*} \textbf{x}^T\textbf{x} &= (A^{-1}b)^TA^{-1}b\\ &=b^TA^{-1}A^{-1}b = 1 \tag{3.3} \end{align*}
Calculating $A^{-1}A^{-1}$ gives:
\[ b^T \begin{bmatrix} 0.0692 & -0.0738\\ -0.0738 & 0.1676 \end{bmatrix} b = 1 \]
Finally letting $b = [b_x, b_y]^T$ will give the equation of a rotated ellipse centered in the origin:
$0.0692b_x^2 +2(-0.0738)b_xb_y + 0.1676b_y^2 = 1 \tag{3.4}$
Plotting the results illustrates how A distorts the unit circle:
 |
Figure 2.
Affects of Linear transformation A on the unit circle
|
 |
Figure 3. Linear transformation A's affect on a set of points
|
Figure 3 shows an important property of the linear transformation. Notice how most of the lines change directions when they go from solid to dashed. You may have also noticed something special. The change in direction between the dashed and solid lines becomes less severe as the points on the ellipse get closer to its principle axis. In fact there are two vectors that don't change directions and both lay on the principle axis:
 |
Figure 4.The two vectors that don't change direction |
Vectors like the two shown in Figure 4 are known as the eigenvectors of $A$. The eigenvectors $v_1,v_2$ of the given $A$ are:
\[ v_1 = \begin{bmatrix} 0.4719\\ -0.8817 \end{bmatrix} , v_2 = \begin{bmatrix} -0.8817\\ -0.4719 \end{bmatrix} \]
Eigenvectors of a linear transformation have the property in which they do not change directions when being transformed, instead they are scaled. Stated mathematically:
$Av_1 = \lambda_1 v_1 \\
Av_2 = \lambda_2 v_2 \tag{3.5}$
Where $\lambda_1$ and $\lambda_2$ denote how much the eigenvectors are scaled by.
Notice that $v_1$ and $v_2$ are both normal vectors. So their magnitudes, or lengths, are both 1. The lengths of the lines illustrated in figure 4 are the same as the major and minor axis of the ellipse. Therefore the two scaling coefficients in (3.5) scale the lengths of $v_1$ and $v_2$ to the lengths of the major and minor axis. $\lambda_1$ and $\lambda_2$ are known as the eigenvalues of $A$. In this example they tell us how long the axes of ellipse are.
(3.5) also gives a relation between $A$'s eigenvectors and eigenvalues. The eigenvectors and eigenvalues give a useful alternate expression for $A$. To see this, let $Q = [v_1 ,v_2]$ be a matrix whose columns are the eigenvectors of $A$. Further let $D_A$ be a diagonal matrix containing the eigenvalues of $A$. Using $Q$, $D_A$ and further manipulation of (3.5) $A$ may be expressed as:
\begin{align*} AQ &= \begin{bmatrix} a_{11} & a_{12}\\ a_{21} & a_{22} \end{bmatrix} \begin{bmatrix} v_{1x} & v_{2x}\\ v_{1y} & v_{2y} \end{bmatrix}\\ &= \begin{bmatrix} a_{11}v_{1x} + a_{12}v_{1y} & a_{11}v_{2x} + a_{12}v_{2y}\\ a_{21}v_{1x} + a_{22}v_{1y} & a_{21}v_{2x} + a_{22}v_{2y} \end{bmatrix} \\ &= \begin{bmatrix} \lambda_1v_{1x} & \lambda_2v_{2x}\\ \lambda_1v_{1y} & \lambda_2v_{2x} \end{bmatrix} = QD_A \Rightarrow \\ AQ &= QD_A \Rightarrow \\ A &= QD_AQ^{-1} \tag{3.6} \end{align*}
(3.6) uses what’s called the eigendecomposition of $A$. Eigendecomposition allows for the powers of $A$ to be expressed as follows:
\begin{align*} A^{-1} &= QD_A^{-1}Q^{-1}\\ A^{1/2} &= QD_A^{1/2}Q^{-1}\\ A^{2} &= QD_A^{2}Q^{-1} \end{align*}
Note: $D_A^{n}$ is computed by applying the power $n$ to each of it's diagonal components. This property will be used shortly, but first it's important to note that not every square matrix has an eigendecomposition. The property that a matrix must have is that it be diagonalizable in order to be factored in the manner shown in (3.6). It will be shown later that we can guarantee the matrices we need to decompose are always diagonalizable. For now it will be assumed that necessary matrices are diagonalizable.
So how do eigenvectors and eigenvalues help in compensating for soft iron errors? Recall that in (3.2) $A$ is what maps the points on the circle to the points on the ellipse. Conversely $A^{-1}$ maps the ellipse into the unit circle. Therefore knowing $A^{-1}$ would allow us to map the output of the mag from an ellipse (affected by soft iron) to a unit circle (calibrated). But we don’t have access to $A^{-1}$ directly. Instead we have to use sections 1 and 2 to estimate the matrix $U = -\frac{1}{w^{\prime}}S$ and then derive $A^{-1}$ from $U$. Eigenvectors and eigenvalues allowed an alternate expression of $A$ in (3.6). This alternate expression is what will be used to derive $A^{-1}$ from $U$ as we'll now see.
$A$ and $U$ can first be related to each other using (3.3):
\begin{align*} b^TA^{-1}A^{-1}b &= b^TUb \Rightarrow\\ U &= A^{-1}A^{-1} \tag{3.7} \end{align*}
$U$'s eigendecomposition is:
$U = Q_UD_UQ_U^{-1}$
Using this, (3.7) and how the powers of $A$ apply to $D_A$, it follows that:
\begin{align*} U &= A^{-1}A^{-1}\\ &=(A^2)^{-1}\\ &= QD_A^{-2}Q^{-1} \Rightarrow\\ UQ &= QD_A^{-2} \Rightarrow \\ U[v_1, v_2] &= [Uv_1, Uv_2] = [\lambda_1^{-2}v_1, \lambda_2^{-2}v_2] \tag{3.8} \end{align*}
Hence (3.8) demonstrates that $U$ and $A$ share the same eigenvectors, further the eigenvalues of $U$ and the eigenvalues of $A$ can be related as:
$D_U^{-1/2} = D_A \tag{3.9}$
We can now solve for $A^{-1}$ using (3.8) and (3.9):
\begin{align*} A^{-1} &= QD_A^{-1}Q^{-1}\\ &= QD_U^{1/2}Q^{-1}\tag{3.10} \end{align*} (3.10) is the result that will be used to calibrate soft iron errors. It shows all that is needed to find $A^{-1}$ are the eigenvalues and eigenvectors of $U$.
Recall that it was stated above that not all square matrix have an eigendecompisition. It was assumed $U$ had one, but is this true in general?
A well documented result of Linear Algebra is if a matrix $U$ is symmetric, then it must be true that:
- The eigenvalues of $U$ are all real
- $U$ is always diagonalizable
\begin{align*} Ax^2 + By^2 + 2Cxy -1 = 0 \Rightarrow\\ \textbf{x}^TU\textbf{x} - 1 = 0 \text{ , } U = \begin{bmatrix} A & C\\ C & B \end{bmatrix} \end{align*}
The same thing occurs in 3D:
\begin{align*} Ax^2 + By^2 + Cz^2 + 2Dxy + 2Exz + 2Fyz -1 = 0 \Rightarrow\\ \textbf{x}^TU\textbf{x} - 1 = 0 \text{ , } U = \begin{bmatrix} A & D & E\\ D & B & F\\ E & F & C \end{bmatrix} \end{align*}
When ellipsoids (and ellipses) are centred at the origin they do not have linear terms, hence they can be expressed by the above forms. Which is why the ellipsoid was translated to the origin at the start of this section, we needed to force it's form to be that of (2.2). Thus the matrix $U$ will always be symmetric, which guarantees that it is diagonalizable.
This example and subsequent findings were in 2D, but everything still applies to 3D using the same methods and analysis. 2D was only chosen for clarity.
Final Algorithm
Now lets tie everything together into an algorithm that will allow calibration of the magnetometer. Given an $n \times 3$ matrix $M$ of magnetometer readings preform the following steps:1. Using magnetometer data construct the $n \times 9$ matrix $D$ from (1.1)
2. Create $\textbf{1}$, an $n \times 1$ vector of ones:
$\textbf{1} = [1 ... 1]^{T}$
3. Find the parameter vector $\textbf{P}$ for the best fit ellipsoid $E_f$:
$\textbf{P} = (D^TD)^{-1}D^T\textbf{1}$
4. Create matrix $S$ and vector $b$ from the parameters in $\textbf{P}$: \[ S = \begin{bmatrix} A &D &E\\ D &B &F\\ E &F &C \end{bmatrix} \textit{ , } b = \begin{bmatrix} G\\ H\\ I \end{bmatrix} \]
5. Find the center $C_0$ of $E_f$ using (2.8):
$C_0 = -S^{-1}b$
6. Calculate $w^{\prime}$ using (3.0):
$w^{\prime} = C_0^TSC_0 + 2C_0^Tb -1$
7. Find the coeficients matrix $U$ describing $E_f$ centered at the origin:
$U = -\frac{1}{w^{\prime}}S$
8. Find the eigenvectors of $U$ and store them columnwise in matrix $Q$:
$Q \gets \textit{eig.vec(}U \textit{)}$
9. Find the eigenvalues of $U$ and store them diagonally in matrix $D_U$:
$D_U \gets \textit{eig.val(}U \textit{)}$
10.Calculate the inverse of the soft iron tranformation matrix $A$:
$A^{-1} = QD_U^{1/2}Q^{-1}$
Magnetometer data can now be calibrated using the following equation:
$M_\text{calib} = A^{-1}(M_\text{raw} - C_0) \tag{4.0} $
To collect proper data sets for $M$ simply spin the magnetometer (or the board it's attached to) around until a satisfactory set of points are obtained. I usually use around 300-500.
As a final note, due to how the above algorithm is constructed, it will transform an ellipsoid into a unit sphere. For whatever reason, if the unit of the magnetometer readings need to be preserved, then the user should pick a value to scale $D_U$ by. Typically this is done by using the minimum principle axis of the ellipsoid. If this is necessary then step 10 on the above Algorithm becomes:
$A^{-1} = Q(r_{\text{min}}D_U^{1/2})Q^{-1} \tag{4.1}$
Where $r_{\text{min}}$ is the minimum principle axis of the fitted ellipsoid. $r_{\text{min}}$ can be found using (3.9) and (3.7). Namely the principle axis of the fitted ellipsoid lies on the diagonal of the matrix $D_A$ as defined in (3.7). However the algorithm gives $D_U$ which can be related to $D_A$ using (3.9) i.e: In other words the inverse square root of each diagonal entry of $D_U$ will give each principle axis of the fitted ellipsoid. From there the minimum one can be chosen as the value which scales $D_U$ in step 10.
This wraps up calibration of hard iron and soft iron distortions. Despite the promise of the method shown in this presentation a potential draw back still exists. That it's potentially inconsistent for very noisy data sets. This inconsistency can lead to very inaccurate estimations in the presence of heavy noise. Next time I will present a solution for this problem using adjusted least square ellipsoid fitting. At the project repository exists a python implementation of the algorithm (EllipsoidFit) and a demo (EFitDemo) that generates an ellipsoid with user defined noise and runs it through the fitting program.
Thanks for reading!
Thanks for sharing.
ReplyDelete