HarmonyFidelisHarmonyFidelis
Connexion
ActualitésGrands projetsActeursAcadémie

1965 : la FFT a rendu l’analyse de Fourier praticable à grande échelle

Publié en avril 1965, l’article de cinq pages de James Cooley et John Tukey a contribué à rendre courant un calcul exigeant : déterminer les composantes fréquentielles d’une suite. Sa transformée de Fourier rapide réorganisait le travail en calculs plus petits et réutilisables. Pour une longueur qui est une puissance de deux, le travail arithmétique croît approximativement comme le nombre de points multiplié par son logarithme, plutôt que comme son carré. Des mathématiciens avaient auparavant découvert des méthodes apparentées. Le changement durable est venu de leur mise en œuvre efficace et de leur diffusion dans le calcul électronique. La rapidité rend de plus grands calculs abordables ; elle n’améliore pas les mesures d’origine et ne supprime pas les limites numériques ou d’échantillonnage.

Source: Cooley and Tukey — Mathematics of Computation, April 1965

1965 : la FFT a rendu l’analyse de Fourier praticable à grande échelle

Couverture : illustration conceptuelle générée par IA d’un signal synthétique et de ses composantes fréquentielles. Elle ne représente ni l’ordinateur d’origine ni des données expérimentales.

Entendre les composantes d’un signal

Imaginez que vous enregistriez un accord musical. Le microphone vous fournit une suite d’amplitudes, alors que votre question porte sur les notes qu’elle contient. L’analyse de Fourier établit un lien mathématique entre ces descriptions. Des questions semblables se posent lorsqu’on étudie des vibrations ou des motifs dans une image.

La transformée de Fourier discrète, ou DFT, prend une suite finie et l’exprime par des composantes fréquentielles discrètes. Une transformée de Fourier rapide, ou FFT, est une manière efficace de calculer cette même transformée. Il ne s’agit pas d’un nouveau type de spectre. La documentation technique de NumPy décrit l’entrée comme un signal dans le domaine temporel et la sortie comme sa représentation dans le domaine fréquentiel. Documentation de la DFT

La difficulté pratique tient aux répétitions. Un calcul direct combine chaque point d’entrée avec chaque fréquence de sortie. Si l’on double le nombre de points, le travail arithmétique dominant devient quatre fois plus grand. Lorsqu’une question scientifique exige des milliers de points, une formule mathématiquement simple peut donc devenir un obstacle au calcul.

Cooley et Tukey ont présenté une factorisation générale et une mise en œuvre commode pour les puissances de deux, avec notamment un moyen de réutiliser l’espace de stockage du tableau d’entrée. Leur article décrit un programme pour IBM 7094 et des temps de calcul, au-delà d’une simple proposition abstraite. Le manuscrit a été reçu le 17 août 1964 et publié dans le numéro d’avril 1965 de Mathematics of Computation. Le jour de publication n’est pas établi ici. Article original, copie universitaire

Une avancée qui s’inscrit dans une longue histoire

L’idée n’est pas apparue de rien en 1965. Les chercheurs en histoire Michael Heideman, Don Johnson et Sidney Burrus ont examiné des textes antérieurs et identifié une décomposition équivalente dans les travaux de Gauss. Ils ont déduit une date de rédaction située vers 1805 ; le manuscrit lui-même n’était pas explicitement daté et a été publié à titre posthume en 1866. Leur étude traite aussi de méthodes ultérieures, dont les travaux de Danielson et Lanczos de 1942 et l’approche de Good par facteurs premiers. Les restrictions de Good diffèrent de celles de la factorisation générale de Cooley–Tukey. Recherche historique et copie d’auteur accessible

Dans leurs souvenirs publiés en 1993, Cooley et Tukey décrivent le rôle de Richard Garwin dans la mise en relation de ces travaux et l’encouragement à leur développement. Ils reconnaissent aussi des découvertes antérieures et expliquent pourquoi les grands calculs électroniques ont rendu cette approche précieuse. Il s’agit d’un récit rétrospectif de participants, avec les limites propres aux souvenirs. Récit de Cooley et Tukey

