The slippery slope of sliding Goertzels

Nearly 30 years after first reading about the sliding-window Goertzel 
algorithm, I finally have cause to try implementing it. :-)

My understanding is you just add one step to the recursive part of the 
usual Goertzel algorithm - you just subtract the input sample from one 
block length back. I tried that, and it basically works. However, there 
are various grubby little details which bother me.

I remember there was some discussion in this group a couple of years 
back when Rick Lyons and Eric Jacobsen wrote an article on this subject, 
and someone was querying the details. However, checking that discussion 
it seemed to fade out without really answering anything.

Because the structure of the Goertzel algorithm is IIR and tuned for 
resonance, it seems a lot like an oscillator. It feels like there should 
be a danger of the output level wandering, in the way that uncontrolled 
IIR oscillators tend to. However, in trying the algorithm with single 
precision floating arithmetic, and certainly loosing some bits along the 
way, it always seems stable. I'm not clear if something I'm missing 
makes it unconditionally stable, or if I am simply not prodding it hard 
enough. The block processing of the standard Goertzel algorithm prevents 
this being an issue there.

If the non-sliding Goertzel algorithm has a block length of N, you clear 
the registers, push N samples of signal through the filter, then a zero, 
then the result is available. It seems with the sliding version I also 
need to push a zero through before calculating the non-recursive half of 
the filter. Can someone confirm that is correct?

The non-sliding Goertzel algorithm only gives optimal results when the 
block length is a whole number of cycles of the filter's centre 
frequency. However, the filter still works well for other block lengths. 
You just get some ripple in the block to block output for a stationary 
input signal. This means if you can be a little tolerant about the 
output level, you can tweak the bandwidth to taste. Typically, this 
seems to be what people do.

With the sliding Goertzel algorithm, it seems that choosing N to be 
anything other than an exact whole number of cycles of the filter's 
centre frequency makes the results fall apart completely. Is there a way 
to tame this behaviour, so the bandwidth of the filter can be tweaked in 
a more flexible manner? It worries me that if things are this sensitive, 
then even when the block length is a multiple of the filter's centre 
frequency, the tiniest error in the frequency dependent coefficient 
could have bad consequences. It seems to work OK in the cases I have 
tried, but can that be guaranteed?

Any help greatly appreciated.

Steve

On Wed, 16 Nov 2005 20:25:52 +0800, Steve Underwood <steveu@dis.org>
wrote:

>Nearly 30 years after first reading about the sliding-window Goertzel 
>algorithm, I finally have cause to try implementing it. :-)
>
  (snipped)
>
>With the sliding Goertzel algorithm, it seems that choosing N to be 
>anything other than an exact whole number of cycles of the filter's 
>centre frequency makes the results fall apart completely. Is there a way 
>to tame this behaviour, so the bandwidth of the filter can be tweaked in 
>a more flexible manner? It worries me that if things are this sensitive, 
>then even when the block length is a multiple of the filter's centre 
>frequency, the tiniest error in the frequency dependent coefficient 
>could have bad consequences. It seems to work OK in the cases I have 
>tried, but can that be guaranteed?
>
>Any help greatly appreciated.
>
>Steve

Hi Steve,
  I'm gettin' ready to take off on a business 
trip tomorrow so I can't spend the time to 
dig through my MATLAB code to help answer all 
your sensible questions.  

I'll send ya a copy of the two articles that 
contained Goertzel & sliding Goertzel info.
Maybe they'll be of some use to you.

[-Rick-] 
 

On Wed, 16 Nov 2005 20:25:52 +0800, Steve Underwood <steveu@dis.org>
wrote:

  (snipped)
>
>Any help greatly appreciated.
>
>Steve

Steve,
  I tried to E-mail you, but the 
<steveu@dis.org> address didn't work.

Gotta run.

[-Rick-]

Steve Underwood wrote:
> Nearly 30 years after first reading about the sliding-window Goertzel 
> algorithm, I finally have cause to try implementing it. :-)
> 
> My understanding is you just add one step to the recursive part of the 
> usual Goertzel algorithm - you just subtract the input sample from one 
> block length back. I tried that, and it basically works. However, there 
> are various grubby little details which bother me.
> 
> I remember there was some discussion in this group a couple of years 
> back when Rick Lyons and Eric Jacobsen wrote an article on this subject, 
> and someone was querying the details. However, checking that discussion 
> it seemed to fade out without really answering anything.
> 
> Because the structure of the Goertzel algorithm is IIR and tuned for 
> resonance, it seems a lot like an oscillator. It feels like there should 
> be a danger of the output level wandering, in the way that uncontrolled 
> IIR oscillators tend to. However, in trying the algorithm with single 
> precision floating arithmetic, and certainly loosing some bits along the 
> way, it always seems stable. I'm not clear if something I'm missing 
> makes it unconditionally stable, or if I am simply not prodding it hard 
> enough. The block processing of the standard Goertzel algorithm prevents 
> this being an issue there.

I haven't applied the Goertzel specifically, but I _have_ applied lots 
of IIR filters.  The Goertzel isn't an oscillator, really, rather it's 
an IIR filter who's poles or exactly on the stability boundary.  As such 
_any_ gain deviation in any of the multiplies will cause the poles to 
deviate from this boundary which will cause the response to either decay 
or grow.

In addition there will be quantization noise which will either add to or 
subtract from the data, most likely in unpredictable ways.
> 
> If the non-sliding Goertzel algorithm has a block length of N, you clear 
> the registers, push N samples of signal through the filter, then a zero, 
> then the result is available. It seems with the sliding version I also 
> need to push a zero through before calculating the non-recursive half of 
> the filter. Can someone confirm that is correct?

I always see the Goertzel expressed as

