La transformation de Fourier rapide (sigle anglais : FFT ou fast Fourier transform) est un algorithme de calcul de la transformation de Fourier discrète (TFD).

Sa complexité varie en O(n log n) avec le nombre n de points, alors que la complexité de l’algorithme « naïf » s'exprime en O(n²). Ainsi, pour n = 1 024, le temps de calcul de l'algorithme rapide peut être 100 fois plus court que le calcul utilisant la formule de définition de la TFD.

C'est en 1965 que James Cooley et John Tukey publient l’article qui « lance » définitivement l'adoption massive de cette méthode en traitement du signal et dans les télécommunications. On a découvert par la suite que l'algorithme avait déjà été imaginé par Carl Friedrich Gauss en 1805, et adapté à plusieurs reprises (notamment par Lanczos en 1942) sous des formes différentes. L'algorithme a aussi été découvert indépendamment par le chanoine Lemaître^[1].

Cet algorithme est couramment utilisé en traitement numérique du signal pour transformer des données discrètes du domaine temporel dans le domaine fréquentiel, en particulier dans les oscilloscopes numériques (les analyseurs de spectre utilisant plutôt des filtres analogiques, plus précis). Son efficacité permet de réaliser des filtrages en modifiant le spectre et en utilisant la transformation inverse (filtre à réponse impulsionnelle finie). Il est également à la base des algorithmes de multiplication rapide (Schönhage et Strassen, 1971), et des techniques de compression numérique ayant mené au format d'image JPEG (1991).

Histoire

Le développement d'algorithmes rapides pour la FFT remonte au travail non publié de Carl Friedrich Gauss en 1805, alors qu'il en avait besoin pour interpoler l'orbite des astéroïdes Pallas et Junon à partir d'observations d'échantillons. Sa méthode était très similaire à celle publiée en 1965 par James Cooley et John Tukey, à qui l'on attribue généralement l'invention de l'algorithme générique moderne de la FFT. Bien que les travaux de Gauss soient antérieurs aux résultats de Joseph Fourier en 1822, il n'a pas analysé le temps de calcul et a finalement utilisé d'autres méthodes pour atteindre son objectif.

Entre 1805 et 1965, certaines versions de la FFT ont été publiées par d'autres auteurs. En 1932, Frank Yates a publié sa version, qui permettait un calcul efficace des transformées de Hadamard et de Walsh. L'algorithme de Yates est toujours utilisé dans le domaine de la conception statistique et de l'analyse des expériences. En 1942, G. C. Danielson et Cornelius Lanczos ont publié leur version permettant de calculer la FFT pour la cristallographie aux rayons X, un domaine où le calcul des transformées de Fourier présentait un formidable goulot d'étranglement. Dans le passé, de nombreuses méthodes s'étaient concentrées sur la réduction de la constante ${\displaystyle \lambda }$ pour laquelle le nombre d'opérations était équivalent à ${\displaystyle \lambda n^{2}}$ en tirant parti des "symétries". Mais Danielson et Lanczos ont réalisé que l'on pouvait utiliser la "périodicité" pour réduire encore cette constante, bien que, comme Gauss, ils n'aient pas analysé que cela conduisait à un algorithme en ${\displaystyle O(n\log(n))}$.

James Cooley et John Tukey ont redécouvert de manière indépendante ces algorithmes antérieurs et ont publié en 1965 une version plus générale de la FFT applicable lorsque ${\displaystyle n}$ est composite et pas nécessairement une puissance de 2, ainsi qu'une analyse de complexité en ${\displaystyle O(n\log(n))}$. Tukey a eu cette idée lors d'une réunion du comité consultatif scientifique du président Kennedy, où l'on discutait de la détection des essais nucléaires de l'Union soviétique en installant des capteurs pour encercler le pays de l'extérieur. Les essais nucléaires avaient en effet lieu sous terre^[2], et la seule manière de les détecter était d'utiliser une transformée de Fourier. En discutant avec Tukey, Richard Garwin a reconnu l'applicabilité générale de l'algorithme non seulement aux problèmes de sécurité nationale, mais aussi à un large éventail de problèmes, dont un qui l'intéressait immédiatement, à savoir la détermination des périodicités des orientations du spin dans un cristal tridimensionnel d'hélium 3^[3]. Comme Tukey ne travaillait pas chez IBM, la brevetabilité de l'idée a été mise en doute et l'algorithme est tombé dans le domaine public, ce qui, grâce à la révolution informatique de la décennie suivante, a fait de la FFT l'un des algorithmes indispensables du traitement numérique du signal.

