This article follows on from the Introduction to the DSP module algorithm at

http://humanise.org/demining/dsp

and the article on the Underlying Maths used by the DSP module at

http://humanise.org/demining/underlying-maths

Averaging / Integration

A simple averaging algorithm that could be used by the DSP module is as follows:

The module averages the results in order to find the average DC bias.  It continues to do this over a fixed number of the previous samples.  The number of samples that is taken each time is such that it normally covers a set discrete number of complete cycles of the input frequency. 
 
This algorithm can be improved in several ways.

1) Weighted Buckets

For simplicity assume we have an input waveform of 10KHz and we are obtaining a result that is proportional to a simple average over the previous second.  We will also assume that we are creating new results at say 10 times a second.  Each tenth of a second we add up all the resultant values across the last tenth of a second (ie, a thousand compete cycles) and we keep a copy of this result in what we will call a bucket.

Without dividing by the number of values this is really integrating rather than averaging.  As the number of data items used for a result is always the same, the result of integrating is directly proportional to the average value so there isn't any real need to perform the division.  In the previous documentation the term averaging was used because it was simpler to understand and corresponded closely to what happens in hardware implementations that simply used low pass filters for averaging.

Assume that at each tenth of a second we want to produce a value proportional to the average over the full last second, not just the last tenth of a second.  Consequently we keep copies of the last 10 of these buckets.  The simplest algorithm we could use is for each tenth of a second the DSP module would create a new bucket and after it creates this bucket it adds this to the addition of the previous 9 buckets and outputs the result as the amplitude for this component.  It would need to do this for the results of both the Inphase and Quadrature components.

We now analyse what happens with this if, as well as the wanted input signal at 10,000Hz, there is additional unwanted signals at say 10,001Hz, 10,001.5Hz and 10,002.5Hz etc contained within the input signal.

When the 10,001 Hz signal is multiplied by the 10,000 Hz reference signal, it will create a 1 Hz signal and a 20,001 Hz signal in the resultant signal.  The 20,001 Hz signal will appear as AC only so isn't of concern and the 1 Hz signal will nicely have a full cycle over our 1 second integration period so this will not be seen in the result.

When we repeat this with the 10,001.5 Hz signal, we get a different result for the 1.5 Hz signal.  As there is one and a half cycles of this within each one second integration, it will effect the result.  That extra half cycle will add a bias.  Similarly with the 2.5 Hz from the 10,002.5 and the 3.5 Hz from the 10,003.5, but as the frequency increases the half cycle involved becomes smaller and the change to the result gets less.

Rather than having specific frequencies that cause errors and other frequencies that don't, we can reduce the error from any individual near frequency by spreading the error effect across the frequency spectrum.  One possible algorithm that does this is as follows:

Each tenth of second the DSP module would create a new bucket and after it creates this bucket it adds this to the addition of the previous 4 buckets.  That is, it has added together the 1st to the 5th bucket.  It then multiplies this total by 10.  It then takes the 6th bucket, multiplies that by 9 and adds it to the total. It then takes the 7th bucket, multiplies that by 8 and adds it to the total. It then takes the 8th bucket, multiplies that by 7 and adds it to the total.  It continues doing this through to the 14th bucket which it just adds to the total.

A weighted bucket algorithm as just described not only results in less maximum error from any individual near frequency but it also gives the result more immediacy while still maximising the advantages of a long integration period.From calculations performed previously we can work out the expected phase (which we will call Ø) between the internal created waveform that we are currently using and the waveform in the input signal. Note that we have a different phase depending on whether we are working on the inphase or quadrature components.

2) Sliding Windows

The effect of an anomaly creating a phase change in the received signal as the coils move towards an anomaly and the reverse phase change as it goes away, is that of changing (increasing and decreasing) the width of the wavelength in the received signal.  This change of wavelength width is only due to the movement of the phase.  If the coils remain stationary over the anomaly, even though there is now a new phase difference, there is no difference in the wavelength.  As we do our averaging on a set number of samples that we assumed would be a discrete number of cycles, the wavelength change will create an error.

The error doesn't diminish if we take our readings over a larger number of cycles.  The increase or decrease in the width just goes up or down in proportion.  Further, this error appears at the time that an anomaly is coming under the coils.  This is exactly the time when we need there to be no errors.

One way to minimise this is to use a sliding window algorithm.  In this algorithm, we use the information on the current phase of the input signal calculated from our previous averaging calculation to position a window across which we take the next set of readings. Each window of samples that we add up is taken such that it is central to the integral number of cycles of the frequency of interest.  Here we assume that waves begin at rising points half way between minima and maxima and end at rising points half way between minima and maxima.  If our waveform is centered on this basis, the error from the beginning will tend to cancel the errors at the end, regardless of whether the wavelength has increased or decreased in length.

