DSP in commercial prospecting and treasure hunting metal detectors ... and GPZ 7000 speculation

[Davor]

Mr Crocker deserves a bit more praise. I like the 3rd (and 5th) removal mechanism. It supports square wave Tx, which means you may design an incredibly simple and effective energy saving Tx that operates at just about any frequency you throw at it, without tuning anywhere. It may enable a completely brain-free and frequency independent implementation of dual frequency method, AND keep the energy saving goodies too.

Now GB.

Ground signal is impractical to extract from the received signal using sumation methods simply because the amplitude of GB is changing with differences in ground proximity while waving around. Therefore in VLF we extract it by means of a GB channel that is orthogonal against the ground response vector. Such a channel that does not contain ground response is used to generate audio. There is not much to it, you simply apply some muffling mechanism to avoid chatters, and do some tone generation, either proportional audio or VCO.
If you only use a GB channel, you get a TR detector that detects all metals, and no ground. Such configuration is the most sensitive VLF scheme.
Tuning GB to properly reject ground is done by bobbing the coil, and thus maximising the amplitude variation of ground response. When GB is correct, there is no response while bobbing.
By applying discrimination you lose a bit of sensitivity.

Practical implementation...
A GB channel in a classic VLF is done by sampling demodulation of input signal at 90° from ground response. You tune it to precisely 90° by adjusting phase of a clock supplied to a sampling demodulator. Resulting voltage has 0V response for any amplitude of ground.
In micros you operate on vector components I and Q, and these are demodulated baseband signals. Each of the components has its own amplitude, and a vector they represent is seen in a Cartesian diagram. Obtaining any phase in between is performed by simple addition/subtraction of the I and Q components.

It goes like this...

Say you have a vector V with amplitude r and angle α
- in Cartesian coordinates it would stand as x= r cos(α), and y= r sin(α)

Say we want to rotate a whole vector by angle β
- we obtain x' and y' as follows: x'= x cos(β) -y sin(β) and y'= x sin(β) + y cos(β)

To obtain ground balance in such constellation we'll use either x' or y' (not both) that has minimum value (preferably 0) while bobbing a coil. To reach at that angle we may either establish the angle β directly using atan2 method as Mr Crocker, or more precisely by successive approximation (PLL/Costas loop style) where noise is less of a factor, and you may substitute sines and cosines with observing zero crossings. The best way of doing this is in a semi-automatic fashion: press button and bob -> find solution -> release button. You may try automating this as well in case you have lots of memory at disposal. It should work reasonably well on most terrains.

There is also a completely different approach using multiple frequencies, and so far amateurs did not tread in that direction. GB is achieved by subtracting (scaled) responses at two (or more) frequencies, and the catch is that practically all ground effects respond equally at different frequencies, while targets respond differently. After subtraction you get only a difference of target responses in a GB channel.

[Qiaozhi]

I thought that ground balance was simply shifting the sample point one the received curve relative to the transmitted curve. In which case you would need a good deal more than 8 values. I would think you would need at least 360 samples, so you could adjust the ground balance to a precise degree, or am I just smoking something

To include ground balance, you would need to take a sample at (or near) the zero-crossing point(s) As far as I can see, Mr. Crocker is taking the first sample at a later time. As you are familiar with the Arduino code, maybe you can determine at what time (relative to the TX signal) that the first sample occurs. There does not appear to be a mechanism whereby the detector user can adjust this sample delay, so it's pretty evident that this functionality has not been included, or that he's even aware of the requirement for GB on a VLF.

In an analog detector (and also in many digital detectors, which do not have direct sampling) the synchronous demodulation is done in hardware. In this case you have a defined sample pulse-width, and the value is stored in a capacitor. With direct sampling, you may be able to get away with a single sample, but subsequent samples would need to be accumulated and integrated over time to remove unwanted noise, as you would normally do in hardware. Probably the best way would be to simply try it. At least, with software, it's relatively simple to test various ideas without having to resort to the soldering iron.

[scrungy_doolittle]

It is actually quite easy. Suppose you have a vector at any arbitrary angle, and you have its Cartesian values x and y. "atan2" returns the angle in radians of such vector with given x and y values. In case your metal detector samples two channels at precisely 90° apart, the signal values in those channels, or bins if you like, are those x and y values. In analogue VLF rigs you have a disadvantage of not having a real atan2 function at hand, so you perform tricks to discriminate, but in a micro, or VDI, you perform atan2.
Bin 0 and 2 are at 90° degrees apart and there you have it.
VDI is some fancy acronym for a gizmo that tells you a degree of a target at hand. Actually atan2 function and a display, nothing more.
Mr. Crocker's device lacks ground balance, and it is not merely a 90 degree thing. People make mistake when they confuse all metal reference with ground balance. Just ask if you want to know more about these.

Dump on me please. I'm all ears to know more about this, as I want to create a open source metal detector based on a micro. While I am a bit math challenged, I have over 30 years experience in the embedded world.. I know about VDI, but Mr Crocker did not know the acronym. So lets talk about how to fit ground balance into code similar to this. And, I would ask, that if anyone creates any code based on these ideas that they post it, either on github or other open source site so all may see it. I take my hat off to Mr. Crocker for having done so.

[scrungy_doolittle]

OK ... so he's actually taking 4 samples for the positive half-cycle, and then repeating the process for the negative half-cycle. Not as I said previously - taking 8 samples per cycle, with 2 samples per quadrant. Which is completely different. Of course, in that case it's possible to accurately extract the amplitude and phase from the RX waveform. But ... there is still no mechanism for ground balancing.

I thought that ground balance was simply shifting the sample point one the received curve relative to the transmitted curve. In which case you would need a good deal more than 8 values. I would think you would need at least 360 samples, so you could adjust the ground balance to a precise degree, or am I just smoking something

[Qiaozhi]

OK ... so he's actually taking 4 samples for the positive half-cycle, and then repeating the process for the negative half-cycle. Not as I said previously - taking 8 samples per cycle, with 2 samples per quadrant. Which is completely different. Of course, in that case it's possible to accurately extract the amplitude and phase from the RX waveform. But ... there is still no mechanism for ground balancing.

[Davor]

It is actually quite easy. Suppose you have a vector at any arbitrary angle, and you have its Cartesian values x and y. "atan2" returns the angle in radians of such vector with given x and y values. In case your metal detector samples two channels at precisely 90° apart, the signal values in those channels, or bins if you like, are those x and y values. In analogue VLF rigs you have a disadvantage of not having a real atan2 function at hand, so you perform tricks to discriminate, but in a micro, or VDI, you perform atan2.
Bin 0 and 2 are at 90° degrees apart and there you have it.
VDI is some fancy acronym for a gizmo that tells you a degree of a target at hand. Actually atan2 function and a display, nothing more.
Mr. Crocker's device lacks ground balance, and it is not merely a 90 degree thing. People make mistake when they confuse all metal reference with ground balance. Just ask if you want to know more about these.

[scrungy_doolittle]

Having looked at the code some more, I think I now understand how this is supposed to work.

Firstly, an array of 4 bins are declared, which are used to accumulate the ADC readings. One cycle is split into 4 equal parts (90 degree quadrants) and are referred to as "phases" (0 to 3). The ADC capture is performed in the ISR (interrupt service routine). Two samples are taken per phase, making a total of 8 readings per cycle. After 1024 sets of readings have occurred, the contents of the bins are deposited into an "averages" array.

Then 4 different "bins" are declared, and the mysterious addition and subtraction process takes place. Initially I assumed that each of the 2 samples per phase were positioned around the zero-crossings and the positive and negative peaks, as this would make most sense. In that case, we would have the following:

phase0 = positive-going zero-crossing -> bins[0] -> averages[0] -> bin0
phase1 = positive peak -> bins[1] -> averages[1] -> bin1
phase2 = negative-going zero-crossing -> bins[2] -> averages[2] -> bin2
phase3 = negative peak -> bins[3] -> averages[3] -> bin3

With the ADC samples positioned around the zero-crossing and the peaks, we can examine the following code snippet:

double bin0 = (averages[0] + halfRoot2 * (averages[1] - averages[3]))/f;
where: halfRoot2 = 0.7071 and f = 200.0

In this instance, the calculation becomes:

bin0 = ((postive-peak - negative-peak) * 0.7071) + positive-going-zero-crossing

f appears to be an empirically defined fudge factor.

Assuming the peaks cancel (not sure why he's bothering to find the RMS value at this stage), bin0 will mainly contain the reading at the positive-going zero-crossing. Likewise, you can analyze the other three phases in the same way, and get the following:

bin0 = positive-going zero-crossing
bin1 = positive peak
bin2 = negative-going zero-crossing
bin3 = negative peak

Next he calculates the absolute amplitude at the zero-crossings:
double amp1 = sqrt((bin0 * bin0) + (bin2 * bin2));

And the absolute amplitude at the peaks:
double amp2 = sqrt((bin1 * bin1) + (bin3 * bin3));

So far so good.

At this point he has all the ingredients to achieve both ground balance and discrimination, together with VDI.
But then things start to go a bit awry.

Instead of calculating the inverse tangent of bin1/bin 0 and the same for bin3/bin2. In other words, calculating the in-phase and quadrature components of the RX signal, he calculates the inverse tangent of bin0/bin2 and likewise bin1/bin3, which really doesn't make any sense at all.

IMHO this shows a basic lack of understanding of how a VLF detector works. In fact, he states:" So we'll say that anything with a phase-shift below -20 deg is non-ferrous". Which indicates that the detector is only designed to give a ferrous/non-ferrous indication.

As Davor correctly said: "..... there are a few details that are missing here, namely ground balance".
Which is a shame, because everything is there to make it happen.

In conclusion, I suspect my initial assumption (that the samples are being taken around the zero-crossings and peaks) is incorrect, and the samples are actually being taken in the middle of each quadrant (phase). In this case, any anomalies in the signal (due to harmonic contamination) will be reduced, but the ability to ground balance and discriminate with VDI will be lost.

According to some comments made concerning this project, it is long defunct, and apparently never made it off the bench.

One thing he pointed out is the additions and subtractions will result in the bins being 0 when no signal is present.
So assume that no signal is 100. a signal comes in with 103 and 97. This will result in the bins having 6's in them.


The adc is biased to approx 1.6 v. I have posted Mr. Crockers reply later on, but here is his explaination of the bins.
This is the biggest question. I want to understand *exactly* what is
going on and why.
I understand the sqrt(bin0^2+bin2^) but why bin 0 and 2 and 1 and 3?
I said earlier that the software was implementing a phase sensitive
detector. More precisely, it is implementing two quadrature phase
sensitive detectors with reference signals 45 degrees apart. Each
quadrature detector itself comprises two phase sensitive detectors with
reference signals 90 degrees apart. The first quadrature detector
produces results in bins 0 and 2, the second produces results in bins 1
and 3.

Consider just one quadrature detector. We know the component of the
received signal at the reference phase, and the component at 90 degrees
to that phase. Assuming the signal is sinusoidal, to find the amplitude
of the actual signal, we need the square root of the sum of the squares
of the two components.

At this point, we could take the output of just one of the quadrature
detectors, and calculate the amplitude and phase from it. But has we
have 2 detectors (which is a consequence of sampling 8 times per cycle
instead of just 4, which is what allowed us to cancel the third
harmonic), it makes sense to use both outputs. So I calculate the signal
amplitude from both phase sensitive detectors, and average them.

I also don't understand how the phase adjustment is working.

The next question concerns:
double phase1 = atan2(bin0, bin2) * radiansToDegrees + 45.0;
double phase2 = atan2(bin1, bin3) * radiansToDegrees;

I know that the arctan will give the phase angle of a vector, but
again why bin 0 and 2 and why the addition of +45 to the resulting
degree measurement.
Bin 0 and 2 because they give the signal amplitudes from samples taken
90 degrees apart. Likewise bins 1 and 3 give the signal amplitudes from
samples taken 90 degrees apart, however those samples are taken 45
degrees later than bins 0 and 2. Hence the +45 so that both phase
results are referenced to the same point in the cycle.
Is this measurement essentially the same as the VDI that takes the X
and Y voltage components, and computes the phase angle? (there is an
arduino VDI posted on the Geotech forums).
Sorry, I don't recognise the acronym VDI.

[scrungy_doolittle]

Errr... no. If averages[0] is taken at 0° (and also at 180°but in counterphase) averages[1] at 45° and so forth, the math is correct. However, it simply can't improve faith of this thing. There are a few details I'm stealing though

In a way there is a similar path to this that could lead to a proper amateur multifrequency rig that could perform well on salt, mineralised, and even viscous grounds. With obscenely simple code and no pompous math. However, 10 bits, and in fact even 12 bits ADC is downright silly.


Here is Mr Crocker's reply to my email.


This is why:
Say you aim for not-too-ambitious 4nV/sqrt(Hz) performance, and your target bandwidth is 25Hz. It gives 4nV x sqrt(25) = 20nV of input referenced noise. Say we put a 100x preamp as well, and we get to 2uV amplified noise level. With such amplification and a good coil, we get, say, 200mV worst case air signal, and S/N ratio of 100000. Throw in a sign bit for negative values, and you get over 17bit resolution. 10bit ADC is simply not in that not-too-ambitious league.

Here is Mr. Crockers reply to my questions.

Dealing with the arduino metal detector, I'm trying to understand
certain parts of the code, what you are doing and most importantly
why. I'm a bit math challenged, but I also
am an embedded firmware engineer with 30+ years experience.
So please make it really clear...

first question is in the ISR . It looks like you are sampling 4
times per half cycle. BUT on the negative half cycle, you are
*subtracting* it from the buffers.
The ADC returns a number from 0-256. Where is it biased? I.e. is the
adc seeing the value as signed?

The hardware (see schematic on github) uses a single rail op amp to
amplify the signal from the receive coil. The output of the op amp is
biased to approximately half the voltage on the 3.3V rail, which is also
used as the ADC voltage reference. So at the zero crossing point, the
ADC reading is approximately half of the 255 maximum ADC reading.

What value does it return for the zero crossing point (degree 0 and
180) of the sine wave.
If you sample at say 45 degrees and 135 degrees, and you get say a +70
for the first one and a -70 for the second, then subtracting the
negative 2nd one from the bucket is the same as
adding it in. Taking the abs and adding it would have been a bit
clearer if that is what is happening.

Taking the abs and adding it would not work. What I have implemented is
a phase sensitive detector in software. Are you familiar with phase
sensitive detection? The idea is to measure a signal of the same
frequency as the reference oscillator, and be insensitive to signals at
other frequencies (including 50/60Hz mains).

If however you would get a 0 for the negative peak of the sine wave,
than all numbers will be positive and you would not be adding things up.
So what is happening here....

if (ctr < 4)
{
*p += (val);
if (*p > 15000) *p = 15000;
}
else
{
*p -= val;
if (*p < -15000) *p = -15000;
}

Suppose the ADC gives a reading of 100 with no signal. The additions and
subtractions will cancel out, so after a whole number of cycles, all the
bins will be zero. It doesn't matter whether the no-signal reading is
100 or some other value around the middle of the ADC range, as long as
it doesn't change too often. Now suppose there is a signal, giving a
reading is 103 on the positive peak of the signal, and 97 on the
negative peak. If one of the 4 samples in the first half cycle catches
the positive peak, and 4 samples later we catch the negative peak, then
for every cycle the corresponding bin will increase by 6.

I'm limiting the value in each bin to +/- 15000 to avoid arithmetic
overflow in the subsequent calculations when the signal is strong.

The next question involves the following:
// Massage the results to eliminate sensitivity to the 3rd harmonic,
and divide by 200
double bin0 = (averages[0] + halfRoot2 * (averages[1] - averages[3]))/f;
double bin1 = (averages[1] + halfRoot2 * (averages[0] + averages[2]))/f;
double bin2 = (averages[2] + halfRoot2 * (averages[1] + averages[3]))/f;
double bin3 = (averages[3] + halfRoot2 * (averages[2] - averages[0]))/f;

What is going on here?
Why the addition of .707 to the raw ADC reading?
Why the division by 200? What is that accomplishing?
Where is the division by 1024 samples to actually compute the averages.
How does this cancel the 3rd harmonic?
and why are some of these additions and some subtractions.

You really need to understand Fourier analysis to appreciate how this
cancels the third harmonic. I can't remember why I divided by 200, I
think it was just to make the displayed numbers more manageable when I
was debugging and writing the 4 bin values to the console. Dividing by
256 might be more efficient, depending on how the floating point maths
library is written. Now that the debugging code has been removed, it
would be simpler to divide the final amplitude by 200 instead, or to
just accept that the final amplitude covers a larger range.

At the start of this calculation, the 4 averages give the amplitude of
the received signal measured at 4 different reference phases, 45 degrees
apart. At the end, the 4 bins also give the amplitude of the received
signal measured at 4 different reference phases 45 degrees apart, but
the 3rd harmonic component of the signal has been cancelled, and the
result has been converted to floating point and scaled down to a more
reasonable range (about -100 to 100 I think).

There is no need to divide by 1024 because we are not interested in the
absolute value of the signal, just its relative value. The code
calculates a scaled average rather than a true average.

This is the biggest question. I want to understand *exactly* what is
going on and why.
I understand the sqrt(bin0^2+bin2^) but why bin 0 and 2 and 1 and 3?

I said earlier that the software was implementing a phase sensitive
detector. More precisely, it is implementing two quadrature phase
sensitive detectors with reference signals 45 degrees apart. Each
quadrature detector itself comprises two phase sensitive detectors with
reference signals 90 degrees apart. The first quadrature detector
produces results in bins 0 and 2, the second produces results in bins 1
and 3.

Consider just one quadrature detector. We know the component of the
received signal at the reference phase, and the component at 90 degrees
to that phase. Assuming the signal is sinusoidal, to find the amplitude
of the actual signal, we need the square root of the sum of the squares
of the two components.

At this point, we could take the output of just one of the quadrature
detectors, and calculate the amplitude and phase from it. But has we
have 2 detectors (which is a consequence of sampling 8 times per cycle
instead of just 4, which is what allowed us to cancel the third
harmonic), it makes sense to use both outputs. So I calculate the signal
amplitude from both phase sensitive detectors, and average them.

I also don't understand how the phase adjustment is working.

The next question concerns:
double phase1 = atan2(bin0, bin2) * radiansToDegrees + 45.0;
double phase2 = atan2(bin1, bin3) * radiansToDegrees;

I know that the arctan will give the phase angle of a vector, but
again why bin 0 and 2 and why the addition of +45 to the resulting
degree measurement.

Bin 0 and 2 because they give the signal amplitudes from samples taken
90 degrees apart. Likewise bins 1 and 3 give the signal amplitudes from
samples taken 90 degrees apart, however those samples are taken 45
degrees later than bins 0 and 2. Hence the +45 so that both phase
results are referenced to the same point in the cycle.

Is this measurement essentially the same as the VDI that takes the X
and Y voltage components, and computes the phase angle? (there is an
arduino VDI posted on the Geotech forums).

Sorry, I don't recognise the acronym VDI.

What I want to do is to modify this code so it works with 6.592 khz
whites detector coils.


I have found a T1 divisor of 243 and a sample divisor of 10 will let
me use a clock frequency that gives me 6.592 khz to generate the
clock, but
That means sampling not at 8 times, but at 10, so rather than 45
degrees it would be 36 degrees. I KNOW from many hours studying
whites schematics, that 90 degree sample is important.

It would not work having 10 samples per cycle. I guess you could use 12
samples (3 quadrature detectors instead of 2), but without doing the
maths I don't now whether you would be able to cancel the third
harmonic. Your detector coils won't mind a slightly different frequency,
so I suggest you use a slightly lower frequency that gives you a
convenient divisor.

I am also working on using an ARM (ST32F411 nucleo board) for a IB and
PI metal detector. In that case, we have a 1 usec sample rate, and I
can digitize the entire waveform using
the DMA feature. Then I can integrate between various samples, but I
need to understand what you are doing here, because that is directly
applicable.

Additionally, I don't see any filters happening here to remove high
frequency and low frequency components, such as 60 cycle hum etc.

The phase sensitive detection makes the system insensitive to
frequencies other than the operating frequency and some of its
harmonics. The amplifier incorporates simple high- and low-pass
filtering. All that is needed is enough to prevent strong signals such
as 50/60Hz mains and long wave radio overloading the amplifier.

[Davor]

Errr... no. If averages[0] is taken at 0° (and also at 180°but in counterphase) averages[1] at 45° and so forth, the math is correct. However, it simply can't improve faith of this thing. There are a few details I'm stealing though

In a way there is a similar path to this that could lead to a proper amateur multifrequency rig that could perform well on salt, mineralised, and even viscous grounds. With obscenely simple code and no pompous math. However, 10 bits, and in fact even 12 bits ADC is downright silly.

This is why:
Say you aim for not-too-ambitious 4nV/sqrt(Hz) performance, and your target bandwidth is 25Hz. It gives 4nV x sqrt(25) = 20nV of input referenced noise. Say we put a 100x preamp as well, and we get to 2uV amplified noise level. With such amplification and a good coil, we get, say, 200mV worst case air signal, and S/N ratio of 100000. Throw in a sign bit for negative values, and you get over 17bit resolution. 10bit ADC is simply not in that not-too-ambitious league.

[Qiaozhi]

..... the samples are 45° apart, so bin 0 and bin 2 are in quadrature, and also bin 1 and bin 3 are in quadrature, but at angle, and atan2 returns the vector angle between them.

The samples are 45° apart, but the quadrants are separated by 90°, which puts them in quadrature. Bins 0 and 2 are actually 180° apart (same for 1 and 3). Remember, there are two samples per quadrant.

Basically, the calculations are a mess.

[Davor]

There is a potential to this, but ...
BTW, the samples are 45° apart, so bin 0 and bin 2 are in quadrature, and also bin 1 and bin 3 are in quadrature, but 45° further, and atan2 returns the vector angle between them.

[Qiaozhi]

Having looked at the code some more, I think I now understand how this is supposed to work.

Firstly, an array of 4 bins are declared, which are used to accumulate the ADC readings. One cycle is split into 4 equal parts (90 degree quadrants) and are referred to as "phases" (0 to 3). The ADC capture is performed in the ISR (interrupt service routine). Two samples are taken per phase, making a total of 8 readings per cycle. After 1024 sets of readings have occurred, the contents of the bins are deposited into an "averages" array.

Then 4 different "bins" are declared, and the mysterious addition and subtraction process takes place. Initially I assumed that each of the 2 samples per phase were positioned around the zero-crossings and the positive and negative peaks, as this would make most sense. In that case, we would have the following:

phase0 = positive-going zero-crossing -> bins[0] -> averages[0] -> bin0
phase1 = positive peak -> bins[1] -> averages[1] -> bin1
phase2 = negative-going zero-crossing -> bins[2] -> averages[2] -> bin2
phase3 = negative peak -> bins[3] -> averages[3] -> bin3

With the ADC samples positioned around the zero-crossing and the peaks, we can examine the following code snippet:

double bin0 = (averages[0] + halfRoot2 * (averages[1] - averages[3]))/f;
where: halfRoot2 = 0.7071 and f = 200.0

In this instance, the calculation becomes:

bin0 = ((postive-peak - negative-peak) * 0.7071) + positive-going-zero-crossing

f appears to be an empirically defined fudge factor.

Assuming the peaks cancel (not sure why he's bothering to find the RMS value at this stage), bin0 will mainly contain the reading at the positive-going zero-crossing. Likewise, you can analyze the other three phases in the same way, and get the following:

bin0 = positive-going zero-crossing
bin1 = positive peak
bin2 = negative-going zero-crossing
bin3 = negative peak

Next he calculates the absolute amplitude at the zero-crossings:
double amp1 = sqrt((bin0 * bin0) + (bin2 * bin2));

And the absolute amplitude at the peaks:
double amp2 = sqrt((bin1 * bin1) + (bin3 * bin3));

So far so good.

At this point he has all the ingredients to achieve both ground balance and discrimination, together with VDI.
But then things start to go a bit awry.

Instead of calculating the inverse tangent of bin1/bin 0 and the same for bin3/bin2. In other words, calculating the in-phase and quadrature components of the RX signal, he calculates the inverse tangent of bin0/bin2 and likewise bin1/bin3, which really doesn't make any sense at all.

IMHO this shows a basic lack of understanding of how a VLF detector works. In fact, he states:" So we'll say that anything with a phase-shift below -20 deg is non-ferrous". Which indicates that the detector is only designed to give a ferrous/non-ferrous indication.

As Davor correctly said: "..... there are a few details that are missing here, namely ground balance".
Which is a shame, because everything is there to make it happen.

In conclusion, I suspect my initial assumption (that the samples are being taken around the zero-crossings and peaks) is incorrect, and the samples are actually being taken in the middle of each quadrant (phase). In this case, any anomalies in the signal (due to harmonic contamination) will be reduced, but the ability to ground balance and discriminate with VDI will be lost.

According to some comments made concerning this project, it is long defunct, and apparently never made it off the bench.

[Davor]

... there are two more of these He says he is removing the 3rd harmonic and dividing by 200. What the heck is this doing

This way of addition of 3 consecutive samples actually removes the 3rd harmonic, and with factors sqrt(0.5) also the 5th as the components of the 3rd and 5th contained in them get mutually subtracted. Higher harmonics are far more volatile with changing ground proximity, so this will give more quiet rig. Even harmonics are removed due to the balanced sampling.

I have no idea why dividing by 200, and it must have something to do with arduino atan2 function that returns arc tangent of y/x. Applying ideal atan2 function would yield same result with and without divide by 200. As I'm not familiar with arduino math, I have no idea what might be the motive there. Most probably atan2 is based on some lookup table and precision may vary.

What would I do different?
First off, I'd avoid micro in the first place, despite the tempting FIR boxcar filtering. Lets pretend I regained my liking for the micros.

First there are a few details that are missing here, namely ground balance. Without it this is merely a nice pinpointer. To add a ground balance, one should play with phases in a way that linear combination of the orthogonal phases (bin0⊥bin2 and bin1⊥bin3) produce zero at ground phase. Simple really, as the routine to establish phase is already implemented. In the end the amplitude (all metal, no ground) is absolute value of the ground balance vector, instead of the sum of squares, as implemented now.

The fun part with this rig is that it is fully static, with simple "calib" function to remove offset. It would be much simpler if there is a simple boxcar highpass added as well, but it would take many more samples to average than just 1024. If there is enough memory to do 1 full second or more, it would be interesting to see perfect FIR slow motion detector with no ringing. Switch on and go - that simple.

So my approach would be to add FIR boxcar high pass in duration of, say 4 seconds, and a button to turn on GB calibration in place of the "calib" function used now. THAT would be something else.

[Daren]

STM32 the is the way to go. STM32F3 or F4 can sample 12 bits at 20nS
You can setup a nice ADC/DMA stright to your buffer then process it. It is absolutely amazing what that chip can do. Plenty of memory, 72MHz, DAC for sound or to drive coil, etc....
Here is a sample plot from the simple PI detector. Signal sampled after the first and only opamp (used to shift and expand signal to the ACD range)
plots of "base line - no target", "steel" and "gold".

After you have the profile then slow/fast integrators are done in software to detect target.
I fried my STM32F3 (12V charged cap plugged into 3.3V ooops); Just got STM32F4 nucleo but did not have time to play yet, the spec are spectacular .

[Davor]

This thing uses 1024 boxcar filtering of 4 phases and thus it increases resolution for additional 5 bits. Simultaneously you have a very neat LPF with very favourable impulse and phase response for a decent detector. I'm not overwhelmed by joy, but this can work.