14

What techniques might one use to estimate the onset time of a sinusoidal tone burst in a noisy signal?

Assume the tone burst has a known fixed frequency (but unknown phase) and a very sharp rise time, and that the goal is to estimate the onset time within better than half the rise time, and/or one period of the frequency of the tone, if possible. How might the estimation techniques change if the S/N ratio is very low (much less than 1)?

Added: Assume the tone burst is of unknown length, but longer than a small multiple of the rise time and the frequency period.

Added: A DFT/FFT shows the very probable existence of a tone. The problem is figuring out exactly precisely where in the FFT window the tone (or perhaps multiple tone bursts of the same frequency) may have started within the FFT window, or determining if the current tone started outside that DFT window, provided I have all that additional time domain data.

Radar pulse detection accuracy is closer to the resolution I need, except I only have an edge, as the tone is of unknown length, and, other than a known rise time, unmodulated. Narrow band pass filters distort the rise time, and thus kill edge arrival estimation resolution.

Peter K.
  • 21,266
  • 9
  • 40
  • 78
hotpaw2
  • 33,409
  • 7
  • 40
  • 88
  • 1
    Can we assume anything about the noise? Is it stationary? Does it follow any sort of a distribution? – Phonon Sep 01 '11 at 19:34
  • 2
    Are false alarms from your detector undesirable? Do you have a specification on the probability of correctly detecting each pulse? This is very similar to (a simplified version of) front-end radar signal processing; locating (possibly-modulated) pulses embedded in noise and estimating their parameters. – Jason R Sep 01 '11 at 19:48
  • Of course, the probability of both false alarms and missed tones should be minimized to whatever the noise level (type and distribution, etc.) might reasonably allow. – hotpaw2 Sep 01 '11 at 21:48
  • @Phonon : Assume stationary white noise. Would 1/f noise make a big difference in the basic estimation algorithm? – hotpaw2 Sep 01 '11 at 21:51
  • 1
    Do you need to do this in real time, or is it an offline analysis? – nibot Sep 03 '11 at 00:08
  • @nibot : As long as the edge estimation has a known latency bound from the data stream, one could do either. – hotpaw2 Sep 03 '11 at 02:47
  • 2
    @hotpaw2: What didn't you like about the [Goertzel algorithm](http://www.eetimes.com/design/embedded/4024443/The-Goertzel-Algorithm) as per [this SO answer](http://stackoverflow.com/questions/4456855/precise-tone-onset-duration-measurement)? – Peter K. Sep 03 '11 at 15:40
  • @Peter - What does the output of this filter have to do with time estimation? – hotpaw2 Sep 03 '11 at 15:45
  • 1
    The Goertzel algorithm is used for tone detection, which seems to be what you are after. The output of the filter is an estimate of the "power" of the signal at the frequency for which it is tuned. Choose a threshold. If the filter output is above this, you have detected a tone. Set your threshold appropriately, and you can detect the onset of the tone earlier (and also be more prone to false alarms). – Peter K. Sep 03 '11 at 15:48
  • @Peter - This seems to have nothing with time. Only power. There is no t in the equation result. – hotpaw2 Sep 03 '11 at 16:22
  • 1
    @hotpaw: If your threshold is $T$ and the filter output is $g(t)$ then the onset time, $t_o$ is that time, $t$ at which $g(t) > T$. Or are you after something else? – Peter K. Sep 03 '11 at 16:44
  • @Peter : The Geortzel article you linked is about a block based algorithm. t disappears in the result of any block based DFT. Perhaps you are thinking of another type of algoritmm or filter where the output depends on t? That's more what I'm looking for. – hotpaw2 Sep 04 '11 at 16:41
  • 1
    @hotpaw2: The Goertzel algorithm isn't block-based. It's basically just a sharp bandpass filter. – Oliver Charlesworth Sep 04 '11 at 17:16
  • @hotpaw2: As Oli says, Goertzel is just a filter. [Try this page](http://www.mstarlabs.com/dsp/goertzel/goertzel.html) if you can't Google "Goertzel filter" and find an explanation that isn't block-based. I assumed that since you were using the FFT that you had to have a block-based approach. I will try to write this up later here, because it seems like it's what you need but you don't seem to grasp it from the snippets I've been posting as comments. – Peter K. Sep 04 '11 at 18:57
  • @Peter : I grasp well that the edge time estimation accuracy of a very narrow band filter is bad. Is there proof that this is anywhere near the best estimate bound possible? – hotpaw2 Sep 04 '11 at 19:25
  • @hotpaw2: If you want proofs, then you'll have to go to another level of detail in your problem description. The Goertzel algorithm can be thought of as an attempt to implement a [matched filter](http://en.wikipedia.org/wiki/Matched_filter), which can be proved to by optimal in terms of DETECTING the tone; I'm not sure about the time-estimate (onset estimate): you'd have to define what you mean by "good"... – Peter K. Sep 04 '11 at 19:50
  • Best case would be 0 rise time at the zero crossing of a known frequency sinusoid. How to estimate the exact time of this 1st zero crossing. Trivial to eye- ball. Now add noise. Now add non- zero rise time Still possible? With what bound given S/N n? – hotpaw2 Sep 04 '11 at 21:07
  • SNR doesn't tell the whole story. Is the duration of the pulse known? How precisely must you discriminate the tone frequency (i.e. do you need to accept tones at frequency $f_1$ and reject tones at a "nearby" frequency $f_2$ )? You talk about detecting the first zero-crossing of the sinusoid. Do you know the pulse's initial phase? That can help also. – Jason R Sep 05 '11 at 00:34
  • The assumption is that I know the exact frequency. Additionally I know the tone is present (say some FFT shows it's presence). So I don't care for any algorithm that tells me the tone occurred unless it can pinpoint the onset time. – hotpaw2 Sep 05 '11 at 07:15
  • @Jason : the answer to your question about phase is in my original question statement. – hotpaw2 Sep 05 '11 at 10:12
  • @Jason : can any radar return processing techniques be used with only one edge (tone of indefinite length)? Say the rise time is known, but not the phase. – hotpaw2 Sep 05 '11 at 15:03
  • @Oli : A non-block based Goertzel is not a filter, but an array of different FIR filters (or varying truncations of a very long IIR) each of whose frequency response and group delay varies enormously with current length or tap number. Sort of hard to use if ones cares about the latency or delay, e.g. when trying to measure delay. – hotpaw2 Sep 05 '11 at 16:54
  • @hotpaw2: It's still a filter (or collection of filters, one for each frequency of interest) though! – Oliver Charlesworth Sep 05 '11 at 17:01
  • @Oli : No, it's a collection of filters, one for each tap or length fraction, with a very different response even at the same design center frequency. – hotpaw2 Sep 05 '11 at 18:29
  • @hotpaw2: the filter responses are all identical, viewed on a linear frequency scale (other than a shift, obviously). – Oliver Charlesworth Sep 05 '11 at 21:21
  • @Oli : The longer you run a time-domain Goertzel algorithm, the longer the equivalent FIR filter. The longer a bandpass FIR filter, the better the adjacent frequency rejection can be. These different length equivalent FIR filters have very different frequency responses with regard to transition and stopband, given the same center frequency. – hotpaw2 Sep 05 '11 at 21:40
  • @Oli : Or are you talking about a sliding window Goertzel (which requires a long delay line for the subtraction, and which may or may not have numerical stability problems)? – hotpaw2 Sep 05 '11 at 21:49

1 Answers1

6

As we've been discussing in the comments, the Goertzel algorithm is the usual way to detect a tone in noise. After the discussion, I'm not sure it's quite what you are after (you want the onset time), but there seemed to be confusion over how the Goertzel algorithm might be applied to your problem, so I thought I'd write it up here.

Goertzel Algorithm

The Goertzel algorithm is good to use if you know the frequency of the tone you are looking for (call it $f_g$), and if you have a reasonable idea of the noise level so that you can select an appropriate detection threshold.

The Goertzel algorithm can be thought of as always calculating the output of ONE FFT bin:

$$y(n) = e^{\jmath 2\pi f_g n} \sum_{k=0}^n x(n) e^{-\jmath 2\pi f_g k}$$

where $f_g$ is the frequency you're looking for.

The Wikipedia page has a better way to calculate this.

Here is a (feeble) Scilab attempt at implementing it:

function [y,resultr,resulti] = goertzel(f_goertzel,x)
realW = 2.0*cos(2.0*%pi*f_goertzel);
imagW = sin(2.0*%pi*f_goertzel);

d1 = 0;
d2 = 0;

for n = 0:length(x)-1,
    y(n+1) = x(n+1) + realW*d1 - d2;
    d2 = d1;
    d1 = y(n+1);
    resultr(n+1) = 0.5*realW*d1 - d2;
    resulti(n+1) = imagW*d1;
end
endfunction

Consider the signal with $f = 0.0239074$ and $\phi = 4.4318752$ :

$$x = \sin(2\pi fn + \phi) + \epsilon(n)$$

where $\epsilon(n)$ is zero-mean, unit variance Gaussian white noise.

In this example, the tone starts one third of the way into the signal at index 1001.

If we run the Goertzel algorithm on it with $f_g = f - 0.001$ then we get the top two traces of the figure.

If we run the Goertzel algorithm on it with $f_g = f$ then we get the bottom two traces of the figure.

The four traces are:

  • $x$ (blue) and $y$ (red) for $f_g = 0.0229074$
  • The resulting $\sqrt{resultr^2 + resulti^2}$
  • $x$ (blue) and $y$ (red) for $f_g = 0.0239074$
  • The resulting $\sqrt{resultr^2 + resulti^2}$ (solid line) and the first result (dashed line).

As you can see, the case where the tone we are interested in is present peaks at about 250. If we set the detection threshold at about half this value (125), then the detection occurs (the square-rooted value is greater than 125) at about index 1450 --- 450 samples after the tone started.

This threshold (125) will not cause a detection in the other case (for this run, anyway), but the maximum value of that output is 115.24, we we cannot reduce the threshold too much without getting a false detection.

Reducing the threshold to 116 will cause detection in the true case (for this run) at index 1401... but we run the risk of more false alarms.

enter image description here

Peter K.
  • 21,266
  • 9
  • 40
  • 78
  • A running Goertzel filter is more suitable if one is only looking for an existence estimate within a fixed length window. A running Goertzel without a loss/decay term changes it's bandwidth over its length, and the narrower bandwidth later in time provides a worsening arrival time estimate, more sensitive to noise and threshold errors. – hotpaw2 Sep 05 '11 at 15:09
  • @hotpaw2: Correct. You can introduce a "forgetting factor" to keep Goertzel running, but otherwise it remembers everything. – Peter K. Sep 05 '11 at 15:27
  • Remembers everything? It's an FIR that can be implemented in recursive form. What have I missed here? – Oliver Charlesworth Sep 05 '11 at 21:28
  • @Oli: If you look at the equation for $y(n)$ above, you'll note that it doesn't end. Yes, it's estimating a (scaled) DFT coefficient, but it's definitely not FIR. – Peter K. Sep 05 '11 at 23:26