Talk:Goertzel algorithm

Straying from topic

The author of the C code was obviously pleased with his efforts and wanted to show it off, but this is not really the place. And the example is far too long for the purpose of demonstrating the algorithm. If used at all, it should demonstrate just the algorithm, not the application, and this is already more appropriately done in the pseudo-code. Not all readers will be C programmers, it should be more generally applicable. The bulk of the topic is a discussion of DTMF which is hardly the point nor the only application. Only a small part of the code is in fact related to the Goertzel algorithm. I'd remove it except that it might look like vandalism. This is not the place to publish school your projects ;) —Preceding unsigned comment added by 217.40.148.115 (talk) 12:19, 26 March 2009 (UTC)

Did the C code work? Because I tried implementing the pseudo-code in C and it does not appear to work. The returned values vary dramatically depending on the value of N samples I give the algorithm. Perhaps there should be some explanation on how to implement it properly? --24.91.253.46 (talk) 06:54, 25 April 2012 (UTC)

Funny

If you enter the word goertzel in Google, this article is the second hit. I started this article and posted the working C code. I had spent years perfecting this code and it was part of a telecom switch. It actually worked and there is the remote chance that if you have used a telephone, you may have exercised this code as you pushed the buttons on your phone. I posted the code hoping to save some fragments of knowledge from a dying industry and in hopes that others smarter than me would make comments on the code and show how such a simple code segment could do such magical work as to extract digits out of signal samples.

Well the code is now gone and is replaced by a mathematical discourse that is quite clear to the small fraction of the population who has survived several years of advanced mathematics and who probably already has several textbooks lying around covered with dust and which contain the same information as presented on this page in its current form.

I also have a stack of 10 or 20 books on signal processing purchased at various times at thrift stores for perhaps a 50 cents or a dollar a book. Which basically tells you the value of mathematical writings of this sort.

To the above poster: If perhaps you would like to make the world a better place, instead of having fun taunting those who have less education than yourself, you might consider sharing some of your knowledge in a form that matters more than those thrift store books. —Preceding unsigned comment added by 70.102.57.178 (talk) 21:49, 28 October 2010 (UTC)

Re: Funny

I am not the poster you are responding to, but I would like to respond to your comment. I sympathize with your conclusion, having struggled through the math in so many Wikipedia articles on DSP. How many times I could have screamed - just give me the code! In the end, however, it appears that the policy of Wikipedia is to just present the high level math. Bit by bit, I am learning to understand that math, though it has taken me years, without a proper academic founding in the subject.

As to the current article, I would like to commend the author(s) for doing an excellent job. After reading many times over several days, I finally understand it in its entirety. The author(s) have provided me with a wealth of insights on the topic.

In conclusion, I would like to ask Wikipedia to consider that pseudocode should be attached to every article on DSP such as this, and the article not be considered complete without it. TropicalCoder (talk) 15:49, 7 July 2015 (UTC)

How does this work?

I see sample code in C but nothing in English that describes the algorithm. --Damian Yerrick (☎) 19:30, 25 March 2006 (UTC)

There have been numerous edits that have fleshed out the article since the tag was added so I removed it. Feel free to add it back, but please leave some specific suggestions as to how the article could be improved. Lunch 02:36, 21 November 2006 (UTC)

I have used the C code and have found that the number of samples used can have a signeficant effect on the accuracy of the filter. I would like to know more about the practical selection of the value for GOERTZEL_N. --R Parker ANU 00:32, 14 February 2007 (UTC)

Problems with the article

Here are some things that should be fixed, but which I don't have time to fix now:

• Goertzel can be used to compute any frequency output of the DTFT (for a finite/windowed signal), not just the discrete "bins" of the DFT.
• The operation counts cited for "the" FFT are for radix-2 Cooley-Tukey only (and are only the leading N log2 N term); it is incorrect to attribute them to "the" FFT algorithm in general (there are many distinct FFT algorithms).
• The question of comparing Goertzel to FFT algorithms is made more complicated than the article suggests by the existence of "pruned" FFT algorithms, which compute M outputs in O(N log M) time instead of Θ(N log N), and may be faster than Goertzel (in my experience) for M as small as 10. See http://www.fftw.org/pruned.html (for example).
• The article intro should clearly state that Goertzel computes M outputs of a DFT (or DTFT) from N samples in Θ(N M) time, the same as evaluating the DFT or DTFT formula directly, and that the advantage of Goertzel over the naive DTFT formula lies in the reduction of the constant coefficient.
• The article should also point out that Goertzel's computational advantages come at a price: it is more susceptible to roundoff error than the DTFT formula, and much more susceptible than a typical FFT (that is, its roundoff errors grow faster with N). (This is well-known in the literature; see this Usenet thread for references, more precise statements, and some example numbers.)

