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 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].
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);
- 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:
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.
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 ?
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)).
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.
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:
