[time-nuts] A simple sampling DMTD

[Jan-Derk Bakker]

Dear Joe,

On Wed, Nov 27, 2019 at 7:22 PM Joseph Gwinn <joegwinn at comcast.net> wrote:

> > [snip]
> > The 1001-point FIR band pass filter is a good reference to get an idea of
> > the best case performance of the system, but it is computationally
> > infeasible to run on an 8-bit processor. While a cheap comb filter can
> take
> > a bite out of the HF noise, canceling offset/drift is harder. Early on I
> > was looking into ways to average all samples of a single period of the
> beat
> > note, but I had trouble finding a closed form solution to the fact that
> the
> > beat note period is never an exact integer multiple of the sampling rate.
> > Through numerical methods I did end up finding a good estimator for the
> > "leakthrough" caused by the fractional part of the beat note period (as a
> > function of the measured period and the phase offset), which was fairly
> > inexpensive to implement. [snip]
>
> In the radar world, the standard solution to the leakthrough problem is
> to batch the data and apply a windowing function to the data in each
> batch.  Typically, the batches overlap such that every sample appears
> in two batches.  The window functions largely eliminate the splice
> error due to the FFT, which is fact splices each batch into a circle,
> causing a discontinuity at the splice.  If the window function is very
> small at the splice, there is little discontinuity power sprayed into
> innocent FFT bins.
>

This is an interesting approach, and one which hadn't occurred to me while
I was working on the offset/drift cancellation. After some consideration I
am afraid it will not help much in this particular case.

What I want is to filter out as much DC/LF energy before running the zero
crossing detector (ZCD) on the rising edge of the sampled sine wave. To do
this right now I also run the ZCD on falling edges, and I remember the
(sub)sample position of the current and previous falling edges. I then
average a full nominal period of the sine wave between these two falling
edges(200 samples at 10Hz beat frequency and 2ksps after the first sinc^2
filter/decimator), to determine the offset at the rising edge in between.
This is indeed mathematically equivalent to calculating bin 0 of a forward
Fourier transform with a rectangular window. Because the actual period is
not an exact integer number of samples, in practice this bin is
contaminated by spectral leakage. Windowing would help here, if it weren't
for the fact that all of my wanted signal energy is in bin 1. As any
windowing function other than the rectangular window trades side-lobe
suppression for main lobe width, the cure would likely be worse than the
disease. This could be ameliorated by averaging over multiple periods
(increasing the transform length), but that would necessarily spread the LF
noise energy over multiple bins, again making it all but impossible to get
optimal rejection.

What I've done instead is make an estimator for the leakage. For a period
of 200 +/- 1 sample, the following has an estimator error of better than
-40dB:

  alpha = 0.5*perOffset^2 + (1.5 - phase)*perOffset - phase
  correction = (1-alpha) * pointN + alpha * pointN1

where

  correction is the estimator (in sample values)
  perOffset is the difference between the measured period (199...201
samples) and the nominal period (200 samples)
  phase is the estimated starting phase of the averaging interval, in
fractional samples (range: [0,1>)
  pointN is the value of the sample immediately after the averaging
interval (ie sample[200], when the averaging interval is
sum(sample[0]...sample[199])
  pointN1 is the value of the second sample after the averaging interval
(ie sample[201]).

This works very well, both on simulated and actual data. I am embarrassed
to admit that I cannot explain *why* it must be this particular formula
(and I'm very open to get hints here).

Sincerely,

JD 'if you can't do the math, run simulations and do curve fits until it
works' B.