Formulation mathématique

Soit ${\displaystyle \lambda _{0},\dots ,\lambda _{n-1}}$des nombres complexes. La transformée de Fourier discrète de ${\displaystyle \lambda _{0},\dots ,\lambda _{n-1}}$ est définie par la formule suivante :

${\displaystyle \sum _{k=0}^{n-1}\lambda _{k}e^{-i2\pi jk/n}}$

pour j = 0, ... , n-1. De manière équivalente, soit ${\displaystyle P=\sum _{k=0}^{n-1}\lambda _{k}X^{k}}$ un polynôme à coefficients complexes. La transformée de Fourier discrète de ${\displaystyle P}$ est définie par la formule suivante :

${\displaystyle (P(\omega ^{j}))_{0\leq j\leq n-1}}$, où ${\displaystyle \omega =e^{-{\frac {2\pi i}{n}}}}$

ou en notation matricielle :

${\displaystyle {\begin{array}{l}{\begin{pmatrix}f_{0}\\f_{1}\\f_{2}\\\vdots \\f_{n-1}\end{pmatrix}}={\begin{pmatrix}1&1&1&\cdots &1\\1&\omega &\omega ^{2}&\cdots &\omega ^{n-1}\\1&\omega ^{2}&\omega ^{4}&\cdots &\omega ^{2(n-1)}\\\vdots &\vdots &\vdots &\ddots &\vdots &\\1&\omega ^{n-1}&\omega ^{2(n-1)}&\cdots &\omega ^{(n-1)^{2}}\end{pmatrix}}\end{array}}{\begin{pmatrix}\lambda _{0}\\\lambda _{1}\\\lambda _{2}\\\vdots \\\lambda _{n-1}\end{pmatrix}}}$, où l'on reconnaît une matrice de Vandermonde ${\displaystyle V(1,\omega ,...,\omega ^{n-1})}$.

Évaluer naïvement ${\displaystyle P}$ coûte (n – 1)² produits complexes et n(n – 1) sommes complexes alors que seuls ${\displaystyle {\frac {n}{2}}\left(\log _{2}(n)-2\right)}$ produits et ${\displaystyle n\log _{2}(n)}$ sommes sont nécessaires avec la version rapide. En général, de tels algorithmes dépendent de la factorisation de ${\displaystyle n}$ mais contrairement à une idée répandue, il y a des transformées de Fourier rapides de complexité ${\displaystyle O(n\log(n))}$ pour tous les ${\displaystyle n}$, même pour les nombres premiers.

Il est possible de générer la transformation inverse de la même manière pour la version rapide. En effet l'inverse de la matrice de Vandermonde ${\displaystyle V(1,\omega ,...,\omega ^{n-1})}$ est simplement ${\displaystyle {\frac {1}{n}}V(1,\omega ^{-1},...,\omega ^{-n+1})}$. Il suffit donc d'appliquer l'algorithme à cet inverse.

Algorithme de Cooley-Tukey

Il s'agit d'un algorithme fréquemment utilisé pour calculer la transformation de Fourier discrète. Il se base sur une approche de type « diviser pour régner » par le biais d'une récursion. Celle-ci subdivise une transformation de Fourier discrète d'une taille composite ${\displaystyle n=n_{1}n_{2}}$ en plusieurs transformées de Fourier discrètes de tailles inférieures ${\displaystyle n_{1}}$ et ${\displaystyle n_{2}}$. Cet algorithme nécessite ${\displaystyle O(n)}$ multiplications par des racines d'unité, plus communément appelés facteurs de rotation. En général, les mises en code essaient d'éviter une récursion pour des questions de performance. Il est aussi possible de mélanger plusieurs types d'algorithme lors des subdivisions.

