Goertzel Algorithm

Contents
Introduction
Recurrence relation to calculate the fixed signal spectral sample
Goertzel Algorithm
Example of Using Goertzel Algorithm
Using Goertzel Algorithm to Decode the DTMF signal
Conclusions
Reference
Introduction
In the previous sections we considered the Discrete Fourier Transform (DFT) and its properties, and algorithms of the Fast Fourier Transform (FFT). We found out that FFT algorithms allow to save compute resources upon calculating the DFT significantly. However sometimes it is necessary to calculate not the full DFT, but only the fixed number of spectral samples. For example, such task rises upon decoding the DTMF signal. In this case using FFT algorithms can be inefficient because it is necessary to calculate the complete DFT to receive one spectral sample.
In this section we will consider the Goertzel algorithm [1] used to calculate the fixed spectral sample set of the Discrete Fourier Transform which is widely used in tasks of decoding the DTMF signal. As it will be shown below, Goertzel algorithm can be realized in the form of the second order IIR-filter.
Recurrence relation to calculate the fixed signal spectral sample
The equation for point Discrete Fourier Transform of the signal , is of the form of:
(1)
Twiddle factors have the following property:
(2)
Thus, multiplying (1) by will not lead to changing the result:
(3)
Break down the sum (3) in:
(4)
Set for the fixed number of the DFT spectral sample. Put outside the brackets in (4) and get:
(5)
Express in terms of . For this purpose put again outside the brackets and get:
(6)
Thus, we get the recurrence relation to calculate at any step in terms of , got at the previous step for the fixed number of the spectral sample :
(7)
where is an input signal sample with the number .
This recurrence relation at the step with the number leads us to , i.e. at the step we will get the spectral sample with the number .
Having analyzed (7), pay attention that the recurrence relation can be interpreted as the difference equation of the first order IIR-filter of the complex factor the block chart of which is presented in Figure 1.
Figure 1. Block Chart of the IIR-Filter Realizing Calculation of the Spectral Sample
with the Number
Set transforms and in terms of and respectively. Then is the transform and the difference equation (7) in the operator form is equal to
(8)
The transfer function of the got IIR-filter is of the form of:
(9)
Thus, we have stated the number of the spectral sample and transformed the DFT equation for one fixed spectral sample to the recurrence relation which allows us to get the required value of the spectral sample at the step . At the same time all intermediate values of the recurrence relation does not interest us, and only causes interest.
Then we interpreted the recurrence relation as the difference equation of the IIR-filter with the transfer function (9). As a result we got the first order filter with the complex factor using of which is per iteration at the output. If is fixed the complex factor value is always constant.
Goertzel Algorithm
To calculate one spectral sample upon using the IIR-filter complex multipliers and adders are required that does not give any advantages in comparison with the direct DFT calculation according to (1).
However computational shortcut can be got upon multiplying a numerator and denominator of the IIR-filter transfer function (9) by :
(10)
Consider the sum in more detail:
(11)
Then taking into account (11), (10) can be presented as follows:
(12)
The transfer function (12) corresponds to the second order IIR-filter structure of which is shown in Figure 2 where .

Figure 2. Structure of the IIR-Filter Realizing Calculation of the Spectral Sample with the Number
Pay attention that is a real coefficient therefore multiplication in the recursive filter branch has become real, and at that, multiplication by cannot be considered. We have one complex coefficient in the transfer function numerator. But we find that intermediate values are not important to us, consequently to multiply by is possible only per iteration when is directly calculated.
Thus, using the second order filter allows to calculate for the fixed with using multipliers, not complex, but real. Also one complex multiplication by per last iteration is required.
Attention! The block chart which differs from the presented one in Figure 2 is provided in Reference [2]. The check up has shown that the second order filter presented in Reference [2] gives the wrong result at the output. The filter presented in Figure 2 gives the exact DFT sample value. You can independently make sure of the filter result validity presented in Figure 2, having executed the following test script goertzel_test.m in the Matlab system or in GNU Octave.
						
        
Warning: include(octave/goertzel_test.m): failed to open stream: Нет такого файла или каталога in /home/host1408394/dsplib.org/htdocs/en/content/goertzel/goertzel.html on line 516

Warning: include(octave/goertzel_test.m): failed to open stream: Нет такого файла или каталога in /home/host1408394/dsplib.org/htdocs/en/content/goertzel/goertzel.html on line 516

Warning: include(): Failed opening 'octave/goertzel_test.m' for inclusion (include_path='.:/usr/local/php/php-5.6/lib/php') in /home/host1408394/dsplib.org/htdocs/en/content/goertzel/goertzel.html on line 516
The wrong structure is given in the package reference MATLAB (though the built-in package function returns the right result), and also in other sources.
Thus, the spectral sample is equal to:
(13)
where are intermediate values which are calculated iteratively:
(14)
Example of Using Goertzel Algorithm
Consider an example of using Goertzel algorithm. Let , the input signal is shown in Figure 3.
Figure 3. Input Signal
Calculate the spectral sample with the number by means of Goertzel algorithm.
Calculate coefficients and beforehand:
(15)
Set as initial calculation specifications.
Calculate , according to (14) iteratively:
(16)
>
Then, according to (13) it is possible to get values for the spectral sample :
(17)
Using Goertzel Algorithm to Decode the DTMF signal
DTMF (Dual-Tone Multi-Frequency) signal represents a dual-tone signal which is used for dialing a number of digital telephone systems.
Each DTMF symbol is encoded with the help of two tones as it is shown in Figure 4. One tone is taken out of the lower group of frequencies, and the second one is taken out of the upper group.
Figure 4. Coding DTMF Symbols
Thus, for DTMF decoding we need to calculate eight spectral DFT samples with indexes which correspond to frequencies of eight tone signals. In case of sampling kHz, for point DPF (), numbers of spectral samples are given in Table 1.

1820222431343842
, (Hz)6977708529411209133614771633
Table 1. Indexes of Spectral Samples Corresponding to DTMF Tones kHz and
Thus, for DTMF decoding we need to calculate only 8 spectral samples instead of calculating point DPF.
The time domain waveform got as a result of DTMF simulating program operation is shown in Figure 5.
Figure 5. Domain Waveform
Spectral samples magnitudes calculated by Goertzel algorithm function for the indexes which correspond to DTMF frequencies (see Table 1) are shown in Figure 6.
Figure 6. DTMF Spectral Samples Magnitudes Calculated by Means of Goertzel Algorithm
Conclusions
In this section we considered Goertzel algorithm used to calculate separate spectral samples by means of the IIR-filter. This algorithm allows to reach high efficiency upon calculating a small amount of spectral samples. At the same time Goertzel algorithm falls short of productivity to algorithms of the Fast Fourier Transform to calculate the full DFT.
Goertzel algorithm finds application upon decoding tone dialing DTMF signals of a digital telephone system in view of analysing eight spectral samples magnitudes.
Reference
[1] Goertzel Gerald An Algorithm for the Evaluation of Finite Trigonometric Series. The American Mathematical Monthly, Vol. 65, No. 1, Jan., 1958, pp. 34-35

[2] Oppenheim Alan V. and Schafer Ronald W. Discrete-Time Signal Processing Second Edition. Prentice-Hall, New Jersey, 1999.

Appendix
GNU Octave (Matlab) and Python simulation scripts:
						
		
						
% N-точечное ДПФ
N = 100;
% Номер спектрального отсчета (бина)
k = 32;

%% Исходный сигнал и расчет его ДПФ    
s = sin(2*pi*(0:N-1)*k/N+pi/6);
X = fft(s);
S_dft = X(k+1); % k+1 потому что Matlab и GNU Octave индексируют массивы начиная с единицы 


%%  Алгоритм Герцеля с DSPLIB.ORG. Формула (12) и рисунок 2.
%
%                W - z^-1
% H(z) = ---------------------------
%          1 - alpha * z^-1 + z^-2
%
%
%	W = W_N^(-k) = exp(2i*pi*k/N).
%	alpha = 2*cos(2*pi*k/N).
%

w = exp(2i*pi*k/N);

%alpha
alpha = 2*cos(2*pi*k/N);

% Выход БИХ-фильтра с DSPLIB.ORG 
b_dsplib = [w, -1, 0];
a_dsplib = [1, -alpha, 1];
X = filter(b_dsplib, a_dsplib, s);
S_dsplib = X(end);


%% Алгоритм Герцеля из книги  
% А. Оппенгейм, Р. Шафер Цифровая обработка сигналов. Москва, Техносфера, 2006.
% страница 634 формула (9.9) рисунок 9.2
%
%              1 - W * z^-1
% H(z) = ---------------------------
%          1 - alpha * z^-1 + z^-2
%
%
%	W = W_N^(k) = exp(-2i*pi*k/N).
%	alpha = 2*cos(2*pi*k/N).
%

w = exp(-2i*pi*k/N);

%alpha
alpha = 2*cos(2*pi*k/N);

b_opp = [1, -w,     0];
a_opp = [1, -alpha, 1];

X = filter(b_opp, a_opp, s);
S_opp = X(end);


%%  PRINT 
fprintf('DFT ....................................%9.4f %+9.4f i\n', real(S_dft), imag(S_dft));
fprintf('dspib.org Algorithm.....................%9.4f %+9.4f i\n', real(S_dsplib), imag(S_dsplib));
fprintf('Oppenheim Algorithm.....................%9.4f %+9.4f i\n', real(S_opp), imag(S_opp));
		


	
					
# -*- coding: utf-8 -*-

import numpy as np
import scipy.signal as signal
						