—Steven G. Johnson 05:45, 3 January 2007 (UTC)

I'm a bit skeptical about the code line of

  power = sprev2*sprev2 + sprev*sprev - coeff*sprev*sprev2 ;

That actually seems like the square of the real component and does not include the imaginary component. Perhaps it should be something more like

 power = sprev2*sprev2 + sprev*sprev - coeff*sprev*sprev2 + wi*sprev2*sprev2 ;

It also does not seem to correctly scale the result by a factor of ( 2 / Nterms). As a stawman, consider replacing that line with

 i = sprev - sprev2 * coef * 2 / Nterms
 q = sprev2 * wi * 2 / Nterms
 power = i*i + q*q

I have not checked that is correct but it would be nice if someone else looked at see if they agree

Cullenfluffyjennings (talk) 23:11, 16 August 2013 (UTC)

This:

  power = sprev2*sprev2 + sprev*sprev - coeff*sprev*sprev2 ;

Is indeed correct. If you simplify this:

  sqrt((cr * sprev - sprev2)^2 + (ci*sprev)^2) 

you get the above answer, obviously with coeff equalling cr*2

The real question is why the imaginary component is:

  ci*sprev

Instead of

  ci*sprev - sprev2

At the moment, I'm getting nothing but strange results no matter how I tweak the algorithm. Oh well.

09:26, 9 August 2014 (UTC)124.190.22.24 (talk) 09:26, 9 August 2014 (UTC)

Algorithm in reverse?

Is there a good external link that could be added (or academic reference) that explains how the algorithm is used in practice in reverse, say to generate DTMF tones? I have some DTMF generation code (the DtmfGenerator code by Plyashkevich Viatcheslav that's floating around in C and Java - don't know where it originally came from) that I think uses Goertzel in reverse for generating as well as recognizing the DTMF signal. But it doesn't explicitly say that, and the math is way over my head to try to figure it out from the code alone - it uses a transform function that goes by the name MPY48SR as part of the algorithm, and looks like it iteratively modifies and shifts a matrix of 6 values. That's all I can see from the code, anyway.. Anyway, anybody know a reference for theory and practice of reversing the algorithm for generation? — Preceding unsigned comment added by Jimw338 (talk • contribs) 20:55, 1 April 2012 (UTC)

Response: Algorithm in reverse

It seems silly to me to employ a Goertzel "in reverse" - whatever that means, when all that is required to generate a sine wave is magnitude * sin(current_phase) for each sample. It doesn't get any simpler than that. I suppose on limited embedded applications one would prefer not to require computing the sine, or even not having to use floating point. If the frequency is an integral of the sampling rate, one could simply save the values for one period in ROM and spit them out - as many periods and a fraction thereof as required at run time. Or save half a period in ROM and change the sign for the second half of the period. Or save a quarter of a period and reverse that as required and add the sign to reconstitute the whole period. TropicalCoder (talk) 15:38, 7 July 2015 (UTC)

Intent to rewrite

After referring to this article to see how well it covered the most important cases of applying the algorithm to real-valued data — and being disappointed — I continued studying why the article presentation seemed uneven. Here is my ill-considered list of problems, inaccuracies, and presentation glitches. My notations, if you wish to track details, are "Paragraph within section" counting parts separated by equation lines as if they were separate paragraphs, and "Sentence" within each paragraph.

Intro

• P1S1 ambiguous grammar.
• P1S2 misses the fundamental relationship that Goertzel algorithm and FFT algorithms are both means of computing terms of the DFT. Discussions should centre on the DFT, since this is what matters the most.
• The introduction misses the important facts that the implementation code is small and well suited to small processors, and friendly to real-valued arithmetic. That might be more important than anything else in the introductory section.

Explanation

