Correcting an Important Goertzel Filter Misconception

Recently I was on the
Signal Processing Stack Exchange web site (a question and answer site for DSP
people) and I read a posted question regarding Goertzel filters [1]. One of the
subscribers posted a reply to the question by pointing interested readers to a
Wikipedia web page discussing Goertzel filters [2]. I noticed the Wiki web site
stated that a Goertzel filter:

That statement can often be misinterpreted. Because of that possible stability misconception, and the fact the literature
of the Goertzel algorithm doesn't discuss Goertzel
filter stability, I decided to write this blog.

which, with
infinite-precision arithmetic, has a single pole located on the
z-plane's unit circle as shown in Figure
1(b).

While the Figure
1(a) resonator can be used to compute individual
N-point DFT spectral samples it has a major flaw. When finite-precision
arithmetic (as in our digital implementations) is used the network's
z-plane pole does not lie exactly on the
unit circle. This means the resonator can be unstable and a computed DFT sample
value will be incorrect. In the literature of DSP the resonator is referred to
as "marginally stable."

We show this in
Figure 2 for various
z-plane poles, in
the frequency range of $ 0 \le 2\pi k/N \le \pi /2 $ radians (zero to one quarter of the input sample rate in Hz), when
5-bit binary arithmetic is used.

which is the z-domain transfer function of the
standard Goertzel filter. The Eq. (4) Goertzel filter's block diagram, shown in
Figure in Figure 3(a), has two conjugate poles and a single pole-canceling zero
on the
z-plane as shown in Figure
3(b).

Misconception:
Because the first ratio in Eq. (2), which is only marginally
stable, times a unity-valued ratio equals Eq. (4), then the
implementation of Eq. (4) must also be only marginally stable.

I'm familiar with this misunderstanding because I was painfully guilty of it myself years ago -- and this is the misconception stated on the Goertzel Wikipedia web site [2].

The Stability of a Goertzel Filter
As it turns out, a Goertzel filter, described by Eq. (4), is guaranteed stable! The filter's poles always lie on the z-plane's unit circle. (Proof of that statement is given in Appendix A.)

So, we can say:

The poles of the Figure 3(a) Goertzel filter always lie
exactly on the z-plane's unit circle making the filter
guaranteed stable regardless of the numerical precision
of its coefficients.

We show this Goertzel filter stability behavior in Figure 4 for various first quadrant z-plane poles for 5-bit binary filter coefficients.

The Limitations of the Goertzel Filter
In finishing this blog I'll mention two important characteristics of a Goertzel filter.

[1] While the Goertzel filter is guaranteed stable, quantized coefficients (a finite number of binary bits) limit the locations where the filter's z-plane poles can be placed on the unit circle. A small number of binary bits used to represent the filter's coefficients limits the precision with which we can specify the filter's resonant frequency.

[2] Quantized coefficients make it difficult to build Goertzel filters that have low resonant frequencies. The lower the desired filter resonent frequency the larger must be the number of coefficient binary bits used.

We illustrate both of those characteristics in Figure 5 and Figure 6 for 5- and 8-bit coefficients.

Conclusion
Here we showed that Goertzel filters will never be unstable due to z-domain pole placement on the z-plane's unit circle. However as pointed out by Ollie Niemitalo, another form of instability is possible with Goertzel filters. A Goertzel filter can become unstable, due to finite precision arithmetic being used for the first feedback accumulator (adder) in Figure 3 to produce the v(n) sequence, for very large values of N. So the user must take care to ensure that, for large N (in the tens of thousands for 32-bit floating point numbers), all binary register word widths can accommodate all necessary arithmetic operations.

APPENDIX A – Proof of Goertzel Filter Guaranteed Stability
Goertzel filters, described by Eq. (4), are guaranteed stable. Proof of that statement was shown to me by the skilled DSP engineer Pavel Rajmic (Brno University of Technology, Czech Republic). Rajmic's straightforward proof of Goertzel filter stability goes like this:

Comments:

Excellent article! This has definitely sharpened my understanding of Geortzel stability. But I still ponder what can be said about "long input sequences". Are they only of concern in the context of the coefficient quantization issues described above. Looking at this thread: http://www.dsprelated.com/showthread/comp.dsp/2367-3.php, there are clearly some relative accuracy issues.

[ - ]

Comment by Rick Lyons●November 12, 2015

Hi Aaron, I just now saw your Comment. I don't understand your question of "what can be said of long sequences?" If you can rephrase your question I'll do my best to answer it. I took a look at the thread you referenced in your URL address. There's no way I have the time to plow through the details of that long thread.Regards,[-Rick-]