C'est en 1965 que James Cooley et John Tukey publient cette méthode^[4],[5] mais il a été découvert par la suite que l'algorithme avait déjà été inventé par Carl Friedrich Gauss en 1805^[6] et adapté à plusieurs reprises sous des formes différentes^[7].

Détail d'un cas classique : ${\displaystyle n_{1}=2}$

L'utilisation la plus classique de l'algorithme de Cooley-Tukey est une division de la transformation en deux parties de taille identique ${\displaystyle n/2}$ et ceci à chaque étape. Cette contrainte limite les tailles possibles, puisque celles-ci doivent être des puissances de deux. Néanmoins on peut traiter n'importe quel polynôme puisqu'un polynôme de degré ${\displaystyle k}$ peut être vu comme un polynôme de ${\displaystyle \mathbb {C} _{2^{n}}[X]}$ où ${\displaystyle k\leq 2^{n}}$.

Principe de l'algorithme

On décompose le polynôme ${\displaystyle P=\sum _{k=0}^{n-1}\lambda _{k}X^{k}}$ en ${\displaystyle P=P_{pair}(X^{2})+X\cdot P_{impair}(X^{2})}$, où ${\displaystyle P_{pair}=\sum _{k=0}^{\frac {n-1}{2}}\lambda _{2k}X^{k}}$ et ${\displaystyle P_{impair}=\sum _{k=0}^{\frac {n-2}{2}}\lambda _{2k+1}X^{k}}$. On applique ensuite l'algorithme en ${\displaystyle P_{pair}}$ et en ${\displaystyle P_{impair}}$, pour en déduire l'évaluation de ${\displaystyle P}$ en ${\displaystyle \omega ^{j},0\leq j\leq n-1}$.

Pseudo-code

Entrée : un polynôme ${\displaystyle P=\sum _{k=0}^{n-1}\lambda _{k}X^{k}}$ et ${\displaystyle \omega }$ une racine primitive ${\displaystyle n}$-ième de l'unité.
Sortie : ${\displaystyle P(\omega ^{0}),...,P(\omega ^{n-1})}$.
fonction FFT(${\displaystyle P,\omega }$) :
|  si ${\displaystyle P}$ est constant :
|  |  retourner ${\displaystyle P(1)}$
|  sinon
|  |  calculer ${\displaystyle P_{pair}}$ et ${\displaystyle P_{impair}}$
|  |  calculer ${\displaystyle P_{pair}(\omega ^{0}),P_{pair}(\omega ^{2}),...,P_{pair}(\omega ^{n-2})}$ et ${\displaystyle P_{impair}(\omega ^{0}),P_{impair}(\omega ^{2}),...,P_{impair}(\omega ^{n-2})}$ via FFT(${\displaystyle P_{pair},\omega ^{2}}$) et FFT(${\displaystyle P_{impair},\omega ^{2}}$) respectivement
|  |  pour k allant de 0 à n-1 :
|  |  |  ${\displaystyle P(\omega ^{k})=P_{pair}(\omega ^{2k})+\omega ^{k}\cdot P_{impair}(\omega ^{2k})}$
|  |  retourner  ${\displaystyle P(\omega ^{0}),...,P(\omega ^{n-1})}$

Cet algorithme a une complexité en ${\displaystyle O(n\log(n))}$.

Schéma explicatif

Schéma correspondant à la FFT (DFT en anglais) d'un polynôme de degré 7.

Dans le schéma ci-contre, on a représenté la FFT d'un polynôme ${\displaystyle P=\sum _{k=0}^{7}x[k]X^{k}}$. On obtient à la fin des valeurs ${\displaystyle X[k]=P(\omega ^{k})}$ Les polynômes intermédiaires obtenus sont ${\displaystyle E}$ (pour even) et ${\displaystyle O}$ (pour odd) qui représentent ${\displaystyle P_{pair}}$ et ${\displaystyle P_{impair}}$.

On obtiendra ensuite des polynômes ${\displaystyle E_{e}}$ et ${\displaystyle E_{o}}$ à partir de ${\displaystyle E}$, ainsi que ${\displaystyle O_{e}}$ et ${\displaystyle O_{o}}$ à partir de ${\displaystyle O}$, et on itérera le procédé. Ainsi ${\displaystyle E_{e}}$ est formé des coefficients ${\displaystyle x[0]}$ et ${\displaystyle x[4]}$, ${\displaystyle E_{o}}$ des coefficients ${\displaystyle x[2]}$ et ${\displaystyle x[6]}$, et de même pour les impairs.

