Background calibration for magnetometers, both gain and offset errors are removed (hard/soft iron).

Dependents:   SML2 optWingforHAPS_Eigen hexaTest_Eigen

Magnetometers are notoriously bad when uncalibrated. There are both gain and offset errors where the offset errors. Gain errors can often be removed by the magnetometer itself, many got options to apply a bias magnetic field of known strength, which can be measured and used to calculate gain errors. However the most significant errors are the offset errors, they can be huge. The offset on an axis can be so large that it will always measure the Earth magnetic field in one direction, which pretty much ruins your day if you want to calculate where the north is.

Just measuring this offset (output its measurement data to your PC, rotate it a bit, see what the offset is) already helps alot. However I consider automatic calibration algorithms nice, how practical they are is a whole different story.

Algorithms

Before I made this I looked a bit around online regarding automatic calibration algorithms for magnetometers. The vast majority are besides very computational intensive also really confusing to me, way too much vector algebra. Then we got those who get the offset simply by averaging their measurements, where each measurement gets a weight factor depending on the amount of rotation measured by the gyroscope. The idea is that if you have for example a quadrocopter it will generally point same amount of time in each direction, and if you use the gyroscope data it will prevent the offset calculations from going horribly wrong when it is pointed one direction for a long time. I am sure it works, but I do not consider it a really neat solution.

Another option is measuring the 'extremes' reported by each axis. If the magnetometer has been in use for a while, and on the x-axis one extreme is -80, and the other one is 120, the offset is obviously 20. And this way you can also immediatly compensate for gain errors. However there is one serious shortcome of this method: if due to whatever reason you make an incorrect measurement of for example 200, your offset and gain estimate is wrong. You can of course lower the impact of this by taking the average of the old value and new value, so one incorrect measurement can only have a limitted impact. However the issue remains that your extremes can only become more 'extreme'. There is no way to use new measurements to lower the value of your calculated extreme value.

Used algorithm

And that is where this algorithm comes into play. In principle it is based on calculating the extreme values of the magnetometer. If a new extreme is found, it will change the old one with a certain gain. However it adds one extra step: it takes into account that when two axes are close to their zero value (compensated for offset), the third one must be at an extreme value. It will then use this data to update the value of its stored extremes: if the measured value is larger than its current offset it will update its positive extreme value, if it is smaller it will update its negative extreme value. That immediatly brings us to the limitation of the algorithm, when the current estimate of the offset is so bad it is outside the actual values the magnetometer can reach, it may not be able to find the actual extremes.

2D example

I hope that was a little bit clear, but it is probably easier to visualize with a 2D example. In the following picture the red circle shows the values our magnetometric sensor can achieve if the world was 2D, with in the center the dot representing its offset. The blue lines connect its current values for its extremes (which are obviously wrong), and the green dot is a new measurement.

/media/uploads/Sissors/mag1.png

Since the y-axis is zero for this measurement, it will know this is a new extreme value. The measured value is smaller than its offset estimate (the point where the lines of the extremes cross), so it is a new negative extreme for the x-axis. So it will moves its negative extreme value a bit to this new measurement (how fast depends on the gain settings), but lets assume it got a whole bunch of measurements on this spot, so its negative extreme is now at this measurement location.

The next figure gives the resulting situation, and after a while the magnetometer gets to a spot where its x-axis value is zero, and it gets a new positive extreme measurement for the y-axis.

/media/uploads/Sissors/mag2.png

The same thing happens again, and it adjusts its y-axis extremes, and also gets a new offset value for this axis. Making the new result:

/media/uploads/Sissors/mag3.png

Now it only needs some measurements to get the gain on its x-axis correct, and it is fully calibrated. When the offset changes for whatever reason, it will automatically use its new measurements to get new estimates of its extremes, and use that to calculate its gain and offset. Besides this algorithm it also directly updates the extreme when it literally finds a new extreme value (still with small steps). So if in the previous example it immediatly gets to the point completely at the leftmost point of the circle it will use this to update the minimum extreme of the x-axis, despite that the y-axis data is not on its offset estimate.

In this example the initial extremes were too small, and an algorithm that only checks if the current measurement value is larger than its stored extreme value would also work fine. However the difference is that this also works if the initial extremes were too large.

Usage of the library

