I'm using an LSM303DLM on a Pololu MinIMU-9 board as part of an AHRS. The software is ported to mbed from Pololu's example code which is in turn ported from ArduPilot 1.5. But that's no what this is about.
This is a saga of one man... with meager math skills... in a desperate fight to calibrate his compass... it's a timeless tale of truimph and defeat...ahem.
Actually, this is mainly a story about soft-iron calibration.
I'm curious to see how accurate the magnetometer can be if properly calibrated since it'll be a primary source of heading reference for my AHRS due to the very short distances, high speeds, and poor GPS reception in the area of operation. So, rather anal proper calibration seems like the right path.
In attempting to calibrate my compass, plain hard iron (offset) and gain calibration didn't help as much as it has on other sensors in other configurations. I wondered why.
I've captured and saved data from the sensor using Hon Bo Xuan's 3D Scatter Processing script (download)
And am using Yury Petrov's Ellipsoid_fit (download) in Octave to attempt to figure out how to calibrate for hard and soft iron offsets. The hard iron offsets are ridiculously easy of course. It's the soft-iron stuff that's frying my feeble mind. :)
Both the ellipsoid fit and a gnuplot rendering of the XY and YZ plots show minimal 'tilt' of the ellipsoid about 2-5 degrees.
The XZ axis is another story. And yes, I've removed nearby ferrous objects to no avail.
Ellipsoid_fit says the eigenvectors (apparently these are the 3 axes of the ellipsoid) reported are:
Ex = -0.471 0.095 0.877
Ey = -0.865 0.145 -0.481
Ez = -0.173 -0.985 0.014
I'm still working out how to figure out the tilt of the x-z ellipse.
Another script, fit_ellipse by Ohad Gal (download) -- it's a 2d fit -- reports approximately -25 degree tilt (for some reason I had to correct x axis mirrored between the two plots; stretched out to resemble the x-z plot above) I don't think the tilt estimate is accurate.
The tilt estimation at least for this latter script is quite sensitive to axis gain. Hm.
I've been trying to wrap my head around the math involved in soft-iron calibration and found a few particularly good articles (several more that are, I think, difficult). I guess the main point is to have a soft-iron compensation matrix involved. Calculating the matrix eludes me as yet.
I've also read a paper or two on swinging a compass (as done in the airline industry) and I'm keeping that approach in my back pocket.
I'd happily ignore soft-iron as I've done in the past but even with offset and gain calibration the sensor is pretty clearly inaccurate at various orientations and I have convinced myself that it's due to soft iron distortion.
If you skip ahead to this reply, it talks about the approach used for soft-iron calibration and shows that, yes indeedy, it seems to be working at least on paper.
So what *should* I be seeing?
Suppose I create a series of samples that fell exactly on a sphere of radius 1. It simulates a perfect magnetometer with zero noise and magnitude 1. Turns out the eigenvectors exx/exy and eyx eyy are all over the place. exz, eyz and ez are all locked in. As I generate more and more points, the eigenvectors converge to ex=[0 0 1] ey=[0 1 0] ez=[1 0 0]. I need 10,000 or more points before it starts to look close which is utterly impractical for manual calibration. So, guess I cannot trust the estimated eigenvectors for this situation.
Next I created a 1000 simulated magnetometer samples with magnitude noise but 0 offset and 0 tilt and ran it through Ellipsoid_fit and got:
offset = [1.6196e-004 -2.5503e-003 -7.0838e-004]
radii = [ 1.00560 1.00101 0.99894 ]
I haven't quantified the relationship between offset/magnitude error and number of points and noise. For now I'll assume that with a few hundred points and typical sensor noise it'll be "close enough".
I wanted to look at magnitude values versus x, y, and z values respectively. With the perfect sphere I get a line from (-x,1) through (x,1) as one might expect. Same for y vs. magnitude and z vs. magnitude. In a perfect sphere, magnitude is 1. As you might imagine, a noisy magnitude produces a noisy plot, not a clean line.
I wonder what my real magnetometer's plots will look like? I did offset correction and scaling based on Ellipsoid_fit. I ran into some problems. The offset and radii were incorrect. The error was particularly bad for the z axis.
Under these conditions the x,y,z vs magnitude plots are a disaster so can't tell what's going on I don't think. More later.
After manually calibrating the real magnetometer for offset and gain as best as I could, I plotted magnitude in 3d.
Then I plotted magnitude versus x and y. Instead of a flat disc we'd get for a sphere, I got this Pringles chip shape.
I believe most of this distortion is due to the soft-iron effect. Thinking about a tilted ellipsoid, it makes sense that magnitude isn't constant.
Of course next is to figure out the best approach for correction. The papers above generally talk about a rotation and using a least squares approach to find the parameters for that rotation. Time to go learn and apply.
I decided to take a look at another magnetometer module: a standalone LSM303DLH. Let's see if it performs any better than the LSM303DLM that's on the MinIMU-9 board.
I collected data and plotted the result in 3d and it looked like it had less soft-iron distortion. I went through a semi-manual process of offset and magnitude calibration. Then plotted the magnitudes and wound up with good results that to me suggest very little soft-iron distorion:
From the Freescale app note "Since the eigenvectors of the matrix A represent the principal axes of the magnetometer measurement ellipsoid (the Principal Axis Theorem in geometry), then constraining the inverse soft-iron matrix W-1 to be symmetric shrinks the ellipsoid into a sphere along the principal axes of the ellipsoid without applying any additional spurious rotation"
Basically what this is saying as far as I am aware is that the soft-iron calibration involves scaling the ellipsoid along its (principal) axes, _not_ rotating as I had mistakenly thought, originally.
So, it looks like I have it.
Using ellipsoid_fit (download) mentioned above, we get back the principal axes of the best fit ellipsoid as well as their radii and the offset of the ellipsoid.
We can create a scaling vector based on the radii of these axes. If R is the vector containing the x, y, and z radii, then set up a scaling matrix in Octave:
[ 1/R(1) 0 0; 0 1/R(2) 0; 0 0 1/R(3) ]
The principal axis is basically just another frame of reference. Our scaling matrix is in that frame of reference. However the magnetometer readings are in the frame of reference of our x, y, and z axes.
So we can create, of all things, a direction cosine matrix that provides a means of converting the scaling matrix into the x, y, and z axes of our magnetometer readings, at which point we multiply those readings and we get what we were hoping for, a nearly perfect sphere of radius 1.
I finished up an Octave script to do all this for me and it then plots out before and after plots and magnitude plots to show it worked.
It's pretty clear that there's some tilted ellipsoid action happening. After running the script we get the following:
Plus, to prove it worked, here's the plot of X, Y, and Z versus magnitude. This is pretty darn close to the ideal (noisy) straight line I posted earlier.
Needless to say I'm quite pleased with the results at least on paper. Plus I learned quite a bit about matrix algebra, transformations, direction cosine matrices, and other goodies along the way. Good fun!
Though I'm confident in the paper results, I plan to test that the resulting compass corrections work on the real world robot and make sure that, yes indeed, the calibration results in sufficient magnetometer accuracy.
I have just started working on to calibrate my magnetometer. I used magnetometer for my heading correction using extended kalman filter. I see that the errors are coming to around 20 deg, and I am very unhappy with this large errors. I considered only the offsets/bias during my simulations. To improve my heading accuracy, I am looking info soft iron and scaling effects.
What is the heading accuracy you have obtained after the corrections? Is your approach working well?
Thanks & Regards,
Thank you for sharing your work. I have some hardware with extreme magnetic errors and also wrote some code using the linked ellipsoid fit in an attempt to make the magnetometer usable. I thought it would be worth sharing these results in this thread.
This first plot shows the uncalibrated data. The grey line is arbitrary rotations and the red, green and blue lines are rotations on a flat level surface around x, y and z (axis pointing down each time). You can see that this is an extreme situation, I found the magnetometer data was useless with hard-iron calibration only.
This next plot shows the measurements calibrated using the ellipsoid method. The ellipsoid fit has left the magnetometer axes at an incorrect orientation. The plot also shows additional RGB vectors drawn to the centroid of each x, y and z rotation dataset. These vectors represent a a non-orthogonal rotation matrix would use to correct for this orientation error.
This final plot shows the calibrated measurements corrected for the orientation error. The x, y and z rotation datasets are now concentric with the principle axes. The plot looks as if a high level of accuracy has been achieved. However, I implemented this within the AHRS hardware and found errors of up to 20 degrees. I wonder if the extreme initial magnetic distortions are too much for this calibration model.
For reference, the calibration parameters for the above data are:
1.2475 -0.0470 0.0429
0.0309 0.8271 -0.0291
-0.1009 -0.6458 0.7490
calibrated = softIronMatrix * uncalibrated - hardIronVector
Interesting, thank for posting this!
20° errors? Wow, that is awful. It seems my soft iron distortion mostly came from the office in which I was working , probably computer / monitor. Ooops. I seemed to get better results outside calibrating several feet off the ground. I haven't looked at this in awhile, having given up magnetometer for my rover few weeks after this due to distortion on the ground,
I seem to recall seeing some weirdness after calibration and thinking if I'd goofed something.
I just found a bug in my code. The above calibration now enables my AHRS device to achieve heading errors < 1 degree!
i'm also using Yury Petrov's Ellipsoid_fit as my magnetometer calibration procedure.
i collect 1000 raw values, and then launch the script.
then i use this formula to get calibrated values, and my magnetometer is calibrated
xt_raw = x_raw - offsetx;
yt_raw = y_raw - offsety;
zt_raw = z_raw - offsetz;
x_calibrated = scalefactor_x * xt_raw + scalefactor_x * yt_raw + scalefactor_x * zt_raw;
y_calibrated = scalefactor_y * xt_raw + scalefactor_y * yt_raw + scalefactor_y * zt_raw;
z_calibrated = scalefactor_z * xt_raw + scalefactor_z * yt_raw + scalefactor_z * zt_raw;
but my question are:
1) is this procedure right?
2) if i install my magnetometer in a device with metal near, can i calibrate again (maybe runtime on atmega) only the scale factor, or should i do a new complete calibration of my device?
I recently tried Yury Petrov's Ellipsoid_fit approach and the method described in this post which essentially implements an approach described in the Freescale documentation found here. Both algorithms do work, but the Freescale approach implementation requires fitting data to a model twice, so it seems like the Ellipsoid_fit approach described here is probably the best way to go.
You do need to be careful with the order of the radii since most math libraries return the eigenvalues and eigenvectors in ascending order. If you set up the scalars as described by Bot Thoughts:
"We can create a scaling vector based on the radii of these axes. If R is the vector containing the x, y, and z radii, then set up a scaling matrix in Octave:
[ 1/R(1) 0 0; 0 1/R(2) 0; 0 0 1/R(3) ]"
You will need to make sure that R(1), R(2) and R(3) are in the right order so they scale the correct radii of the ellipsoid. The order you need probably won't be ascending order.
If you apply the radii in the wrong order, you get something that looks like this image: the blue sphere is the ideal fit, the red ellipsoid is the offset input data, and the green ellipsoid is the output data that has been corrected for offset and skew. However, the radii scalar was applied to the green data in the wrong order and it is slightly mis-fit.
Once I had figured out how to solve the radii ordering problem, I got a really nice fit from my simulated data. The red sphere has radii of 1.4, 1.3, and 1.2 and an offset of 2,2,2 from the center 0,0,0. The blue sphere has a radius of 1 and a center of 0,0,0. The green sphere is the red ellipsoid after being scaled and centered. It has a radius of 1 and a center of 0,0,0, so it worked perfectly.
In the real world, I was working with an Android device and you can see that gathering enough points to generate a good data set can be a pain.
Assuming you spend some time getting a good data set to fit to the model, you can get a very accurate compass reading from a magnetic field that was unusable before calibration. However, it is very sensitive to changes in location, so it you move it even a little bit, you need to re-calibrate.
today i've post a python script to help magnetometer calibration, http://davidegironi.blogspot.it/2013/01/magnetometer-calibration-he...
it just helps collection of raw values.
the ellipsoit fit is Yury Petrov script, your java implementation can be usefull! :)