Algorithme de Goertzel

L' algorithme de Goertzel est une technique de traitement numérique du signal (DSP) permettant d'évaluer efficacement les termes individuels de la transformée de Fourier discrète (DFT). Il est utile dans certaines applications pratiques, telles que la reconnaissance des tonalités de signalisation multifréquence à double tonalité (DTMF) produites par les boutons-poussoirs du clavier d'un téléphone analogique traditionnel . L'algorithme a été décrit pour la première fois par Gerald Goertzel en 1958.

Comme la DFT, l'algorithme de Goertzel analyse une composante de fréquence sélectionnable à partir d'un signal discret . Contrairement aux calculs DFT directs, l'algorithme de Goertzel applique un seul coefficient à valeur réelle à chaque itération, en utilisant l'arithmétique à valeur réelle pour les séquences d'entrée à valeur réelle. Pour couvrir un spectre complet (sauf lors de l'utilisation d'un flux continu de données où les coefficients sont réutilisés pour les calculs ultérieurs, ce qui a une complexité de calcul équivalente à celle de la DFT glissante ), l'algorithme de Goertzel a un ordre de complexité plus élevé que les algorithmes de transformée de Fourier rapide (FFT), mais pour calculer un petit nombre de composantes de fréquence sélectionnées, il est plus efficace numériquement. La structure simple de l'algorithme de Goertzel le rend bien adapté aux petits processeurs et aux applications embarquées.

L'algorithme de Goertzel peut également être utilisé « à l'envers » comme fonction de synthèse sinusoïdale, ce qui ne nécessite qu'une multiplication et une soustraction par échantillon généré.

Le calcul principal de l'algorithme de Goertzel se présente sous la forme d'un filtre numérique , raison pour laquelle l'algorithme est souvent appelé filtre de Goertzel . Le filtre opère sur une séquence d'entrée en cascade de deux étapes avec un paramètre , donnant la fréquence à analyser, normalisée en radians par échantillon.

La première étape calcule une séquence intermédiaire, :

La deuxième étape applique le filtre suivant à , produisant une séquence de sortie :

Le premier étage de filtrage peut être observé comme étant un filtre RII du second ordre avec une structure à forme directe . Cette structure particulière a la propriété que ses variables d'état internes sont égales aux valeurs de sortie passées de cet étage. Les valeurs d'entrée pour sont présumées toutes égales à 0. Pour établir l'état initial du filtre afin que l'évaluation puisse commencer à l'échantillon , les états du filtre se voient attribuer des valeurs initiales . Pour éviter les risques d'aliasing , la fréquence est souvent limitée à la plage de 0 à π (voir le théorème d'échantillonnage de Nyquist-Shannon ) ; l'utilisation d'une valeur en dehors de cette plage n'est pas dénuée de sens, mais équivaut à l'utilisation d'une fréquence aliasée à l'intérieur de cette plage, puisque la fonction exponentielle est périodique avec une période de 2π dans .

Le filtre de deuxième étape peut être considéré comme un filtre FIR , puisque ses calculs n'utilisent aucune de ses sorties passées.

Les méthodes de transformation en Z peuvent être appliquées pour étudier les propriétés de la cascade de filtres. La transformation en Z du premier étage de filtre donnée dans l'équation (1) est

La transformée Z du deuxième étage de filtrage donnée dans l'équation (2) est

La fonction de transfert combinée de la cascade des deux étages de filtrage est alors

Cela peut être transformé en une séquence de domaine temporel équivalente, et les termes déroulés jusqu'au premier terme d'entrée à l'index :

On peut observer que les pôles de la transformée en Z du filtre sont situés à et , sur un cercle de rayon unitaire centré sur l'origine du plan complexe de la transformée en Z. Cette propriété indique que le processus de filtrage est marginalement stable et vulnérable à l'accumulation d'erreurs numériques lorsqu'il est calculé à l'aide d'une arithmétique de faible précision et de longues séquences d'entrée. Une version numériquement stable a été proposée par Christian Reinsch .

Pour le cas important du calcul d’un terme DFT, les restrictions spéciales suivantes sont appliquées.

• Le filtrage se termine à l'index , où est le nombre de termes dans la séquence d'entrée de la DFT.
• Les fréquences choisies pour l'analyse de Goertzel sont limitées à la forme spéciale

• Le numéro d'index indiquant le « bac de fréquence » de la DFT est sélectionné dans l'ensemble des numéros d'index

En effectuant ces substitutions dans l'équation (6) et en observant que le terme , l'équation (6) prend alors la forme suivante :

Nous pouvons observer que le côté droit de l'équation (9) est extrêmement similaire à la formule de définition du terme DFT , le terme DFT pour le numéro d'indice , mais pas exactement le même. La sommation illustrée dans l'équation (9) nécessite des termes d'entrée, mais seuls les termes d'entrée sont disponibles lors de l'évaluation d'une DFT. Un expédient simple mais peu élégant consiste à étendre la séquence d'entrée avec une valeur artificielle supplémentaire . Nous pouvons voir à partir de l'équation (9) que l'effet mathématique sur le résultat final est le même que celui de la suppression du terme de la sommation, ce qui permet d'obtenir la valeur DFT souhaitée.

Il existe cependant une approche plus élégante qui évite le passage de filtre supplémentaire. À partir de l'équation (1), nous pouvons noter que lorsque le terme d'entrée étendu est utilisé dans l'étape finale,

Ainsi, l'algorithme peut être complété comme suit :

• terminer le filtre IIR après avoir traité le terme d'entrée ,
• appliquer l'équation (10) pour construire à partir des sorties précédentes et ,
• appliquer l'équation (2) avec la valeur calculée et avec celle produite par le calcul direct final du filtre.

Les deux dernières opérations mathématiques sont simplifiées en les combinant algébriquement :

Notez que l'arrêt des mises à jour du filtre à terme et l'application immédiate de l'équation (2) plutôt que de l'équation (11) manquent les mises à jour de l'état final du filtre, ce qui donne un résultat avec une phase incorrecte.

La structure de filtrage particulière choisie pour l'algorithme de Goertzel est la clé de l'efficacité de ses calculs DFT. Nous pouvons observer qu'une seule valeur de sortie est utilisée pour calculer la DFT, de sorte que les calculs pour tous les autres termes de sortie sont omis. Étant donné que le filtre FIR n'est pas calculé, les calculs de l'étape IIR , etc. peuvent être abandonnés immédiatement après la mise à jour de l'état interne de la première étape.

Cela semble laisser un paradoxe : pour compléter l'algorithme, l'étape de filtrage FIR doit être évaluée une fois en utilisant les deux dernières sorties de l'étape de filtrage IIR, tandis que pour l'efficacité du calcul, l'itération du filtre IIR rejette ses valeurs de sortie. C'est là que les propriétés de la structure de filtre à forme directe sont appliquées. Les deux variables d'état internes du filtre IIR fournissent les deux dernières valeurs de la sortie du filtre IIR, qui sont les termes requis pour évaluer l'étape de filtrage FIR.

En examinant l'équation (6), un dernier filtre IIR passe pour calculer le terme en utilisant une valeur d'entrée supplémentaire et applique un multiplicateur complexe de magnitude 1 au terme précédent . Par conséquent, et représentent une puissance de signal équivalente. Il est tout aussi valable d'appliquer l'équation (11) et de calculer la puissance du signal à partir du terme ou d'appliquer l'équation (2) et de calculer la puissance du signal à partir du terme . Les deux cas conduisent à l'expression suivante pour la puissance du signal représentée par le terme DFT :

Dans le pseudo-code ci-dessous, les données d'entrée à valeur réelle sont stockées dans le tableau x et les variables sprevstockent sprev2temporairement l'historique de sortie du filtre IIR. Ntermsest le nombre d'échantillons dans le tableau et Ktermcorrespond à la fréquence d'intérêt, multipliée par la période d'échantillonnage.

Nterms définis ici Kterm sélectionné ici ω = 2 × π × Kterme / Ntermes ; coeff := 2 × cos(ω) sprev := 0 sprev2 := 0 pour chaque index n dans la plage 0 à Nterms-1, faites s := x[n] + coeff × sprev - sprev2 sprev2 := sprev sprev := s fin puissance := sprev 2 + sprev2 2 - (coeff × sprev × sprev2) 

Il est possible d'organiser les calculs de manière à ce que les échantillons entrants soient livrés individuellement à un objet logiciel qui maintient l'état du filtre entre les mises à jour, le résultat de puissance final étant accessible une fois l'autre traitement terminé.

Le cas des données d'entrée à valeur réelle se présente fréquemment, en particulier dans les systèmes embarqués où les flux d'entrée résultent de mesures directes de processus physiques. Lorsque les données d'entrée sont à valeur réelle, les variables d'état internes du filtre peuvent sprevégalement sprev2être observées comme étant à valeur réelle. Par conséquent, aucune arithmétique complexe n'est requise dans la première étape IIR. L'optimisation pour l'arithmétique à valeur réelle est généralement aussi simple que l'application de types de données à valeur réelle appropriés pour les variables.

Une fois les calculs utilisant le terme d'entrée et les itérations de filtre terminés, l'équation (11) doit être appliquée pour évaluer le terme DFT. Le calcul final utilise une arithmétique à valeurs complexes, mais celle-ci peut être convertie en arithmétique à valeurs réelles en séparant les termes réels et imaginaires :

Par rapport à l'application du spectre de puissance, la seule différence réside dans le calcul utilisé pour terminer :

(Mêmes calculs de filtre IIR que dans l'implémentation de la puissance du signal) XKreal = sprev * cr - sprev2; XKimag = sprev * ci; 

Cette application nécessite la même évaluation du terme DFT , comme indiqué dans la section précédente, en utilisant un flux d'entrée à valeur réelle ou à valeur complexe. La phase du signal peut alors être évaluée comme

prendre les précautions appropriées pour les singularités, le quadrant, etc. lors du calcul de la fonction tangente inverse.

Les signaux complexes se décomposant linéairement en parties réelles et imaginaires, l'algorithme de Goertzel peut être calculé en arithmétique réelle séparément sur la séquence de parties réelles, ce qui donne , et sur la séquence de parties imaginaires, ce qui donne . Après cela, les deux résultats partiels à valeurs complexes peuvent être recombinés :

• Selon la théorie de la complexité computationnelle , le calcul d'un ensemble de termes DFT à l'aide d'applications de l'algorithme de Goertzel sur un ensemble de données avec des valeurs avec un « coût par opération » de présente une complexité de .

Pour calculer un seul bin DFT pour une séquence d'entrée complexe de longueur , l'algorithme de Goertzel nécessite des multiplications et des additions/soustractions dans la boucle, ainsi que 4 multiplications et 4 additions/soustractions finales, pour un total de multiplications et d'additions/soustractions. Cette opération est répétée pour chacune des fréquences.

• En revanche, l’utilisation d’une FFT sur un ensemble de données contenant des valeurs présente une complexité .

Ceci est plus difficile à appliquer directement car cela dépend de l'algorithme FFT utilisé, mais un exemple typique est une FFT radix-2, qui nécessite des multiplications et des additions/soustractions par bac DFT , pour chacun des bacs.

Dans les expressions d'ordre de complexité, lorsque le nombre de termes calculés est inférieur à , l'avantage de l'algorithme de Goertzel est évident. Mais comme le code FFT est relativement complexe, le facteur « coût par unité de travail » est souvent plus élevé pour une FFT, et l'avantage pratique favorise l'algorithme de Goertzel même pour plusieurs fois plus grand que .

En règle générale, pour déterminer si une FFT à base 2 ou un algorithme de Goertzel est plus efficace, ajustez le nombre de termes dans l'ensemble de données vers le haut jusqu'à la puissance exacte la plus proche de 2, en appelant cela , et l'algorithme de Goertzel sera probablement plus rapide si

Les implémentations FFT et les plates-formes de traitement ont un impact significatif sur les performances relatives. Certaines implémentations FFT effectuent des calculs internes de nombres complexes pour générer des coefficients à la volée, augmentant ainsi considérablement leur « coût K par unité de travail ». Les algorithmes FFT et DFT peuvent utiliser des tables de valeurs de coefficients précalculées pour une meilleure efficacité numérique, mais cela nécessite davantage d'accès aux valeurs de coefficients mises en mémoire tampon dans la mémoire externe, ce qui peut conduire à une augmentation des conflits de cache qui contrecarre une partie de l'avantage numérique.

Les deux algorithmes gagnent environ un facteur 2 en efficacité lorsqu'ils utilisent des données d'entrée à valeur réelle plutôt qu'à valeur complexe. Cependant, ces gains sont naturels pour l'algorithme de Goertzel mais ne seront pas obtenus pour la FFT sans utiliser certaines variantes d'algorithmes spécialisées dans la transformation de données à valeur réelle .