If the phase has changed on the input signal between the previous window and the window previous to that by delta_Ø, then we can assume that this is due to the speed of the coils across the ground. As coils have momentum when moving and take time to slow down or speed up, the best assumption is that the phase will continue changing for this window and the next in a similar manner. Consequently, a reasonable assumption is that the same delta_Ø change in phase will occur for this window.

The change in frequency (delta_f) on the input waveform can be calculated as follows:

delta_f = delta_Ø (degrees) / (360 * window_time (seconds))

This is the complete change in frequency from the reference frequency to the current frequency seen on the input. We are assuming no significant change from window to window.

It should also be noted that the delta_Ø for each window is also the change in length (in degrees) for the window. Consequently delta_Ø / (360 * f) is the change in time for each window. In this case we are calculating the window to window change in time on the input waveform in seconds.

The centering is done on the waveform of interest in the resultant waveform after the input waveform has been multiplied by our reference or quadrature reference waveform. Assume our reference frequency is f and signal of interest within the input signal is f + delta_f. The resultant frequency of the waveform that we are performing our averaging on has a frequency of (f + f + delta_f).

If the signal of interest within the input signal has an amplitude of A then the waveform being produced is now as follows:

- 1/2 A cos((360(2f + delta_f)t) + Ø)

We can start and finish our sliding window at positive to negative zero crossing rather than negative to positive. This positive to negative crossing may be easier due to the negative at the beginning of the above specified waveform produced.

Just as the Ø phase shift in the input signal already has delta_Ø included and that Ø is calculated as the average phase of the input waveform across the whole window, the Ø phase shift in the resultant waveform is the average phase shift for this resultant waveform across the window. The time to start each sliding window is simply the following number of seconds after the reference waveform is at the zero phase position.

Ø / (360 * (2f+delta_f))

3) Nulling Intermodulation

There will also be a delta_f frequency waveform present in the result. If we can calculate where this is then we can work out what effect this has on the result and subtract it out of the result.

From our original equation in the mathematics section, where we have an input waveform of Amplitude A and frequency f + delta_f and multiply this by our reference waveform, then:

(A sin((360(f + delta_f)t) + Ø)) * sin(360 f t)
produces a waveform - 1/2 A cos((360(2f + delta_f)t) + Ø)
plus a second waveform 1/2 A cos((360 delta_f t) + Ø)

Note that in the above equations, Ø is the currently expected phase in the input, not the standard phase expected in the absence of an anomaly. I'm calling the second waveform produced an intermodulation frequency. Note that while we're using Ø as part of the specification of this very slow long intermodulation waveform, it represents a long period of time while in the first produced waveform it only represents a short period of time.

Given the length of the window we need to work out the integral for the 1/2 A cos((360 delta_f t) + Ø) waveform across the current window and subtract that from the other window integral. With sliding windows, since the window that we are calculating doesn't start until after Ø / (360 * (2f+delta_f)) seconds, this integral shouldn't start until this time point also and it should end at a window time after this.

4) Minimising Width of Frequency Range Based on Speed of Coils Across the Ground.

In order to extract the smallest possible harmonics from noise in the input signal and to precisely measure its amplitude and phase, we need to sample over as long a period as possible.  Taking our average over a long period of time is the equivalent of narrowing the frequency range that we are measuring.  Narrowing the frequency range too much isn't practicable when the coil is moving over the ground at speed.  The effect of an anomaly creating a phase change as the anomaly comes towards the coils and the reverse phase change as it goes away (just as it is moving) is that of changing the frequency up and down.  If this new frequency is outside the narrow range of frequencies being extracted then we will lose our information.

Once an operator recognises that s/he is over an anomaly, the operator normally holds the detector still in order to read from the display.  If the detector can detect the speed over the ground, it should be able to adjust the algorithm so that it can use this period where the detector is still, to take a better reading by increasing the period of time that the reading is taken over. 

If the algorithm maintains buffers containing the previous buckets over a long period of time, then these buckets are instantly available from the time its recognised that the coils are stationary.

Note that there are alternatives for this in user procedures.  We could ask that whenever the user recognises an anomaly, the user stops the machine, pushes a button to indicate that this procedure is being done and leaves the machine to stand for several seconds.  The detector would give a beep when a set time was up.  This would work better when the detector's weight is carried on wheels (as some machines are) behind the coils so that it can come to a proper stand still.