Les flèches vers les valeurs ${\displaystyle X[k]}$ signifient que la valeur de départ de la flèche influence la valeur de ${\displaystyle X[k]}$ : par exemple la valeur de ${\displaystyle X[1]}$ dépend de celles de ${\displaystyle E[1]}$ et ${\displaystyle O[1]}$.

Autres algorithmes

Il existe d'autres algorithmes qui permettent de calculer la transformée de Fourier discrète. Pour une taille n = n₁n₂, avec des nombres premiers entre eux n₁ et n₂, il est possible d'utiliser l'algorithme PFA (Good-Thomas) basé sur le théorème des restes chinois. Le PFA est similaire à celui de Cooley-Tukey.

L'algorithme de Rader-Brenner^[8] est aussi une variante de Cooley-Tukey avec des facteurs de rotation purement imaginaires qui améliorent les performances en réduisant le nombre de multiplications mais au détriment de la stabilité numérique et une augmentation du nombre d'additions. Les algorithmes qui procèdent aussi par des factorisations successives sont celui de Bruun (en) et l'algorithme QFT. Les versions originales travaillent sur des fenêtres dont la taille est une puissance de 2 mais il est possible de les adapter pour une taille quelconque. L'algorithme de Bruun considère la transformée de Fourier rapide comme une factorisation récursive du polynôme zⁿ – 1 en des polynômes à coefficients réels de la forme z^m – 1 et z^2m + az^m + 1.

L'algorithme de Winograd factorise zⁿ – 1 en un polynôme cyclotomique, dont les coefficients sont souvent –1, 0 ou 1, ce qui réduit le nombre de multiplications. On peut voir cet algorithme comme la version optimale en termes de multiplications. Winograd a montré que la transformée de Fourier discrète peut être calculée avec seulement O(n) multiplications, et ce minorant est atteint pour les tailles qui sont des puissances de 2. Toutefois, des additions supplémentaires sont nécessaires, ce qui peut être pénalisant sur les processeurs modernes comportant des unités arithmétiques performantes.

L'algorithme de Rader (en) est quant à lui destiné aux fenêtres dont la taille est un nombre premier. Il profite de l'existence d'une génératrice pour le groupe multiplicatif modulo n. La transformation discrète dont la taille est un nombre premier s'exprime ainsi comme une convolution cyclique de taille n – 1. On peut ensuite la calculer par une paire de transformation de Fourier rapide.

Finalement, un autre algorithme destiné aux transformations avec des tailles qui sont des nombres premiers est due à Leo Bluestein, en 1968. On l'appelle plus souvent l'algorithme chirp Z (en). Ici encore, la transformation est vue comme une convolution dont la taille est identique à la fenêtre originale. On utilise à cet effet l'identité de polarisation : jk = –(j – k)²/2 + j²/2 + k²/2.

Algorithmes spécialisés dans le traitement de données réelles ou/et symétriques

Dans beaucoup d'applications, les données en entrée de la transformation discrète de Fourier sont uniquement des nombres réels. Dans ce cas, les sorties satisfont la symétrie suivante : ${\displaystyle f_{n-j}=f_{j}^{*},}$ et des algorithmes efficaces ont été conçus pour cette situation, par exemple celui de Sorensen en 1987^[9]. Une approche possible consiste à prendre un algorithme classique comme celui de Cooley-Tukey et à enlever les parties inutiles dans le calcul. Cela se traduit par un gain de 50 % en mémoire et en vitesse de calcul. Alternativement, il est possible d'exprimer une transformation discrète sur des nombres réels (avec une taille paire) en une transformation avec des nombres complexes mais dont la taille a été divisée par deux (les parties imaginaires sont les éléments impairs et les parties réelles sont les éléments pairs) suivie d'un décodage dont la complexité est de l'ordre de O(n) opérations.

