 ## Goertzel Algorithm and C Implementation Using the Octave GNU Tool

### 2. Introduction

##### The FFT Algorithm requires: ### 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: ##### 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 HzThen 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

### 3. Windows Overview

##### 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: • ##### 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 windowIf 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

##### Figure 4-1: Magnitude Response of the Discrete Fourier Transform (DFT) Filter with the Hamming Window. N=50, Fs=12kHz, Fk=3.12kHz

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

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,…,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,…,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