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² à $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 :

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}$ :

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 :

void fft_rec(int n, complex a[], complex omega, complex fft_a[]) {
  if (n == 1) {
    fft_a[0] = a[0];
  } else {
    int m = n/2;
    complex
      *b = malloc(m * sizeof(complex)),        /* complex b[m]; */
      *c = malloc(m * sizeof(complex)),
      *fft_b = malloc(m * sizeof(complex)),
      *fft_c = malloc(m * sizeof(complex)),
      omega2 = complexe_mult(omega,omega),
      alpha = {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] = complexe_add(fft_b[i], 
                              complexe_mult(alpha, fft_c[i]));
      alpha = complexe_mult(alpha, omega);
    }
    for (i=0; i<m; i++) {
      fft_a[m+i] = complexe_add(fft_b[i], 
                                complexe_mult(alpha, fft_c[i]));
      alpha = complexe_mult(alpha, omega);
    }
    free(b);
    free(c);
    free(fft_b);
    free(fft_c);
  }
}

On remarquera qu'en l'absence de tableaux dynamiques, les arguments des appels récursifs sont des adresses de blocs alloués dynamiquement, et que ces blocs sont explicitement dés-alloués.

  

Figure 35: Arbre des appels de la FFT sur 8 points

La figure 35 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 complex a[], 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_... ...q - \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(complex a[], int j, int k, complex alpha) {
  complex u = a[j];
  complex v = complexe_mult(alpha, a[k]);

  a[j] = complexe_add(u,v);
  a[k] = complexe_sub(u,v);
}

void fft(int h, const complex a[], 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;
    complex alpha = { 1, 0 };
    complex omega_m;

    omega_m.re = cos(2*M_PI/m);
    omega_m.im = 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 = complexe_mult(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³=8, on a N=257, et $\omega=2^2=4$ ; un n-uplet et sa transformée sont :