Magnetometer coil equalization

Question

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.

Answer 1

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.

Answer 2

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')
solution = pickle.load(f)
f.close()
print(ccode(solution[c0], c0))
print(ccode(solution[c1], c1))
print(ccode(1 - c0 - c1, c2))
print(ccode(solution[d0], d0))
print(ccode(solution[d1], d1))
print(ccode(solution[d2], d2))

It may be better to solve the system iteratively, which is something to investigate.