On pensait que les transformations avec des nombres réels pouvaient être plus efficacement calculées via une transformation de Hartley discrète (en) mais il a été prouvé par la suite qu'une transformation de Fourier discrète modifiée pouvait être plus efficace que la même transformation de Hartley. L'algorithme de Bruun était un candidat pour ces transformations mais il n'a pas eu la popularité escomptée.

Il existe encore d'autres variantes pour les cas où les données sont symétriques (c’est-à-dire des fonctions paires ou impaires) avec un gain supplémentaire de 50%. Dans ce cas, on utilise une transformée en cosinus discrète.

Problèmes numériques et approximations

Par leur nature analytique, tous les algorithmes proposés ci-dessus calculent la transformée sans aucune erreur. Toutefois, il existe des algorithmes qui peuvent s'accommoder d'une marge d'erreur pour accélérer les calculs. En 1999, Edelman et al. proposent une approximation à la transformée de Fourier rapide^[10]. Cette version est destinée à une mise en œuvre en parallèle. Une approximation basée sur les ondelettes est proposée en 1996 par Guo et Burrus^[11] et tient compte de la répartition dans les entrées/sorties. Un autre algorithme a encore été proposé par Shentov et al. en 1995^[12]. Seul l'algorithme de Edelman fonctionne bien avec n'importe quel type de données, il profite de la redondance dans la matrice de Fourier plutôt que de la redondance dans les données initiales.

Toutefois, même les algorithmes décrits de manière analytique présentent des erreurs lorsqu'ils sont codés avec des virgules flottantes dont la précision est limitée. L'erreur est cependant limitée. Une borne supérieure d'erreur relative pour Cooley-Tukey est donnée par ${\displaystyle O(\epsilon \log(n))}$ alors qu'elle est de ${\displaystyle O(\epsilon n^{3/2})}$ pour la formulation triviale de la transformée de Fourier discrète^[13]. Le terme ${\displaystyle \epsilon }$ représente ici la précision relative en virgule flottante. En fait, l'erreur quadratique moyenne est encore plus limitée avec seulement ${\displaystyle O(\epsilon {\sqrt {\log(n)}})}$ pour Cooley-Tukey et ${\displaystyle O(\epsilon {\sqrt {n}})}$ pour la version triviale. Il ne faut malgré tout pas oublier que la stabilité peut être perturbée par les différents facteurs de rotation qui interviennent dans les calculs. Un manque de précision sur les fonctions trigonométriques peut fortement augmenter l'erreur. L'algorithme de Rader est par exemple nettement moins stable que celui de Cooley-Tukey en cas d'erreurs prononcées.

Avec une arithmétique en virgule fixe, les erreurs s'accumulent encore plus vite. Avec Cooley-Tukey, l'augmentation de l'erreur quadratique moyenne est de l'ordre de ${\displaystyle O(\epsilon {\sqrt {n}})}$. De plus, il faut tenir compte de la magnitude des variables lors des différentes étapes de l'algorithme.

Il est possible de vérifier la validité de l'algorithme grâce à une procédure qui vise à déterminer la linéarité et d'autres caractéristiques de la transformation sur des entrées aléatoires^[14].

Exemple d'application : la multiplication de polynômes

Quand on veut multiplier deux polynômes représentés par la liste de leurs coefficients, la multiplication se fait en ${\displaystyle O(n^{2})}$. En revanche, si on les représente par une liste ${\displaystyle (x_{k},P(x_{k}))_{0\leq k\leq n-1}}$ d'évaluations en des points tous distincts, la multiplication se fait en ${\displaystyle O(n)}$ (en supposant que l'on a assez de valeurs pour avoir ${\displaystyle P}$ défini de manière unique). On a donc intérêt à chercher un algorithme permettant de passer rapidement d'une de ces représentations à l'autre.

Or la FFT permet l'évaluation d'un polynôme en ${\displaystyle O(n\log(n))}$ et, inversement, son interpolation en ${\displaystyle n+1}$ points (comme explicité dans la partie « Formulation mathématique »). Ainsi, la multiplication de deux polynômes représentés par la liste de leurs coefficients s'effectue en ${\displaystyle O(n\log(n))}$.

Notes et références

Voir aussi