Accelerometer Calibration III: Improving Accuracy With Least-Squares and the Gauss-Newton Method

This is the third post in a series.  

  1. Introduction
  2. Simple Methods
  3. Least-Squares and Gauss Newton
  4. Streaming Gauss-Newton on an ATMEGA
  5. Error Analysis
  6. ?

In the last post we looked at some simple calibration methods and they didn’t seem to get us results that were better than the ones we got by just trusting the datasheet.  There were two obvious problems:

  1. Our calibration methods assumed we had placed our accelerometer in perfect axial alignment each time.  We were, um, less than perfect.
  2. We ignored the fact that there is noise in the sensor readings.  Eyeballing it, this alone would give us about ~1% error.

We could address problem (1) by finding ways to place our device very carefully before taking measurements. We could address problem (2) by averaging our samples (and assuming that the average was a good estimate of the real values).

Instead, we will use some Mathematics to develop a calibration method that is robust to noise and is not dependent on careful sensor placement. This will give more accurate results even though it requires less care by the user.

Least-Squares

Here is the idea.  We’ve been assuming that there are just six numbers we need to find to get a good calibration: m_x, m_y, m_z, \delta_x, \delta_y, \text{ and } \delta_z.  (See the last post for details and definitions.  Also note that if there is any skew we’ll need a more complicated model.)  Now matter how we place our sensor before we take readings, we assume it is still and is detecting exactly 1G of acceleration in some direction.

So if we read value x on the X-pin, y on the Y-pin, and z on the Z-pin we can look at the acceleration vector \mathbf{a} = (a_x, a_y, a_z) given by

a_x = \frac{x-m_x}{\delta_x} \quad a_y = \frac{y-m_y}{\delta_y} \quad a_z = \frac{z - m_z}{\delta_z}

If our measurements have no noise and our parameters are correct, the length of this vector should be exactly 1 no matter which way it is pointing.  That is

a_x^2 + a_y^2 + a_z^2 = 1

Or

\left(\frac{x-m_x}{\delta_x}\right)^2 + \left(\frac{y-m_y}{\delta_y}\right)^2 + \left(\frac{z-m_z}{\delta_z}\right)^2 = 1

To calibrate our model we just need to find parameters m_x, m_y, m_z, \delta_x, \delta_y, \text{ and } \delta_z so that this equation is always true. Unfortunately

  • With only one observation we get one equation in 6 unknowns.  There will be infinitely many solutions in general.  How do we pick one?
  • With exactly six distinct observations  we have 6 equations in 6 unknowns and  should get a unique solution.  The observations should be as different as possible so the solution isn’t error prone.  But we know there is noise in our readings to this won’t be exact.
  • With more than six observations we capture information about the noise in the data but we get more than 6 equations in 6 unknowns and cannot expect to find a solution.

It seems that the right thing to do is work with many observations.  We won’t be able to solve our equations, but we can choose parameters that make the errors as small as possible.

Let’s presume we used the circuit we built last time, collected data for the 6-point calibration, and ended up with N=100 readings.  Write x_i, y_i, and z_i for the reading on the X-pin, Y-pin, and Z-pin  in the i-th sample.  Then the error we get for sample i is

\epsilon_i = \left(\frac{x_i-m_x}{\delta_x}\right)^2 + \left(\frac{y_i-m_y}{\delta_y}\right)^2 + \left(\frac{z_i-m_z}{\delta_z}\right)^2 - 1

Now we want to find parameters so that all of the \epsilon_i are small.  We expect some errors so we don’t want to over-penalize these, but large errors suggest significant problems and should be penalized heavily.  It turns out that if we look at square errors, \epsilon_i^2 , they have just this sort of property and make the algorithms that follow simpler.   So we will try to pick parameters so that the sum of square errors is small:

\sum_{i=1}^{N} \epsilon_i^2 \text{ is as small as possible}

Why square errors? In statistics we always look at square errors first, even if it isn’t always right.  In this case I believe using square errors are both reasonable and convenient even if there isn’t a clear reason we shouldn’t use some other convex function of the error.  And you’ll see we get pretty good results.

But I digress.

Now we can formally state the problem we are trying to solve.

Least Squares Problem. Given samples (x_i, y_i, z_i) for i = 1, \dots, N, find parameters m_x, m_y, m_z, \delta_x, \delta_y, \delta_z such that the penalty function

\sum_{i=1}^N \epsilon_i^2 = \sum_{i=1}^{N}\left[ \left(\frac{x_i=m_x}{\delta_x} \right)^2 + \left(\frac{y_i-m_y}{\delta_y} \right)^2 + \left(\frac{z_i-m_z}{\delta_z} \right)^2 - 1\right]^2 

is as small as possible.

This is a classic nonlinear least-squares problem and it can be solved numerically using the Gauss-Newton method.  One day I may write more about these things — it is a beautiful algorithm and we’ll need to adjust it to work with the limited memory of an ATMEGA chip — but for today I will just offer you my code.

Performing Gauss-Newton to Solve the Least-squares Problem