y_n = 2 cos(th) y_{n-1} - y_{n-2} + x_n.                    (1)

This not only makes the algorithm more sensitive to gain and data 
quantization than necessary, it also makes it a bit obscure to pull out 
the "actual" value of the sine wave if it happens to be at a null when 
you finish your data.  For resonant filters I prefer to use a structure like

[ y1_n ]   [ cos(th)   sin(th) ] [ y1_{n-1} ]   [ 1 ]
[      ] = [                   ] [          ] + [   ] x_n.   (2)
[ y2_n ]   [-sin(th)   cos(th) ] [ y2_{n-1} ]   [ 0 ]

For a sinusoid at x_n = cos(th * n) the y1 and y2 will contain inphase 
and quadrature parts of the response, respectively.  So the strength and 
phase of the incoming signal can be "read" directly from y = y1 + i y2 
(or the other way around -- I'm doing this from memory, it's your job to 
check the math).

An additional advantage (and the reason I use this for resonant filters) 
is that the amount of precision needed to attain a given accuracy will 
drop by roughly a factor of 2.  This means that if you needed double 
precision to implement (1) adequately you only need single precision to 
implement (2).  Unless your processor has such large gonads that it 
takes double-precision floating point in stride this can be a big advantage.

> 
> The non-sliding Goertzel algorithm only gives optimal results when the 
> block length is a whole number of cycles of the filter's centre 
> frequency. However, the filter still works well for other block lengths. 
> You just get some ripple in the block to block output for a stationary 
> input signal. This means if you can be a little tolerant about the 
> output level, you can tweak the bandwidth to taste. Typically, this 
> seems to be what people do.
> 
> With the sliding Goertzel algorithm, it seems that choosing N to be 
> anything other than an exact whole number of cycles of the filter's 
> centre frequency makes the results fall apart completely. Is there a way 
> to tame this behaviour, so the bandwidth of the filter can be tweaked in 
> a more flexible manner? 

With either form you could predict the value of the states resulting 
from an input some number of samples earlier.  With the technique you 
state and equation (1) you would use an even number of cycles, which 
would guarantee that for nd = even number of cycles y_{n-1} would be 
equal to (x_{n-nd} + y_{n-nd-1})/(2 cos(th)), so subtracting out 
x_{n-nd} would give you the state before you added the thing in in the 
first place.  Similar results would hold for (2), with only y1_n being 
affected by x_{n-nd}.

If you wanted to set the bandwidth to some non-integer number of cycles 
then you would have to modify both y_n and y_{n-1} in the case of 
equation (1), or both y1_n and y2_n in the case of equation (2).  Just 
calculate the response of y_n and y_{n-1} at n = nd for an impulse at n 
= 0; that's your gains for subtracting out the effect.  For equation (2) 
these gains will be (I think, check) cos(th * nd) and sin(th * nd).

> It worries me that if things are this sensitive, 
> then even when the block length is a multiple of the filter's centre 
> frequency, the tiniest error in the frequency dependent coefficient 
> could have bad consequences. It seems to work OK in the cases I have 
> tried, but can that be guaranteed?
> 
I don't think it can be guaranteed, no.  In fact I _do_ think that "no" 
can be guaranteed.  Because the results of sine and cosine operations 
are irrational numbers there will always be dangling bits, so at some 
point the error will build up.  Unless you can guarantee that your 
system will only run for to some definite maximum amount of time before 
the filters are cleared you cannot guarantee that the inevitable error 
build-up won't mess things up.

I think you could sidestep many of the precision issues by multiplying 
the incoming signal by e^(i th n) and doing a running average in the 
resulting complex number.  You'll have some precision issues in the 
computation of e^(i th n) = cos(th n) + i sin(th n), but if you store 
the _result_ to subtract out of your accumulator, and if you used 
fixed-point math, you'll be guaranteed of a bit-for-bit match.

If you must use a sliding-Goertzel-like-algorithm I suggest you add some 
damping.  Use the transfer function

                whatever
H(z) = -------------------------
        z^2 - 2 d cos(th) z + d^2

with d ever so slightly less than 1.  I'd suggest you use the form of 
(2) but that's your choice.  For the subtraction you'll definitely need 
to calculate what to take away from your states; at best you can arrange 
this to be zero on an even sample boundary and just subtract some d^nd * 
x_{n-nd}.

Using a damped form will give you a filter that is not a perfect sinc 
function centered on your frequency of interest; it'll be close to a 
sinc with a little bit of oddball phase shift and imperfect nulls -- but 
if you're just looking for a bandpass filter it'll do.

> Any help greatly appreciated.
> 
> Steve

Sorry I didn't respond to this earlier -- I was busy last week when I 
saw this & put it aside "for later".  I didn't get to it until Rick's 
replies brought it back to the front on my newsreader.

One final comment:

Whenever I see someone propose a sliding average I have to ask "is that 
really the low-pass filter you want to implement?".  So I'll ask the 
equivalent question here:  Is this really the bandpass filter you want 
to implement?

You should be able to implement a 4th-order IIR bandpass filter with 
only a bit more processor time and much less memory than a sliding 
Goertzel, without the worries about the metastable filter core running 
away with your accuracy.  True, the filter will have a different shape, 
most particularly in that it won't have the nulls that a sliding 
Goertzel does, but do you really need the shape?

Similarly, you should be able to implement a FIR bandpass filter.  With 
a DSP it should even be without a great sacrifice of processor time 
unless it's quite long -- in which case you can demodulate, decimate, 
and filter with something shorter.

If it were me  I'd look to alternatives to make sure that the sliding 
Goertzel was right for my application, and even then I'd strongly 
consider using a leaky sliding Goertzel to keep things safe.

-- 

Tim Wescott
Wescott Design Services
http://www.wescottdesign.com