# Goertzel Algorithm for a Non-integer Frequency Index

If you've read about the Goertzel algorithm, you know it's typically presented as an efficient way to compute an individual *k*th bin result of an *N*-point discrete Fourier transform (DFT). The integer-valued frequency index *k* is in the range of zero to *N*-1 and the standard block diagram for the Goertzel algorithm is shown in Figure 1. For example, if you want to efficiently compute just the 17th DFT bin result (output sample X_{17}) of a 64-point DFT you set integer frequency index *k* = 17 and *N* = 64, and implement the Goertzel algorithm.

I won't go into any details of how the Goertzel algorithm works because there's plenty of information on the Internet, as well as in most DSP textbooks, describing that algorithm. However, this blog shows how to implement the Goertzel algorithm for any non-integer value of frequency index *k* in the range 0 ≤ *k* < *N*.

__Using a Non-integer Frequency Index__

A straightforward generalized Goertzel algorithm, derived and presented by Pavel Rajmic and Petr Sysel in reference [1], gives us the opportunity to have non-integer values for the frequency index *k*. This yields the flexible situation where the center frequency of the spectral analysis need not be restricted to be an integer multiple of the DFT's fundamental frequency of *f*_{s}/*N*, where *f*_{s} is the *x*(*n*) input signal's sample rate in Hz. I use the word "flexible" because, for example, Matlab's goertzel() command requires that frequency index *k* must be an integer [2].

My interpretation of the Rajmic and Sysel algorithm is shown in Figure 2, where the frequency index *k* can be any real value in the range 0 ≤ *k* < *N*.

To implement this non-integer *k* Goertzel algorithm, we perform the processing shown in Figure 3.

Angle variables α and β (both measured in radians) are defined as:

α = 2π*k/N*, and (1)

