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