Algoritmo di Goertzel

L'algoritmo di Goertzel è una tecnica di digital signal processing (DSP) utilizzata per identificare le diverse componenti in frequenza di un segnale. La più generale Fast Fourier transform (FFT) considera la completa banda passante del segnale; l'algoritmo di Goertzel si riferisce invece ad alcuni punti specifici predeterminati.

Premessa modifica

Dato il segnale tempo-discreto x[n] causale e di lunghezza N, dal quale si desidera estrarre una o più componenti in frequenza, si consideri la sequenza complessa:

y[n]=x[n]+y[n-1]e^{j\omega }

anch'essa causale di lunghezza N, e dipendente dal parametro $\omega$ . L'algoritmo di Goertzel sfrutta il fatto che l'ultimo campione (ossia l'(N-1)-esimo) di questa sequenza coincide, a meno di un riscalaggio, con la trasformata di Fourier tempo-discreta (DTFT) di x[n], calcolata nella pulsazione angolare $\omega$ (rad). Per poter verificare ciò, è sufficiente passare alla funzione di trasferimento del sistema IIR che calcola y[n]:

H(z)={\frac {Y(z)}{X(z)}}={\frac {1}{1-e^{j\omega }z^{-1}}}

Svolgendo il rapporto tra polinomi, si ottiene:

H(z)=1+e^{j\omega }z^{-1}+e^{j\omega 2}z^{-2}+...=\sum _{k=0}^{+\infty }e^{j\omega k}z^{-k}

per cui:

y[n]=\sum _{k=0}^{+\infty }e^{j\omega k}x[n-k]=\sum _{m=n}^{-\infty }e^{j\omega (n-m)}x[m]=e^{j\omega n}\sum _{m=-\infty }^{n}e^{-j\omega m}x[m]=e^{j\omega n}\sum _{m=0}^{n}x[m]e^{-j\omega m}

ove nel terzultimo passaggio è stato operato il cambio di variabile m = n-k. La sommatoria è inoltre stata troncata sull'intervallo (0;n) essendo x[n] una sequenza nulla per n < 0. L'ultimo campione risulta quindi:

y[N-1]=G(\omega )=e^{j\omega (N-1)}\sum _{m=0}^{N-1}x[m]e^{-j\omega m}

cioè la trasformata di Fourier tempo-discreta $X(\omega )$ del segnale x[n] (calcolata nella pulsazione $\omega$ ) moltiplicata per la costante complessa $e^{j\omega (N-1)}$ .

Ponendo $X(\omega )=e^{-j\omega (N-1)}G(\omega )$ , il problema che l'algoritmo di Goertzel affronta è il calcolo di $|X(\omega )|^{2}$ passando attraverso un certo numero di operazioni puramente reali. Notare che, per semplice definizione di trasformata, non è possibile ottenere con operazioni reali la quantità complessa $X(\omega )$ ; pur tuttavia, l'algoritmo di Goertzel mira all'ottenimento del modulo, che spesso riveste maggior importanza pratica. Ad esempio, la rilevazione sul segnale x[n] della presenza di un tono DTMF può essere compiuta andando ad esaminare l'ampiezza delle componenti spettrali alle due frequenze di tono, non considerando minimamente la loro fase. Gli algoritmi di calcolo della DFT e della DTFT richiedono invece il passaggio attraverso calcoli non reali anche nel caso sia richiesto il solo modulo delle componenti, pertanto da questo punto di vista sono più svantaggiosi rispetto all'algoritmo di Goertzel.

Spiegazione dell'algoritmo modifica

Per poter calcolare $X(\omega )=e^{-j\omega (N-1)}G(\omega )$ , occorre prima di tutto ottenere la sequenza y[n] partendo dall'ingresso x[n]. L'algoritmo di Goertzel adempie a ciò interpretando la trasformata Z riportata in precedenza come la cascata tra le funzioni di trasferimento di due sistemi discreti: un filtro IIR del secondo ordine (a coefficienti puramente reali), ed un filtro FIR del primo ordine (a coefficienti complessi). Infatti si ha:

H(z)={\frac {1}{1-e^{j\omega }z^{-1}}}={\frac {1-e^{-j\omega }z^{-1}}{(1-e^{j\omega }z^{-1})(1-e^{-j\omega }z^{-1})}}={\frac {1}{1-2\cos(\omega )z^{-1}+z^{-2}}}(1-e^{-j\omega }z^{-1})={\frac {S(z)}{X(z)}}{\frac {Y(z)}{S(z)}}

L'uscita s[n] del filtro IIR è facilmente calcolabile per via ricorsiva come:

s[n]=x[n]+2\cos(\omega )s[n-1]-s[n-2]

La sequenza y[n] espressa in funzione di s[n] risulta invece:

y[n]=s[n]-e^{-j\omega }s[n-1]

A questo punto è sufficiente calcolare:

X(\omega )=e^{-j\omega (N-1)}y[N-1]=e^{-j\omega (N-1)}(s[N-1]-e^{-j\omega }s[N-2])

Si noti come in realtà non occorre affatto tutta la sequenza y[n], ma solo l'ultimo suo campione. Per ottenerlo, sono a loro volta necessari gli ultimi due samples di s[n], sequenza che comunque dovrà essere calcolata per intero (dipendendo ciascun campione dai due precedenti). Si adottano come condizioni iniziali s[-1] = s[-2] = 0.

Nel caso comune, piuttosto che la DTFT si preferisce lavorare con la trasformata discreta di Fourier (essendo questa per l'appunto una sequenza); la pulsazione $\omega$ d'interesse è ottenuta in questo caso discretizzando l'asse delle frequenze sull'intervallo $(0;2\pi )$ :

\omega \to \omega _{k}={\frac {2\pi }{N}}k,k=0..N-1

Adottando questa sostituzione, si ottiene:

N={\frac {2\pi k}{\omega _{k}}}

e^{-j\omega (N-1)}\to e^{-j\omega _{k}(N-1)}=e^{-j\omega _{k}({\frac {2\pi k}{\omega _{k}}}-1)}=e^{-j2\pi k}e^{j\omega _{k}}=e^{j\omega _{k}}

X(\omega _{k})=e^{j\omega _{k}}(s[N-1]-e^{-j\omega _{k}}s[N-2])=e^{j\omega _{k}}s[N-1]-s[N-2]

Essendo infine s[n] una sequenza reale (se x[n] è reale), il complesso coniugato di $X(\omega _{k})$ risulta banalmente:

X^{*}(\omega _{k})=e^{-j\omega _{k}}s[N-1]-s[N-2]

Pertanto, è possibile ottenere l'intensità della componente di pulsazione $\omega _{k}$ come:

|X(\omega _{k})|^{2}=X(\omega _{k})X^{*}(\omega _{k})=s^{2}[N-1]+s^{2}[N-2]-2\cos(\omega _{k})s[N-1]s[N-2]

cioè con una sequenza di operazioni reali. Ripetendo il calcolo di s[n] e |X|² per ulteriori altre pulsazioni, è possibile ottenere lo spettro di ampiezza del segnale x[n], sebbene l'algoritmo di Goertzel sia maggiormente impiegato per l'estrazione di singole componenti a frequenze predeterminate.

Segue un esempio di programma che riconosce la presenza di toni DTMF in un segnale digitale (ottenuto ad esempio per campionamento e quantizzazione), sfruttando l'algoritmo di Goertzel per l'individuazione delle componenti spettrali (DFT bins).

Codice d'esempio per un decodificatore Goertzel modifica

#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
 *----------------------------------------------------------------------------
 * Qui è dove guardiamo nei bin e decidiamo se abbiamo un segnale valido.
 */
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"}};

    /* Trova il maggiore nel gruppo delle righe. */
    row = 0;
    maxval = 0.0;
    for ( i=0; i<4; i++ )
    {
        if ( r[i] > maxval )
        {
            maxval = r[i];
            row = i;
        }
    }

    /* Trova il maggiore nel gruppo delle colonne. */
    col = 4;
    maxval = 0.0;
    for ( i=4; i<8; i++ )
    {
        if ( r[i] > maxval )
        {
            maxval = r[i];
            col = i;
        }
    }

    /* Controlla per energia minima */

    if ( r[row] < 4.0e5 )   /* 2.0e5 ... 1.0e8 non è significativo */
    {
        /* energia non abbastanza elevata */
    }
    else if ( r[col] < 4.0e5 )
    {
        /* energia non abbastanza elevata */
    }
    else
    {
        see_digit = TRUE;

        /* Controllo di inversione
         * CEPT => inversione < 6dB
         * AT&T => inversione diretta < 4dB e inversione inversa < 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] )
        {
            /* Inversione diretta */
            max_index = col;
            if ( r[row] < (r[col] * 0.398) )    /* inversione > 4dB, errore */
                see_digit = FALSE;
        }
        else /* if ( r[row] > r[col] ) */
        {
            /* Inversione inversa */
            max_index = row;
            if ( r[col] < (r[row] * 0.158) )    /* inversione > 8db, errore */
                see_digit = FALSE;
        }

        /* Controllo del rapporto segnale/rumore
         * AT&T afferma che il rumore deve essere inferiore di 16dB rispetto al segnale.
         * Qui contiamo il numero di segnali sopra la soglia e
         * dovrebbero essercene solo due.
         */
        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;

    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;
    }
}

Collegamenti esterni modifica

Portale Informatica: accedi alle voci di Wikipedia che trattano di informatica