La transformée de Fourier rapide

La FFT, ou transformée de Fourier rapide, est un des quelques algorithmes dont la publication a provoqué une véritable révolution dans le champ technique. Généralement associé aux noms de J.W. Cooley et J.W. Tuckey qui l'ont publié en 1965, cet algorithme de calcul de la transformée de Fourier discrète avait été maintes fois << redécouvert >> depuis Gauss, notamment par Danielson et Lanczos en 1942. La FFT permet de ramener le calcul de la transformée de Fourier discrète de $N^2$ à $N\log N$ opérations ; cette réduction de complexité suffit à faire passer d'impossibles à facilement résolubles nombre de problèmes.

La transformée de Fourier discrète d'un $n$-uplet de nombres complexes $a = (a_0, \ldots, a_{n-1})$ est le $n$-uplet $\hat{a} = (\hat{a}_0, \ldots, \hat{a}_{n-1})$ défini par

$$\hat{a}_q = \sum_{p=0}^{n-1} a_p \omega_n^{pq}$$

où $\omega_n = e^{2\pi i/n}$. Notons que $\hat{a}_{q+n} = \hat{a}_q$ et que pour $n=1$, ceci se réduit en $\hat{a} = a$. La transformée inverse se calcule ainsi :

$$a_p = \frac{1}{n} \sum_{q=0}^{n-1} \hat{a}_q \omega_n^{-pq}$$

Supposons que $n=2m$ et séparons les éléments de $a$ d'indices pairs de ceux d'indices impairs :

$$\begin{eqnarray*} a^{[0]} &= &(a_0, a_2, \ldots, a_{n-2}) \\ a^{[1]} &= &(a_1, a_3, \ldots, a_{n-1}) \end{eqnarray*}$$

Alors, en séparant la somme en deux termes regroupant d'une part les $p=2j$ et d'autre part les $p=2j+1$, et en utilisant l'égalité $\omega_{2m}^{2jq} = \omega_{m}^{jq}$ :

