# 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.

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 References
 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

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 March 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 March 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 February 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 February 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 March 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 March 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  and go through its mathematics to see why that final phase shift is necessary.

[ - ]
Comment by July 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.

[ - ]
Comment by July 7, 2023 