# Magnetometer coil equalization

I have a search coil magnetometer with two coils, one at each end of a ferromagnetic core of high magnetic permeability: Figure 1. Double search coil magnetometer. (Modified from a drawing CC BY-SA 3.0 Christophe coillot.)

Each coil converts change in the magnetic field into a voltage signal. The two signals are digitized. I'm only interested in local anomalies in the magnetic field as they appear in the difference signal between the two coils when the magnetometer is moved through the anomalies. The anomalies are in Earth's magnetic field, created by ferromagnetic objects buried in the ground. So this is a kind of metal detector.

The problem is that the coils are not aligned perfectly symmetrically and are not built perfectly identically. The difference signal is therefore polluted by imperfectly cancelled noise. The noise originates primarily from accidental shaking of the magnetometer that changes its orientation with respect to Earth's magnetic field.

I can do calibration by going to nature and lifting up the magnetometer in the air, away from stray magnetic fields and magnetic anomalies, and shaking it, producing pure noise that ideally would be identical at both coils, but is not, as described above. During such calibration, I have tried adaptively changing the gain of one of the coil signals to match that of the other, which works, but only to an extent. I think I would get better results by matching the frequency responses as well. How can I do this? Figure 2. Possible signal flow diagrams with equalization of A) just one of the channels to match the other, B) both channels to match a common target response.

I don't need to flatten the frequency responses, as long as they are equal, because I'm not trying to properly quantify the anomalies, just detect them. As far as I understand, the equivalent circuit of each search coil is: Figure 3. LTspice equivalent circuit of a search coil. Voltage from the voltage source is proportional to magnetic flux time derivative. The component values are very approximate, but probably close enough to reveal the general shape of the frequency response. Figure 4. LTspice simulated frequency response of the search coil.

Gain must remain as one of the variables to equalize, because some of the gain difference comes from the digitization. I can set the sampling frequency as high as needed. I'd rather not do the calibration in a test bench with a generated artificial magnetic field, but adaptively in the field.

• "The noise originates primarily from accidental shaking of the magnetometer that changes its orientation with respect to Earth's magnetic field." Is this due to mechanical movement or gradients in the earth's magnetic field? Also, as I read the application, you are sweeping the search coil back and forth and looking for peaks in output as the wand is moved. Is the wand hand operated or on a cart? What are typical signal levels? – rrogers Jun 14 '17 at 19:12
• @rrogers Mechanical movement. It's hand-operated. Also the detector frequently bumps into vegetation. I'm not sure about the signal levels. – Olli Niemitalo Jun 14 '17 at 19:30
• It would be most helpful if you monitored the voltages produced by both coils in a controlled situation: say place a piece of rebar on the floor and record the raw data for 5 sweeps. And then remove the rebar and try to do a similar operation. Post! Since it's mechanical the raw data is important. Not knowing any more I would suggest using the sum of the two voltages as a gain control (roughly because there is a problem of division by zero. And looking for a null in difference. Not a real null but midway between positive and negative extremes. – rrogers Jun 14 '17 at 20:29
• can't you use an LMS or NLMS to adapt one coil to the other? with no difference stimuli, then the error signal ($e[n]$) should go to zero. at some point, the adaptation gain should get very low so that when there is a legit difference signal, you will see it and it won't be adapted out. – robert bristow-johnson Jun 15 '17 at 16:15
• concept is pretty simple, but it's FIR, not IIR. if you want IIR, i think it's called RLS filter, but i have forgotten exactly how that one works. – robert bristow-johnson Jun 15 '17 at 17:44

I would equalize both channels. The most robust method is to cut both back to a lower bandwidth. Then you only need to equalize the gains.
Assuming that your scanning by hand a 1-2 second sweep rate seems reasonable. Unless your imaging, you can narrow the sweep down to smaller distances when a signal needs investigation.
A 10-50Hz bandwidth should accommodate human response times. In this case you can use simple "RC" filtering to a fixed standard on both samples.
You mention construction variations; so first, you need to buffer yourself away from the coils. Then apply a fixed RC filter. Since you have isolated yourself both in time and impedance from the coil variations you have complete control of the balance between the two sides.
Of course, this might all be done in software but make sure your A/D has sufficient resolution (time and amplitude) and compliance.

I've made some progress with plan B: Figure 1. Plan B: both coils equalized.

Each Equalizer is a 3-tap FIR filter, which should be sufficient to cancel the apparent two poles of the coil (see Fig. 4 of the question). The impulse response of Equalizer 0 is $[c_2, c_1, c_0]$, and the impulse response of Equalizer 1 is $[d_2, d_1, d_0]$ Equalizer 0 is constrained to have a 0 Hz response of 1 by $c_0 + c_1 + c_2 = 1$. I think it is a reasonable normalization that avoids setting all coefficients to zero. Coil 0 samples are denoted by $x[k]$ and Coil 1 samples by $y[k].$ The equalized Coil 0 samples are denoted by $x_\text{eq}[k] = c_0x[k] + c_1x[k+1] + c_2x[k+2],$ and the equalized Coil 1 samples are denoted by $y_\text{eq}[k] = d_0y[k] + d_1y[k+1] + d_2y[k+2].$ The Difference between the equalized samples is $x_\text{eq}[k] - y_\text{eq}[k].$ For pure noise input we want to find the coefficients that minimize the sum of squared differences, equivalent to minimizing the mean square power of the difference signal. Sum of squared differences is:

