Discrete Fourier transform with second-order Goertzel algorithm

The Goertzel algorithm implements the discrete Fourier transform *X*(*k*) as the convolution of an *N*-point input *x*(*n*), *n* = 0, 1, …, *N* – 1, with the impulse response

$${h}_{k}(n)={e}^{-j2\pi k}\text{\hspace{0.17em}}{e}^{j2\pi kn/N}\text{\hspace{0.17em}}u(n)\equiv {e}^{-j2\pi k}\text{\hspace{0.17em}}{W}_{N}^{-kn}\text{\hspace{0.17em}}u(n),$$

where *u*(*n*), the unit step sequence, is 1 for *n* ≥ 0 and 0 otherwise. *k* does not have to be an integer. At a
frequency *f* =
*k**f*_{s}/*N*, where *f*_{s} is the sample rate, the transform has a value

$$X(k)={{y}_{k}(n)|}_{n=N},$$

where

$${y}_{k}(n)={\displaystyle \sum _{m=0}^{N}x(m)\text{\hspace{0.17em}}{h}_{k}(n-m)}$$

and *x*(*N*) = 0. The Z-transform of the impulse response is

$${H}_{k}(z)=\frac{\left(1-{W}_{N}^{k}{z}^{-1}\right){e}^{-j2\pi k}}{1-2\mathrm{cos}\left(\frac{2\pi k}{N}\right)\text{\hspace{0.17em}}{z}^{-1}+{z}^{-2}},$$

with this direct form II implementation:

Compare the output of `goertzel`

to the result of a direct
implementation of the Goertzel algorithm. For the input signal, use a chirp sampled at 50 Hz
for 10 seconds and embedded in white Gaussian noise. The chirp's frequency increases linearly
from 15 Hz to 20 Hz during the measurement. Compute the discrete Fourier transform at a
frequency that is not an integer multiple of *f*_{s}/*N*. When calling `goertzel`

, keep in mind that MATLAB^{®} vectors run from 1 to *N* instead of from 0 to *N* – 1. The results agree to high
precision.

fs = 50; t = 0:1/fs:10-1/fs; N = length(t); xn = chirp(t,15,t(end),20)+randn(1,N)/100; f0 = 17.36; k = N*f0/fs; ykn = filter([1 -exp(-2j*pi*k/N)],[1 -2*cos(2*pi*k/N) 1],[xn 0]); Xk = exp(-2j*pi*k)*ykn(end) dft = goertzel(xn,k+1) df = abs(Xk-dft)

df = 4.3634e-12

You can also compute the DFT with:

`fft`

: less efficient than the Goertzel algorithm when you only need the DFT at a few frequencies.`fft`

is more efficient than`goertzel`

when you need to evaluate the transform at more than log_{2}*N*frequencies, where*N*is the length of the input signal.`czt`

:`czt`

calculates the chirp Z-transform of an input signal on a circular or spiral contour and includes the DFT as a special case.

[1] Burrus, C. Sidney, and Thomas W.
Parks. *DFT/FFT and Convolution Algorithms: Theory and Implementation*.
New York: John Wiley & Sons, 1985.

[2] Proakis, John G., and Dimitris G.
Manolakis. *Digital Signal Processing: Principles, Algorithms, and
Applications*. 3rd Edition. Upper Saddle River, NJ: Prentice Hall,
1996.

[3] Sysel, Petr, and Pavel Rajmic.
“Goertzel Algorithm Generalized to Non-Integer Multiples of Fundamental Frequency.”
*EURASIP Journal on Advances in Signal Processing*. Vol. 2012, Number
1, December 2012, pp. 56-1–56-8. https://doi.org/10.1186/1687-6180-2012-56.