Documentation

This is machine translation

Translated by Microsoft
Mouse over text to see original. Click the button below to return to the English verison of the page.

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

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

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

collapse 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=ej2π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.

Introduced before R2006a

Was this topic helpful?