Goertzel Algorithm and C Implementation Using the Octave GNU Tool
1. Abbreviation
DSP – Digital Signal Processing
DFT – Discrete Fourier Transform
IDFT – Inverse Discrete Fourier Transform
FFT – Fast Fourier Transform
FIR – Finite Impulse Response
IIR – Infinite Impulse Response
a + jb – complex number where
2. Introduction
Often in Digital Signal Processing applications we need to evaluate the Discrete Fourier Transform (DFT) for only one relative frequency ω=2πk/N using the formula:
where n = 0…N1 signal index and k = 0…N1 frequency index
In this case, the relative value of the frequency ω = 2nk/N is related to the real frequency Ω by the formula:
where Ts – sampling time
Then:
where Fs=1/Ts sampling frequency
Let’s estimate the number of operations for formula (1) if the k is fixed. It requires (N1) complex multiplications and (N1) complex additions. The complex multiplication:
(a + jb) (c + jd) = ac – bd + j(ad + bc)
requires 4 multiplications for real numbers and complex addition:
a + jb + c + jd = a+c + j(b + d)
requires 2 additions of real numbers.
Then formula (1) requires:
Multiplications = 4(N1) (3)
Additions = 2(N1) (4)
The FFT Algorithm requires:
See the “Fast Fourier Transform (FFT) and C Implementation Using the Octave GNU Tool“
Formulas (3) – (6) are showed that using the fast Fourier transform (FFT) is not optimal for calculating the spectrum for a single 2nk/N frequency. Usually a special filter is developed for the task.
Note:
In the following I will only take the number of multiplications to estimate the complexity of the algorithms
2. DFT with the FIR Filter in Recursive Form
The article “FIR Filter Implementation Using Octave GNU Tool and C Language” is described the FIR Filter in Recursive Form. The FIR realization is used by the formula (7):
The formula (7) for calculating the spectrum for only one frequency ω = 2nk/N will take the form:
Note: The factor (1/N) is not used for the spectrum calculation.
Formula (8) can be got by another way directly from (1). For it the property of the complex exponential is used:
Using (9) and (1) can be written:
where
We have got a FIR filter with a impulse response h(n), which calculates the spectrum value for the frequency 2nk/N for the last N samples of the input signal. The FIR filter has the H (z):
Next using the wellknown formula for the sum of the geometric progression:
Then (12) can be written in the final form:
So we got the formula (8) again. This filter is implemented:
See the Figure 21:
Figure 21: Discrete Fourier Transform with FIR Recursive Filter
Example:
N = 50, Fs=12000Hz=12kHz, k=13 => Fk = Fs k / N = 12000Hz*13/50 = 3120 Hz
Then the DFT can be calculated for the frequency 3120Hz using a filter:
The magnitude response of the filter is showed Figure 22 for this example.
Figure 22: Magnitude Response of the Discrete Fourier Transform (DFT) Filter with the Rectangular Window. N=50, Fs=12kHz, Fk=3.12kHz
Сonclusions:

FIR filter in recursive form (14) amplifies the input signal by a N factor for 2nk/N frequency. In our example (Figure 22) the magnitude response peak at the 3.12 kHz frequency has:

FIR filter in recursive form (14) has a pole on the unit circle exp (j2nk/N). The pole is compensated by zero from numerator
therefore the filter is stable. Really calculations have a limited accuracy and mistakes will be accumulated in the history y(n1). As result the filter will have an infinite impulse response or the filter will become unstable after some time. Thus the filter requires periodic reinitializations: clear the history y(n1) 
Filter (14) requires to save the input history: x(n1),…,x(nN) for N values and output history: y(n1) for one value. See the Figure 21

The algorithm complexity of the filter is determined by 4N multiplications (or N complex multiplications). The first DFT value is got after the first N cycles. Then after each next filter cycle we will get the DFT for the previous N samples of the input signal

The DFT calculation is done only for a small fragment of the input signal, which results in Gibbs phenomenon. Spectrum is got ripples (distortions) in the frequency response. See the Figure 22: main lobe peak 34dB (f=3120Hz) and many sidelobes with highlevel peaks (up to 20dB). To minimize the effect is used special smooth windows for the input signal. The windows: Bartlett (triangular), Blackman, Hamming, Hanning, Kaiser, Lanczos, Turkey are used in DSP. I use mostly Hamming or Kaiser windows in my practice.
The recursive FIR filter (14) is used the rectangular window and Gibbs phenomenon is maximal in the case
3. Windows Overview
The signal fragment can be written:
where x(n) full input signal, w(n) window with the length N: n=0…N1
If the signal x(n) has the spectrum X(exp(jω)), and the window w (n) has the spectrum W( exp(jω)) then the signal fragment
has the spectrum as convolution X and W:
By analyzing formula (18) it is possible to determine the requirements for the window spectrum in order to minimize the Gibbs phenomenon:

The entire spectrum of the window should be concentrated in a narrow frequency band around the maximum at frequency f = 0

Minimize the number of side lobes and their peaks
Thus choosing the optimal window is a mathematical task. There are many designed windows. Let’s look at a few of them.

Rectangular window. See the Figure 31 and Figure 32. Many side lobes are ocured on the spectrum. The maximum side lobe reaches the 13dB level. The width of the main lobe is 4π/(N+1). This is the worst window. The formulas (1), (15) is used the rectangular window
Figure 31: Rectangular Window N = 50
Figure 32: Spectrum of the Rectangular Window N = 50

Bartlett or triangular window. See the Figure 33 and Figure 34. This is the simplest window which has a triangular shape. The window is easily implemented in practice even when the system has limitations in computing capabilities. The maximum side lobe reaches the 26dB level. The width of the main lobe is 8π/N. Octave GNU Tool has a function: bartlett(N), which calculates the window with the N length
Figure 33: Bartlett or Triangular Window N = 50
Figure 34: Spectrum of the Bartlett Window N = 50

Hamming window. See the Figure 35 and Figure 36. The maximum side lobe reaches the 42dB level. The width of the main lobe is 8π/N. The Octave GNU Tool has a function: hamming(N), which calculates the window with the N length
Figure 35: Hamming Window N = 50
Figure 36: Spectrum of the Hamming Window N = 50

Kaiser window. This is a universal window with the extra parameter β. By selecting this option we can get windows with any characteristics:
β = 0 according to Rectangular window
β = 1.33 according to Bartlett window
β = 3.86 according to Hanning window
β = 4.86 according to Hamming window
If the β is increased then peaks of the side lobes are decreased (good result), but the width of the main lobe is increased (bad result). Thus compromise shall be found a between reducing the side lobes and limiting the width of the main lobe.
The Figure 37 and Figure 38 are shown an Kaiser window example with β = 3.86. The peak side lobe reaches the 30dB level and the width of the main lobe is 8π/N. This window corresponds to the characteristics of the Hanning window.
The Octave GNU Tool has a function: kaiser(N, β), which calculates the window with the N length.
Figure 37: Kaiser Window β = 3.86 and N = 50
Figure 38: Spectrum of the Kaiser Window β = 3.86 and N = 50
4. DFT with the FIR Filter Using Smoothing Window
Consider the DFT using FIR filter with impulse response (11) and one Window from previous section. Then the impulse response can be written:
where n=0…N1, w(n) – Using Window
Let’s use the example from section 2:
N = 50, Fs=12kHz, k=13 => Fk = 3120 Hz
Take for example a Hamming window
Then Magnitude Response of the FIR filter (19) using Hamming window
Figure 41: Magnitude Response of the Discrete Fourier Transform (DFT) Filter with the Hamming Window. N=50, Fs=12kHz, Fk=3.12kHz
Compare the Figure 22 and Figure 41:

The level of the side lobes of the filter has significantly decreased:
Figure 22: Rectangular Window +20dB
Figure 41: Hamming Window 13dB
Thus 33dB improvement (in 45x)!!! 
At the same time, the width of the main lobe (unfortunately) is increased by about 2 times, which also corresponds to the theory.
Note:
You can implement FIR filter (19) using my article “FIR Filter Implementation Using Octave GNU Tool and C Language”
5. Goertzel Algorithm
There are N input samples of the real signal: x (0), x (1),…, x (N1) and it is necessary to get the DFT for the frequency 2nk/N with a minimum calculations. In the case the formula (14) can be rewritten in more simple form:
Multiply the numerator and denominator by
Then:
Formula (22) is described the Goertzel algorithm, which can be implemented as shown in Figure 51:
Figure 51: Goertzel Filter. Direct Form 2 canonic
The implementation of the Goertzel algorithm:
Obviously the DFT value is got on the Nth cycle (n=N1) of the filter: y(N1)=X(k)
Thus the values: y(0), y(1),…, y(N2) are not used. It means that the complex multiplication by
in (24) should be done only at the last Nth (n=N1) step.
The Goertzel algorithm is required N real multiplications by 2cos (2nk/N) in formula (23) and one complex multiplication in the last step n=(N1) in the formula (24). Thus (N+4) real multiplications are required:
Multiplications = N + 4 (25)
The Goertzel algorithm (23), (24) requires almost 4 times less multiplications compared to formula (1). See the (3) and (25).
Steps of the Goertzel algorithm:

Measurement of the input signal: x (0), x (1),…, x (N1)

The formula (17) is calculated: x(n) w(n) using one of the smoothing window
Note: If a rectangular window is used then the item is not executed 
The values: z(0), z(1),…, z(N1) are calculated using (23)