Le tournant durable a été la diffusion, dans la pratique, d’un calcul réutilisable. Le récit de Princeton consacré au jalon reconnu par l’IEEE le relie au calcul scientifique, au traitement du signal, à l’imagerie médicale et à la transmission de données. Ces applications dépendent aussi d’instruments, de modèles physiques et d’autres algorithmes ; la FFT a apporté une puissante brique de calcul. Cette reconnaissance ultérieure étaye l’importance historique, sans faire d’une commémoration récente la date de découverte. Princeton et le jalon IEEE

Quatre points : comprendre la réutilisation avant la formule générale

Prenons la suite sans dimension [1, 2, 3, 4]. Cet exemple volontairement petit est notre calcul pédagogique, pas une série de données de l’article de 1965. Notons i l’unité imaginaire, avec i² = −1. La DFT sans facteur d’échelle et à exposant négatif donne quatre sorties :

Indice de fréquenceSortie
010
1−2 + 2i
2−2
3−2 − 2i

La première sortie additionne les quatre entrées. Les autres combinent des poids positifs, négatifs et imaginaires. Une sortie complexe contient deux composantes ; son module et sa phase décrivent ensemble une composante fréquentielle.

Séparons maintenant les entrées d’indices pairs [1, 3] des entrées d’indices impairs [2, 4]. Pour chaque paire, il suffit de calculer une somme et une différence :

  • Paire d’indices pairs : E = [4, −2].
  • Paire d’indices impairs : O = [6, −2].

Combinons ces résultats plus courts. Le premier papillon produit 4 + 6 = 10 et 4 − 6 = −2. Le second commence par faire tourner le résultat impair en le multipliant par −i, ce qui donne (−i)(−2) = 2i. Ses deux sorties sont −2 + 2i et −2 − 2i.

« Papillon » désigne cette opération appariée de somme et de différence. Ses lignes croisées servent à représenter l’organisation du calcul, pas une connexion physique entre des ondes. Les deux sorties réutilisent les mêmes résultats intermédiaires. Pour une taille plus grande, chaque transformée plus courte peut à nouveau être divisée : la réutilisation forme alors une hiérarchie, au lieu d’un raccourci ponctuel.

De l’exemple au changement d’échelle

Pour N valeurs d’entrée, définissons notre convention par

Xk=∑n=0N−1xne−2πink/N.X_k = \sum_{n=0}^{N-1} x_n e^{-2\pi i nk/N}.Xk​=n=0∑N−1​xn​e−2πink/N.

Ici, n indexe les entrées, k indexe les sorties, et tous deux vont de zéro à N − 1. Les facteurs sont des rotations complexes de module unité. Notre transformée directe n’a pas de facteur de normalisation ; une transformée inverse utilise l’exposant positif et un facteur 1/N. D’autres conventions existent. L’article original part d’une somme de Fourier à exposant positif : il faut donc comparer explicitement les signes et la normalisation. Convention moderne

Pour N pair, séparons la somme entre les indices d’entrée pairs et impairs. Notons Eₖ et Oₖ leurs transformées de longueur N/2, et posons W_N = exp(−2πi/N). Pour 0 ≤ k < N/2, la reconstruction s’écrit

Xk=Ek+WNkOk,Xk+N/2=Ek−WNkOk.X_k = E_k + W_N^k O_k, \qquad X_{k+N/2} = E_k - W_N^k O_k.Xk​=Ek​+WNk​Ok​,Xk+N/2​=Ek​−WNk​Ok​.

La seconde égalité utilise le changement de signe de la rotation à mi-parcours du cercle. Nous calculons donc deux transformées plus courtes et effectuons un nombre de combinaisons proportionnel à N. En répétant cette division pour N = 2ᵐ, on obtient m = log₂N niveaux.

Chaque niveau demande une quantité de travail proportionnelle à la longueur de la suite. Le total croît donc comme N log₂N. La notation O(N log N) décrit cette croissance, pas un décompte exact des instructions du processeur ou des secondes. La construction pédagogique en base 2 est limitée aux puissances de deux ; les factorisations générales de Cooley–Tukey peuvent utiliser d’autres longueurs composées.

Sur un ordinateur réel, la rapidité dépend aussi de l’organisation des données, des transferts en mémoire et de la mise en œuvre. Un nombre plus faible d’opérations arithmétiques ne suffit pas, à lui seul, à établir un gain précis sur le temps écoulé avec votre appareil.

Ce qu’un calcul plus rapide ne corrige pas

Si des mesures régulièrement espacées ont un intervalle Δt en secondes, la fréquence d’échantillonnage est de 1/Δt hertz et l’écart entre les cases fréquentielles de 1/(NΔt) hertz. L’interprétation des indices de sortie les plus élevés exige de tenir compte de l’ordre choisi pour les fréquences positives et négatives. Conventions fréquentielles

Un échantillonnage fini limite toujours ce que l’on peut déduire. Ici, t est le temps en secondes, f la fréquence du signal en hertz, et fₛ = 1/Δt la fréquence d’échantillonnage en hertz. Par exemple, les échantillons de cos(2πft) aux instants t = n/fₛ restent identiques si l’on remplace f par f + fₛ : la phase ajoutée vaut 2πn. Aucun algorithme plus rapide ne peut distinguer ces deux signaux à partir de ces seuls échantillons. Le bruit et les modèles de mesure inadéquats restent eux aussi des problèmes de mesure.

L’arithmétique numérique introduit une autre limite. Les facteurs de rotation complexes nécessitent généralement des approximations numériques, et changer l’ordre des opérations peut modifier les arrondis. L’analyse de l’exactitude publiée par FFTW souligne l’importance de facteurs de rotation précis et décrit des différences entre mises en œuvre et environnements. L’équivalence mathématique ne promet pas des configurations de bits identiques en virgule flottante. Analyse de l’exactitude

Vérifier le calcul pédagogique

Notre contrôle a utilisé CPython 3.12.10 et uniquement la bibliothèque standard. Les racines du cas à quatre points [1, −i, −1, i] sont représentées exactement : l’égalité exacte convient donc à ces petites entrées entières. Ce critère ne doit pas être transposé mécaniquement à des transformées plus longues en virgule flottante.

Les calculs élémentaires suivants ont été exécutés dans notre contrôle enregistré :

PYTHON
ROOTS = (1, -1j, -1, 1j)

def direct(x, roots=ROOTS):
    return [
        sum(x[n] * roots[(n*k) % 4] for n in range(4))
        for k in range(4)
    ]

def butterfly(x):
    even = (x[0] + x[2], x[0] - x[2])
    odd = (x[1] + x[3], x[1] - x[3])
    return [
        even[0] + odd[0], even[1] - 1j*odd[1],
        even[0] - odd[0], even[1] + 1j*odd[1],
    ]

Comparez les deux fonctions sur [1, 2, 3, 4], une impulsion [1, 0, 0, 0], une constante [1, 1, 1, 1], des zéros et [1+i, −2, 3−i, 2i]. Les cinq résultats concordaient exactement. Le premier correspondait aussi au tableau ci-dessus. Inverser les signes imaginaires des racines a modifié le spectre du cas asymétrique, et le contrôle l’a rejeté comme une convention différente.

Cela valide les cas pédagogiques fixés. Cela ne prouve pas que toute mise en œuvre est correcte et ne réplique pas de manière indépendante les performances historiques. L’article original ne fournit ni le programme complet, ni les entrées du banc de mesure, ni l’ensemble des conditions de chronométrage nécessaires à une telle reproduction.

Article rédigé et traduit par un système d’IA à partir des sources citées, avec des contrôles automatisés selon la méthode éditoriale News. Le calcul pédagogique a été exécuté comme décrit ; il ne constitue ni une évaluation par les pairs, ni une validation scientifique humaine, ni une reproduction du programme historique.