La transformée de Fourier rapide

Une transformée de Fourier discrète de N échantillons correspond à N sommes de N termes. Divisez cette somme en sommes plus petites et le même résultat coûte beaucoup moins cher. Observez cette décomposition dans le diagramme en papillon, ou telle que Gauss l'a effectuée, sur douze valeurs concernant l'astéroïde Pallas.

Exemple résolu — La transformée de Fourier rapide

Combien de fois moins	5,3

Ce sont les valeurs affichées par le calculateur à l’ouverture de la page. Modifiez une entrée ci-dessus pour les recalculer.

Chargement de la simulation interactive...

Leçon

La théorie — La transformée de Fourier rapide

La transformée de Fourier rapide (FFT) désigne tout algorithme qui calcule la transformée de Fourier discrète X[k] = Σ x[n] W^kn, W = e^−2πi/N, de manière exacte aux arrondis près, en environ N log N opérations au lieu de N². La forme de Cooley-Tukey décompose N = N₁N₂ : N₂ transformées entrelacées de longueur N₁, dont chaque sortie k₁ de la transformée n₂ est multipliée par son propre facteur de rotation W^n₂k₁, suivies de N₁ transformées de longueur N₂ à travers elles. Appliquée de façon répétée à N = 2^m, elle devient la FFT en base 2 : m étapes de N/2 papillons.

Ce que signifie chaque symbole

N: le nombre d'échantillons, ainsi que de fréquences calculées
W: la racine de l'unité e^−2πi/N ; ses puissances sont les facteurs de rotation
m: log₂N, le nombre d'étapes d'une FFT en base 2

Gauss a noté la décomposition de la FFT environ 160 ans avant Cooley et Tukey 🖖

Cliquez sur Gauss, 3 périodes de 4. Ces douze nombres proviennent du manuscrit de Gauss : la déclinaison de la limite nord de la bande de ciel dans laquelle l'astéroïde Pallas peut errer, tous les 30° d'ascension droite, en minutes d'arc, copiée d'après le journal de von Zach. Il cherchait la série trigonométrique qui les traverse. Au lieu d'une unique somme à douze termes, il a pris une valeur sur trois (0°, 90°, 180°, 270° ; puis 30°, 120°, … ; puis 60°, …), a ajusté chaque ensemble de quatre indépendamment, puis a combiné les trois ajustements courts. Les courbes en pointillés représentent ces ajustements. Elles divergent : quatre points ne peuvent pas percevoir les oscillations plus rapides. Chaque ensemble obtient donc un résultat différent, ce qui l'a poussé à chercher la combinaison. Le tableau sous les résultats reprend son calcul : la première période donne 779,5, −415,0, −726,0 et 43,5, soit les nombres de son manuscrit. La série combinée correspond bien à la version directe à douze termes. Il a procédé d'une deuxième manière, avec quatre périodes de trois, et a obtenu la même série. L'arithmétique des indices correspond à celle que Cooley et Tukey ont publiée en 1965, comme l'ont montré Heideman, Johnson et Burrus en 1984. Ce que le manuscrit omet compte tout autant. Il ne calcule jamais le travail économisé. Là où Gauss affirme que la décomposition est plus rapide, il en laisse la preuve à la pratique. Praxis tentantem docebit, celui qui essaie le verra. Ses exemples développés comportent deux facteurs, 12 = 3 × 4 et 36 = 6 × 6. L'article 27 se contente d'indiquer que la décomposition s'étend à davantage de facteurs. Le texte n'a pas été publié de son vivant. Il est paru en latin dans ses œuvres complètes en 1866 et ne porte aucune date. L'éditeur des Werke le qualifie de remaniement d'un brouillon apparemment commencé en octobre 1805. Heideman, Johnson et Burrus le situent en octobre ou novembre 1805 à partir de plusieurs indices. Parmi eux figure une excentricité de Junon, 0,254236, que le manuscrit utilise et qui a été imprimée en mai 1805.

L'économie n'est pas un facteur fixe : elle croît avec N 🖖

Cliquez sur 1 024 points. La somme directe multiplie chaque échantillon par chaque facteur de rotation : N² = 1 048 576 multiplications. La FFT exécute log₂N = 10 étapes de N/2 = 512 papillons avec une multiplication chacun : 5 120. C'est 204,8 fois moins, et le rapport vaut 2N/log₂N. Il continue donc de croître. Cliquez sur Un million de points : 1 099 511 627 776 multiplications contre 10 485 760, un rapport de 104 857,6. Doublez N et le travail de la somme directe est multiplié par quatre. Celui de la FFT augmente de 2,2 fois pour un millier de points, et tend vers le double à mesure que vous progressez. À 8 points, l'économie n'est que de 5,3 fois. C'est pourquoi le diagramme ressemble à un câblage complexe pour un gain minime. Les deux décomptes incluent les multiplications par W⁰ = 1, tout comme le fait le N² de Cooley et Tukey. Ils ont compté les paires de multiplication et d'addition pour obtenir un nombre inférieur à 2N log₂N. Dans les deux cas, le rapport croît comme N/log N. Aucun ordinateur ne modifie ces décomptes ni leur rapport. Une machine plus rapide modifie les temps d'exécution, et pas nécessairement du même facteur pour les deux. Le bouton situé sous les tableaux chronomètre les deux méthodes dans votre navigateur. Le rapport mesuré augmente également avec N, bien qu'il n'atteigne pas ces valeurs. Le trafic mémoire et les additions coûtent du temps, or le décompte les ignore.

Pour de grands N, la méthode rapide est aussi la plus exacte 🖖

Le tableau d'exactitude exécute les deux méthodes en arrondissant chaque opération à la simple précision, soit environ sept chiffres significatifs, face à un résultat en double précision. À 1 024 points, la somme directe s'écarte de 5,75×10⁻⁷ de la taille de la réponse, contre 1,23×10⁻⁷ pour la FFT. C'est environ 4,7 fois moins. Passer de 16 points à 1 024, soit 64 fois plus, a multiplié l'erreur de la somme directe par 7,3 et celle de la FFT par 1,4. Chaque sortie de la somme directe est le total cumulé de N produits arrondis. Ses erreurs d'arrondi s'accumulent donc sur 1 024 additions. Chaque sortie de la FFT ne traverse que log₂N = 10 papillons. À 16 points l'écart est négligeable, et sur cette entrée la FFT se révèle même légèrement moins bonne, avec 8,74×10⁻⁸ contre 7,91×10⁻⁸. Les deux méthodes obtiennent chaque facteur de rotation de la même façon : calculé en double précision puis arrondi une seule fois. Le tableau compare ainsi les deux algorithmes et non deux manières de calculer un sinus.

Problème entièrement résolu

L'autre exemple de Gauss : Junon, 36 valeurs 4 étapes
L'article 41 du manuscrit de Gauss interpole l'équation du centre de Junon à partir de 36 valeurs, une tous les 10°, décomposées en 6 × 6. Selon le décompte de Cooley et Tukey, combien d'opérations la somme directe exige-t-elle, combien la décomposition de Gauss en économise-t-elle, et une décomposition plus fine ferait-elle mieux ? Cliquez d'abord sur Gauss, 3 périodes de 4 : la page effectue le même décompte pour douze valeurs.

Ouvrez ce problème dans l'outil
1. La somme directe, une opération par échantillon et par fréquence :
  
  $N^{2} = 3 6^{2} = 1, 296$
2. La décomposition de Gauss, six périodes de six : le décompte en deux étapes de Cooley et Tukey équivaut à N fois la somme des facteurs.
  
  $T = N (N_{1} + N_{2}) = 36 \times (6 + 6) = 432$
3. La décomposition poussée jusqu'au bout, 36 = 2 × 2 × 3 × 3 :
  
  $T = 36 \times (2 + 2 + 3 + 3) = 360$
4. La règle de Cooley et Tukey : un facteur r = st avec s, t > 1 présente toujours s + t < r, sauf si s = t = 2 :
  
  $s + t < s t (s, t > 1, (s, t) \neq = (2, 2)); 6 = 2 \times 3 : 2 + 3 = 5 < 6$
Réponse

1 296 opérations pour la somme directe, 432 avec la décomposition 6 × 6 de Gauss, 360 en décomposant jusqu'au bout. La décomposition de Gauss accomplit déjà trois fois moins de travail, et la décomposition la plus fine économise encore un sixième du reste. Les douze valeurs de la page suivent la même logique : 144 en somme directe et 84 avec la décomposition, ou 12 × (2 + 2 + 3) = 84 à nouveau jusqu'au bout, car diviser un facteur de 4 en 2 × 2 coûte autant.

$1, 296 \to 432 \to 360; 1 2^{2} = 144 \to 12 \times (4 + 3) = 12 \times (2 + 2 + 3) = 84$

Parcours

La transformée de Fourier en pratique

Références (3)

Gauss's manuscript on trigonometric interpolation, printed after his death. Article 28 is the page's Gauss mode: twelve declinations of the northern limit of Pallas's zodiac at every 30° of right ascension, split into three periods of four and into four periods of three. Article 27 says the split extends to three or more factors and that practice will teach whoever tries it ("praxis tentantem docebit"); article 41 works Juno's equation of the centre from 36 values with the eccentricity 0.254236. The editor's note calls the text a reworking of a draft apparently begun in October 1805: C. F. Gauss, Theoria interpolationis methodo nova tractata (Nachlass). Werke, Band III, Königliche Gesellschaft der Wissenschaften, Göttingen, 1866; second printing 1876, where the treatise begins on p. 265 and arts. 25–28 are on pp. 303–310.
The paper that showed Gauss's split is the Cooley–Tukey FFT: the change of index n = N₂n₁ + n₂, k = k₁ + N₁k₂ written as Gauss's procedure, his examples (N = 12 as 4 × 3 and 3 × 4, N = 36 as 6 × 6), the dating to October–November 1805 from the Juno eccentricity and other clues, and what the manuscript does not contain: "Gauss did not, however, go on to quantify the computational requirements of his method": M. T. Heideman, D. H. Johnson and C. S. Burrus, Gauss and the history of the fast Fourier transform. IEEE ASSP Magazine 1(4), 14–21 (October 1984). doi:10.1109/MASSP.1984.1162257
The paper that made the FFT standard: the split N = r₁r₂ with T = N(r₁ + r₂) operations, its repetition to T = N(r₁ + … + r_m), the radix-2 form with bit-reversed order, and the count of N² operations for the direct sum against fewer than 2N log₂N, where an operation is a complex multiplication followed by a complex addition: J. W. Cooley and J. W. Tukey, An algorithm for the machine calculation of complex Fourier series. Mathematics of Computation 19(90), 297–301 (1965). doi:10.1090/S0025-5718-1965-0178586-1

Par l’équipe 42Math · Vérifié sur des sources publiées. Corrections bienvenues.