$$\begin{eqnarray*} \hat{a}_q &= &\sum_{p=0}^{n-1} a_p \omega_n^{pq} \qquad (0\le... ...a_{m}^{jq} \\ &= &\hat{a}^{[0]}_q + \omega_n^q \hat{a}^{[1]}_q \end{eqnarray*}$$

Ainsi, pour calculer la transformée de Fourier sur $n$ points, il suffit de calculer deux transformées de Fourier sur $n/2$ points, de faire $n$ multiplications et $n$ additions. Si $n$ est une puissance de 2, cette décomposition peut s'appliquer récursivement, ce qui donne pour équation de complexité : $T(n) = 2T(n/2) + \Theta(n)$. Il en résulte une complexité en $\Theta(n\log n)$ opérations.

La définition récursive de la FFT résulte directement de cette décomposition ; elle utilise le type double_complex qui est déclaré dans <complex> :

#include <complex>

void fft_rec(int n,
             double_complex a[],
             double_complex omega,
             double_complex fft_a[])
{
  if (n == 1) {
    fft_a[0] = a[0];
  } else {
    int m = n/2;
    double_complex
      *b = new double_complex[m],
      *c = new double_complex[m],
      *fft_b = new double_complex[m],
      *fft_c = new double_complex[m],
      omega2 = omega*omega,
      alpha = double_complex(1,0);
    int i;
    
    for (i=0; i<m; i++) {
      b[i] = a[2*i];
      c[i] = a[2*i+1];
    }
    fft_rec(m, b, omega2, fft_b);
    fft_rec(m, c, omega2, fft_c);
    for (i=0; i<m; i++) {
      fft_a[i] = fft_b[i] + alpha*fft_c[i];
      alpha *= omega;
    }
    for (i=0; i<m; i++) {
      fft_a[m+i] = fft_b[i]+alpha*fft_c[i];
          alpha *= omega;
    }
    delete[] b;
    delete[] c;
    delete[] fft_b;
    delete[] fft_c;
  }
}

On remarquera que les arguments des appels récursifs sont des adresses de tableaux dynamiques alloués sur le tas, et que ces tableaux dynamiques sont explicitement dés-alloués.

La figure 33 représente l'arbre des appels de ce calcul récursif pour $n=8$ points : aux feuilles, la FFT sur 1 point, qui est l'identité. L'ordre dans lequel les $a_p$ sont rangés dans les feuilles, de gauche à droite, est << l'inverse binaire >> de $p$. La fonction bit_rev(h,p) renverse la représentation binaire de $n$, pour $0\le p\le 2^h-1$ ; par exemple, $\mathtt{bit\_rev}(3,3) =\mathtt{bit\_rev}(3, \mathtt{011}_2) = \mathtt{110}_2 = 6$ ; on utilise les opérateurs binaires <<, >>, &, | et ^, plus efficaces pour opérer au niveau des bits d'un entier :

int bit_rev(int h, int p) {
  int y = 0;
  int n = 1 << h;        // n = 2 puissance h       
  int r;

  while (h>0) {
    h--;
    n >>= 1;             // n = n/2 
    r = p&1;             // r = p%2 
    p >>= 1;             // p = p/2 
    if (r == 1) {
      y |= n;            // y = y+n 
    }
  }
  return y;
}

On va donc commencer par permuter le $n$-uplet $a$ pour ranger ses éléments dans cet ordre :

void bit_rev_array(int h, const double_complex a[], double_complex b[]) {
  int p;

  for (p=0; p < 1<<h; p++) {
    b[bit_rev(h,p)] = a[p];
  }
}

Cela va permettre de construire un algorithme itératif, qui parcourt l'arbre des feuilles vers la racine.

En utilisant la périodicité de la transformée de Fourier, et l'égalité $\omega_n^m = -1$, on peut réécrire la dernière expression de $\hat{a}$ :

$$\begin{array}{rcl} \hat{a}_q &= &\hat{a}^{[0]}_q + \omega_n... ...- \omega_n^q \hat{a}^{[1]}_q \end{array}\qquad (0\le q\le m-1)$$

Cette transformation, de la forme $(x,y) \mapsto (x+\alpha y, x-\alpha y)$, est appelée un papillon (en anglais, butterfly) ; on l'appliquera successivement avec $\alpha = 1$, $\omega$, $\omega^2$, ..., $\omega^{m-1}$. La fonction butterfly s'écrit :

void butterfly(double_complex a[], int j, int k, double_complex alpha) {
  double_complex u = a[j];
  double_complex v = alpha * a[k];

  a[j] = u+v;
  a[k] = u-v;
}

void fft(int h, double_complex a[], double_complex fft_a[]) {
  const int n = 1 << h;       // n = 2 puissance h 
  int l, i, j;

  bit_rev_array(h, a, fft_a);
  for (l=1; l<=h; l++) {
    int m = 1<<l;
    double_complex alpha = 
      double_complex(1,0);
    double_complex omega_m = 
      double_complex(cos(2*M_PI/m), sin(2*M_PI/m));
    
    for (i=0; i<m/2; i++) {
      for (j=i; j<n; j = j+m) {
        butterfly(fft_a, j, j+m/2, alpha);
      }
      alpha *= omega_m;
    }
  }
}

Il existe aussi une FFT en arithmétique modulo $N$, où $N = 2^{tn/2} +1$, $t$ est un entier quelconque ; on pose $\omega = 2^t$, qui est une racine $n$-ème principale de l'unité dans $\mathbf{Z}/N\mathbf{Z}$ ; on a toujours $n=2^h$. Par exemple, quand $t=2$ et $n=2^3=8$, on a $N=257$, et $\omega=2^2=4$ ; un $n$-uplet et sa transformée sont :

$$\begin{eqnarray*} a &= &(255, 8, 0, 226, 37, 240, 3, 0) \\ \hat{a} &= &(255, 85, 127, 200, 78, 255, 194, 75) \end{eqnarray*}$$