Documentation

goertzel

Discrete Fourier transform with second-order Goertzel algorithm

Syntax

dft_data = goertzel(data)
dft_data = goertzel(data,freq_indices)
dft_data = goertzel(data,freq_indices,dim)

Description

dft_data = goertzel(data) returns the discrete Fourier transform (DFT) of the input data, data, using a second-order Goertzel algorithm. If data is a matrix, goertzel computes the DFT of each column separately. You can specify data types as double or single precision.

dft_data = goertzel(data,freq_indices) returns the DFT for the frequency indices freq_indices. The values of freq_indices can be any whole number or fraction.

dft_data = goertzel(data,freq_indices,dim) computes the DFT of the matrix data along the dimension dim.

Examples

collapse all

Estimate the frequencies of the tone generated by pressing the 1 button on a telephone keypad.

The number 1 produces a tone with frequencies 697 and 1209 Hz. Generate 205 samples of the tone with a sample rate of 8 kHz.

Fs = 8000;
N = 205;
lo = sin(2*pi*697*(0:N-1)/Fs);
hi = sin(2*pi*1209*(0:N-1)/Fs);
data = lo + hi;

Use the Goertzel algorithm to compute the DFT of the tone. Choose the indices corresponding to the frequencies used to generate the numbers 0 through 9.

f = [697 770 852 941 1209 1336 1477];
freq_indices = round(f/Fs*N) + 1;
dft_data = goertzel(data,freq_indices);

Plot the DFT magnitude.

stem(f,abs(dft_data))

ax = gca;
ax.XTick = f;
xlabel('Frequency (Hz)')
title('DFT Magnitude') Generate a noisy cosine with frequency components at 1.24 kHz, 1.26 kHz, and 10 kHz. Specify a sample rate of 32 kHz. Reset the random number generator for reproducible results.

rng default

Fs = 32e3;
t = 0:1/Fs:2.96;
x = cos(2*pi*t*10e3) + cos(2*pi*t*1.24e3) + cos(2*pi*t*1.26e3) ...
+ randn(size(t));

Generate the frequency vector. Use the Goertzel algorithm to compute the DFT. Restrict the range of frequencies to between 1.2 and 1.3 kHz.

N = (length(x)+1)/2;
f = (Fs/2)/N*(0:N-1);
indxs = find(f>1.2e3 & f<1.3e3);
X = goertzel(x,indxs);

Plot the mean squared spectrum expressed in decibels.

plot(f(indxs)/1e3,mag2db(abs(X)/length(X)))

title('Mean Squared Spectrum')
xlabel('Frequency (kHz)')
ylabel('Power (dB)')
grid Algorithms

The Goertzel algorithm implements the DFT as a recursive difference equation. To establish this difference equation, express the DFT as the convolution of an N-point input, x(n), with the impulse response$h\left(n\right)={W}_{N}^{-kn}u\left(n\right)$, where ${W}_{N}^{-kn}={e}^{-j2\pi k/N}$ and u(n) is the unit step sequence.

The Z-transform of the impulse response is

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

The direct form II implementation is: Alternatives

You can also compute the DFT with:

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

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

References

 Proakis, John G., and Dimitris G. Manolakis. Digital Signal Processing: Principles, Algorithms, and Applications. 3rd Edition. Upper Saddle River, NJ: Prentice Hall, 1996, pp. 480–481.

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