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.

dft_data = goertzel(data,freq_indices) returns the DFT for the frequency indices freq_indices.

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

Examples

expand all

Estimate Telephone Keypad Frequencies

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')

Resolve Frequency Components of a Noisy Tone

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

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.

More About

expand all

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 responseh(n)=WNknu(n), where WNkn=ei2πk/N and u(n) is the unit step sequence.

The Z-transform of the impulse response is

H(z)=1WNkz112cos(2πk/N)z1+z2.

The direct form II implementation is:

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.

Was this topic helpful?