The library does not use any knowledge about how large the signal is supposed to be, it does not matter if they got a size of +/- 0.001, or +/- 10,000. This results in a problem since it really has no clue about the initial estimates it should start with, and for this reason the first 100 measurement points it only stores the extreme values found, and its output is identical to its input. For good results you want to rotate the sensor alot during these initial points, but really it is only useful when running it the first time to get a global estimate of the extremes.

For normal usage it is best to use the setExtremes function in the beginning to supply the filter with its starting points. When setExtremes is called it will disable those first 100 measurement points, and will function normally from the start. It is probably ideal to at the end of your program let the mbed write the latest extreme values of the filter (getExtremes) to a non-volatile memory, such as the localFileSystem. Then at the start of your program you get these values again, so it will always use the most accurate initial conditions.

The output of the filter is the magnetic vector, compensated for gain and offset, normalized to a length of one. Generally the length of the magnetic vector is not used for anything, so this should not be a problem. Do not this normalization is not absolute: it just uses its gain and offset estimates to normalize to a length of roughly one (in practise a bit more than one), but it does not actually normalize the vector to have exactly a length of one.

Because the filter only uses fairly simple calculations: multiplications, additions and comparisons mostly, and not computational intensive stuff such a sine and square root calculations, it is pretty lightweight. I think execution takes around 10us, but I will update that value when I check that again.

Committer:
Sissors
Date:
Wed Jul 11 12:26:43 2012 +0000
Revision:
0:dc85e939b642
Child:
1:8620dbd88cd4
[mbed] converted /Quadrocopter/CalibrateMagneto

Who changed what in which revision?