I use the following workflow to calibrate an ADXL335 using the Gauss-Newton method.  I do things on windows 7, but it is all implemented with Python and Octave – free tools that are available many platforms.  Small adjustments may be needed.  (You can download Python here and download octave here.)

  1. Pick a directory to store data and code.  I’ll call this C:\your\calibration\directory\.
  2. Collect data using the pushbutton circuit described in the last post.  Put the sensor in a variety of positions (the 6 positions you used for 6-point calibration would be a good start).   Store the serial output in a text file.  I’ll assume it is stored in C:\your\calibration\directory\calData.txt.
  3. download the Python script octaveprep.py and place it in C:\your\calibration\directory\, then run it, passing the data file name as an argument.  At the windows command line (already in the right directory) I typeC:\Python27\python.exe octaveprep.py calData.txt

    This will write a new file, an octave script called calTest_octave.m.  This script will be used to perform the real data analysis.

  4.   Download the file gaussnewton.mand place it in  C:\your\calibration\directory\. Open octave and load the script by typingsource “C:\\your\\calibration\\directory\\calData_octave.m”

    This will compute a 6-number array beta which holds all of our calibration parameters.  The correspondence is as follows:

    m_x = beta(1)
    m_y = beta(2)
    m_z = beta(3)
    1/\delta_x = beta(4)
    1/\delta_y = beta(5)
    1/\delta_z = beta(6)

When I run this workflow on the data I collected for the 6-point calibration in the last post I get the following:
beta =  5.1448e+002  5.0231e+002  5.1665e+002  9.5485e-003  9.4631e-003  9.6936e-003

which means that m_x = 514.48, m_y = 502.31, m_z = 516.65 and \delta_x = 1/0.0095485 = 104.73, \delta_y = 105.67, \delta_z = 103.16.  Sticking these parameters in the unit conversion loop I used to evaluate the 6-point calibration method in the last post, I get readings

  X-up     0.9980144   0.0254738   0.0615826
  X-down  -0.9976190   0.0254738   0.0131148
  Y-up     0.0145204   1.0001717  -0.0062723
  Y-down  -0.0332220  -0.9965395   0.0131148
  Z-up    -0.0618675   0.0538631   0.9921641
  Z-down   0.0336174  -0.0123785  -1.0047088

We did not get a bunch of beautiful readings like “0 0 1”, “0 0 -1”, etc. because we still had sloppy placement going in.  But if you compute the magnitudes of each of these vectors you will see that they range from 0.99555 to 1.00535.  They are almost all within one half of one percent of the perfect 1G reading!

If you have trouble getting these scripts to work for you please leave a comment and let me know.

Now we have a calibration that uses exactly the same measurements that our 6-point calibration system did but that has very accurate results.  If we automated the workflow described above, we’d have an easy accurate calibration method.

Better Accuracy Through Better Samples

This is starting to look pretty good, but we can still do better. For the last example I just used data collected from six basic positions: Z-up, Z-down, Y-up, Y-down, X-up, and X-down.  We can get better results if we sample different directions, but wich new directions should we add?

You can think of these positions we already sampled as the ones we would get if we stuck our accelerometer on each of the faces of a cube.  We want to get new directions that are “as different as possible”. To do this, just chop off each of the corners of the cube to get 8 new faces and set the accelerometer on each of these new faces too.

I don’t really use a truncated cube to take measurements.  Actually I tape my accelerometer to a ruler and clamp it in a third hand from my soldering bench to position it securely.

After 30 minutes of fiddling, I got a stream of slightly noisy data from each of these 14 positions. I calibrate and get the following parameters:

beta =  5.1409e+002  5.0231e+002  5.1662e+002  9.5632e-003  9.4572e-003  9.6777e-003

These are very close to the parameters we got from the 6-position data.  In fact, on those 6 positions these parameters do not give better results (they can’t because we had chosen the best parameters for those 6 positions).  But for the new positions and for positions scattered around the sphere we see significantly better accuracy.  I’ll show you the details when we talk about error analysis later.

Six-point Calibration Overall Grades:  Accuracy – Great Easiness – OK

What’s Next

We could go through this process every time we get a new sensor, store the parameters in EEPROM, and carry on. But I want to be able to calibrate a device on the slopes or on the lift and I won’t be toting around my laptop and Octave installation.  If you looked through some of that code you would see that we were dealing with some very large matrices, ones that would quickly overrun the memory on an Arduino. In the next post we’ll see how to overcome these problems and implement the Gauss-Newton method directly on the Arduino.

Advertisements
This entry was posted in Electronics and tagged , , , , , , . Bookmark the permalink.

