Goertzel Algorithm

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 $N-$point Discrete Fourier Transform of the signal $s(n)$, $n = 0 \ldots N-1$ is of the form of:
$$ S(k) = \sum_{n = 0}^{N-1} s(n) \cdot W_N^{n \cdot k}; $$ $$ W_N^{n \cdot k} = \exp \left( - j \cdot \frac{2 \pi}{N} \cdot n \cdot k \right); \textrm{ } k = 0 \ldots N-1. $$
Twiddle factors $W_N^{n \cdot k}$ have the following property:
$$ W_N^{- k \cdot N} = \exp \left( j \cdot \frac{2 \pi}{N} \cdot N \cdot n \right) = 1; \textrm{ } k = 0 \ldots N-1. $$
Thus, multiplying (1) by $W_N^{- k \cdot N}$ will not lead to changing the result:
$$ S(k) = W_N^{-k \cdot N} \sum_{n = 0}^{N-1} s(n) \cdot W_N^{n \cdot k} = \sum_{n = 0}^{N-1} s(n) \cdot W_N^{k \cdot (n-N)}. $$
Break down the sum (3) in:
$$ S(k) =s(0) \cdot W_N^{-k \cdot N} + s(1) \cdot W_N^{-k \cdot (N-1)} + s(2) \cdot W_N^{-k \cdot (N-2)} + \ldots $$ $$ \ldots + s(N-2) \cdot W_N^{-k \cdot 2} + s(N-1) \cdot W_N^{-k}. $$
Set $S(k) = y_{N-1}(k)$ for the fixed number $k$ of the DFT spectral sample. Put $W_N^{-k}$ outside the brackets in (4) and get:
$$ y_{N-1}(k) = W_N^{-k} \cdot \Big( s(N-1) + \underbrace{W_N^{-k} \cdot s(N-2) + \ldots + s(0) \cdot W_N^{-k \cdot (N-1)} }_{=y_{N-2}(k)}\Big) = $$ $$ = W_N^{-k} \cdot \big( s(N-1) + y_{N-2}(k)\big) . $$
Express $y_{N-2}(k)$ in terms of $y_{N-3}(k)$. For this purpose put again $W_N^{-k}$ outside the brackets and get:
$$ y_{N-2}(k) = W_N^{-k} \cdot s(N-2) + W_N^{-2\cdot k} \cdot s(N-3) + \ldots + s(0) \cdot W_N^{-k \cdot (N-1)} = $$ $$ = W_N^{-k} \cdot \Big( s(N-2) + \underbrace{W_N^{-k} \cdot \ s(N-3) + y_{N-2}(k) \ldots + s(0) \cdot W_N^{-k \cdot (N-2)}}_{=y_{N-3}(k)} \Big) . $$
Thus, we get the recurrence relation to calculate $y_{r}(k)$ at any step $r = 0 \ldots N-1$ in terms of $y_{r-1}(k)$, got at the previous step for the fixed number of the spectral sample $k$:
Figure 1. Block Chart of the IIR-Filter Realizing
Calculation of the Spectral Sample
with the Number $k$
$$ y_r(k) = W_N^{-k} \cdot \Big( s(r) + y_{r-1}(k) \Big), $$
where $s(r)$ is an input signal sample with the number $r$.
This recurrence relation at the step with the number $r = N-1$ leads us to $S(k) = y_{N-1}(k)$, i.e. at the $N-1$ step we will get the spectral sample with the number $k$.
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 $W_N^{-k}$ the block chart of which is presented in Figure 1.
Set $z-$transforms $x(r)$ and $y_r(k)$ in terms of $X(z)$ and $Y(z)$ respectively. Then $z^{-1} \cdot Y(z) $ is the $z-$transform $y_{r-1}(k)$ and the difference equation (7) in the operator form is equal to
$$ Y(z) = W_N^{-k} \cdot \left( X(z) + z^{-1} \cdot Y(z) \right). $$
The transfer function of the got IIR-filter is of the form of:
$$ H(z) = \frac{Y(z)}{X(z)} = \frac{W_N^{-k}}{1-W_N^{-k} \cdot z^{-1}}. $$
Thus, we have stated the number of the spectral sample $k$ and transformed the DFT equation for one fixed spectral sample to the recurrence relation which allows us to get the required value $k$ of the spectral sample $S(k)$ at the step $N-1$. At the same time all intermediate values of the recurrence relation does not interest us, and only $y_{N-1}(k) = S(k)$ 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 $W_N^{-k}$ using of which is $S(k)$ per $N-1$ iteration at the output. If $k$ is fixed the complex factor value $W_N^{-k}$ is always constant.
Goertzel Algorithm
To calculate one spectral sample $S(k)$ upon using the IIR-filter $N$ 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 $1-W_N^k \cdot z^{-1}$:
$$ H(z) = \frac{W_N^{-k}}{1-W_N^{-k} \cdot z^{-1}} = \frac{W_N^{-k} \cdot \left(1-W_N^k \cdot z^{-1} \right)}{ \left(1-W_N^{-k} \cdot z^{-1} \right) \cdot \left(1-W_N^k \cdot z^{-1} \right)} = \frac{W_N^{-k} - z^{-1} }{1- z^{-1} \cdot \left( W_N^{-k} + W_N^{k} \right) + z^{-2}}. $$
Consider the sum in more detail:
$$ W_N^{-k} + W_N^{k} = \exp \left( j \cdot \frac{2\pi}{N} \cdot k \right) + \exp \left( -j \cdot \frac{2\pi}{N} \cdot k \right) = 2 \cdot \cos \left( \frac{2\pi}{N} \cdot k \right). $$
Then taking into account (11), (10) can be presented as follows:
$$ H(z) = \frac{W_N^{-k} - z^{-1} }{1- 2 \cdot \cos \left( \frac{2\pi}{N} \cdot k \right) \cdot z^{-1} + z^{-2}}. $$
The transfer function (12) corresponds to the second order IIR-filter structure of which is shown in Figure 2 where $\alpha = 2 \cdot \cos \left( 2\pi \cdot \frac{k}{N}\right)$.
Figure 2. Structure of the IIR-Filter Realizing
Calculation of the Spectral Sample with the Number $k$
Pay attention that $\alpha = 2 \cdot \cos \left( 2\pi \cdot \frac{k}{N}\right)$ is a real coefficient therefore multiplication in the recursive filter branch has become real, and at that, multiplication by $-1$ cannot be considered. We have one complex coefficient $W_N^{-k}$ in the transfer function numerator. But we find that intermediate values $y_r(k)$ are not important to us, consequently to multiply by $W_N^{-k}$ is possible only per $N-1$ iteration when $S(k)$ is directly calculated.
Thus, using the second order filter allows to calculate $S(k)$ for the fixed $k$ with using $N$ multipliers, not complex, but real. Also one complex multiplication by $W_N^{-k}$ 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.
% N-point DFT
N = 100;
% Number of the spectral sample (bean)
k = 32;

%% Pickup signal and calculation of its DFT    
s = sin(2*pi*(0:N-1)*k/N+pi/6);
X = fft(s);
S_dft = X(k+1); % k+1 because Matlab and GNU Octave index arrays beginning with 1 

%%  Goertzel Algorithm by means of DSPLIB.ORG. Equation (12) and Figure 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 = 2*cos(2*pi*k/N);

% IIR-filter output by means of DSPLIB.ORG 
b_dsplib = [w, -1, 0];
a_dsplib = [1, -alpha, 1];
X = filter(b_dsplib, a_dsplib, s);
S_dsplib = X(end);

%% Goertzel Algorithm by Oppenheim  Alan V [2] pages 633-634
% Figure 9.1, equation (9.9)
%              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 = 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('Pickup DFT ...........................%9.4f %+9.4f i\n', real(S_dft), imag(S_dft));
fprintf('Algorithm dspib.org......................%9.4f %+9.4f i\n', real(S_dsplib), imag(S_dsplib));
fprintf('By Оппенгейма.....................%9.4f %+9.4f i\n', real(S_opp), imag(S_opp));
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 $S(k)$ is equal to:
$$ S(k) = y_{N-1}(k) = W_N^{-k} \cdot v(N-1) - v(N -2), $$
where $v(r)$ are intermediate values which are calculated iteratively:
$$ v(r) = s(r) + 2 \cdot \cos \left(2\pi\cdot\frac{k}{N} \right) \cdot v(r-1) - v(r-2). $$
Example of Using Goertzel Algorithm
Consider an example of using Goertzel algorithm. Let $N=8$, the input signal $s(n)$ is shown in Figure 3.
Figure 3. Input Signal
Calculate the spectral sample $S(1)$ with the number $k=1$ by means of Goertzel algorithm.
Calculate coefficients $\alpha$ and $W_N^{-k}$ beforehand:
$$ \alpha = 2 \cdot \cos \left(2\pi\cdot\frac{1}{8} \right) = 1.4142; $$ $$ W_N^{-k} = W_8^{-1} = \exp\left( j \cdot \frac{2\pi}{8} \right) = 0.7071 + j \cdot 0.7071. $$
Set $v(-2) = v(-1) = 0$ as initial calculation specifications.
Calculate $v(r)$, $r = 0 \ldots N-1$ according to (14) iteratively:
\begin{align} v(0) =& s(0) + \alpha \cdot v(-1) -v(-2) = 3; \\ v(1) =& s(1) + \alpha \cdot v(0) - v(-1)= 6.2426;\\ v(2) =& s(2) + \alpha \cdot v(1) -v(0) = 6.8284; \\ v(3) =& s(3) + \alpha \cdot v(2) - v(1)= 2.4142;\\ v(4) =& s(4) + \alpha \cdot v(3) -v(2) = -2.4142; \\ v(5) =& s(5) + \alpha \cdot v(4) - v(3)= -7.8284;\\ v(6) =& s(6) + \alpha \cdot v(5) -v(4) = -11.6569; \\ v(7) =& s(7) + \alpha \cdot v(6) - v(5)= -10.6569. \end{align}
Then, according to (13) it is possible to get values for the spectral sample $S(1)$:
$$ S(1) = W_8^{-1} \cdot v(7) - v(6) = 4.1213 - j \cdot 7.5355. $$
Using Goertzel Algorithm to Decode the DTMF signal
Figure 4. Coding DTMF Symbols
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.
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 $F_{\textrm{s}} = 8$ kHz, for $205-$point DPF ($N = 205$), numbers of spectral samples $k$ are given in Table 1.

$f(k)$, (Hz)6977708529411209133614771633
Table 1. Indexes of Spectral Samples $k$ Corresponding to DTMF Tones $F_{\textrm{s}} = 8$ kHz and $N = 205$
Thus, for DTMF decoding we need to calculate only 8 spectral samples instead of calculating $205-$point DPF.
Simulating DTMF signals and calculating spectral samples were made with using DSPL library. The source program code is given in the DSPL library examples section.
The time domain waveform got as a result of simulating program operation is shown in Figure 5.
Figure 5. Domain Waveform
Spectral samples magnitudes calculated by means of dspl_goertzel 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
Users of Matlab and GNU Octave can implement goertzel_dtmf.m script to plot Figures 5 and 6.
clear all;
close all;

% Sampling frequency
Fs = 8000;

% Pickup signal vector size
N  = 205;

% DTMF frequency vector
frq = [697, 770, 852, 941, 1209, 1336, 1477, 1633];

% Spectral sample DTMF index vector (indexation begins with 0)
ind = [18,  20,  22,  24,  31,   34,   38,   42];

% Sampling moments of the pickup DTMF signal
t = (0:N-1) / Fs;

for(k = 1:4)  
  for(m = 1:4)
    % DTMF signal
    s = sin(2*pi*frq(k)*t) + sin(2*pi*frq(m+4)*t);
    % Spectral sample vector
    S = goertzel(s, ind);    
    % Domain waveform
    figure(1); subplot(4,4,(k-1)*4+m); 
    plot(t,s); axis([0, 0.025, -2, 2]);
    % Spectral samples magnitudes
    figure(2); subplot(4,4,(k-1)*4+m); 
    stem(frq,abs(S), '.'); axis([500, 1800, 0, 130]);

function S = goertzel(s, ind)
% S = goertzel(s, ind)
% Function returns DFT samples indexes ind, calculated by Goertzel algorithm
% for input signal s.
% Input parameters
%  s   - input signal vector [N x 1];
%  ind - DFT samples indexes vector (indexation starts from 0)
% Ouptut parameters
%  S   - DFT samples corresponds to indexes ind
% Author: Sergey Bakhurin (dsplib.org)

  N = length(s);
  for k = 1:length(ind);
    w = exp(2i*pi*(ind(k)-1)/N);
    alpha = 2*cos(2*pi*(ind(k)-1)/N);
    % second-order IIR 
    b = [w, -1, 0];
    a = [1, -alpha, 1];
    X = filter(b, a, s);
    S(k) = X(end);

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. The result of simulating DTMF signal with using dspl_goertzel function which is a part of DSPL library, and also m-scripts of Matlab and GNU Octave is given in this section.

You can provide your any feedbacks, questions and wishes in the message board or guestbook.
[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.

Select spelling error with your mouse and press Система Orphus