Goertzel Algorithm and C Implementation Using the Octave GNU Tool

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…N-1 signal index and k = 0…N-1 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 (N-1) complex multiplications and (N-1) 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(N-1) (3)

Additions = 2(N-1) (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 well-known 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 2-1:

Figure 2-1: 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 2-2 for this example.

Figure 2-2: 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 2-2) 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(n-1). 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(n-1)
  

• Filter (14) requires to save the input history: x(n-1),…,x(n-N) for N values and output history: y(n-1) for one value. See the Figure 2-1

• 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 2-2: main lobe peak 34dB (f=3120Hz) and many sidelobes with high-level 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…N-1

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:

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

2. 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 3-1 and Figure 3-2. 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 3-1: Rectangular Window N = 50

 

Figure 3-2: Spectrum of the Rectangular Window N = 50

 

• Bartlett or triangular window. See the Figure 3-3 and Figure 3-4. 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 3-3: Bartlett or Triangular Window N = 50

 

Figure 3-4: Spectrum of the Bartlett Window N = 50

 

• Hamming window. See the Figure 3-5 and Figure 3-6. 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 3-5: Hamming Window N = 50

 

Figure 3-6: 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 3-7 and Figure 3-8 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 3-7: Kaiser Window β = 3.86 and N = 50

 

Figure 3-8: 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…N-1, 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 4-1: Magnitude Response of the Discrete Fourier Transform (DFT) Filter with the Hamming Window. N=50, Fs=12kHz, Fk=3.12kHz

 

Compare the Figure 2-2 and Figure 4-1:

• The level of the side lobes of the filter has significantly decreased:
  Figure 2-2: Rectangular Window +20dB
  Figure 4-1: 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 (N-1) 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 5-1:

Figure 5-1: Goertzel Filter. Direct Form 2 canonic

 

The implementation of the Goertzel algorithm:

Obviously the DFT value is got on the N-th cycle (n=N-1) of the filter: y(N-1)=X(k)
Thus the values: y(0), y(1),…, y(N-2) are not used. It means that the complex multiplication by

in (24) should be done only at the last N-th (n=N-1) 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=(N-1) 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:

1. Measurement of the input signal: x (0), x (1),…, x (N-1)

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

3. The values: z(0), z(1),…, z(N-1) are calculated using (23)

4. Calculate the spectrum value using the formula (26):
   
   

6. Conclusions

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

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

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

4. 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:

1. DFT FIR in Recursive Form:
   – Magnitude Response Plot;

2. DFT FIR using smoothing Window
   – Magnitude Response Plot;

3. Smoothing Windows: Bartlett, Hamming, Kaiser:
   – Generate the window array;
   – Magnitude Response Plot;

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

1. 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 -32767

   typedef 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 = {…};

2. Main export functions:
   

   /* Fixed Point Arithmetic */
   /*
   Fixed Point Goertzel Algorithm
   z(n)=x(n)+z(n-1)cos(2*pi*k/N)-z(n-2)
   y(N-1)=z(N-1)+z(N-2)*(-exp(-j2*pi*k/N))
   Input
   sint16* ptr_buffer – pointer on the input buffer: x[0],…,x[N-1]
   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(n-1)cos(2*pi*k/N)-z(n-2)
   y(N-1)=z(N-1)+z(N-2)*(-exp(-j2*pi*k/N))
   Input
   sint16* ptr_buffer – pointer on the input buffer: x[0],…,x[N-1]
   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 “Discrete-Time Signal Processing“, published by Pearson Education, Inc, publishing as Prentice Hall, 1999

[2] L.Rabiner, B.Gold “Theory and application of digital signal processing“, Prentice-Hall, Ing Englewood Cliffs, New Jercy 1975

[3] John G. Proakis, Dimitris K. Manolakis, “Digital Signal Processing”, Fourth Edition, Pearson Education Limited 2014, ISBN 10: 1-292-02573-5, ISBN 13: 978-1-292-02573-5