Goertzel Alogorithm Stability in Low SNR

On Nov 15, 3:43 pm, "superlou" <lsim...@safeflight.com> wrote:
> >The gain for a Goertzel is the same as for a DFT (ignoring numerical
> >issues). So if you are feeding an input whose frequency matches that
> >of the Goertzel, then your output will by scaled up by N/2 (amplitude)
> >or (N/2)^2 for energy, where N is the number of samples. I don't think
> >you've told us how many samples you are feeding into the thing.
>
> >So expect your output to be 32767*N/2 for the expected amplitude.
>
> >Since you didn't give the number of samples, we can back calculate to
> >see how many your have
>
> >(if 156dB is amplitude)
> >156 = 10 log(32767*N/2) ->  N = 2.4E11 samples if 156dB is your
> >amplitude
>
> >(if 156dB is energy)
> >156 = 20 log(32767*N/2) -> N = 3851 samples
>
> >Are you using about 3850 samples?  If so then your numbers seem about
> >right.
>
> >Why don't you let us know: sample rate, number of samples, amplitude
> >of signal
>
> >Then we can check to see if the algo's output makes sense.
>
> >Clay
>
> The output of the Goertzel algorithm is in terms of amplitude squared, and
> so 10*log(output^2) = 20*log(output).  If I had received 157dB output
> it would have indicated 4321 data points and the output I gave was simply
> truncated down.  Altering the program to give me another two significant
> figure, I get 156.32dB which leads to 3996 points.  Since I'm using a 4000
> sample block of data for the Goertzel, this seems ok to me.
>
> The signal amplitude was 32767 (range [-32767 32767]) and the sampling
> rate is 4000Hz.  In order to provide more than one energy reading per
> second, I am shifting in data from a 40 point input buffer and
> recalculating the energy after each shift in.  The input buffer size in my
> testing has ranged from about 40 to 400 in order to free up processing time
> for filtering and other processing, but the goertzel output remains the
> same (just gets calculated less often).
>
> The interesting thing I am noticing is that when I run the algorithm with
> the 32767 ampl. wave, the 60Hz energy level is a constant output while
> another Goertzel algorithm with a target frequency at 50Hz is oscillating
> between 0 and 12 dB (with an occasional -112).  When I run the exact same
> test with an input 1/2 the size (16382 ampl.), the 60Hz energy retreats to
> 150dB (expected 6.1 dB decrease for 1/2 peak to peak range) but the 50Hz
> energy is constant at -108dB (no oscillation).  Decreasing the amplitude
> by another factor of 2 sees the 60Hz energy to 144dB and 50Hz energy to
> -114dB.  It looks almost like something is saturating in the algorithm
> with the full scale wave, but all the data types are double and I have not
> had any luck finding it yet.  Being fairly inexperience with DSP in general
> and the Goertzel algorithm in particular, I was wondering if anyone had run
> across a similar situation.
>
> Thanks,
> Louis- Hide quoted text -
>
> - Show quoted text -

Hello Louis,

I think your code is running fine now. When you input 60Hz and you see
the expected value coming from the 60 Hz detector and a nearly zero
value coming from the 50Hz detector - this is what you should get. The
confusion you have stems from the dB representation of the numbers. If
your 50 Hz level is over 100dB below that of the 60 Hz level, you are
essentially getting zero of the 50Hz.  You know doubles have
essentially 16 decimal places in their mantissas. So when you subtract
two of them with similar scaling, a lot of the digits go away.

So I'm not suprised that you are seeing some non-linear stuff because
of numerical noise. But your detector is working with well over a
100dB range. How much do you need? Try normalizing your numbers to
32767*(4000/2) and displaying the numbers in a linear scale. I think
this will put stuff into perspective. If you recall the limit of
log(x) as x->0 is -infinity, so your near zero outputs on a log scale
will be all over the place, but in contrast to a full scale output,
they indicate essentially zero out.

IHTH,
Clay

superlou wrote:

   ...

> The signal amplitude was 32767 (range [-32767 32767]) and the sampling
> rate is 4000Hz.  In order to provide more than one energy reading per
> second, I am shifting in data from a 40 point input buffer and
> recalculating the energy after each shift in.  The input buffer size in my
> testing has ranged from about 40 to 400 in order to free up processing time
> for filtering and other processing, but the goertzel output remains the
> same (just gets calculated less often).
> 
> The interesting thing I am noticing is that when I run the algorithm with
> the 32767 ampl. wave, the 60Hz energy level is a constant output while
> another Goertzel algorithm with a target frequency at 50Hz is oscillating
> between 0 and 12 dB (with an occasional -112).  When I run the exact same
> test with an input 1/2 the size (16382 ampl.), the 60Hz energy retreats to
> 150dB (expected 6.1 dB decrease for 1/2 peak to peak range) but the 50Hz
> energy is constant at -108dB (no oscillation).  Decreasing the amplitude
> by another factor of 2 sees the 60Hz energy to 144dB and 50Hz energy to
> -114dB.  It looks almost like something is saturating in the algorithm
> with the full scale wave, but all the data types are double and I have not
> had any luck finding it yet.  Being fairly inexperience with DSP in general
> and the Goertzel algorithm in particular, I was wondering if anyone had run
> across a similar situation.

