La transformée de Fourier rapide

L'algorithme de FFT

Je n'ai pas l'intention, dans cette page, d'exposer les fondements de la transformée de Fourier. Je vous renvoie à votre cours d'analyse ou de traitement du signal favori. Si besoin, on trouve sur le web des exposés plus ou moins détaillés.
Je vous recommande également la lecture du bouquin de A. Garcia (Numerical Methods for Physics) qui présente une introduction pratique à la transformée de Fourier discrète (TFD). Ce bouquin a été ma source d'inspiration pour l'algorithme de FFT proposé ci-dessous et pour le programme exemple.
La TFD est un outil puissant de traitement du signal, mais qui présente l'inconvénient majeur de nécessiter une quantité de calculs astronomique. Plus précisément, le nombre d'opérations pour calculer une TFD sur n points est de l'ordre de O(n²). Imaginez ce que cela peut donner lorsqu'on calcule, en temps réel de préférence, une TFD sur 1024 points. Bref, pour les ordinateurs du milieu du XXeme siècle, c'était un peu juste....
En 1965, John Tukey et James Cooley ont inventé (à partir de travaux précédents datant des années 40) un algorithme particulier dit "FFT" pour Fast Fourier Transformation. Cet algorithme doit également beaucoup à un lemme très astucieux de Danielson et Lanczos, lemme qui impose que le nombre de points traités soit une puissance de 2. Pour l'anecdote, John Tukey, mort en 2000, est aussi l'inventeur du mot "software".
L'algorithme FFT réduit le nombre d'opérations nécessaires pour calculer une TFD à un ordre O (nlog₂ n)), ce qui est un progrès considérable. On peut maintenant calculer une TFD sur 512 ou 1024 points, en temps réel, sur un PC ou même un microcontrôleur standard (plutôt un DSP, mais bon...). C'est cet algorithme que je vous propose d'implémenter sur votre PC pour nos petites manipulations de physique numérique.

Attention: l'usage de la TFD et de l'algorithme FFT implique le respect d'un certains nombres de règles particulières, sans quoi les résultats obtenus sont... curieux! J'en aborderai quelques un plus bas et d'autres dans les expériences numériques que je vous proposerai.

Implémentation de l'algorithme FFT

On trouve dans la littérature un grand nombre d'implémentations de l'algorithme de FFT. La plus répandue étant sans doute celle des "Numerical Recipes" en C ou en FORTRAN. A.Garcia dans son ouvrage déjà cité, en propose une en C++. Et signalons que les langages Matlab et Scilab proposent une fonction fft "built-in".
Ces algorithmes sont basés sur la même technique, qui repose sur la décomposition du calcul de la TFD sur n points en calcul d'un certain nombre de TFD sur deux points puis en une combinaison de ces TFD élémentaires selon une méthode appelée poétiquement "papillon" ou "butterfly" selon la langue, après réarrangement particulier de l'ordre des termes...
Il existe d'autres méthodes qui font appel à la récursivité, mais que je préfère éviter: je me méfie comme d'une peste de la récursivité, qui rend très difficile le debugage des programmes.
L'algorithme que je vous propose ci-dessous est un mixte de l'algorithme de A.Garcia et de celui des Numerical Recipes. Je vous le propose en FORTRAN, car je trouve que l'écriture est plus élégante qu'en C, en grande partie parce que le FORTRAN possède un type COMPLEX natif. Ce qui donne:

C Calcul d'une transformee de Fourier rapide FFT sur un

C tableau de complexes. Inspire de A.Garcia (Numerical

C Methods for Physics).

SUBROUTINE FFT(A, N)

IMPLICIT NONE

C Description de l'interface de la subroutine

INTEGER N

! Nombre de points du vecteur complexe d'entree (doit etre une puissance de 2)

COMPLEX A(N)

! Vecteur complexe d'entree

C En entree le tableau de complexes donnee contient les donnees

C En sortie le tableau de complexes donnee contient les transformees

C Declaration des constantes et parametres

REAL PI

PARAMETER (PI = 3.141592654)

C Declaration des variables

INTEGER half_N, Nm1,i,j,k, m,ki,ki1,ip

COMPLEX T, U, W

REAL angle

C Initisalisation des variables

half_N = N/2

Nm1 = N-1

m = INT(LOG(REAL(N))/LOG(2.0) + 0.5)

! N = 2**m

C Bit-scramble du vecteur d'entree

j = 1

DO i=1,Nm1

IF (i.LT.j) THEN ! on swappe A(i) et A(j)

T = A(j)

A(j) = A(i)

A(i) = T

ENDIF

k = half_N

DO WHILE (k.LT.j)

j = j-k

k = k/2

ENDDO

j = j+k

ENDDO

C Calcul de la transformee (Butterfly operation)

DO k=1,m

ki = 2**k

ki1 = ki/2

U = (1.0, 0.0)

angle = PI/REAL(ki1)

W = CMPLX(COS(angle),-SIN(angle))

DO j=1,ki1

DO i=j,N, ki

ip = i+ki1

T = A(ip)*U

A(ip) = A(i)-T

A(i) = A(i)+T

ENDDO

U = U*W

ENDDO

RETURN

END

Quelques commentaires:

L'usage de variables complexes simplifie largement la programmation. Comparez la même implémentation en C (dans les Numerical Recipes in C par exemple)...
Le tableau A en entrée de la routine est de type COMPLEX.
Il est indispensable que le nombre de points N soit une puissance de 2. Cela est imposé par l'algorithme de FFT.

Le fichier source de la routine FFT est disponible en FORTRAN 77 dans FFT.for

Exemple d'application de la FFT

Pour illustrer l'usage de la routine FFT, je vous propose un programme de test ultra classique : le calcul du spectre d'un signal périodique. L'analyse spectrale est en effet un des domaines principaux de l'usage de la FFT. Le signal choisi est trivial : sin(2*π*f + φ).
Il s'agit d'abord de constituer un vecteur réel y de n points, qui stocke sous forme discrète chaque valeur y_i du signal, i variant de 1 à n. n est bien sur une puissance de 2 (dans mon exemple 1024).
Pour calculer la FFT, on charge ce vecteur y dans un vecteur complexe yc, en chargeant y_i dans la partie réelle yc_i et en déclarant nulle la partie imaginaire de yc_i.
Finalement, je calcule le spectre de fréquences du signal en remarquant que la valeur du spectre au point i est la norme du yc_i.
Voici le code du programme TESTFFT.for:

C Programme de test de FFT

C Dominique Lefebvre avril 2006

PROGRAM TESTFFT

IMPLICIT NONE

EXTERNAL FFT

C Declaration des constantes

INTEGER n

REAL Fech, DEUXPI

PARAMETER (n = 1024,DEUXPI=6.283185307)

C Declaration des variables

INTEGER i

REAL frequence, phi

REAL t(n), y(n), f(n), Spectre(n)

COMPLEX yc(n)

C Saisie de la valeur de la frequence

WRITE(*,'(a,$)')'Frequence du signal (Hz): '

READ(*,*)frequence

C Saisie de la valeur d'echantillonnage

WRITE(*,'(a,$)')'Le critere de Nyquist impose que la frequence'

WRITE(*,*)'d''echantillonnage soit > 2*frequence du signal'

WRITE(*,*)''

WRITE(*,'(a,$)')'Frequence d''echantillonnage (Hz): '

READ(*,*)Fech

C Saisie de la valeur de la phase

WRITE(*,'(a,$)')'Phase (rd): '

READ(*,*)phi

C Verification du critere de Nyquist

IF (Fech .GT. 2*frequence)THEN

C Calcul du vecteur signal

DO i=1,n

t(i) = (i-1)/Fech

y(i) = SIN(DEUXPI*t(i)*frequence + phi)

f(i) = Fech*(i-1)/n

ENDDO

C Calcul de la transformee de Fourier du signal

DO i=1,n

yc(i) = CMPLX(y(i), 0.0)

ENDDO

CALL FFT(yc, n)

C Calcul du spectre du signal

DO i=1,n

Spectre(i) = REAL(yc(i))**2 + AIMAG(yc(i))**2

ENDDO

C Stockage des donnees dans un fichier

OPEN(10,FILE = 'Spectre.dat')

DO i=1,n

WRITE(10,*) t(i),y(i),f(i),REAL(yc(i)),AIMAG(yc(i)), * Spectre(i)

ENDDO

CLOSE(10)

ELSE

WRITE(*,*)'Le critere de Nyquist n''est pas respecte'

ENDIF

STOP

END

Quelques commentaires:

Sur le plan de la programmation pure, ce programme n'appelle pas de commentaire particulier.
Sa structure est simple:
- lecture des paramètres du signal d'entrée (fréquence du signal, phase) et de la fréquence d'échantillonnage.
- constitution du vecteur de signal d'entrée.
- calcul de la FFT après avoir chargé le signal réel dans un vecteur complexe.
- calcul du spectre.
- enregistrement des résultats, pour tracé.
Le spectre calculé n'est pas normalisé. Dans nos manips futures, nous normaliserons le spectre, mais dans ce programme exemple, cela n'a pas beaucoup d'importance.
Vous avez sans doute noté le contrôle fait sur la fréquence d'échantillonnage. Vous vous souvenez bien sur du théorème de Shannon, duquel Nyquist a tiré un critère d'échantillonnage : "Un signal analogique x(t) ayant des fréquences plus petites que Fmax peut être reconstruit exactement à partir de ses échantillons x[n]=x(nTs) si les échantillons sont produits à un taux Fs=1/Ts qui est plus grand que 2Fmax". On contrôle donc que la fréquence d'échantillonnage est bien au moins deux fois plus grande que la fréquence du signal."

Pour faire fonctionner ce programme, il faut créer un projet EssaiFFT avec le programme VFORT (voir Introduction au FORTRAN), puis créer les 2 fichiers sources FFT et TestFFT. Après exécution, on peut visualiser les résultats en utilisant les fichiers de commandes GnuPlot TraceFFT.p et TraceSpectre.p
Voici les graphes obtenus avec les paramètres suivants:

fréquence du signal : 0.2 Hz
phase : 0
fréquence d'échantillonnage : 1 Hz
nombre de points de la FFT : 1024

Tracé de la transformée de Fourier:

Tracé du spectre de puissance non normalisé: