Goertzel algorithm好文收集

最新推荐文章于 2021-03-29 22:37:53 发布

xiaomantoubaba

最新推荐文章于 2021-03-29 22:37:53 发布

阅读量3k

点赞数

文章标签： algorithm processing fft signal components constants

本文链接：https://blog.csdn.net/sshcx/article/details/2757655

版权

Goertzel algorithm

From Wikipedia, the free encyclopedia

Jump to: navigation, search

The Goertzel algorithm is a digital signal processing (DSP) technique for identifying frequency components of a signal, published by Dr. Gerald Goertzel in 1958. While the general Fast Fourier transform (FFT) algorithm computes evenly across the bandwidth of the incoming signal, the Goertzel algorithm looks at specific, predetermined frequency.

A practical application of this algorithm is that of recognizing the tones produced by the buttons pushed on a telephone keypad. This application is illustrated by an implementation of the algorithm in C as shown below to produce a DTMF tone detector.

Contents

[hide]

<script type=text/javascript> // </script> [edit] Explanation of algorithm

The Goertzel algorithm computes a sequence, s(n), given an input sequence, x(n), as

s(n) = x(n) + 2cos(2πω)s(n − 1) − s(n − 2)

where s( − 2) = s( − 1) = 0 and ω is some frequency of interest, in cycles per sample, which should be less than 1/2. This effectively implements a second-order IIR filter with poles at e ^{+ 2πiω} and e ^{− 2πiω}, and requires only one multiplication (assuming 2cos(2πω) is pre-computed), one addition and one subtraction per input sample. For real inputs, these operations are real.

The Z transform of this process is

Applying an additional, FIR, transform of the form

will give an overall transform of

The time-domain equivalent of this overall transform is

which, when x(n) = 0 for all n < 0, becomes

or, the equation for the (n + 1)-sample DFT of x, evaluated for ω and multiplied by the scale factor e ^{+ 2πiωn}.

Notice that applying the additional transform Y(z)/S(z) only requires the last two samples of the s sequence. Consequently, upon processing N samples x(0)...x(N − 1), the last two samples from the s sequence can be used to compute the value of a DFT bin which corresponds to the chosen frequency ω as

X(ω) = y(N − 1)e ^{− 2πiω(N − 1)} = (s(N − 1) − e ^{− 2πiω}s(N − 2))e ^{− 2πiω(N − 1)}

For the special case often found when computing DFT bins, where ωN = k for some integer, k, this simplifies to

X(ω) = (s(N − 1) − e ^{− 2πiω}s(N − 2))e ^{+ 2πiω} = e ^{+ 2πiω}s(N − 1) − s(N − 2)

In either case, the corresponding power can be computed using the same cosine term required to compute s as

X(ω)X'(ω) = s(N − 2)² + s(N − 1)² − 2cos(2πω)s(N − 2)s(N − 1)

When implemented in a general-purpose processor, values for s(n − 1) and s(n − 2) can be retained in variables and new values of s can be shifted through as they are computed, assuming that only the final two values of the s sequence are required. The code may then be as follows:

s_prev = 0

s_prev2 = 0

coeff = 2*cos(2*PI*normalized_frequency);

for each sample, x[n],

s = x[n] + coeff*s_prev - s_prev2;

s_prev2 = s_prev;

s_prev = s;

end

power = s_prev2*s_prev2 + s_prev*s_prev - coeff*s_prev2*s_prev;

[edit] Computational complexity

In order to compute a single DFT bin for a complex sequence of length N, this algorithm requires 2N multiplies and 4N add/subtract operations within the loop, as well as 4 multiplies and 4 add/subtract operations to compute X(ω), for a total of 2N+4 multiplies and 4N+4 add/subtract operations (for real sequences, the required operations are half that amount). In contrast, the Fast Fourier transform (FFT) requires 2log₂N multiplies and 3log₂N add/subtract operations per DFT bin, but must compute all N bins simultaneously (similar optimizations are available to halve the number of operations in an FFT when the input sequence is real).

When the number of desired DFT bins, M, is small (e.g., when detecting DTMF tones), it is computationally advantageous to implement the Goertzel algorithm, rather than the FFT. Approximately, this occurs when

or if, for some reason, N is not an integral power of 2 while you stick to a simple FFT algorithm like the 2-radix Cooley-Tukey FFT algorithm, and zero-padding the samples out to an integral power of 2 would violate

Moreover, the Goertzel algorithm can be computed as samples come in, and the FFT algorithm may require a large table of N pre-computed sines and cosines in order to be efficient.

If multiplications are not considered as difficult as additions, or vice versa, the 5/6 ratio can shift between anything from 3/4 (additions dominate) to 1/1 (multiplications dominate).

[edit] Practical considerations

The term DTMF or Dual-Tone Multi Frequency is the official name of the tones generated from a telephone keypad. (AT&T used the trademark "Touch-Tone" for its DTMF dialing service.^[1]) The original keypads were mechanical switches triggering RC controlled oscillators.^{[citation needed]} The digit detectors were also tuned circuits. The interest in decoding DTMF is high because of the large numbers of phones generating these types of tones.

At present, DTMF detectors are most often implemented as numerical algorithms on either general purpose computers or on fast digital signal processors. The algorithm shown below is an example of such a detector.

However, this algorithm needs an additional post-processing step to completely implement a functional DTMF tone detector. DTMF tone bursts can be as short as 50 milli-seconds or as long as several seconds. The tone burst can have noise or dropouts within it which must be ignored. The Goertzel algorithm produces multiple outputs; a post-processing step needs to smooth these outputs into one output per tone burst.

One additional problem is that the algorithm will sometimes produce spurious outputs because of a window period that is not completely filled with samples. Imagine a DTMF tone burst and then imagine the window superimposed over this tone burst. Obviously, the detector is running at a fixed rate and the tone burst is not guaranteed to arrive aligned with the timing of the detector. So some window intervals on the leading and trailing edges of the tone burst will not be entirely filled with valid tone samples. Worse, RC-based tone generators will often have voltage sag/surge related anomalies at the leading and trailing edges of the tone burst. These also can contribute to spurious outputs.

It is highly likely that this detector will report false or incorrect results at the leading and trailing edges of the tone burst due to a lack of sufficient valid samples within the window. In addition, the tone detector must be able to tolerate tone dropouts within the tone burst and these can produce additional false reports due to the same windowing effects.

The post-processing system can be implemented as a statistical aggregator which will examine a series of outputs of the algorithm below. There should be a counter for each possible output. These all start out at zero. The detector starts producing outputs and depending on the output, the appropriate counter is incremented. Finally, the detector stops generating outputs for long enough that the tone burst can be considered to be over. The counter with the highest value wins and should be considered to be the DTMF digit signaled by the tone burst.

While it is true that there are eight possible frequencies in a DTMF tone, the algorithm as originally entered on this page was computing a few more frequencies so as to help reject false tones (talkoff). Notice the peak tone counter loop. This checks to see that only two tones are active. If more than this are found then the tone is rejected.

[edit] Sample code for a DTMF detector

#define SAMPLING_RATE 8000

#define MAX_BINS 8

#define GOERTZEL_N 92

int sample_count;

double q1[ MAX_BINS ];

double q2[ MAX_BINS ];

double r[ MAX_BINS ];

* coef = 2.0 * cos( (2.0 * PI * k) / (float)GOERTZEL_N)) ;

* Where k = (int) (0.5 + ((float)GOERTZEL_N * target_freq) / SAMPLING_RATE));

* More simply: coef = 2.0 * cos( (2.0 * PI * target_freq) / SAMPLING_RATE );

double freqs[ MAX_BINS] =

{

697,

770,

852,

941,

1209,

1336,

1477,

1633

};

double coefs[ MAX_BINS ] ;

/*----------------------------------------------------------------------------

* calc_coeffs

*----------------------------------------------------------------------------

* This is where we calculate the correct co-efficients.

void calc_coeffs()

{

int n;

for(n = 0; n < MAX_BINS; n++)

{

coefs[n] = 2.0 * cos(2.0 * 3.141592654 * freqs[n] / SAMPLING_RATE);

}

/*----------------------------------------------------------------------------

* post_testing

*----------------------------------------------------------------------------

* This is where we look at the bins and decide if we have a valid signal.

void post_testing()

{

int row, col, see_digit;

int peak_count, max_index;

double maxval, t;

int i;

char * row_col_ascii_codes[4][4] = {

{"1", "2", "3", "A"},

{"4", "5", "6", "B"},

{"7", "8", "9", "C"},

{"*", "0", "#", "D"}};

/* Find the largest in the row group. */

row = 0;

maxval = 0.0;

for ( i=0; i<4; i++ )

{

if ( r[i] > maxval )

{

maxval = r[i];

row = i;

}

/* Find the largest in the column group. */

col = 4;

maxval = 0.0;

for ( i=4; i<8; i++ )

{

if ( r[i] > maxval )

{

maxval = r[i];

col = i;

}

/* Check for minimum energy */

if ( r[row] < 4.0e5 ) /* 2.0e5 ... 1.0e8 no change */

{

/* energy not high enough */

}

else if ( r[col] < 4.0e5 )

{

/* energy not high enough */

}

else

{

see_digit = TRUE;

/* Twist check

* CEPT => twist < 6dB

* AT&T => forward twist < 4dB and reverse twist < 8dB

* -ndB < 10 log10( v1 / v2 ), where v1 < v2

* -4dB < 10 log10( v1 / v2 )

* -0.4 < log10( v1 / v2 )

* 0.398 < v1 / v2

* 0.398 * v2 < v1

if ( r[col] > r[row] )

{

/* Normal twist */

max_index = col;

if ( r[row] < (r[col] * 0.398) ) /* twist > 4dB, error */

see_digit = FALSE;

}

else /* if ( r[row] > r[col] ) */

{

/* Reverse twist */

max_index = row;

if ( r[col] < (r[row] * 0.158) ) /* twist > 8db, error */

see_digit = FALSE;

}

/* Signal to noise test

* AT&T states that the noise must be 16dB down from the signal.

* Here we count the number of signals above the threshold and

* there ought to be only two.

if ( r[max_index] > 1.0e9 )

t = r[max_index] * 0.158;

else

t = r[max_index] * 0.010;

peak_count = 0;

for ( i=0; i<8; i++ )

{

if ( r[i] > t )

peak_count++;

}

if ( peak_count > 2 )

see_digit = FALSE;

if ( see_digit )

{

printf( "%s", row_col_ascii_codes[row][col-4] );

fflush(stdout);

}

/*----------------------------------------------------------------------------

* goertzel

*----------------------------------------------------------------------------

void goertzel( int sample )

{

double q0;

ui32 i;

sample_count++;

/*q1[0] = q2[0] = 0;*/

for ( i=0; i<MAX_BINS; i++ )

{

q0 = coefs[i] * q1[i] - q2[i] + sample;

q2[i] = q1[i];

q1[i] = q0;

}

if (sample_count == GOERTZEL_N)

{

for ( i=0; i<MAX_BINS; i++ )

{

r[i] = (q1[i] * q1[i]) + (q2[i] * q2[i]) - (coefs[i] * q1[i] * q2[i]);

q1[i] = 0.0;

q2[i] = 0.0;

}

post_testing();

sample_count = 0;

}

[edit] References

^ USPTO trademark entry for Touch-Tone, retrieved on March 29, 2007.

[edit] External links

Modified Goertzel Algorithm in DTMF Detection Using the TMS320C80 DSP

The Goertzel Algorithm

The Goertzel algorithm can perform tone detection using much less CPU horsepower than the Fast Fourier Transform, but many engineers have never heard of it. This article attempts to change that.

Most engineers are familiar with the Fast Fourier Transform (FFT) and would have little trouble using a "canned" FFT routine to detect one or more tones in an audio signal. What many don't know, however, is that if you only need to detect a few frequencies, a much faster method is available. It's called the Goertzel algorithm.

Tone detection

Many applications require tone detection, such as:

DTMF (touch tone) decoding
Call progress (dial tone, busy, and so on) decoding
Frequency response measurements (sending a tone while simultaneously reading back the result)-if you do this for a range of frequencies, the resulting frequency response curve can be informative. For example, the frequency response curve of a telephone line tells you if any load coils (inductors) are present on that line.

Although dedicated ICs exist for the applications above, implementing these functions in software costs less. Unfortunately, many embedded systems don't have the horsepower to perform continuous real-time FFTs. That's where the Goertzel algorithm comes in.

In this article, I describe what I call a basic Goertzel and an optimized Goertzel.

The basic Goertzel gives you real and imaginary frequency components as a regular Discrete Fourier Transform (DFT) or FFT would. If you need them, magnitude and phase can then be computed from the real/imaginary pair.

The optimized Goertzel is even faster (and simpler) than the basic Goertzel, but doesn't give you the real and imaginary frequency components. Instead, it gives you the relative magnitude squared. You can take the square root of this result to get the relative magnitude (if needed), but there's no way to obtain the phase.

In this short article, I won't try to explain the theoretical background of the algorithm. I do give some links at the end where you can find more detailed explanations. I can tell you that the algorithm works well, having used it in all of the tone detection applications previously listed (and others).

A basic Goertzel

First a quick overview of the algorithm: some intermediate processing is done in every sample. The actual tone detection occurs every Nth sample. (I'll talk more about N in a minute.)

As with the FFT, you work with blocks of samples. However, that doesn't mean you have to process the data in blocks. The numerical processing is short enough to be done in the very interrupt service routine (ISR) that is gathering the samples (if you're getting an interrupt per sample). Or, if you're getting buffers of samples, you can go ahead and process them a batch at a time.

Before you can do the actual Goertzel, you must do some preliminary calculations:

Decide on the sampling rate.
Choose the block size, N.
Precompute one cosine and one sine term.
Precompute one coefficient.

These can all be precomputed once and then hardcoded in your program, saving RAM and ROM space; or you can compute them on-the-fly.

Sampling rate

Your sampling rate may already be determined by the application. For example, in telecom applications, it's common to use a sampling rate of 8kHz (8,000 samples per second). Alternatively, your analog-to-digital converter (or CODEC) may be running from an external clock or crystal over which you have no control.

If you can choose the sampling rate, the usual Nyquist rules apply: the sampling rate will have to be at least twice your highest frequency of interest. I say "at least" because if you are detecting multiple frequencies, it's possible that an even higher sampling frequency will give better results. What you really want is for every frequency of interest to be an integer factor of the sampling rate.

Block size

Goertzel block size N is like the number of points in an equivalent FFT. It controls the frequency resolution (also called bin width). For example, if your sampling rate is 8kHz and N is 100 samples, then your bin width is 80Hz.

This would steer you towards making N as high as possible, to get the highest frequency resolution. The catch is that the higher N gets, the longer it takes to detect each tone, simply because you have to wait longer for all the samples to come in. For example, at 8kHz sampling, it will take 100ms for 800 samples to be accumulated. If you're trying to detect tones of short duration, you will have to use compatible values of N.

The third factor influencing your choice of N is the relationship between the sampling rate and the target frequencies. Ideally you want the frequencies to be centered in their respective bins. In other words, you want the target frequencies to be integer multiples of sample_rate/N.

The good news is that, unlike the FFT, N doesn't have to be a power of two.

Precomputed constants

Once you've selected your sampling rate and block size, it's a simple five-step process to compute the constants you'll need during processing:

w = (2*π/N)*k
cosine = cos w
sine = sin w
coeff = 2 * cosine

For the per-sample processing you're going to need three variables. Let's call them Q₀, Q₁, and Q₂.

Q₁ is just the value of Q₀ last time. Q₂ is just the value of Q₀ two times ago (or Q₁ last time).

Q₁ and Q₂ must be initialized to zero at the beginning of each block of samples. For every sample, you need to run the following three equations:

Q₀ = coeff * Q₁ - Q₂ + sample
Q₂ = Q₁
Q₁ = Q₀

After running the per-sample equations N times, it's time to see if the tone is present or not.

real = (Q₁ - Q₂ * cosine)
imag = (Q₂ * sine)
magnitude² = real² + imag²

A simple threshold test of the magnitude will tell you if the tone was present or not. Reset Q₂ and Q₁ to zero and start the next block.

An optimized Goertzel

The optimized Goertzel requires less computation than the basic one, at the expense of phase information.

The per-sample processing is the same, but the end of block processing is different. Instead of computing real and imaginary components, and then converting those into relative magnitude squared, you directly compute the following:

magnitude² = Q₁² + Q₂²-Q₁*Q₂*coeff

This is the form of Goertzel I've used most often, and it was the first one I learned about.

Pulling it all together

Listing 1 shows a short demo program I wrote to enable you to test-drive the algorithm. The code was written and tested using the freely available DJGPP C/C++ compiler. You can modify the #defines near the top of the file to try out different values of N, sampling_rate, and target_frequency.

The program does two demonstrations. In the first one, both forms of the Goertzel algorithm are used to compute relative magnitude squared and relative magnitude for three different synthesized signals: one below the target_frequency, one equal to the target_frequency, and one above the target_frequency.

You'll notice that the results are nearly identical, regardless of which form of the Goertzel algorithm is used. You'll also notice that the detector values peak near the target frequency.

In the second demonstration, a simulated frequency sweep is run, and the results of just the basic Goertzel are shown. Again, you'll notice a clear peak in the detector output near the target frequency. Figure 1 shows the output of the code in Listing 1.

Figure 1: Demo output

For SAMPLING_RATE = 8000.000000 N = 205 and FREQUENCY = 941.000000, k = 24 and coeff = 1.482867

For test frequency 691.000000:

real = -360.392059 imag = -45.871609

Relative magnitude squared = 131986.640625

Relative magnitude = 363.299652

Relative magnitude squared = 131986.640625

Relative magnitude = 363.299652

For test frequency 941.000000:

real = -3727.528076 imag = -9286.238281

Relative magnitude squared = 100128688.000000

Relative magnitude = 10006.432617

Relative magnitude squared = 100128680.000000

Relative magnitude = 10006.431641

For test frequency 1191.000000:

real = 424.038116 imag = -346.308716

Relative magnitude squared = 299738.062500

Relative magnitude = 547.483398

Relative magnitude squared = 299738.062500

Relative magnitude = 547.483398

Freq=	641.0	rel mag^2=	146697.87500	rel mag=	383.01160
Freq=	656.0	rel mag^2=	63684.62109	rel mag=	252.35812
Freq=	671.0	rel mag^2=	96753.92188	rel mag=	311.05292
Freq=	686.0	rel mag^2=	166669.90625	rel mag=	408.25226
Freq=	701.0	rel mag^2=	5414.02002	rel mag=	73.58002
Freq=	716.0	rel mag^2=	258318.37500	rel mag=	508.25031
Freq=	731.0	rel mag^2=	178329.68750	rel mag=	422.29099
Freq=	746.0	rel mag^2=	71271.88281	rel mag=	266.96796
Freq=	761.0	rel mag^2=	437814.90625	rel mag=	661.67584
Freq=	776.0	rel mag^2=	81901.81250	rel mag=	286.18494
Freq=	791.0	rel mag^2=	468060.31250	rel mag=	684.14935
Freq=	806.0	rel mag^2=	623345.56250	rel mag=	789.52234
Freq=	821.0	rel mag^2=	18701.58984	rel mag=	136.75375
Freq=	836.0	rel mag^2=	1434181.62500	rel mag=	1197.57324
Freq=	851.0	rel mag^2=	694211.75000	rel mag=	833.19373
Freq=	866.0	rel mag^2=	1120359.50000	rel mag=	1058.47034
Freq=	881.0	rel mag^2=	4626623.00000	rel mag=	2150.95874
Freq=	896.0	rel mag^2=	160420.43750	rel mag=	400.52521
Freq=	911.0	rel mag^2=	19374364.00000	rel mag=	4401.63184
Freq=	926.0	rel mag^2=	81229848.00000	rel mag=	9012.76074
Freq=	941.0	rel mag^2=	100128688.00000	rel mag=	10006.43262
Freq=	956.0	rel mag^2=	43694608.00000	rel mag=	6610.18994
Freq=	971.0	rel mag^2=	1793435.75000	rel mag=	1339.19226
Freq=	986.0	rel mag^2=	3519388.50000	rel mag=	1876.00330
Freq=	1001.0	rel mag^2=	3318844.50000	rel mag=	1821.76965
Freq=	1016.0	rel mag^2=	27707.98828	rel mag=	166.45717
Freq=	1031.0	rel mag^2=	1749922.62500	rel mag=	1322.84644
Freq=	1046.0	rel mag^2=	478859.28125	rel mag=	691.99658
Freq=	1061.0	rel mag^2=	284255.81250	rel mag=	533.15643
Freq=	1076.0	rel mag^2=	898392.93750	rel mag=	947.83594
Freq=	1091.0	rel mag^2=	11303.36035	rel mag=	106.31726
Freq=	1106.0	rel mag^2=	420975.65625	rel mag=	648.82635
Freq=	1121.0	rel mag^2=	325753.78125	rel mag=	570.74841
Freq=	1136.0	rel mag^2=	36595.78906	rel mag=	191.30026
Freq=	1151.0	rel mag^2=	410926.06250	rel mag=	641.03516
Freq=	1166.0	rel mag^2=	45246.58594	rel mag=	212.71245
Freq=	1181.0	rel mag^2=	119967.59375	rel mag=	346.36337
Freq=	1196.0	rel mag^2=	250361.39062	rel mag=	500.36127
Freq=	1211.0	rel mag^2=	1758.44263	rel mag=	41.93379
Freq=	1226.0	rel mag^2=	190195.57812	rel mag=	436.11417
Freq=	1241.0	rel mag^2=	74192.23438	rel mag=	272.38251

To avoid false detections in production code, you will probably want to qualify the raw detector readings by using a debouncing technique, such as requiring multiple detections in a row before reporting a tone's presence to the user.

As you can see, the Goertzel algorithm deserves to be added to your signal processing toolbox.

Kevin Banks has been developing embedded systems for 19 years, as a consultant and as an employee at companies including SCI, TxPORT, DiscoveryCom, and Nokia. Currently he is back in consulting mode. He can be reached at kbanks@hiwaay.net.

References

Here are some links to other resources you may find useful:

www.numerix-dsp.com/goertzel.html

ptolemy.eecs.berkeley.edu/papers/96/dtmf_ict/www/node3.html

www.analogdevices.com/library/dspManuals/Using_ADSP-2100_Vol1_books.html

www.analogdevices.com/library/dspManuals/pdf/Volume1/Chapter_14.pdf

www.analogdevices.com/library/dspManuals/Using_ADSP-2100_Vol2_books.html

www.analogdevices.com/library/dspManuals/pdf/2100Chapter_8.pdf

www.delorie.com/djgpp/

xiaomantoubaba

关注

0
点赞
踩
5

收藏

觉得还不错? 一键收藏
0
评论
Goertzel algorithm好文收集

Goertzel algorithmFrom Wikipedia, the free encyclopediaJump to: navigation, searchThe Goertzel algorithm is a digital signal processing (DSP) technique for identifying frequency components of a
复制链接

扫一扫