Arduino, Wire, and I2C part 5b: Data Analysis and gnuplot

I've been using a feature of gnuplot lately that until yesterday I hadn't really given much thought to its potential. Gnuplot has a curve fitting tool which can be used to come up with those parameters to the sinusoidal functions in the previous post. It's not magic in that it won't necessarily always give you the right answer every time, but it can help with cleaning up estimates. Here's an example. First, define the function you're trying to fit to. In this case, I'm trying to fit to a simple sinusoidal function:

gnuplot> f(x)=a2*sin(omega*x+phi)+d2

This is the function described in the previous post, in gnuplot syntax. Now, you can try and fit the data to the function. The total data set isn't a consistent wave form so I'm going to reduce the data set to points between 150 and 250.

gnuplot> fit [150:250] f(x) "spin.dat" using 0:7 via a2,omega,phi,d2

The above instructs gnuplot to execute an iterative process, attempting to derive the values for a2, omega, phi, and d2 in the function defined as f(x). Gnuplot will output the results in various stages of the process. The most critical results are the final ones of course, which I reprint here:

Final set of parameters            Asymptotic Standard Error
=======================            ==========================

a2              = 0.774445         +/- 11.02        (1423%)
omega           = 0.973404         +/- 0.4104       (42.16%)
phi             = 6.881            +/- 82.27        (1196%)
d2              = -66.2102         +/- 6.535        (9.871%)

You may notice two things:

1. The equation parameters are quite different from those that I posted previously, and
2. The standard error percentage of the parameters is quite high (even the lowest, 9% is pretty significant).

You can "prime" the fit algorithm by seeding the parameter to reasonable approximations, i.e. give the fit a good starting point (there's probably a mathematical term for this but I'm not a mathematician). Using the parameters guessed at in the last post:

gnuplot> a2=93.5
gnuplot> d2=-75.6
gnuplot> omega=1/4.6
gnuplot> phi=-6
gnuplot> fit [150:250] f(x) "spin.dat" using 0:7 via a2,omega,phi,d2

Which ultimately resulted in:

Final set of parameters            Asymptotic Standard Error
=======================            ==========================

a2              = 90.8588          +/- 1.735        (1.909%)
omega           = 0.220797         +/- 0.0006131    (0.2777%)
phi             = -8.0615          +/- 0.1239       (1.537%)
d2              = -74.2817         +/- 1.185        (1.596%)

gnuplot> plot "spin.dat" using 7 axes x1y1 t "mag y" with lines, a2*sin(omega*x+phi)+d2

Which looks like a fairly good fit, with less than 2% error. A better fit could be made using equations with more independent variables, however, gnuplot's curve fit feature is limited to 5, of which the above already was using 4. I don't think there is much chance of refining this estimation any further using the fit feature of gnuplot.

Arduino, Wire, and I2C part 5: Data Analysis and Modeling

Now that I've managed to reduce noise levels in the sensors to a mostly manageable level, it's time to record some data and try and model the board. It's worth reiterating what it is I'm attempting to achieve here: a hand-held underwater metal detector. Magnetometers can be a great way to achieve this, but this particular device is measuring the strength of the magnetic field.

I'm fairly confident that linear movement of the magnetometer is not going to significantly change the sensor measurements. This is because Earth's magnetic "lines of force" don't change particularly rapidly. You would probably have to move hundreds of miles to be able to see any significant change in the magnetometer's measurements, and can you really say you've moved in a "straight line" by then? More than likely you'll have been moving in an arc relative to gravitational center of the earth, which would be even less likely to show a significant change in measurement. Therefore, I'm going to ignore ADXL measurements.

The drawing above shows the board and the axes of measurement for each device. Unfortunately, sparkfun mounted the magnetometer such that it is 90° out of phase with the other devices. This means that the positive X axis of the accelerometer and gyro is the positive Y axis of the magnetometer, and the positive Y axis of the accelerometer and gyro is the negative X axis of the magnetometer. An inconvenience, but hopefully just a minor one.

Data capture was achieved using a common spice rack. I wanted something that didn't have much, if any, metal in it, while having sufficient space to throw the entire test apparatus on. The sensor board is roughly centered (it's a little hard to see given that the tape holding it down is a similar color to the spice rack, but you can see the connector to it above the battery pack), with the power pack and the Arduino UNO, XBee shield and XBee around the edges. Data was transferred via the XBee radio to a PC running Linux. Data loss made it necessary to make the data processing software a bit more robust.