• P1 Presentation is fragmented. The algorithm produces y(n) given an input sequence x(n). The s(n) terms are only intermediate. The entire filter should be presented up front, not just the IIR stage.
• P2S3 It is false that omega parameter "should be" less than 1/2. That is entirely dependent on the application. While there are reasons related to sampling theory and the Nyquist limit for paying attention to this, the algorithm works just fine thank you for other values, even negative ones. Case in point: half of the calculations for a complex-valued DFT violate the stated "should be" rule.
• P3 "For real inputs" is extraneous to the development of the algorithm and belongs later in the discussion of real-valued implementations. Unfortunately, that topic is essentially missing from the rest of the article.
• P4 "Applying an additional FIR transform" — this will no longer seem plucked out of the air after P1 is improved. Then both of the Z transforms can be introduced together and dispensed with at one time.
• P8 "...or, the equation for the (n+1) sample DFT sequence of x." This, unforunately, is seriously wrong. A simple change of notation that will make this clear. Replace w with K/N, where K is the DFT term index, and where N is the number of terms in the DFT input sequence. This will reveal (by comparing with the DFT definition) that this is NOT the expression for a N+1 term DFT, because then the frequency terms would have the form K/(N+1).The upper limit on the summation is not consistent with the definition of a DFT. This can become a serious problem if the algorithm is not terminated correctly.
• P9 "which corresponds to the chosen frequency as..." Unfortunately, the equation isn't quite correct. Comparing to the immediately preceding equation, the value of y(n) should be used, otherwise the algorithm has not run to completion. However, there are conditions under which using the next-to-the-last term is valid... not discussed.
• P11 "power does not depend on the phase factor" True... but if the presentation is done in a cleaner way, you don't need to think about it, because at least for DFT calculations, all of the extra phase terms have reduced to constant 1 and drop out, which is not obvious.
• Though the theory is there, this article does not cover the important application details of applying the Goertzel algorithm to compute the value of a DFT term, in the same way that it does signal power. This is somewhat tricky (see P8 above). Wikipedia could do a great service by pointing out where implementations can easily go wrong.
• In general, there is considerable notational ambiguity between n, the current term that the filter is calculating, and N, the particular term at which the filter finishes calculating.

Implementation

• P1 "When implemented in a general purpose processor..." Unclear. If you don't save the filter states one way or another, the filtering won't work on any kind of processor.
• The fact that the FIR filter is not evaluated until the last step is assumed but never mentioned specifically. This is actually a major feature of the algorithm, contributing to is efficiency, and this omission is a strange oversight.
• The particular choice of a direct-form II IIR filter for the first stage is an important feature of the algorithm, because in general, the internal state representation of an IIR filter is not unique, and would not necessarily preserve the required history terms for the FIR filter to use. Not stating this explicitly is an oversight.

Complexity

