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:
Mon Jul 16 18:09:30 2012 +0000
Revision:
1:8620dbd88cd4
Parent:
0:dc85e939b642
v1.0

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