UserRevisionLine numberNew contents of line
Sissors 0:dc85e939b642 1 //TODO: ABLE TO DISABLE CALIBRATION UPDATE MORE EFFICIENT
Sissors 0:dc85e939b642 2 //TODO: CHANGE AVERAGE DUE TO NEW EXTREME FOUND
Sissors 0:dc85e939b642 3 #include "CalibrateMagneto.h"
Sissors 0:dc85e939b642 4
Sissors 0:dc85e939b642 5 CalibrateMagneto::CalibrateMagneto( void ) {
Sissors 0:dc85e939b642 6 measurementNumber = 0;
Sissors 0:dc85e939b642 7 this->setAbsGain(0.2);
Sissors 0:dc85e939b642 8 this->setRelGain(0.1);
Sissors 0:dc85e939b642 9 };
Sissors 0:dc85e939b642 10
Sissors 0:dc85e939b642 11
Sissors 0:dc85e939b642 12 void CalibrateMagneto::setAbsGain( float gain ) {
Sissors 0:dc85e939b642 13 if (gain>=0)
Sissors 0:dc85e939b642 14 abs_gain=gain;
Sissors 0:dc85e939b642 15 }
Sissors 0:dc85e939b642 16
Sissors 0:dc85e939b642 17 void CalibrateMagneto::setRelGain( float gain ) {
Sissors 0:dc85e939b642 18 if (gain>=0 && gain<=1)
Sissors 0:dc85e939b642 19 rel_gain=gain;
Sissors 0:dc85e939b642 20 }
Sissors 0:dc85e939b642 21
Sissors 0:dc85e939b642 22 float CalibrateMagneto::getAbsGain( void ) {
Sissors 0:dc85e939b642 23 return abs_gain;
Sissors 0:dc85e939b642 24 }
Sissors 0:dc85e939b642 25
Sissors 0:dc85e939b642 26 float CalibrateMagneto::getRelGain( void ) {
Sissors 0:dc85e939b642 27 return rel_gain;
Sissors 0:dc85e939b642 28 }
Sissors 0:dc85e939b642 29
Sissors 0:dc85e939b642 30 void CalibrateMagneto::setExtremes(float *min, float *max) {
Sissors 0:dc85e939b642 31 x_min=min[0];
Sissors 0:dc85e939b642 32 y_min=min[1];
Sissors 0:dc85e939b642 33 z_min=min[2];
Sissors 0:dc85e939b642 34 x_max=max[0];
Sissors 0:dc85e939b642 35 y_max=max[1];
Sissors 0:dc85e939b642 36 z_max=max[2];
Sissors 0:dc85e939b642 37
Sissors 0:dc85e939b642 38 //Disable initial points
Sissors 0:dc85e939b642 39 measurementNumber = INITIAL_POINTS+1;
Sissors 0:dc85e939b642 40 }
Sissors 0:dc85e939b642 41
Sissors 0:dc85e939b642 42 void CalibrateMagneto::getExtremes(float *min, float *max) {
Sissors 0:dc85e939b642 43 min[0]=x_min;
Sissors 0:dc85e939b642 44 min[1]=y_min;
Sissors 0:dc85e939b642 45 min[2]=z_min;
Sissors 0:dc85e939b642 46 max[0]=x_max;
Sissors 0:dc85e939b642 47 max[1]=y_max;
Sissors 0:dc85e939b642 48 max[2]=z_max;
Sissors 0:dc85e939b642 49
Sissors 0:dc85e939b642 50 }
Sissors 0:dc85e939b642 51
Sissors 0:dc85e939b642 52 void CalibrateMagneto::run(float *input, float *output) {
Sissors 0:dc85e939b642 53 float x_avg, y_avg, z_avg, x_range, y_range, z_range;
Sissors 0:dc85e939b642 54 bool x_zero, y_zero, z_zero;
Sissors 0:dc85e939b642 55 float temp;
Sissors 0:dc85e939b642 56
Sissors 0:dc85e939b642 57
Sissors 0:dc85e939b642 58 if (measurementNumber==0) { //Initial measurement just presets everything
Sissors 0:dc85e939b642 59 x_min=input[0];
Sissors 0:dc85e939b642 60 x_max=input[0];
Sissors 0:dc85e939b642 61 y_min=input[1];
Sissors 0:dc85e939b642 62 y_max=input[1];
Sissors 0:dc85e939b642 63 z_min=input[2];
Sissors 0:dc85e939b642 64 z_max=input[2];
Sissors 0:dc85e939b642 65 measurementNumber++;
Sissors 0:dc85e939b642 66 output[0] = input[0];
Sissors 0:dc85e939b642 67 output[1] = input[1];
Sissors 0:dc85e939b642 68 output[2] = input[2];
Sissors 0:dc85e939b642 69 } else {
Sissors 0:dc85e939b642 70 x_avg = (x_min+x_max)/2;
Sissors 0:dc85e939b642 71 y_avg = (y_min+y_max)/2;
Sissors 0:dc85e939b642 72 z_avg = (z_min+z_max)/2;
Sissors 0:dc85e939b642 73
Sissors 0:dc85e939b642 74
Sissors 0:dc85e939b642 75 if (measurementNumber < INITIAL_POINTS) { //No abs gain used, no rel gain used, no range used
Sissors 0:dc85e939b642 76 if (input[0]<x_min)
Sissors 0:dc85e939b642 77 x_min = input[0];
Sissors 0:dc85e939b642 78 else if (input[0]>x_max)
Sissors 0:dc85e939b642 79 x_max = input[0];
Sissors 0:dc85e939b642 80
Sissors 0:dc85e939b642 81 if (input[1]<y_min)
Sissors 0:dc85e939b642 82 y_min = input[1];
Sissors 0:dc85e939b642 83 else if (input[1]>y_max)
Sissors 0:dc85e939b642 84 y_max = input[1];
Sissors 0:dc85e939b642 85
Sissors 0:dc85e939b642 86 if (input[2]<z_min)
Sissors 0:dc85e939b642 87 z_min = input[2];
Sissors 0:dc85e939b642 88 else if (input[2]>z_max)
Sissors 0:dc85e939b642 89 z_max = input[2];
Sissors 0:dc85e939b642 90
Sissors 0:dc85e939b642 91 //Return inputs, since filter is still bad
Sissors 0:dc85e939b642 92 output[0] = input[0];
Sissors 0:dc85e939b642 93 output[1] = input[1];
Sissors 0:dc85e939b642 94 output[2] = input[2];
Sissors 0:dc85e939b642 95
Sissors 0:dc85e939b642 96 measurementNumber++;
Sissors 0:dc85e939b642 97 } else { //Now we should have enough data
Sissors 0:dc85e939b642 98 x_range=x_max-x_min;
Sissors 0:dc85e939b642 99 y_range=y_max-y_min;
Sissors 0:dc85e939b642 100 z_range=z_max-z_min;
Sissors 0:dc85e939b642 101
Sissors 0:dc85e939b642 102 //If measurement is a new extreme:
Sissors 0:dc85e939b642 103 if (input[0]<x_min) {
Sissors 0:dc85e939b642 104 temp = rel_gain*(x_min - input[0]);
Sissors 0:dc85e939b642 105 if (temp > abs_gain*(x_max-x_min))
Sissors 0:dc85e939b642 106 temp = abs_gain*(x_max-x_min);
Sissors 0:dc85e939b642 107 x_min = x_min - temp;
Sissors 0:dc85e939b642 108 } else if (input[0]>x_max) {
Sissors 0:dc85e939b642 109 temp = rel_gain*(input[0] - x_max);
Sissors 0:dc85e939b642 110 if (temp > abs_gain*(x_max-x_min))
Sissors 0:dc85e939b642 111 temp = abs_gain*(x_max-x_min);
Sissors 0:dc85e939b642 112 x_max = x_max + temp;
Sissors 0:dc85e939b642 113 }
Sissors 0:dc85e939b642 114
Sissors 0:dc85e939b642 115 if (input[1]<y_min) {
Sissors 0:dc85e939b642 116 temp = rel_gain*(y_min - input[1]);
Sissors 0:dc85e939b642 117 if (temp > abs_gain*(y_max-y_min))
Sissors 0:dc85e939b642 118 temp = abs_gain*(y_max-y_min);
Sissors 0:dc85e939b642 119 y_min = y_min - temp;
Sissors 0:dc85e939b642 120 } else if (input[1]>y_max) {
Sissors 0:dc85e939b642 121 temp = rel_gain*(input[1]-y_max);
Sissors 0:dc85e939b642 122 if (temp > abs_gain*(y_max-y_min))
Sissors 0:dc85e939b642 123 temp = abs_gain*(y_max-y_min);
Sissors 0:dc85e939b642 124 y_max = y_max + temp;
Sissors 0:dc85e939b642 125 }
Sissors 0:dc85e939b642 126
Sissors 0:dc85e939b642 127 if (input[2]<z_min) {
Sissors 0:dc85e939b642 128 temp = rel_gain*(z_min - input[2]);
Sissors 0:dc85e939b642 129 if (temp > abs_gain*(z_max-z_min))
Sissors 0:dc85e939b642 130 temp = abs_gain*(z_max-z_min);
Sissors 0:dc85e939b642 131 z_min = z_min - temp;
Sissors 0:dc85e939b642 132 } else if (input[2]>z_max) {
Sissors 0:dc85e939b642 133 temp = rel_gain*(input[2]-z_max );
Sissors 0:dc85e939b642 134 if (temp > abs_gain*(z_max-z_min))
Sissors 0:dc85e939b642 135 temp = abs_gain*(z_max-z_min);
Sissors 0:dc85e939b642 136 z_max = z_max + temp;
Sissors 0:dc85e939b642 137 }
Sissors 0:dc85e939b642 138
Sissors 0:dc85e939b642 139 //If measurement is near the zero of the other two axes
Sissors 0:dc85e939b642 140 x_zero=false;
Sissors 0:dc85e939b642 141 y_zero=false;
Sissors 0:dc85e939b642 142 z_zero=false;
Sissors 0:dc85e939b642 143 if (abs( input[0] - x_avg ) < (x_range * ZERO_RANGE))
Sissors 0:dc85e939b642 144 x_zero=true;
Sissors 0:dc85e939b642 145 if (abs( input[1] - y_avg ) < (y_range * ZERO_RANGE))
Sissors 0:dc85e939b642 146 y_zero=true;
Sissors 0:dc85e939b642 147 if (abs( input[2] - z_avg ) < (z_range * ZERO_RANGE))
Sissors 0:dc85e939b642 148 z_zero=true;
Sissors 0:dc85e939b642 149
Sissors 0:dc85e939b642 150
Sissors 0:dc85e939b642 151 if (x_zero && y_zero) {
Sissors 0:dc85e939b642 152 if (input[2]>z_avg) {
Sissors 0:dc85e939b642 153 temp = rel_gain*(input[2] - z_max);
Sissors 0:dc85e939b642 154 if (abs(temp) > abs_gain*(z_max-z_min))
Sissors 0:dc85e939b642 155 temp = sign(temp)*abs_gain*(z_max-z_min);
Sissors 0:dc85e939b642 156 z_max = z_max + temp;
Sissors 0:dc85e939b642 157 }
Sissors 0:dc85e939b642 158 if (input[2]<z_avg) {
Sissors 0:dc85e939b642 159 temp = rel_gain*(z_min - input[2]);
Sissors 0:dc85e939b642 160 if (abs(temp) > abs_gain*(z_max-z_min))
Sissors 0:dc85e939b642 161 temp = sign(temp)*abs_gain*(z_max-z_min);
Sissors 0:dc85e939b642 162 z_min = z_min - temp;
Sissors 0:dc85e939b642 163 }
Sissors 0:dc85e939b642 164 }
Sissors 0:dc85e939b642 165
Sissors 0:dc85e939b642 166 if (x_zero && z_zero) {
Sissors 0:dc85e939b642 167 if (input[1]>y_avg) {
Sissors 0:dc85e939b642 168 temp = rel_gain*(input[1] - y_max);
Sissors 0:dc85e939b642 169 if (abs(temp) > abs_gain*(y_max-y_min))
Sissors 0:dc85e939b642 170 temp = sign(temp)*abs_gain*(y_max-y_min);
Sissors 0:dc85e939b642 171 y_max = y_max + temp;
Sissors 0:dc85e939b642 172 }
Sissors 0:dc85e939b642 173 if (input[1]<y_avg) {
Sissors 0:dc85e939b642 174 temp = rel_gain*(y_min - input[1]);
Sissors 0:dc85e939b642 175 if (abs(temp) > abs_gain*(y_max-y_min))
Sissors 0:dc85e939b642 176 temp = sign(temp)*abs_gain*(y_max-y_min);
Sissors 0:dc85e939b642 177 y_min = y_min - temp;
Sissors 0:dc85e939b642 178 }
Sissors 0:dc85e939b642 179 }
Sissors 0:dc85e939b642 180
Sissors 0:dc85e939b642 181 if (y_zero && z_zero) {
Sissors 0:dc85e939b642 182 if (input[0]>x_avg) {
Sissors 0:dc85e939b642 183 temp = rel_gain*(input[0] - x_max);
Sissors 0:dc85e939b642 184 if (abs(temp) > abs_gain*(x_max-x_min))
Sissors 0:dc85e939b642 185 temp = sign(temp)*abs_gain*(x_max-x_min);
Sissors 0:dc85e939b642 186 x_max = x_max + temp;
Sissors 0:dc85e939b642 187 }
Sissors 0:dc85e939b642 188 if (input[0]<x_avg) {
Sissors 0:dc85e939b642 189 temp = rel_gain*(x_min - input[0]);
Sissors 0:dc85e939b642 190 if (abs(temp) > abs_gain*(x_max-x_min))
Sissors 0:dc85e939b642 191 temp = sign(temp)*abs_gain*(x_max-x_min);
Sissors 0:dc85e939b642 192 x_min = x_min - temp;
Sissors 0:dc85e939b642 193 }
Sissors 0:dc85e939b642 194 }
Sissors 0:dc85e939b642 195
Sissors 0:dc85e939b642 196 //And now the actual filter part:
Sissors 0:dc85e939b642 197 output[0] = 2* (input[0] - x_avg)/x_range;
Sissors 0:dc85e939b642 198 output[1] = 2* (input[1] - y_avg)/y_range;
Sissors 0:dc85e939b642 199 output[2] = 2* (input[2] - z_avg)/z_range;
Sissors 0:dc85e939b642 200
Sissors 0:dc85e939b642 201
Sissors 0:dc85e939b642 202 }
Sissors 0:dc85e939b642 203 }
Sissors 0:dc85e939b642 204 }
Sissors 0:dc85e939b642 205
Sissors 0:dc85e939b642 206
Sissors 0:dc85e939b642 207
Sissors 0:dc85e939b642 208 float CalibrateMagneto::sign(float input) {
Sissors 0:dc85e939b642 209 if (input<0)
Sissors 0:dc85e939b642 210 return -1.0;
Sissors 0:dc85e939b642 211 else
Sissors 0:dc85e939b642 212 return 1.0;
Sissors 0:dc85e939b642 213 }