$$\sum_{k=0}^{N-1}\left(x_\text{eq}[k] - y_\text{eq}[k]\right)^2\\ = \sum_{k=0}^{N-1}\big(c_0x[k] + c_1x[k+1] + (1 - c_0 - c_1) x[k+2] - d_0y[k] - d_1y[k+1] - d_2y[k+2]\big)^2\\ = c_0^2\sum_{k=0}^{N-1}x^2[k] + c_0c_1\sum_{k=0}^{N-1}x[k]x[k+1]+\\ \ldots\\ +d_1d_2\sum_{k=0}^{N-1}y[k+1]y[k+2] - d_2^2\sum_{k=0}^{N-1}y^2[k+2]$$

Because the complete expanded equation is very long, the above shows only the first two and the last two terms. Minimization of the sum of squared differences is achieved by finding the set of coefficient values for which the partial derivatives of the sum of squared differences, with respect to each coefficient, is zero. This gives a system of equations to solve, written in code for easy copy-pasting:

c0*(x0x0 - 2*x0x2 + x2x2) + c1*(x0x1 - x1x2 - x0x2 + x2x2) - d0*(x0y0 - x2y0) - d1*(x0y1 - x2y1) - d2*(x0y2 - x2y2) + x0x2 - x2x2 = 0
c0*(x2x2 - x0x2 + x0x1 - x1x2) + c1*(x1x1 - 2*x1x2 + x2x2) - d0*(x1y0 - x2y0) - d1*(x1y1 - x2y1) - d2*(x1y2 - x2y2) + x1x2 - x2x2 = 0
- c0*(x0y0 - x2y0) - c1*(x1y0 - x2y0) + d0*y0y0 + d1*y0y1 + d2*y0y2 - x2y0 = 0
- c0*(x0y1 - x2y1) - c1*(x1y1 - x2y1) + d0*y0y1 + d1*y1y1 + d2*y1y2 - x2y1 = 0
- c0*(x0y2 - x2y2) - c1*(x1y2 - x2y2) + d0*y0y2 + d1*y1y2 + d2*y2y2 - x2y2 = 0


The left side of each equation in the system is half the partial derivative of the sum of squares of differences with respect to coefficients $c_0, c_1, d_0, d_1, d_2,$ respectively. The 4-letter variables represent dot products between time-shifted signals. For example x0y2 with $i=0$ and $j=2$ represents dot product $\sum_{k=0}^{N-1}x[k+i]y[k+j]$. The dot products are calculated from the noise training signal.

Unfortunately the exact solution is a set of formulas for the coefficients that are each a lengthy degree 5 rational function of the dot products. This gives problems with numerical precision when evaluating them using floating point arithmetic. The exact solution can be found by this Python script that takes about an hour to execute:

from sympy import solve, symbols, cse
from sympy.printing import ccode
import pickle
c0, c1, c2, d0, d1, d2 = symbols('c0,c1,c2,d0,d1,d2')
x0x0, x0x1, x0x2, x0y0, x0y1, x0y2, x1x1, x1x2, x1y0, x1y1, x1y2, x2x2, x2y0, x2y1, x2y2, y0y0, y0y1, y0y2, y1y1, y1y2, y2y2 = symbols('x0x0,x0x1,x0x2,x0y0,x0y1,x0y2,x1x1,x1x2,x1y0,x1y1,x1y2,x2x2,x2y0,x2y1,x2y2,y0y0,y0y1,y0y2,y1y1,y1y2,y2y2')
solution = solve([c0*(x0x0 - 2*x0x2 + x2x2) + c1*(x0x1 - x1x2 - x0x2 + x2x2) - d0*(x0y0 - x2y0) - d1*(x0y1 - x2y1) - d2*(x0y2 - x2y2) + x0x2 - x2x2, c0*(x2x2 - x0x2 + x0x1 - x1x2) + c1*(x1x1 - 2*x1x2 + x2x2) - d0*(x1y0 - x2y0) - d1*(x1y1 - x2y1) - d2*(x1y2 - x2y2) + x1x2 - x2x2, - c0*(x0y0 - x2y0) - c1*(x1y0 - x2y0) + d0*y0y0 + d1*y0y1 + d2*y0y2 - x2y0, - c0*(x0y1 - x2y1) - c1*(x1y1 - x2y1) + d0*y0y1 + d1*y1y1 + d2*y1y2 - x2y1, - c0*(x0y2 - x2y2) - c1*(x1y2 - x2y2) + d0*y0y2 + d1*y1y2 + d2*y2y2 - x2y2], c0, c1, d0, d1, d2)
f = open('solution.txt', 'w')
pickle.dump(solution, f)
f.close()


After solving, the solution can be loaded from solution.txt and emitted as C code by:

from sympy import symbols
from sympy.printing import ccode
import pickle
c0, c1, c2, d0, d1, d2 = symbols('c0,c1,c2,d0,d1,d2')
x0x0, x0x1, x0x2, x0y0, x0y1, x0y2, x1x1, x1x2, x1y0, x1y1, x1y2, x2x2, x2y0, x2y1, x2y2, y0y0, y0y1, y0y2, y1y1, y1y2, y2y2 = symbols('x0x0,x0x1,x0x2,x0y0,x0y1,x0y2,x1x1,x1x2,x1y0,x1y1,x1y2,x2x2,x2y0,x2y1,x2y2,y0y0,y0y1,y0y2,y1y1,y1y2,y2y2')
f = open('solution.txt', 'r')