• I note that some of Steven Johnson's points were not addressed and should be.
• My practical experience is that on some processors, the math is cheap but operations in external memory (i.e. large code, lookup tables, things that don't fit into registers and level-1 cache) suffer from a large latency penalty. For this reason, the complexity order predictions can underrate the Goertzel algorithm performance compared to an FFT.
• "For real-valued input sequences, the amount of computation is halved in both cases." Well, this is approximately true, but only for the case of FFT variants specialized for real-valued data.
• "must compute all N bins in parallel" Order of execution does not seem relevant to a complexity analysis.

Overall I am proposing a rewrite that will cover all of the issues above while discarding very little from what is currently in the article. Not ready to post, I need some numerical validation first. So past contributors, do not panic. Your contributions should still be there, but maybe located in a different place, maybe with some different dressing. Mostly, a much greater emphasis on more specific and orderly presentation with respect to practical calculations of DFT terms. I have some concerns about overall bulk. ParaTechNoid (talk) 00:45, 29 October 2012 (UTC)

Posted

I have posted the rewrite, which I think offers an improved balance between theoretical accuracy and practical usefulness. But then again, nothing is perfect. The things that some might find controversial:

• "Your final equation for calculating a DFT term do not agree with reference blah blah which says..." Well, references already cited in the article didn't agree, so there was no way to please everybody. Maybe the choice boils down to what is more important: respecting the authority of references, or computing a correct result. Before complaining, you might want to download a copy of a test script I found (http://www.mstarlabs.com/kb/files/GoertFix.m , probably not suitable as an article reference), and compare whether your favourite implementation (and this Wiki article!) works or doesn't.
• In the Complexity section, I shifted emphasis from multiply-counting toward complexity-theoretical. Perhaps this section should be reduced. To calibrate complexity expressions enough to make them useful, you need to run practical benchmarks; but after you have the benchmarks, complexity theory doesn't add much.

ParaTechNoid (talk) 21:05, 18 November 2012 (UTC)

Confusion between omega and f

All the equations in this page use 2*pi*omega as the argument of trigonometric functions. This is confusing at best. They should either be written as omega (which I recommend) or 2*pi*f. Larry Doolittle (talk) 18:14, 17 October 2013 (UTC)

Citations

I do not have a signal processing background but whoever sprayed the text with Citation needed flags was obviously not qualified to do so. There is a childish level of mathematical knowledge assumed that does not need citation. The most egregious example is the citation request for equation 6 which is simply derivable from equation 5 which has no similar request. Equation 2 is in context obviously required to be this form for the algorithm to make sense and the initialisation of the state vector to zero is in some sense arbritary but if the algorithm is to be used for a DFT is obviously required to have an unbiased result. Anyone to whom this is not obvious should not be commenting. It would not make sense to require a citation for a statement concerning real numbers like 1 + 1 = 2. 109.239.88.226 (talk) 11:21, 7 March 2014 (UTC)AJ

Seconding the above. This "citation needed" business is plain silly and induces facepalm. There is a T101 in your kitchen (talk) 17:41, 24 November 2014 (UTC)

Thirding the above. My exact same thought (including the example!). Can we remove those "citation needed" flags now? — Preceding unsigned comment added by 213.123.251.51 (talk) 09:45, 10 February 2015 (UTC)

• Yes, the cn tags are over done, but the section misuses big-O notation. I'd trim some content (no K in big-O notation) and drop the cn tags on reasonable statements. Glrx (talk) 06:39, 14 February 2015 (UTC)

• I have removed the first 3 of the 4 "citation needed" flags, as I have spent considerable time studying this page and agree with the comments above that these should be removed. One remains in the Z-transform section. Somebody else can decide about that.TropicalCoder (talk) 00:59, 8 July 2015 (UTC)

• I took out all the rest of the citation needed tags, and replaced them with one 'refimprove' tag for that section. –jacobolus (t) 05:42, 13 October 2015 (UTC)

Comparing the Goertzel with the DFT

I have added the section "Comparing the Goertzel with the DFT". This also demonstrates obtaining the phase and magnitudes as we would expect to receive them from a DFT. Though an embedded implementation of the Goertzel will likely have no interest in this kind of information, one certainly needs it to verify their code. My C++ implementation of the pseudo code listing that I linked to in the External Links section has been rigorously tested with dozens of frequencies, each with 16 different initial phase offsets. Worst case errors recorded were in the range of 250 PPM for phase and 30 PPM for magnitude in the case of the Goertzel, whereas it was only 5 PPM for phase and 8 PPM for magnitude for the DFT. Where the number of samples in a period was a whole number, the errors were pretty close to 0 for either algorithm. TropicalCoder (talk) 00:34, 11 July 2015 (UTC)

Incomplete

The article, as written, fails to address the most important point: connecting the algorithm back to its intended purpose. The article never relates the results of the transform y back to the DFT terms that are being sought. Nor does it explain how the final summation is to be used. As a result, a reader unfamiliar with the topic is unable to get almost any useful information from the article as written. Can someone familiar with the topic add this important information? — Preceding unsigned comment added by Wmagro (talk • contribs) 13:09, 31 October 2015 (UTC)

Directed at User: Glrx

I see it was tagged a possible conflict of interest. How did you come to that conclusion? Lakeweb.net is a personal website and all the content is open source. Please respond. http://lakeweb.net/MSP430/PrivateLine/Goertzel_Algorithm.html Thanks, Dan (talk) 19:42, 12 July 2016 (UTC)

Advertising your personal project on WP is WP:COI. In addition, the web page says little (it points to WP as a source). It is not a reliable source. Glrx (talk) 20:14, 12 July 2016 (UTC)

I am not offering something about myself, I am not advertising. The two other links do not offer any complete working example. In the VS solution there are three working and very usable projects. You should take a look at it.Dan (talk) 23:16, 12 July 2016 (UTC)

That really doesn't look like an appropriate EL; just a personal project page. Dicklyon (talk) 02:45, 13 July 2016 (UTC)

Ok, If it is about how it 'looks' rather than content... But as far as I know it is the only source of a complete working set of code to test the algorithm, develop the target code, and generate a target header of coefficients. I won't bother you with this anymore. Dan (talk) 16:23, 13 July 2016 (UTC)

Reuse of index variable in DFT section for a different job is confusing

The section "DFT computations" reuses $k$ for the "bin number" of the discretized frequency, swapping in $n$ without warning for the resulting sequence's sequence index in equation (9). This is tremendously confusing if you're following along with the derivation, as it leads you to believe that the frequency bin number should be indexed along with the sum, as seen in this Math.SE question.

I propose using $l$ or some other letter instead of $k$. An edit that does this looks like the below:

— Preceding unsigned comment added by 107.3.128.7 (talk • contribs) 05:06, 13 October 2016 (UTC)

DFT computations

For the important case of computing a DFT term, the following special restrictions are applied.

• The filtering terminates at index $n=N$, where $N$ is the number of terms in the input sequence of the DFT.
• The frequencies chosen for the Goertzel analysis are restricted to the special form

$\omega _{0}=2\pi {\frac {l}{N}}.$ (7)

• The index number $l$ indicating the selected "frequency bin" of the DFT can be any of the set of index numbers over which the sequence is nonzero (the $k$ term from equation (6)):

$l\in \{0,1,2,...,N-1\}.$ (8)

Making these substitutions into equation (6) and observing that the term $e^{+j2\pi l}=1$, equation (6) then takes the following form:

$y[N]=\sum _{n=0}^{N}x[n]e^{-j2\pi {\frac {nl}{N}}}.$ (9)

— Preceding unsigned comment added by 107.3.128.7 (talk • contribs) 05:06, 13 October 2016 (UTC)

Appalling bad article

This article is an example of everything which is wrong with Wikipedia.

Instead of giving a simple explanation and a block diagram, it launches immediately into mind-boggling maths, which will only serve to repel most readers.

Presumably this happens is because the authors want to show how clever they are, but in actual fact have no real grasp on the subject at all.

The Goertzel algorithm is an evolution of the age-old Super-Regenerative radio receiver, as invented by Armstrong back in 1913.

see https://en.wikipedia.org/wiki/Regenerative_circuit

It uses an Oscillator which is keyed on for a short time by a Quench signal. When the Oscillator first starts, the Q of the tuned circuit causes the amplitude to build slowly.

However, if there is a signal present (at the right frequency) the output builds more quickly. If the Oscillator is stopped before it saturates, the amplitude remaining is proportional to the incoming frequency.

In the Goertzel algorithm, exactly the same thing happens, except it uses a simple IIR oscillator.

— Preceding unsigned comment added by Gutta Percha (talk • contribs) 07:45, 10 March 2017 (UTC)

What a clever relationship you have noticed and pointed out! If you can find a source about that, do add it to the article. Actually, though, the Goertzel algorithm's filter is infinite-Q, or conditionally stable, while the filters used in super-regeneration are actually unstable (unlike simple regenerative, which is finite Q and stable). So the Goertzel algorithm can be viewed as being on the exact border between regenerative and super-regenerative. These do serve a similar purpose, as you point out. You can probably find a second-order IIR filter diagram and specialize it for this article if you like, to show that your motivations are indeed more pure than all us other wikipedia editors, who just like to show off and make things look complicated when they're not. Dicklyon (talk) 18:15, 19 March 2017 (UTC)

Download to PDF drops equations

Just an observation from a basic user of Wikipedia. The PDF converter doesn't pick up (many of) the equations in the article, rendering blank spaces in the PDF instead of the equations.

Maybe the format of the equations could be changed to something the PDF converter recognizes?

Best regards to the community. — Preceding unsigned comment added by 130.65.41.187 (talk) 00:04, 11 April 2017 (UTC)

Application section ci is not used for imaginary part

In the Application section value ci is set however it is not used in the pseudcode offered. The ci value is used to compute the imaginary part which is not used in the example shown. Either the full algorithm should be shown else the ci line should be removed. SoftwareThing (talk) 22:22, 7 September 2018 (UTC)

equation given in 'the algorithm' has unidentified term

what is the j in the exponent term?

it says bla bla to the power of -jw0. i understand the w0 part, but what is j? — Preceding unsigned comment added by 101.98.178.115 (talk) 10:32, 24 July 2019 (UTC)

j is the square root of −1; it is the same as i in typical mathematics.. It is a common variable in electrical engineering and is often not identified. Functions such as sin and cos are also infrequently explained. Glrx (talk) 19:08, 13 August 2019 (UTC)

I came here to ask the same question: what does j mean? I now understand that it means the same "square root of -1" that I label i. It would be nice if the article stated this clearly. Maybe this can be addressed as part of rewriting to fix all the other difficulties which others have noted above. --Jdlh | Talk 18:45, 14 May 2020 (UTC)

Stability

Rick Lyons (author of "Understanding DSP") wrote a blog post to specifically to address a perceived mistake in this very wikipedia article, regarding the stability of goertzel filters: https://www.dsprelated.com/showarticle/796.php — Preceding unsigned comment added by 46.151.207.2 (talk) 12:55, 25 July 2023 (UTC)