# N-точечное ДПФ
N = 100
# Номер спектрального отсчета (бина)
k = 32

# Исходный сигнал и расчет его ДПФ 
t = np.linspace(0, N, N, endpoint = False, dtype = 'float64') 
s = np.sin(2*np.pi*t*k/N+np.pi/6)
X = np.fft.fft(s)
S = X[k]

"""
Алгоритм Герцеля с DSPLIB.ORG. Формула (12) и рисунок 2.

               W - z^-1
H(z) = ---------------------------
         1 - alpha * z^-1 + z^-2


W = W_N^(-k) = exp(2i*pi*k/N).
alpha = 2*cos(2*pi*k/N).
"""

w = np.exp(2j * np.pi * k/N);

#alpha
alpha = 2*np.cos(2*np.pi*k/N);

# Выход БИХ-фильтра с DSPLIB.ORG 
b0 = np.array([w, -1, 0],     dtype = 'complex128')
a0 = np.array([1, -alpha, 1], dtype = 'complex128')
X = signal.lfilter(b0, a0, s)
Y = X[N-1]

"""
Алгоритм Герцеля из книги  
А. Оппенгейм, Р. Шафер Цифровая обработка сигналов. Москва, Техносфера, 2006.
страница 634 формула (9.9) рисунок 9.2

             1 - W * z^-1
H(z) = ---------------------------
         1 - alpha * z^-1 + z^-2

W = W_N^(k) = exp(-2i*pi*k/N).
alpha = 2*cos(2*pi*k/N).
"""

w = np.exp(-2j * np.pi * k/N);

#alpha
alpha = 2*np.cos(2*np.pi*k/N);

b1 = np.array([1, -w,     0],     dtype = 'complex128')
a1 = np.array([1, -alpha, 1],     dtype = 'complex128')

X = signal.lfilter(b1, a1, s)
Z = X[N-1]


# PRINT 
print('FFT......................', S, 'модуль', np.abs(S))
print('DSPLIB.ORG...............', Y, 'модуль', np.abs(Y))
print('А. Оппенгейм.............', Z, 'модуль', np.abs(Z))

	
					
# -*- coding: utf-8 -*-
"""
Created on Sat Aug 26 15:19:12 2017

@author: sergey
"""

import numpy as np 
import scipy.signal as signal 


def goertzel(s, ind):
    """
    Функция возвращает занчения спектральных отсчетов ДПФ,
    рассчитанные с использованием алгоримта Герцеля для
    входного сигнала, заданного параметром s.
    
    Входные параметры
        s   - вектор значений входного сигнала
        ind - вектор индексов отсчетов ДПФ (индексация начинается с 0)
    
    Выходные параметры
        S   - Вектор спектральных отсчетов ДПФ, соответсвующих индексам ind
    
    Author: Sergey Bakhurin (dsplib.org)
    """
    N = len(s)
    NF = float(N)
    S = np.zeros(len(ind), dtype='complex128')
    for k in range(len(ind)):
        w = np.exp(2j * np.pi * float(ind[k]) / NF)
        alpha = 2.0 * np.cos(2.0 * np.pi * float(ind[k])/NF)
        b = np.array([w, -1.0, 0.0])
        a = np.array([1.0, -alpha, 1.0])
        X = signal.lfilter(b, a, s)
        S[k] = X[N-1]
    return S


# Частота дискретизации DTMF
Fs = 8000.0

# Количество отсчетов DTMF сигнала
N  = 205

# Вектор частот DTMF тонов
frq = np.array([697, 770, 852, 941, 1209, 1336, 1477, 1633], dtype = 'float64')

# Вектор индексов спектральных отсчетов 
ind = np.array([18,  20,  22,  24,  31,   34,   38,   42], dtype ='int32')

# Символы DTMF 
smb = "123a456b789cs0rd"

# Вектор моментов дискретизации DTMF
t = np.linspace(0, N, N, endpoint = False, dtype = 'float64') / Fs

# Цикл моделирования сигналов  DTMF и расчета значений спектральных отсчетов
for k in range(4):
    for m in range(4):
        
        # DTMF сигнал символа smb[k*4+m]
        s = np.sin(2.0 * np.pi * frq[k]   * t) + np.sin(2.0 * np.pi * frq[m+4] * t)
        
        # Модули спектральных отсчетов с индексами ind
        S = np.abs(goertzel(s, ind))
        
        # Сохранение временной осциллограммы в текстовый файл
        fname = "dat/dtmf_sym_%c_time.csv" % smb[k*4+m]
        np.savetxt(fname, np.transpose([t, s]), fmt="%+.9e")
        
        # Сохранение амплитуд спектральных отсчетовы в текстовый файл
        fname = "dat/dtmf_sym_%c_freq.csv" % smb[k*4+m]
        np.savetxt(fname, np.transpose([frq, S]), fmt="%+.9e")