DSPRelated.com
Blogs

Goertzel Algorithm for a Non-integer Frequency Index

Rick LyonsOctober 7, 201325 comments

If you've read about the Goertzel algorithm, you know it's typically presented as an efficient way to compute an individual kth 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 X17) 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 fs/N, where fs 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].

This article is available in PDF format for easy printing

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 w1 and w2 to zero in Figure 3.

If k is an integer the Figure 3 processing computes the kth 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 kfs/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 w1 and w2 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);

[ - ]
Comment by RobinLeueringFebruary 28, 2021

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!

[ - ]
Comment by Rick LyonsFebruary 28, 2021

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.


[ - ]
Comment by atulsMarch 28, 2019
HI this is very nice work to increase the speed of processing. I just need to know that I want use it in one conference paper, and don't know how to refer it in reference section. Can you please give me any proper link if you have mentioned the same work anywhere.

Thanks

Atul

[ - ]
Comment by Rick LyonsMarch 28, 2019

Hi Atuls.

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

[ - ]
Comment by nilbhatMarch 22, 2023
Hi Rick,  How do I intuitively understand the the one extra complex multiplication at the end. Is it like interpolating?

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?

[ - ]
Comment by Rick LyonsMarch 24, 2023

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.

[ - ]
Comment by lferreyroJuly 4, 2023

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! 

[ - ]
Comment by Rick LyonsJuly 7, 2023

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.

[ - ]
Comment by iggeApril 4, 2024

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

[ - ]
Comment by Rick LyonsApril 7, 2024

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));



[ - ]
Comment by iggeApril 7, 2024

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?

[ - ]
Comment by Rick LyonsApril 8, 2024

Hello igge.

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

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

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

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

Your second four-sample input sine wave sequence is:

    x2(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 x2(n) also has a "peak amplitude" of 1. (Look at the "1.0" amplitude factor at the beginning of the above x2(n) equation.)

igge, how are you measuring, estimating, the peak amplitude of x2(n)?

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

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

[ - ]
Comment by iggeApril 8, 2024

Thanks Rick, yes, amplitude 1, freq 1.01Hz, phase 0 => 4 samples x2(n) - so we know frequency, amplitude, starting angle and so I do expect that RGOE[x2(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 ?

[ - ]
Comment by Rick LyonsApril 8, 2024

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 x2(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)).

[ - ]
Comment by iggeApril 9, 2024

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!

[ - ]
Comment by Rick LyonsApril 10, 2024

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.

[ - ]
Comment by iggeApril 11, 2024

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

[ - ]
Comment by Rick LyonsApril 11, 2024

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.


[ - ]
Comment by iggeApril 15, 2024

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"

[ - ]
Comment by iggeApril 10, 2024

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

[ - ]
Comment by Rick LyonsApril 13, 2024

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.


[ - ]
Comment by iggeApril 24, 2024

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

[ - ]
Comment by SinuhetJune 27, 2024
Dear Mr Lyons

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

[ - ]
Comment by Rick LyonsJune 28, 2024

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?


[ - ]
Comment by Rick LyonsJune 28, 2024

Hello Sinuhet.

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

sinuhet spectrum-i_50018.jpg

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

sinuhet spectrum-i_20125.jpg

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: