Jump to content

Goertzel algorithm

From Wikipedia, the free encyclopedia
This is an old revision of this page, as edited by 209.210.211.18 (talk) at 01:25, 29 June 2005. The present address (URL) is a permanent link to this revision, which may differ significantly from the current revision.
(diff) ← Previous revision | Latest revision (diff) | Newer revision → (diff)

The Goertzel algorithm is DSP technique for identifying frequency components of a signal. The general FFT algorithm computes evenly across the bandwidth of the incoming signal. The Goertzel algorithm looks at specific, predetermined points.

Sample Code for a DTMF Detector


/* Stuff for goertzel algorithm */
#define CHANNEL             0
#define MAX_BINS            12
#define GOERTZEL_N          92

int         sample_count;
double      q1[ MAX_BINS ];
double      q2[ MAX_BINS ];
double      r[ MAX_BINS ];
double      coefs[ MAX_BINS*2 ] = {
1.7088388,
1.6339398,
1.5514226,
1.4616719,

1.1533606,
1.0391679,
0.7968022,
0.5395935,

-0.6697592,
-1.0391679,
-1.3651063,
-1.7088388};


/*----------------------------------------------------------------------------
 *  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( "%c", row_col_ascii_codes[row][col-4] );
            fflush(stdout);
        }
    }
}



/*----------------------------------------------------------------------------
 *  goertzel
 *----------------------------------------------------------------------------
 */
void
goertzel( int sample )
{
double      q0;
ui32        i;

    if ( sample_count < GOERTZEL_N )
    {
        sample_count++;
        for ( i=0; i<MAX_BINS; i++ )
        {
            q0 = coefs[i] * q1[i] - q2[i] + sample;
            q2[i] = q1[i];
            q1[i] = q0;
        }
    }
    else
    {
        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;
    }
}