[ - ]

Comment by Aaron45●December 5, 2015

No problem, Rick. It took me a while to see your response. Since then it looks like a fruitful discussion took place. I think Olli was on to what I was asking. Since then, I've also read alot more about the lack of numerical stability of Goertzel and long input sequences, taking the Gentleman article (http://comjnl.oxfordjournals.org/content/12/2/160.full.pdf) as a base and following citations.

I suspect that most people wouldn't be terribly concerned from a practical standpoint given modest precision and input lengths. But in my case I was dealing with very long input sequences (up to 1e9 samples!) For a while I was looking at pre-decimating with a CIC+FIR as a mitigation. But after further consideration, I realized that the Reincsh modification would always work in my application. The reason is that for very large input sequences, my frequency is well away from pi/2 (where Goertzel performs numerically better than Reinsch). When the frequency is closer to pi/2, the sample size is much smaller where rounding error accumulation is negligible.

[ - ]

Comment by Olli Niemitalo●November 12, 2015

The Goertzel filter is not numerically stable. It is not unstable due to quantization of coefficients, as you have shown. It is unstable because of accumulation of rounding errors. Let's say the resonator's coefficient for the previous output has a bit depth of 10 bits. Multiplication of an N bit number by an M bit number gives an N+M bit number. So on the 10th iteration the resonator's output needs already a bit depth of about 100 bits. We can't store that and begin to throw away bits and accumulate rounding error.

Hi Ollie. I want to understand what you're saying here. Can you tell me, what does the phrase "resonator's coefficient for the previous output" mean. I'm trying to figure out what that phrase means relative to my above Figure 3(a).Regards,[-Rick-]

[ - ]

Comment by Olli Niemitalo●November 13, 2015

Hi Rick! I mean the coefficient 2cos(2?k/N) of the recursive part of the Goertzel filter: out[t] = in[t] + 2cos(2?k/N)out[t-1] - out[t-2].

[ - ]

Comment by Rick Lyons●November 14, 2015

Hi Ollie, Ah, OK. You're bringing up both a good and an important implementation topic here. Yes, multiplication of a 10 bit number by a 10 bit number can produce a 20-bit number. But after the 10th Goertzel iteration the accumulated v(n) result may not require a full 100-bit word depth. Accumulator register overflow is not guaranteed. The required v(n) accumulator word depth depends on the magnitudes of the numerical values being summed. When implementing the Goertzel algorithm using fixed-point numbers, the user must do one of two things. (1), make sure the v(n) accumulator can accommodate the number of iterations, N, and the magnitudes of the input signal samples. Or (2), truncate the v(n) accumulator output after each accumulation. The following web page discusses the truncation scheme.http://www.embedded.com/design/configurable-systems/4008873/2/PRODUCT-HOW-TO-Efficient-Fixed-Point-Implementation-of-the-Goertzel-Algorithm-on-a-Blackfin-DSP

[ - ]

Comment by Olli Niemitalo●November 14, 2015

I did not mean overflow due to input values. Even a fixed-point implementation of the Goertzel algorithm is fed just an unit impulse followed by zeros, the bit string of the fractional part of v(n) will grow in length every iteration, by about the length of the bit string of the coefficient. So (1) is hardly realistic. I tested a 32-bit floating point implementation using the round-to-nearest mode, k=1, N=50, and got an evolution of the magnitude of y(n) that looks a lot like random walk. That is because the rounding error is effectively a pseudorandom sequence of ups and downs. Eventually it enters a loop in its finite states. See http://imgur.com/yMQl2Ex Horizontal axis: iteration count, vertical axis: magnitude of y(n).

[ - ]

Comment by Rick Lyons●November 20, 2015

Hi Olli, OK, I'm back on duty. Your plotted output curve is interesting. Can you tell me, what was the input signal you used for our testing? Also, I don't know what "round-to-nearest mode" means with regard to floating point numbers.[-Rick-]

[ - ]

Comment by Olli Niemitalo●November 25, 2015

Hi Rick, the input signal was just a unit impulse followed by zeros. Floating point rounding mode determines how a result that cannot be exactly represented is rounded. Possible modes for IEEE binary floating point are round to nearest (ties to even), round towards negative infinity, round towards zero, round towards positive infinity. See: https://en.wikipedia.org/wiki/IEEE_floating_point#Rounding_rules

I will paste my code below. It can be compiled with G++, a C++ compiler. If you direct the program's output to a file output.txt, you can plot it in Octave using x = csvread("output.txt"); plot(x(:, 1), x(:, 2));