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 , 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 .
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" (), 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
 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.
 A quote from Maverick (Tom Cruise) in the 1986 movie Top Gun. http://www.youtube.com/watch?v=CUpwLhZh66A
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);
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
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.
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?
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  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.
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: