TangenteX.com

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(n2). 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 (nlog2 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 en FORTRAN

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
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
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)
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
ENDDO
RETURN
END

Quelques commentaires:

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

Exemple d'application de la FFT en FORTRAN

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\pi f + \phi) \).
Il s'agit d'abord de constituer un vecteur réel y de n points, qui stocke sous forme discrète chaque valeur yi 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 yi dans la partie réelle yci et en déclarant nulle la partie imaginaire de yci.
Finalement, je calcule le spectre de fréquences du signal en remarquant que la valeur du spectre au point i est la norme du yci.
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
Spectre(i)
WRITE(10,*) t(i),y(i),f(i),REAL(yc(i)),AIMAG(yc(i)), * ENDDO
ENDDO
CLOSE(10)
ELSE
WRITE(*,*)'Le critere de Nyquist n''est pas respecte'
ENDIF
STOP
END

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:

Tracé de la transformée de Fourier:

TraceFFT

Tracé du spectre de puissance non normalisé:

Trace spectre

Contenu et design par Dominique Lefebvre - www.tangenteX.com avril 2013 Licence Creative Commons Contact :