Calculate the spectrum value using the formula (26):
6. Conclusions

Using filters (14) or (19) can be calculated the DFT for a single frequency. It requires 4N real multiplications

The Goertzel (23), (24) algorithm requires only N+4 real multiplications. It is almost 4 times less as in the (1), (14), (19).

Sometimes the DFT for M frequencies shall be calculated. Which algorithm is optimal to use? FFT or Goertzel? Comparing (5) and (25), we can write down the condition for choosing the Goertzel algorithm:
M x Goertzel_Multiplications < FFT_Multiplications
Then
Thus the formula (27) is the criterion for choosing the Goertzel algorithm. If condition (27) is not met then it is more efficient to use FFT 
Smoothing window (Hamming or Kaiser or …) is effective to use for calculating the DFT
Note: these windows are still used to design FIR filters using IIR prototypes to limit the infinite impulse response.
7. Octave GNU file Goertzel_Filter.m
The Goertzel_Filter.m script prepare/support:

DFT FIR in Recursive Form:
– Magnitude Response Plot; 
DFT FIR using smoothing Window
– Magnitude Response Plot; 
Smoothing Windows: Bartlett, Hamming, Kaiser:
– Generate the window array;
– Magnitude Response Plot; 
Generate the help arrays for the C Code
– Smoothing windows;
– exp(j2πk/N);
– cos(2πk/N);
– Test signals;
8. C Language GoertzelFilter.c/h
The GoertzelFilter.c/h is supported two variants: fixed point arithmetic and float arithmetic

Defines, structures and constants:
/* DFT Point Number */
#define DFT_POINTS 50u /* N = 50 *//* Float Arithmetic */
typedef struct float_goertzel_filter_s {
uint8 isWinActive; /* FALSE – none window is used, TRUE – window is used */
float * ptr_history;
flt_goertzel_const_t const * ptr_goertzel_const;
float const * ptr_window;
} flt_goertzel_filter_t;const flt_goertzel_filter_t flt_goertzel_filter = {…};
/* Fixed Point Arithmetic */
#define LOG2_INTEGER_FACTOR 14u
#define INTEGER_FACTOR (1u<<LOG2_INTEGER_FACTOR)
/* SINT16: Max and Min Values */
#define SINT16_MAX_POSITIVE_NMB 32767
#define SINT16_MIN_NEGATIVE_NMB 32767typedef struct integer_goertzel_filter_s {
uint8 isWinActive; /* FALSE – none window is used, TRUE – window is used */
sint16 * ptr_history;
int_goertzel_const_t const * ptr_goertzel_const;
sint16 const * ptr_window;
} int_goertzel_filter_t;const int_goertzel_filter_t int_goertzel_filter = {…};

Main export functions:
/* Fixed Point Arithmetic */
/*
Fixed Point Goertzel Algorithm
z(n)=x(n)+z(n1)cos(2*pi*k/N)z(n2)
y(N1)=z(N1)+z(N2)*(exp(j2*pi*k/N))
Input
sint16* ptr_buffer – pointer on the input buffer: x[0],…,x[N1]
int_goertzel_filter_t* ptr_goertzel_filter – pointer on the Goertzel parameters
Return
int_complex_nmb_t – DFT value
*/
int_complex_nmb_t int_goertzel_algorithm(sint16* ptr_buffer, int_goertzel_filter_t const* ptr_goertzel_filter);/* Float Arithmetic */
/*
Float Goertzel Algorithm
z(n)=x(n)+z(n1)cos(2*pi*k/N)z(n2)
y(N1)=z(N1)+z(N2)*(exp(j2*pi*k/N))
Input
sint16* ptr_buffer – pointer on the input buffer: x[0],…,x[N1]
flt_goertzel_filter_t* ptr_goertzel_filter – pointer on the Goertzel parameters
Return
flt_complex_nmb_t – DFT value
*/
flt_complex_nmb_t flt_goertzel_algorithm(sint16* ptr_buffer, flt_goertzel_filter_t const* ptr_goertzel_filter);int main(void) – only test
9. Download the Goertzel_Filter.m and GoertzelFilter.c/h
You can download the files:
Goertzel_Filter.m
GoertzelFilter.c/h
with the button:
10. Literature / References
[1] A.Oppenheim, R.Schafer “DiscreteTime Signal Processing“, published by Pearson Education, Inc, publishing as Prentice Hall, 1999
[2] L.Rabiner, B.Gold “Theory and application of digital signal processing“, PrenticeHall, Ing Englewood Cliffs, New Jercy 1975
[3] John G. Proakis, Dimitris K. Manolakis, “Digital Signal Processing”, Fourth Edition, Pearson Education Limited 2014, ISBN 10: 1292025735, ISBN 13: 9781292025735