Have you verified that no intermediate calculations overflow or saturate 
at any frequency?

Jerry
-- 
Engineering is the art of making what you want from things you can get.
�����������������������������������������������������������������������

On Nov 15, 12:43 pm, "superlou" <lsim...@safeflight.com> wrote:
> >The gain for a Goertzel is the same as for a DFT (ignoring numerical
> >issues). So if you are feeding an input whose frequency matches that
> >of the Goertzel, then your output will by scaled up by N/2 (amplitude)
> >or (N/2)^2 for energy, where N is the number of samples. I don't think
> >you've told us how many samples you are feeding into the thing.
>
> >So expect your output to be 32767*N/2 for the expected amplitude.
>
> >Since you didn't give the number of samples, we can back calculate to
> >see how many your have
>
> >(if 156dB is amplitude)
> >156 = 10 log(32767*N/2) ->  N = 2.4E11 samples if 156dB is your
> >amplitude
>
> >(if 156dB is energy)
> >156 = 20 log(32767*N/2) -> N = 3851 samples
>
> >Are you using about 3850 samples?  If so then your numbers seem about
> >right.
>
> >Why don't you let us know: sample rate, number of samples, amplitude
> >of signal
>
> >Then we can check to see if the algo's output makes sense.
>
> >Clay
>
> The output of the Goertzel algorithm is in terms of amplitude squared, and
> so 10*log(output^2) = 20*log(output).  If I had received 157dB output
> it would have indicated 4321 data points and the output I gave was simply
> truncated down.  Altering the program to give me another two significant
> figure, I get 156.32dB which leads to 3996 points.  Since I'm using a 4000
> sample block of data for the Goertzel, this seems ok to me.
>
> The signal amplitude was 32767 (range [-32767 32767]) and the sampling
> rate is 4000Hz.  In order to provide more than one energy reading per
> second, I am shifting in data from a 40 point input buffer and
> recalculating the energy after each shift in.  The input buffer size in my
> testing has ranged from about 40 to 400 in order to free up processing time
> for filtering and other processing, but the goertzel output remains the
> same (just gets calculated less often).
>
> The interesting thing I am noticing is that when I run the algorithm with
> the 32767 ampl. wave, the 60Hz energy level is a constant output while
> another Goertzel algorithm with a target frequency at 50Hz is oscillating
> between 0 and 12 dB (with an occasional -112).

Note that 50 Hz and 60 Hz sinusoids will be orthogonal
every 0.1 seconds of window length (every 400 samples
in your example).  They will not be orthogonal over other
durations of filtering, thus the Goertzel filter output
will not be zero for those durations.  If you want only
orthogonal correlations, then you can try running two
Goertzel filters starting 400 samples apart and subtract.




IMHO. YMMV.
--
rhn A.T nicholson d.0.t C-o-M

After a lot of tedious code stepping, I'm fairly certain nothing is
overflowing at reasonable frequencies.  The 0-12 dB oscillation with 3dB
mean is something I should have expected, and is insignificant since it
pertains to a signal of manitude 0-30 while the majority are around
1x10^15 or much larger.

The initial issue had to do with the Goertzel output being noisy with a
mean of about 90dB and range of 20dB up and down.  As the signal
increased, the mean and variance first decreased, and then the mean began
increasing as expected.  Since I am performing detection based on
thresholding against ambient noise, this initial decrease limits the
sensitivity of the device.

By simulating a 60Hz input in software plus additive noise (approx.
white), I have been able to approximate the 90dB out with the 20dB range. 
However, as I raise the magnitude of the 60Hz signal, the mean rises
monotonically without the initial drop (once it overcomes the noise).  It
is looking less and less like a software issue and I may have to go back
to hunting on the hardware side of things.

>
>Note that 50 Hz and 60 Hz sinusoids will be orthogonal
>every 0.1 seconds of window length (every 400 samples
>in your example).  They will not be orthogonal over other
>durations of filtering, thus the Goertzel filter output
>will not be zero for those durations.  If you want only
>orthogonal correlations, then you can try running two
>Goertzel filters starting 400 samples apart and subtract.
>

I'm not certain how the orthogonality of the two is applied since the 50Hz
and 60Hz levels are being analyzed with two seperate implementations of the
goertzel and returning indepedant single values.  In this case are you
talking about using a goertzel alorithm's states to provide outputs of a
filter and then subtracting a two time (sample) delayed results?  I
haven't had any luck finding an explanation of this technique and would
appreciate if you could point me towards an example or literature.

Thank you,
Louis