Data capture was successful, and the graph to the right is a plot of the gyro measurements. As you can see in the graph, the rotation was primarily occurring around the Z axis, which matches the diagram of the axes and the picture of the apparatus. What you can't really see very well in the graph is that the Y and X axes were also showing measurable rotation. Those measurements show the same oscillation that can be seen on the Z-axis measurements. This is always going to be the case in any real-world measurements. You're simply never going to be able to get all the axes completely lined up and the axis of rotation perfectly aligned with the Z-axis (or whichever axis you're picking as "up"). You might get a lot closer than what I achieved with some really expensive lab equipment. The Z rotation is negative, indicating that the rotation was in a counter-clockwise direction around the Z axis.

Before looking at the magnetometer data, it's time to take a quick excursion back (for me at least) to Euclid and his geometry, though for the sake of simplicity, I'm going to stick with a 2-dimensional geometry for now. When translating between polar coordinates (angle and radius) and Euclidian coordinates, the following formulae apply:

$\begin{array}{rcl}x_{1} & = & r \cos \alpha,\\y_{1} & = & r \sin \alpha.\end{array}$

where r is the radius and α is the angle. In Euclidian geometry, your axes are always orthogonal, meaning they're at angles of 90° to each other. As such, the trigonometric operations sine and cosine are also 90° out of phase with each other. This phase offset shows up clearly in the recorded data. Given that the experiment was intended to rotate the board at a reasonably constant rate around the Z axis, the X and Y magnetometer measurements will look like sinusoidal waves 90° out of phase with each other. The following two graphs demonstrate. The first graph is the actual data, while the second graph is a simulation of the data set using only cos and sin functions, which are scaled to the appropriate mean and amplitude.

Recorded Measurements

Simulated with Trig functions

The "simulation" in the second plot is of the following functions in standard form:

$\begin{array}{rcl}y_{1}(t) & = & A_{1}\cdot \cos (\omega t + \varphi ) + D_{1},\\y_{2}(t) & = & A_{2}\cdot \sin (\omega t + \varphi ) + D_{2},\\A_{1} & = & 82,\\D_{1} & = & 15.73,\\A_{2} & = & 93.5,\\D_{2} & = & -75.6,\\\omega & = & \frac{1}{4.6},\\\varphi & = & -6.\end{array}$

Note that the frequency (ω) and phase (φ) are the same for both equations. Only the amplitude (A) and center amplitude (D) are different.

One important thing to note (in fact, pretty much the whole point of this post) is that ω is almost certainly a function of the Z-measurement of the gyro. In the simplified environment of a perfect rotation about a co-axial Z axis for the two sensors, this might be expressed as:

$\omega = \gamma z_{g}$

where γ is some constant. The reason for this is that the gyro is measuring some Δα, that is, the rotation rate around a given axis. This is, in fact, the frequency ω. It only needs some scale factor applied to it to match the measurements to the model. The value of that constant scale factor can probably be derived from the data, but I probably need to take a few more sample runs at various speeds of rotation before I can feel comfortable quantifying it.

More to come...

Arduino, Wire, and I2C Part 4: Noise and analysis

These weekend I spent some time recording sample data sets from the sparkfun SEN-10724 and looking at it. I'd say I was analyzing the data, except that might give the impression that I actually have a good idea what I'm doing.

The recording of data, I actually set up using the sparkfun "ethernet pro" board, which is basically an Arduino with built-in ethernet and seemingly poorly designed voltage regulation (it gets pretty hot if you power it with the intended voltage levels). I felt more comfortable using the ethernet rather than serial-over-USB since I would be sending more data than was reasonable for the 9600bps that the Arduino UNO is fixed at (over USB - the firmware for the atmega8u2 that handles the USB-to-serial interface is programmed such that it will only operate at 9600bps).

After finding a few errors with decoding the data (mostly in the PC side, the microcontroller was programmed to send the raw measurements as-is), I found that the ADXL345 accelerometer has noise levels way out of spec. I've plotted the measurements and attached images of said plots at the end of this post. The XYZ measurements of the ADXL345 all had significant levels of noise (RMSD 5-12 LSBs) and the magnetometer had noise on the X-axis only (RMSD 9 LSBs). During the test data recording, the sensor board was in a breadboard sitting on my desk with no significant sources of vibration.

One concern I had is that the schematic for the SEN-10724 had two .1μF capacitors "near" the ADXL345. The datasheet calls for a single .1μF capacitor between ground and V_DDI/O, and for a 1μF tantalum capacitor at V_S, with an optional 10μF tantalum capacitor in parallel.

If the circuit instead only has a .1μF capacitor at V_S, that might explain why the measurements look so noisy. In-circuit measurement of the capacitors is impossible, and I don't have a stockpile of SMD capacitors (strangely enough, after the 3 separate LED cube builds) to replace it with. Time to email customer support.

Arduino, Wire, and I2C

In previous postings, I described connecting a sparkfun.com sensor board to the Arduino and getting magnetometer readings from it.  In this post, I intend to go into a little more detail about how I²C works, while providing similar functionality using the accelerometer on the same board.

The first thing to understand is that I²C is a completely 8-bit bus protocol, which means that all addresses and all data are 8-bit quantities (0-255).  The bus consists of a single "master" with up to 112 unique "slave" devices.  The limit of 112 is due to the fact that a device address is actually a 7-bit quantity (more on this later), and 16 addresses are reserved.  Multiple masters are apparently possible, but aren't relevant for this project.  Slave devices can't talk to other slave devices.  Off-the-shelf devices, like the three on the SEN-10724 board, have pre-assigned addresses managed by a central authority.

Devices have a 7-bit address.  This 7-bit address is stored in the 7 MSB (most significant bits) of an 8-bit quantity when sequenced on the bus, with the least significant bit indicating whether the operation is a read or a write (1 or set being "read").  As an example, the HMC5883L datasheet lists 8-bit addresses 0x3D for read, and 0x3C for write.  Convert these to binary and you'll find a 7-bit address of 0x1E followed by a binary 1 for read or 0 for write.  Some manufacturers will list only a 7-bit address in their data sheet, some will list 8-bit addresses, some will do both.

The datasheet for the Analog Devices ADXL345 lists both - 0x1D is the 7-bit address, which results in 0x3A for writes and 0x3B for reads.  However, it also supports an alternate addressing of 0x53/0xA6/0xA7 if pin 12 is tied to ground. The schematic for the sparkfun board indicates that pin 12 is indeed tied to ground, so let's assume that we're going to be using alternate addressing for that device.

Each I²C device will also have a number of registers that are readable and/or writable.  Table 16 on page 14 of the ADXL345 datasheet lists this device's register map.  For this project, we're most interested in registers 0x32-0x37.  If you look at the table, you'll see how manufacturers design for 16-bit quantities, which is (at least in this case and in the case of the HMC5883L) to have adjacent 8-bit registers containing the most significant and least significant bytes of a 16-bit quantity.

Let's start with the X axis on the accelerometer, which according to the board's silkscreen is the side-to-side direction (along the short edge of the board).  The following Arduino sketch implements a simple program that prints out the X axis acceleration once a second.

#include <Wire.h>
// Wire library uses 7-bit addresses and automatically
// sets the r/w bit
static const uint8_t ADXL345 = 0x53;
static const uint8_t DATAX0 = 0x32;

void setup()
{
  Wire.begin();
  Serial.begin(9600);
}

void loop()
{
  // Put the read address on the bus
  Wire.beginTransmission(ADXL345);
  Wire.send(DATAX0);
  Wire.endTransmission();
  // Request data starting with the X register
  Wire.beginTransmission(ADXL345);
  // 2 bytes gets us the LSB and MSB of the X data
  Wire.requestFrom(ADXL345, (uint8_t)2);
  // store the data in a signed integer quantity
  int16_t x;
  // pointer to use to store the data
  byte *p = (byte*)&x;
  // wait for 2 bytes to be available
  while (Wire.available() < 2) {}
  *p = Wire.receive();
  p++; // advance the pointer to the next 8 bits
  *p = Wire.receive();
  Wire.endTransmission();
  Serial.println(x); // finally, print the quantity
  delay(1000); // delay a bit to avoid flooding
}

If you upload this sketch to your Arduino, with the IMU connected as described in the earlier post, you'll probably get a sequence of zeroes in the serial monitor. If, for example, you tried to change the address in the above code to the primary (0x1D), you probably wouldn't see anything at all. So why is it only showing zeroes, even if I shake the heck out of the board? My guess is that the device doesn't start up in a mode where continuous measurements are taken (which is how the HMC5883L operates), and that to get these measurements, other registers must be manipulated...

To be continued...

Magnetometer Metal Detection

Part of what motivated me to start playing with the magnetometer was an episode of Deep Sea Detectives I had recently watched on Netflix.  I believe it was the "Damn the Torpedoes" episode, where they used a magnetometer to find a ship that presumably had been buried by silt in a hurricane.  The magnetometer used in DSD was a cesium magnetometer, which is more sensitive, but probably a lot more expensive, even if it were a simple matter for an ordinary person to get their hands on cesium.  The HMC5883L that I'm using has a resolution of 2 milli Gauss, or 200nT (nanotesla).  A cesium magnetometer is probably several orders of magnitude more sensitive.

The plot shows "detection" of metal.  In this case, I moved the remains of my old internal Zip drive (which I had disassembled after viewing a "maker" youtube video along those lines) between the sensor and ground (as in the stuff beneath your feet, not the electrical sense in this case).  The off-scale measurements in the chart are probably an implementation error in the HMC5883L library mentioned in my previous post.  The maximum range of the measurements is -2047 - 2048, so I suspect that the values are being incorrectly converted from the 12-bit measurements to a 16-bit quantity.

The following list, which I unabashedly stole from geometrics.com, should give you an idea of how that rates in a real situation where a marine magnetometer is employed:

Object	flux density	measurement range
Ship 1000 tons	0.5 to 1 nT	800 ft (244 m)
Anchor 20 tons	0.8 to 1.25 nT	400 ft (120 m)
Automobile	1 to 2 nT	100 ft (30 m)
Light Aircraft	0.5 to 2 nT	40 ft (12 m)
Pipeline (12 inch)	1 to 2 nT	200 ft (60 m)
Pipeline (6 inch)	1 to 2 nT	100 ft (30 m )
100 kg of iron	1 to 2 nT	50 ft (15 m)
100 lbs of iron	0.5 to 1 nT	30 ft (9 m)
10 lbs of iron	0.5 to 1 nT	20 ft (6 m)
1 lb of iron	0.5 to 1 nT	10 ft (3 m)
Screwdriver 5 inch	0.5 to 2 nT	12 ft (4 m)
1000 lb bomb	1 to 5 nT	100 ft (30 m)
500 lb bomb	0.5 to 5 nT	50 ft (16 m )
Grenade	0.5 to 2 nT	10 ft (3 m )
20 mm shell	0.5 to 2 nT	5 ft (1.8 m)

Of course, their product is trying to do a high-level survey, where the magnetometer is tethered to a ship and is well into the water column.  In my application, I intend to have a hand-held device that I'm moving along the sea/lake bed, so I'm not going to be trying to detect a screwdriver from 4 meters, I'm going to be trying to detect a screwdriver from, maybe 1m.

I have no idea how much one would have to pay for the geometrics magnetometer, but I figure it's one of those situations where "if you have to ask...".

There is still plenty to do to even approach a usable metal detector.  Unlike traditional inductive detectors, magnetometers are measuring minute changes in the magnetic field flux density (no, really, I'm not trying to go all technobabble).  Inductive detectors can detect non-ferrous metals, but have a much lower sensitivity and therefore lower range than a magnetometer.  Obviously, then, the magnetometer is restricted to ferrous (iron-based) metals but can have a significantly higher range of detection.  The other issue, and probably the one more significant, is that any rotation of the sensor along any axis will show up in the data.  The trick will be to distinguish between movement of the detector and actual metal.

"9" degrees of freedom IMU on Arduino

This past weekend I finally got around to playing with the "9 degrees of freedom" sensor board (SEN-10724) I'd purchased from sparkfun.com.  I'm not sure if it makes sense to call it that, but what you get is a 3-axis gyro, a 3-axis accelerometer and a 3-axis magnetometer.  The gyros give you information about rotational acceleration, the accelerometer measures linear acceleration, and the magnetometer works like a 3-axis compass.

The circuit off to the side there shows how I managed to get it hooked up and talking through my Arduino UNO.  The short of it is this: it's a 3.3V device (well, 3 devices) and it needs pull-up resistors (I used 4.7K) on the data and clock lines in order to function.  It's wired up to the analog inputs 4 and 5 on the Arduino board, which is what the provided "Wire" library uses to talk to devices like this, that use the I2C protocol.

Honestly, I'm a little bit unclear as to how it's actually able to talk to the sensor stick/IMU (inertial measurement unit) without doing level conversion between the 5V Arduino and the 3.3V IMU, but it does seem to work reliably in this configuration.  If I were doing something more significant (a production board, for example) I'd probably be a bit more careful about matching the signal levels.

The Wire library for Arduino is just a basic library for talking to I2C devices.  Getting into the specific interfaces is another matter, though for a quick start, I used the HMC588L compass library provided by Love Electronics in the UK.  This was enough to get me started and verify that I was able to communicate with the magnetometer.  There's a pretty decent tutorial on that page on how to use the library, though the sample works pretty well as-is.  If you have trouble compiling the sample code, I found that for some reason it has a period (".") at the very beginning of the file - remove that period to make it compile.

Quick update - looking at the example .pde file in linux using hd (hex dump), it turned out there were three unprintable characters at the beginning of the file.  The easiest way I found to fix the problem was to remove any odd characters before the /* at the start of the file, then save it as a new sketch.  That new .pde file can then be used to replace the original example .pde file in the library, or you can just use the fixed version in your sketchbook.