18 Responses to Accelerometer Calibration III: Improving Accuracy With Least-Squares and the Gauss-Newton Method

  1. Great tutorial. I’ll be trying to use your method shortly. Looking forward to the next parts.

  2. kipod says:

    Good and interesting Post.
    One question: When tilting the accelerometer to 45 deg (for example) in any axis, does the total reading in 3 axis is still 1?

    • If the accelerometer is laying “flat” so that we get a reading of 1 in the z-axis and 0 in the x and y axes then any rotation about the z-axis leaves the reading unchanged. In other words we can spin it flat on the table. However, rotations about the x axis will change both y and z axis readings. Rotations about the y axis change both x and z axis readings.

      Hopefully I’ll get the next installment out soon, or at least just post the code. My time has all disappeared futzing with packaging devices for wireless outdoor use — which should be fun when it’s done — but it can get frustrating.

      • kipod says:

        Rotation around Z axis will result exactly what you describe. however when I rotate my accelerometer around X or Y, I get g values of ~0.6g for each axis.
        The total g from all axis is greater then 1 and I did not figure it out yet.
        Did you had similar results after your calibrations?

  3. kipod says:

    OK, the total acceleration is sqrt(Ax^2 + Ay^2 + Az^2).
    Another reason to get the calibration right 🙂

    • Yep, here are the readings I’m getting. They don’t up perfectly on the axes because I soldered the thing in crookedly. I guess I’m trying to use math to compensate for my lousy dexterity 🙂

      Anyway, laying flat on the table I get
      x= -0.0368 y= 0.16039 z= 0.99891
      Magnitude = \sqrt{x^2 + y^2 + z^2} = 1.01237

      Tilted on the side:
      x= -0.0276 y= -0.99397 z= 0.12783
      Magnitude = \sqrt{x^2 + y^2 + z^2} = 1.000254

      And about halfway between the two:
      x= -0.0368 y= -0.71930 z= 0.73312
      Magnitude = \sqrt{x^2 + y^2 + z^2} = 1.027722

      This was with done with Gauss-Newton but only using six data points for input. Taking 14-20 readings to calibrate gets better accuracy (especially on the “tilted” readings). Still, this is good enough for a few of my applications.

  4. Pingback: Connecting to Sparkfun’s 9DOF “Sensor Stick”: I2C access to ADXL345, ITG-3200, and HMC5843 | Chionotech

  5. nick says:

    why Ax^2 + Ay^2 + Az^2 = 1 ????

    Btw excellent posts…

    • Daniel says:

      I’m going to go ahead and answer this.

      Applying a sensitivity scaling (deltas) transforms the raw accelerometer output to g (gravitational acceleration) units. If no external acceleration is applied to the accelerometer (i.e. when it is still), it will only detect the gravitational acceleration, which measured in g’s is 1.

      The formula is merely Phytagoras theorem squared in 3D linear space. To get the resulting acceleration for all three directions, you would apply the square root to it. Since sqrt(1) = 1, the 1 remains intact.

  6. Pingback: SensorLib: Using Calibration | Chionotech

  7. Pingback: Implementing the Gauss-Newton Algorithm for Sphere Fitting (2 of 3) | Chionotech

  8. mini says:

    Hi , Gr8!!!! tutorial for a starter .
    I’m new to python , where the filename has to be replaced in the code …
    passing as argument , means???

  9. Arpit Kumar says:

    Hi,

    I am using a MPU6050 IMU which is not an analog based sensor but uses I2C communication to get the IMU data. Using the datasheet based values for bias and scale, I got the accelerometer calibrated but it is not good enough. So, I tried to use your code to make the calibration of my accelerometer more accurate. And on running the python and octave scripts successfully got these beta values:
    5.1921e+02 5.1921e+02 5.1921e+02 1.1173e-03 1.1131e-03 1.1057e-03

    Obviously these are not giving me the correct calibrated acceleration values. I see you have used some constants [512,512,512,0.0095,0.0095,0.0095] somehow in your code. Probably this offset is caused due to that (My sensor was already giving data approximately within -1 to +1 but the values were not perfectly correct). As per my understanding for finding the bias and scale values for an acclerometer it should be independent of whether my sensor is analog or digital. Can you please help me what’s wrong I am doing here.

    Thanks a lot.

    Arpit

  10. I also met the same trouble with Arpit Kumar, here is the beta results I get from running python and octave code:

    5.1922e+002 5.1922e+002 5.1922e+002 1.1130e-003 1.1138e-003 1.1092e-003

    Compare to yours I realize that my beta(4), (5) and (6) are smaller for about 10 times. Could you please explain this for me?

  11. trungoc says:

    Hi,

    I am also using ADXL335 with Arduino UNO and I also need to calibrate my data. I have followed your instruction from the start to here, but when I run your Python and Octave code, I meet the same trouble with Arpit Kumar. My beta values I have got:

    5.1922e+002 5.1922e+002 5.1922e+002 1.1130e-003 1.1138e-003 1.1092e-003

    I realized that my beta(4), (5) and (6) values are much smaller than yours (about 10 times) and cause my acceleration values smaller than it should be for about 10 times (for example, instead of
    0.10 0.10 1.00 now they turn to 0.01 0.01 0.10). Could you please explain me about this difference between your result and mine? Even if I use the result of 6-point calibration you posted in previous post, the beta values are still not the same with ones you posted here.

    Thanks a lot,

    trungoc

  12. Cam says:

    If you fit 6 measurements to a 6 parameter ellipsoid function, shouldn’t the unique solution also have zero error? Why aren’t the mapped values all identically 1.0 norm? I am wondering if this error comes from the iterative gauss-newton. It is possible to fit an ellipsoid to minimize the functional error analytically, and I believe this method should have zero error fitting 6 points.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s