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 à 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 -uplet de nombres complexes est le -uplet défini par
où . Notons que et que pour , ceci se réduit en . La transformée inverse se calcule ainsi : Supposons que et séparons les éléments de d'indices pairs de ceux d'indices impairs :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.
La figure représente l'arbre des appels de ce calcul
récursif pour points : aux feuilles, la FFT sur 1 point, qui est
l'identité. L'ordre dans lequel les sont rangés dans les
feuilles, de gauche à droite, est << l'inverse binaire >> de . La
fonction bit_rev(h,p) renverse la représentation binaire de
, pour ; par exemple, ; 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 -uplet 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é , on peut réécrire la dernière expression de :
Cette transformation, de la forme , est appelée un papillon (en anglais, butterfly ) ; on l'appliquera successivement avec , , ,..., . 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 , où , est un entier quelconque ; on pose , qui est
une racine -ème principale de l'unité dans ;
on a toujours . Par exemple, quand et , on a
, et ; un -uplet et sa transformée sont :