β = 2π*k*(*N*-1)/*N*. (2)

Before the processing begins, we initialize the delay line contents by setting *w*1 and *w*2 to zero in Figure 3.

If *k* is an integer the Figure 3 processing computes the *k*th DFT bin result of an *N* point DFT. But it's worth noting that when *k* is a non-integer the processing in Figure 3 computes the discrete time Fourier transform (DTFT) of input *x*(*n*) for the frequency value of 2π*k/N* radians/sample, or equivalently *kf*_{s}/*N* Hz.

__Reducing Computational Workload__

If you "feel the need for speed" ([3]), there's a way to reduce the computational workload of the final complex-valued operations. Assuming the input sequence *x*(*n*) is real-valued and the –*e*^{-jα} and *e*^{-jβ} coefficients are in rectangular form, implementing the final complex-valued computations requires six real multiply and three real addition operations. I was able to reduce that computational workload by formulating the network shown in Figure 4 where *w*1 and *w*2 are nodes in the Goertzel delay line.

The coefficient values in Figure 4, precomputed and stored in memory prior to processing, are:

*a* = cos(β),

*b* = -sin(β),

*c* = sin(α)*sin(β) -cos(α)*cos(β), and

*d* = sin(2π*k*).

Implementing the real-valued computations in Figure 4 requires only four real multiply and two real addition operations.

Matlab code is provided below to demonstrate my interpretation of the 'non-integer *k*' Goertzel algorithm. Pavel Rajmic's Matlab code is available at: http://www.mathworks.com/matlabcentral/fileexchange/35103-generalized-goertzel-algorithm/content/goertzel_general_shortened.m

**References**

[1] P. Sysel and P. Rajmic, "Goertzel Algorithm Generalized to Non-integer Multiples of Fundamental Frequency", EURASIP Journal on Advances in Signal Processing, 2012, 2012:56.

[2] http://www-rohan.sdsu.edu/doc/matlab/toolbox/signal/goertzel.html

[3] A quote from Maverick (Tom Cruise) in the 1986 movie Top Gun. http://www.youtube.com/watch?v=CUpwLhZh66A

__Matlab Code__

function Xk = goertzel_non_integer_k(x,k) % Computes an N-point DFT coefficient for a % real-valued input sequence 'x' where the center % frequency of the DFT bin is k/N radians/sample. % N is the length of the input sequence 'x'. % Positive-valued frequency index 'k' need not be % an integer but must be in the range of 0 –to- N-1. % [Richard Lyons, Oct. 2013] N = length(x); Alpha = 2*pi*k/N; Beta = 2*pi*k*(N-1)/N; % Precompute network coefficients Two_cos_Alpha = 2*cos(Alpha); a = cos(Beta); b = -sin(Beta); c = sin(Alpha)*sin(Beta) -cos(Alpha)*cos(Beta); d = sin(2*pi*k); % Init. delay line contents w1 = 0; w2 = 0; for n = 1:N % Start the N-sample feedback looping w0 = x(n) + Two_cos_Alpha*w1 -w2; % Delay line data shifting w2 = w1; w1 = w0; end Xk = w1*a + w2*c + j*(w1*b +w2*d);

- Comments
- Write a Comment Select to add a comment

Dear Mr Lyons,

thanks for your great work. I implemented your code in python by using numpy. I was able to shorten the code a little bit by using an addition theorem, where

```
c = sin(Alpha)*sin(Beta)-cos(Alpha)*cos(Beta);
c = -cos(Alpha+Beta);
c = -cos(2*pi*k);
c = -cos(Alpha*N);
```

My whole code is

def goertzel(x,k): N = np.size(x) A = 2*np.pi*k/N B = 2*np.pi*k*(N-1)/N f = 2*np.cos(A) y1 = 0 y2 = 0 for n in range(N): y0 = x[n]+y1*f-y2 y2 = y1 y1 = y0 return 2*np.hypot(y1*np.cos(B)-y2*np.cos(A*N),y2*np.sin(A*N)-y1*np.sin(B))/N

I know I only calc. the magnitude instead of the whole complex number with phase argument.

Hope it helps someone else. Take care and stay healthy!

Hello Robin Leuering. I'm not familiar with numpy or python, but thanks for posting your code. Your code may indeed help other people here on this web site.

Thanks

Atul

Thank you be "being professional." I have not mentioned this work anywhere other than in this blog. If I had to *reference* my blog here, I'd use the following:

[?] R. Lyons, "Goertzel Algorithm for a Non-integer Frequency Index", dsprelated.com website blog. Available online: https://www.dsprelated.com/showarticle/495.php#comments

To ask it in other way. what is standard Goertzel computing when k is fractional? and what is this extra step of complex multiplication is doing?

Hi nilbhat.

That final complex multiply is not any kind of interpolation. It's some sort of phase shift. It shifts the phase of the complex-valued output of the rightmost adder. I don't have the time right now to reread Reference [1] and go through its mathematics to see why that final phase shift is necessary.

Dear Mr. Lyons,

This is a very nice and smooth application of the Goertzel Filter. One thing I last to understand is if the equalization by N is required afterwards. I mean, devide Xk by N.

Doing a DFT requires to normalize by the size of the window where the DFT was applied, isn't it necessary in this algorithm you wrote for solving the Goertzel Filter?. For recovering the amplitude (in the time domain) of the input signal, of course.

Thanks in advance!

Hello Iferreyro. Forgive me for my delayed reply to your question.

In using the DFT some people divide two times a single DFT output sample's magnitude (2|Xk|) by 'N' to obtain an estimate of the peak amplitude of an input sinusoid whose frequency corresponds to a given frequency index 'k'.

If you wish you can divide two times |Xk| by 'N' but such an estimate of the peak amplitude of an input sinusoid whose frequency corresponds to a given non-integer frequency index 'k' is only reasonably accurate if N is large, say N > 200.

QUICKLY: Goertzel Algorithm (GOE) for a Non-integer Frequency Index (RGOE)* DOES NOT WORK in term to get real amplitude*. Orig article is "Goertzel algorithm generalized to non-integer multiples of fundamental frequency"

DETAILS: __17th DFT bin result __(output sample X17) of a 64-point DFT in this article .. means "__integer ____Frequency Index__" and thus "**non-integer Frequency Index**" means e.g. **17.01 Hz** ... right? I am asking for such obvious things only because results from algo are not correct - so I want to make sure that e.g. 1Hz means integer and 1.01Hz is non-integer index (or multiplyer) as is explained for example here: https://www.researchgate.net/figure/Illustration-o...

Lets do quick test with single signal freq __f=17Hz__ (GOE:__standard Goertzel algo__) or **f=17.01Hz** (RGOE: Sysel-Rajmic-Goertzel algo) of a 64 samples // consider magnitude =1, phase =0 //

1/ GOE@64samples: __Standard Goertzel__ "integer Frequency Index"__f=17Hz__ return magnitude 1! NICE - as expected (changing phase has not effect at magnitude). at signal f=17.01Hz magnitude out is 0.994856 (int freq. detection = default) or 0.990096 (non-int freq detection = just for test)

2/ RGOE@64samples: **Sysel-****Rajmic-Goertzel **"**non-integer Frequency Index**" **f=17.01Hz** return magnitude 1.00004 ... whhhaatt ? expected is amplitude (magnitude) 1 .. not this magnitudes variations (also phase changes affect magnitude significantly - expected is immunity to this angle phase variations) .. so *it seems algo does not work well in term of meaning "integer vs non-integer Frequency Index***"**

REQUEST @4samples: Please help on this how to detecting non-integer freq. (try to explain numerically at 4 samples, f=1.01, magnitude=1, phase =0) .. samples =[ 0,0.99988,-0.03141,-0.99889] (FYI: second sample asin(0.03141) means 180+1.8deg thus 360+3.6deg thus f=1+0.01Hz). IF you use basic goniometry you can calc from this samples correctly f1.01, magnitude 1, phase 0 and immunity to phase changes .. but .. with discussed **Sysel**-**Rajmic-Goertzel** algo .. it failed. its pitty as it will be very usefull supplement

Hello igge.

In your third paragraph I do not understand what "consider magnitude =1, phase =0" means. Are those words describing a 64-sample input sequence to the non-integer DFT? If those words describe an x(n) input sequence, what are the amplitude values of the 64 input samples??? Your words "magnitude" and "phase" really do confuse me.

In any case igge, create an input 'x' sequence of length N = 64. Select a value for 'k'. Apply that input to my blog's MATLAB code and compute a result. Then compare your result to the "true non-integer DFT result" (in MATLAB code) of:

True_DFT_at_index_k = sum(x.*exp(-j*2*pi*n*k/N));

Thanks Rick, yes first test was at 64 samples. then i tried make it as simple as possible (4 samples is math limit here) to see how the RGOE algo works. magnitude = amplitude in this example. //GOE (Standard Goertzel int freq index //RGOE (Sysel-Rajmic Gortzel for non-int freq index)

Lets consider 4 samples single freq. 1 (integer) or 1.01Hz (non-int multiply fundamental freq) ...

inp x(__1Hz__) = [ 0,1,0,-1]

=> magnitude DFT@__1Hz__ = 1 (correct) //DFT=[0j,2j,-0j,(-0-2j)]

=> magnitude GOE@__1Hz__, RGOE@__1Hz__ = 1 (correct)

inp x(**1.01Hz**) = [ 0,0.99988,-0.03141,-0.99889]

=> magnitude DFT@1Hz = 0.999507 //DFT=[(-0.0304+0j),(0.0314+1.9988j),(-0.0324-0j),(0.0314-1.9988j)]

=> magnitude GOE@__1Hz__, RGOE@__1Hz__ = 0.999507 (just test 1Hz GOE,RGOE target at single signal 1.01Hz magnitude 1)

=> magnitude GOE@**1.01Hz **(NA, just test), RGOE@**1.01Hz** (<<) = 0.999383 (**wrong** magnitude! expected is 1, also not immunity to phase changes) << *this is proof that RGOE algo DOES NOT work* .. you could found 5 others RGOE algo similarly implemened ... the question is ... is this idea/algo: wrong/--? or correct/wrong?

Hello igge.

When *N* = 4 and *n* = [0, 1, 2, 3], your first four-sample input sine wave sequence is:

* x*1(*n*) = 1.0*sin(2*pi**n**1/*N*) = [0, 1, 0, -1].

The word "magnitude" should not be used to describe *x*1(*n*). The word "magnitude" should only be used when discussing complex-valued sequences. The input *x*1(*n*) has a "peak amplitude" of 1.

Without knowing the time period between *x*1(*n*)'s samples we do not know what is the frequency of *x*1(*n*). All that we can say about *x*1(*n*)'s frequency is that its frequency is 0.25 cycles/sample.

Your second four-sample input sine wave sequence is:

* x*2(*n*) = 1.0*sin(2*pi**n**1.01/*N*) = [0, 0.99988, -0.03141, -0.99889].

It is important for you to know that input *x*2(*n*) also has a "peak amplitude" of 1. (Look at the "1.0" amplitude factor at the beginning of the above *x*2(*n*) equation.)

igge, how are you measuring, estimating, the peak amplitude of *x*2(*n*)?

In any case, when k = 1.01, the RGOE@k=1.01 value for *x*2(*n*) is equal to:

X(k=1.01) = -0.0314 - 1.9985i.

Thanks Rick, yes, amplitude 1, freq 1.01Hz, phase 0 => 4 samples *x*2(*n*) - so we know frequency, amplitude, starting angle and so I do expect that RGOE[*x*2(*n*)]@k=1.01Hz returns values which gives me amplitude 1

(from you) RGOE@k=1.01 => X(k=1.01) = -0.0314 - 1.9985i. =>

amplitude RGOE@k=1.01 => SQRT(0.0314^2 + 1.9985i^2) = 1.9987 / 2 = 0.99937 (expected is 1 for phase 0, also try another phase to generate different wrong ampitude values) .. any idea ?

Hello igge.

OK, after all this confusion I now see what you're doing and why you used the word "magnitude" so often. The magnitude of RGOE@k=1 or RGOE@k=1.01 is NEVER equal to 1. You've been "normalizing" RGOE@k=1 and RGOE@k=1.01 by dividing them by N/2.

When the input sine wave has a peak amplitude of 'A' the magnitude of RGOE@k=1 (my variable 'X(k=1)') is ALWAYS equal to AN/2.

Now, by the very nature of the DFT for non-integer values of frequency index 'k' the magnitude of RGOE@k=1.01 (my variable 'X(k=1.01)') will fluctuate by a small amount. And the amount of magnitude fluctuations depend on both the initial phase shift of the input *and* the value of 'N'. For examples, the normalized magnitudes of the signal *x*2(*n*) = 1.0*sin(2*pi**n**1.01/*N*) are:

N |RGOE@k=1| |RGOE@k=1.01|

4 1 0.99938

8 1 0.99182

16 1 0.99041

32 1 0.99013

64 1 0.99009

128 1 0.99010

256 1 0.99011

512 1 0.99012

1024 1 0.99012

In any case, the Sysel-Rajmic Gortzel algorithm for non-integer 'k' DOES work properly! That's because the algorithm produces results equal to the correct non-integer DFT results of:

True_DFT_at_index_k = sum(x.*exp(-j*2*pi*n*k/N)).

Hello Rick - you hit exactly the point of discussion. So my and my coleague's expectations was wrong: we **expected to get exact amplitude** (and/or phase) of the signal frequency such is **1.01Hz**. Based on this expectation the RGOE algorithm does NOT works - does not fulfill request to get exactly the amplitude at non-integer freq. **Result from this discussion/investigation is ... RGOE/GOE does NOT detect exactly "amplitude" and "phase" at non-integer fundamental freq***.*

ARTICLE: Suggested RGOE at k=1.01Hz return Re+Im of "side-lobes" of the DFT freq.-spectrum .. Thus __algorithm RGOE works to identify DFT lobes__! so the Generaliation means RGOE can calc any part in the spectrum

//Note that applied GOE doing usually k=Int(k=1.01)=1 or Round() so when you remove Int() from GOE and use real k=1.01Hz then you will also get then same result like RGOE does and/or zoomed/zeropaddded DFT. Thus it seems to me that the improvement RGOE vs GOE is about 1/saving 1 loop iteration and finalize it by extra code after loop and .. 2/RGOE use real k=1.01 instead of k=Int(1.01) (... which both "inventions" some embedded engineers use anyway before RGOE release .. and from this fact i guess we have this discussion as we probably expected more then detecting DFT lobes).

so maybe its good to differentiate naming GOE/Int() and RGOE/non-Int() as GOEK (int index k only) vs GOEG (general, real: int or non-int index k)

Thanks again Rick for support!

Hello igge.

I did not understand much of your April 9th "Reply."

However, I can assure you that the Non-Integer Goertzel algorithm in my blog (Sysel-Rajmic Gortzel) does, indeed, compute correct results for non-integer 'k' DFTs.

Hello Rick - its just conclusion for me and engineers they asked me to clarify why in this example RGOE does not return real amplitude and phase. Reason is that **RGOE does DFT for non-integer k (lobes)**, ... thus

__expectation failed__

__(expectation was wrong) to calculate with__

__RGOE__(at k=1.01Hz)

__exactly__"

__amplitude__and

__phase__" from given sample stream x(k=

**1.01Hz, ammplitude=1, phase=0**) = [ 0,0.99988,-0.03141,-0.99889] ... as RGOE does DFT for "non-int k" correctly - thus we do not get exactly real input amplitude and phase.

(from you) RGOE@k=1.01 => X(k=1.01) = -0.0314 - 1.9985i.

amplitude RGOE@k=1.01 => SQRT(0.0314^2 + 1.9985i^2) = 1.9987 / 2 = 0.99937 ... __failed to get amplitde 1__

igge

Hello igge.

I am not sure if you understand what we have done here. But I will say this:

1] The traditional (standard) "Non-Integer k" discrete Fourier transform (DFT) does not compute correct frequency domain magnitude results.

2] The 'Non-Integer k" Goertzel algorithm computes a single-bin DFT magnitude result exactly equal to what is computed by the traditional DFT.

3] Therefore the 'Non-Integer k" Goertzel algorithm does not compute correct frequency domain magnitude result.

Hello Rick, your explanation is correct - and so its clear why RGOE@k=1.01Hz returns correct DFT Re+Im output (apmplitude != 1) which does not reflect "input reality: amplitude = 1 at freq 1.01Hz and phase 0"

Hello Rick, justt a tip: it looks that there is another GOE improvement called JCO-Goertzel (GOEJ) pdf

I would appreciate qualitative comparison between GOE (GOEK, GOEG) and the GOEJ :)

igge

Hello igge.

I had much trouble trying to read and understand that paper.

I don't know if their network (algorithm), given in their Figure 3, computes a correct single-bin DFT result or not. But I can say that the Figure 3 network is definitely NOT more computationally efficient than the traditional Goertzel algorithm.

it will be nice to have algorithm if this JCO GOE to compare it with GOE

I have a synthetic composite of 900 data points consisting of 4 sine waves with period of 150, 200, 240 and 270 data points. The amplitude of each since wave is the same as its period (150/150, 200/200 etc.). If I apply a standard Goertzel algorithm, it does not distinguish the periods of 240 and 270 from each other as the both periods are "too close to each other" while I have relatively small sample with just 900 data points. Both 240 and 270 periods are displayed on spectrogram as one peak with incorrect amplitude and incorrect phase.

Am I correct to assume if I apply instead correctly the generalized Goertzel algorithm or your version of it as you outlined above, I will be able to distinguish all 4 periods with correct phase and amplitude despite having just 900 data points?

Please accept my apologies if my question is stupid, but my background is not math (I am physician) so before I consider using the generalized Goertzel or your modification of it, I want to be absolutely sure I understand its limitations.

Thank you

Sinuhet

Hello Sinuhet.

I am assuming that your software modeling is using floating point arithmetic. (That is a critically important assumption.)

My floating point Goertzel modeling results *__for single sine wave__* input signals (N = 900) are:

Periods k Goertzel result Goertzel magnitude

150 6 6.7500e+04 -2.4316e-08i 6.7500e+04

200 4.5 9.0000e+04 -5.4802e-08i 9.0000e+04

240 3.75 1.0812e+05 -4.5826e+03i 1.0822e+05

270 3.333 1.1909e+05 -4.4086e+03i 1.1917e+05

How do my Goertzel results compare to your Goertzel results?

Hello Sinuhet.

Your four sinewaves are so close in frequency that when N = 900 (samples) you have the following spectrum for your composite signal:

If you collect 7200 samples of your composite signal you have the following spectrum for your composite signal:

In the above second (N = 7200) spectrum you can clearly see four individual spectral components. When a '*k*' value is an integer then that spectral component exhibits no spectral leakage. For the *Period* = 270 sine wave the '*k*' value is not an integer so that sine wave's spectrum suffers from spectral leakage.

To post reply to a comment, click on the 'reply' button attached to each comment. To post a new comment (not a reply to a comment) check out the 'Write a Comment' tab at the top of the comments.

Please login (on the right) if you already have an account on this platform.

Otherwise, please use this form to register (free) an join one of the largest online community for Electrical/Embedded/DSP/FPGA/ML engineers: