La Transformée de Fourier Rapide (FFT)

La rapidité d'exécution du calcul est essentielle lorsqu'il s'agit de traiter un signal en temps réel. Elle déterminera la fréquence maximale qui puisse être analysée. Repartons de la formulation complexe :

Posons

Comme nous l'avons vu et expérimenté précédemment, l'application directe de cette formulation consiste à faire évoluer les indices n et k de 0 à N ce qui conduit à effectuer opérations ce qui représente vite, même avec des processeurs rapides, un temps de calcul important (plusieurs ms sur un microcontrôleur tournant à 20MHz, limitant la fréquence d'analyse à quelques dizaines de Hz). Mais on peut réduire ce nombre d'opérations en appliquant quelques astuces algorithmiques :

Posons :

1Valeurs remarquables de :

1.1

Démonstration. de (2) :

1.2 pour

Rappels :

Tableau 1.

Démonstration. de (3):

1.3

Démonstration. de (4) :

1.4

Démonstration. de (5) :

1.5

Démonstration. de (6) :

2Algorithme de Transformée rapide

2.1Séparation en deux moitiés :

Séparons les échantillons d'ordre pair et impair, en tenant compte de ce que nous avons vu plus haut, il vient:

Nous voyons que la transformée de Fourier de taille N est devenue la somme de deux transformées de Fourier de taille N/2. Si le calcul de départ prenait opérations, maintenant il ne nécessite plus que soit environ la moitié.

Si N est une puissance de 2, on peut répéter le fractionnement jusqu'à l'obtention d'une somme de base 2, l'ensemble des calculs ne nécessitant plus que opérations ce qui peut être 100 fois moins lourd.

Nous remarquons également que le facteur se retrouve à l'identique dans les termes de rang pair et impair de même ordre et ne sera donc calculé qu'une fois sur deux.

Le facteur qui ne dépend pas de (ni de ) (raison pour laquelle on l'a « sorti » de la somme ) figurant dans les termes de rang impair, pourra bénéficier de la propriété(4).

L'algorithme le plus célèbre qui en découle est celui inventé en 1805 par Carl Friedrich Gauss pour calculer les trajectoires des astéroïdes Pallas et Juno, puis réinventé en 1965 par James Cooley (IBM) et John Tukey (Princeton) qui l'adaptèrent au calcul sur ordinateur et lui donnèrent leurs noms.

2.1.1transformée de Fourier des échantillons de rang pair :

est la transformée de Fourier d'ordre des échantillons de rang pair.

(F comme Fourier, indice P comme Pair, et indice 2 comme ordre N/2)

Remarque, d'après (5):

Démonstration. de (9) :

2.1.2transformée de Fourier des échantillons de rang impair :

est la transformée de Fourier d'ordre des échantillons de rang impair

remarque, toujours d'après (5) :

La transformée de Fourier complète devient :

Nous pouvons maintenant décrire l'algorithme de Cooley-Tukey qui concrétise les équations qui précèdent. Toutefois toutes les descriptions que j'ai pu voir (en français) de cet algorithme (agrémentés de croquis de « papillons ») m'ont paru claires comme du jus de chique. Je vais donc essayer de rendre la chose limpide et j'avoue que je ne sois par sûr d'y parvenir ! Je n'ai moi-même compris la chose qu'en lisant des travaux parus en anglais.

Je pense que le plus simple consiste à se fixer une valeur de k petite, 8 par exemple, d'écrire dans un tableau les 8 lignes de calculs correspondantes, puis de voir les simplifications qu'on peut y apporter.

soit donc N=8, et N/2=4.

Appliquons le résultat connu (9) : ainsi que (11) :

Les coefficients qui multiplient les seconds termes sont appelés en anglais « twiddle factor » ce qui est assez cocasse surtout si vous le traduisez mot à mot... Nous allons nous en occuper maintenant.

reprenons (14)

avec n = 1, il vient:

Nous pouvons alors diminuer à nouveau le nombre de termes différents: (on fait sauter à )

N'apparaissent plus que les transformées de Fourier d'ordre N/2 à et à soit quatre TF de rang pair, quatre TF de rang impair et quatre « facteurs de tripatouilages » (), ces derniers pouvant être pré-calculés et mémorisés.

Récapitulons afin d'être clair:

c'est une des composantes de la transformée de Fourier que l'on cherche à obtenir. C'est une raie dans le domaine fréquenciel.

c'est une transformée partielle d'indice k calculée sur les échantillons pairs du signal temporel.

Elle nécessite de calculer une somme basée sur tous les échantillons pairs (balayage par ).

c'est une transformée partielle d'indice k calculée sur les échantillons impairs du signal temporel.

Elle nécessite de calculer une somme basée sur tous les échantillons impairs.

Les transformées partielles et sont obtenues en prenant successivevent et en balayant toutes les valeurs de (donc tous les échantillons) pour chacune de ces valeurs de .

Donc le calcul de chaque valeur de nécessite d'utiliser tous les échantillons (pairs et impairs) du signal au même titre que les transformées partielles.

2.2Représentations graphiques

Voici une représentation graphique de ces calculs. Les moulinettes qui brassent tous les échantillons se situent dans les gros pavés oranges.

Remarquons que les mécanismes de ces deux moulinettes sont identiques, la seule différence se trouve dans les échantillons qu'on leur donne à traiter, pairs pour celle du haut, impairs pour celle du bas.

Remarquons également que celle du haut traite un nombre pair d'échantillons (quatre), on peut donc la remplacer par deux moulinettes traitant chacune la moitié de ces échantillons séparés en échantillons pairs et impairs, leurs sorties une fois multiplexés suivant le même schéma que celui utilisé sur la figure actuelle.

Et la même remarque vaut pour la moulinette du bas.

Allez, hop au boulot, traçons ce nouveau schéma qui devrait nous permettre de diviser par deux le nombre d'opérations mathématiques à effectuer.

Oui mais pourquoi les échantillons qui entrent dans les boites à gauche ne sont plus disposés dans le même ordre? 0, 4, 2, 6 par exemple pour les quatre du haut au lieu de l'ordre naturel des entiers pairs, à savoir 0, 2, 4, 6 ???

Et bien c'est parce qu'il faut renuméroter les échantillons du premier diagramme :

le premier échantillon qui entrait dans la boite était le numéro 0

le second était le numéro 2

le troisième était le numéro 4

le quatrième était le numéro 6.

Ces quatre échantillons étant la totalité de ce qui entre dans la TF, on les renumérote par ordre d'arrivée, puis on prend les deux pairs pour la petite moulinette du haut puis les deux impairs pour la petite moulinette suivante. Les deux pairs par ordre d'arrivée ce sont les num 0 et 4, les deux impairs ce sont les num 2 et 6.

Voilà, je pense que la seule difficulté pour appréhender la recette réside dans ce point. Le reste c'est de la cuisine ordinaire.

Nous pouvons aller encore plus loin en fractionnant les TFR d'ordre N/4 en deux TFR d'ordre N/8 et ne traitant plus chacune q'un seul échantillon.

Qu'avons-nous gagné à fractionner ainsi les transformées de Fourier ? On gagne le fait que les sommes qui brassent les échantillons s'effectuent avec des boucles bien plus petites, ce qui conduit à au total au lieu de opérations (plus les quelques multiplications par les facteurs de tripatouillages...)

Notons qu'avec un N bien plus grand, (mais toujours égal à une puissance de 2) on également multiplier le fractionnement jusqu'à l'obtention de cellules de base, ce qui divise d'autant le nombre d'opérations. Le diagramme devient difficile à tracer, toutefois l'ordinateur n'aura aucun mal à le suivre avec des procédures qui peuvent être récursives, mais pas obligatoirement. Ainsi pour N=1024 on divisera le nombre d'opérations nécessaires (par la partie gauche du diagramme) par... 1024. Mais au total, en comptant les opérations de multiplications par les Wn et les additions, le gain sera moins important.

Remarquons que pour N=8, la TF paire d'ordre N/8 devient:

et la la TF impaire d'ordre N/8 devient:

La somme se réduit ainsi à sa plus simple expression ! (On pouvait s'y attendre, les « moulinettes à échantillons » n'ayant plus alors qu'un seul échantillon à mouliner !! )

Ce qui nous donne au final le diagramme suivant :

... dans lequel les grandes sommes ont totalement disparu.

Rappel :

Les multiplications par les facteurs de tripatouillages (« twiddle factor ») sont donc des multiplications de nombres complexes avec partie réelle et partie imaginaire. On voit qu'il définissent une rotation discrète dans le plan complexe.

Ces facteurs ne seront calculés qu'une fois et gardés en mémoire dans un tableau.

Dans le programme en langage C que je vais vous proposer prochainement, j'ai créé une classe (un objet) « Complexe » pourvue d'un opérateur de multiplication interne. Pour le portage sur microcontrôleur, il sera sans doute plus pertinent de décomposer les multiplications complexes en multiplications réelles effectuées par des fonctions.

Remarque 1. (pouvant être utilisée lors de l'implantation sur microcontrôleur)

Les échantillons sont présentés en entrée suivant l'ordre 0,4,2,6,1,5,3,7. Représentons ces nombres dans un tableau avec leur écriture en binaire puis ces valeurs binaires lues de droite à gauche (inversion des bits):

Nous voyons qu'il faut inverser les bits des nombres entiers naturels écrits en binaire pour obtenir les numéros des échantillons dans l'ordre qui convient au traitement.

Remarque 2. Concernant les « twiddle factor » :

Voici un petit tableau qui repésente les valeurs de twiddle factors pour une TFR avec 16 échantillons (c'est l'équivalent des diagrammes vus plus haut, pour 16 lignes au lieu de 8) :

La seconde colonne indique les numéros d'échantillons à traiter après mise dans l'ordre « bits reverse »

Les colonnes qui suivent donnent les valeurs « k » pour les à appliquer.

On voit vite la logique des opérations: en partant de la droite vers la gauche, il faut prendre les entiers naturels (n), puis une valeur sur deux (soit 2n), puis de nouveau une valeur sur deux (2x2 n)... etc... En numérotant « c » les colonnes de 0 à 3 en partant de la gauche, la valeur des « k » pour la colonne « c » est donc:

avec n = entiers naturels 1, 2, 3...

Et il faut prendre des groupes de N/2 puis N/4 puis N/8... k c'est à dire un nombre de k.

Il fallait préciser tout cela car notre programme traitera 1024 échantillons en dix étapes.

Voici le coeur de ce programme en C++ et Qt4, la fonction « calcul_fourier() » qui parcourt les dix étapes, traite pour chacune de ces étapes les 1024 lignes d'échantillons, c'est a dire réalise les multiplications complexes par les twiddles factors ainsi que les additions en suivant les croisements en « papillons » vus précédemment. Comme vous pouvez le constater cette fonction comme tous le reste de mon programme est commentée en français, c'est pas du luxe !

Et voici une saisie d'écran de son exécution sous Linux Mint:

J'ai écrit une page présentant ce programme, avec accès à la totalité du code source :

L'étape suivante consistera à implanter tout ça sur un microcontrôleur ATmega avec acquisition du signal en temps réel...




Silicium628
22178