TangenteX.com

L'algorithme de FFT

Pourquoi un algorithme rapide pour le calcul de la TFD

Dans la page de TangenteX consacrée à une introduction de la transformée de Fourier discrète (TFD), j'ai eu l'occasion d'aborder les aspects mathématiques de la TFD.
La TFD est un outil puissant d'analyse harmonique, mais qui présente l'inconvénient majeur de nécessiter une quantité de calculs astronomique. Plus précisément, la complexité algorithmique du calcul d'une TFD sur N points est de N2. Imaginez ce que cela peut donner lorsqu'on calcule en temps réel, une TFD sur 1024 points. Pour les ordinateurs du milieu du XXeme siècle, c'était un peu juste....
En 1965, James Cooley et John Tukey ont inventé, à partir de travaux précédents de Danielson et Lanczos datant de 1942, un algorithme particulier dit "FFT" pour Fast Fourier Transformation. Cet algorithme doit beaucoup à un lemme très astucieux de Danielson et Lanczos. Pour l'anecdote, John Tukey, mort en 2000, est aussi l'inventeur du mot "software".
L'algorithme FFT réduit la complexité algorithmique du calcul d'une TFD à Nlog2N, ce qui est un progrès considérable. On peut maintenant calculer une TFD sur 512 ou 1024 points sur un PC, un DSP (Digital Signal Processor)  ou même un microcontrôleur standard. C'est cet algorithme que je vous propose d'implémenter sous différentes formes pour nos besoins.

Comment fonctionne l'algorithme de FFT

L'algorithme est basé 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.

La transformée de Fourier discrète est définie par, en laissant tomber les facteurs constants, en particulier le \( \dfrac{1}{\sqrt{2\pi}} \) :

\( \displaystyle \widehat{c_k} =   \sum_{n=0}^{N-1} y_n e^{-i\dfrac{2\pi kn}{N}} \)

que je vais écrire différemment en posant \( W = e^{-i\dfrac{2\pi n}{N}} \) :

\( \displaystyle \widehat{c_k} =   \sum_{n=0}^{N-1} y_n W^{nk} \)

En 1942, Danielson et Lanczos montrent dans un lemme, qu'il est possible d'écrire la définition de la TFD de longueur N sous la forme de la somme de deux TFD de longueur N/2, l'une étant constituée par les éléments pairs de la TFD originale et l'autre par les éléments impairs, soit pour reprendre leur formulation originale :

\( \displaystyle \widehat{F_k} = F_k^{even} + W^k F_k^{odd} \)

\( F_k \) étant notre \( \widehat{c_k} \), et \( W^k \) le terme défini ci-dessus.

Cette décomposition fonctionne parce que la fonction exponentielle est très particulière, qui permet d'observer l'égalité suivante sur le terme \(W^k \) :

\( W_N^{-2kn} = W_{N/2}^{-kn} \).

Appliquons maintenant la même décomposition en deux sommes de longueur N/4 sur les sommes \( \displaystyle F_k^{even} \) et \( \displaystyle F_k^{odd} \). Et ainsi de suite jusqu'à ce que cela ne soit plus possible, c'est à dire que \( \dfrac{N}{2^p} = 1 \). Vous voyez maintenant pourquoi il est vivement conseillé d'utiliser un nombre de points N qui soit une puissance de 2.

Dernière chose : les résultats obtenus à l'issue de ces calculs ne sont pas naturellement ordonnés. On dit qu'ils présentent un "entrelacement en temps" ou effet "butterfly". Il faut donc les réordonner. C'est ici qu'intervient la procédure de renversement des bits, qui permet de remettre les termes dans le bon ordre. La page wikipedia consacrée à l'algorithme de Cooley et Tukey présente assez complétement la chose.

Différentes implémentations de la 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. Cependant, la plupart des langages de programmation scientifiques modernes comme Python, MatLab, Maple ou Scilab proposent une fonction fft "built-in".

Tous les programmes et les scripts cités dans cette page sont téléchargeables dans la bibliothèque de codes de TangenteX.

Pour présenter ces différentes implémentations en Python et FORTRAN, je vais travailler avec le même signal qui m'a servi pour étudier la TFD, en changeant toutefois le nombre de points, que je fixerai ici à 32, une puissance de 2. Je vous rappelle ici ses caractéristiques :

# paramètres du signal
L = 2.0             # durée du signal
fd = 1.0/L          # domaine fréquentiel du signal
a0 = 1.             # amplitude du signal
f0 = 1.0            # fréquence du signal
# paramètres de l'échantillonnage
N = 32              # nombre de points d'échantillonnage
dt = L/N            # intervalle temporel entre deux points
fe = 1./dt          # fréquence d'échantillonnage
# construction du signal
t = linspace(0.0,L,N)
s = a0*sin(2*pi*f0*t)

Vous pouvez modifier ces données, en prenant toutefois garde à respecter la condition de Shannon-Nyquist dans le choix de la fréquence d'échantillonnage par rapport à la fréquence du signal. Attention aussi au choix de la durée de fenêtrage du signal L, qui doit être très grand par rapport à la période du signal.

Dans toutes les implémentations ci-dessous, j'utilise une fenêtre d'acquisition rectangulaire. J'aborderai la notion importante de fenêtrage dans une autre page. Le fenêtrage rectangulaire est intuitif. Lorsque j'effectue une acquisition sur un signal f(t), l'acquisition débute à un instant t0 et s'achève à l'instant t1. Je découpe une tranche du signal réel, d'une durée L = t1 - t0, et je considère que mon signal est L-périodique, c'est à dire que le bout de signal se répète de la même manière avant t0 et après t1. Pratiquement, on peut imaginer que l'on ferme un interrupteur à t0, en laissant passer la partie du signal que je mesure, puis j'ouvre l'interrupteur en t1. Mathématiquement, on peut exprimer la portion de signal s(t) obtenu par la mesure de f(t), en utilisant la fonction échelon d'Heaviside H, comme \( s(t) = H(t - t_0)H(t_1 - t)f(t) \).

Cette manière de procéder est artificielle et provoque des effets indésirables sur le calcul du spectre, que nous aborderons dans une autre page.

L'algorithme FFT en FORTRAN

Commençons par examiner une implémentation de FFT destinée à être employée avec des langages qui ne disposent pas d'une fft "built-in", comme C ou FORTRAN. Avec ces langages, il est aussi possible, et peut-être plus rapide, d'utiliser les routines fft de la librairie GSL. Vous trouverez un exemple d'utilisation avec C dans la page consacrée à cette librairie.

A propos de FORTRAN, vous constaterez que j'ai rajeuni les meubles ! Je travaille désormais en FORTRAN 90 en utilisant Eclipse/Photran comme environnement de développement et le compilateur/linker gfortran, c'est à dire le GNU Fortran qui est intégré à l'ensemble des compilateurs GCC. Le F90 a le bon goût d'être normalisé et compatible ascendant avec les versions 2000, 2003 et 2008 du FORTRAN.

L'algorithme FFT en FORTRAN

Cette implémentation est basée sur la récursivité de l'algorithme de Danielson-Lanczos, que j'ai évoqué plus haut. F90 permet de déclarer des subroutines recursives, propriété que je vais exploiter ici. Voici le code :

  recursive subroutine FFT(signal, coeff, N)

    implicit none

    double complex, parameter :: I = cmplx(0.0d0,1.0d0)

    integer, intent(in) :: N
    double complex, intent(in) :: signal(0:N-1)
    double complex, intent(inout) :: coeff(0:N-1)
    double complex :: Fk_even(0:N/2-1)
    double complex :: Fk_odd(0:N/2-1)
    double complex :: W
    integer :: k

    if (N == 1) then
       coeff(0) = signal(0)
       return
    else
       call fft(signal(0:N-1:2), Fk_even, N/2)
       call fft(signal(1:N-1:2), Fk_odd, N/2)
       W = exp(-2*pi*I/N)
       do k = 0, N/2-1
          coeff(k) = Fk_even(k) + W**k * Fk_odd(k)
          coeff(k + N/2) = Fk_even(k) - W**k * Fk_odd(k)
       enddo
       return
    endif
  end subroutine FFT

Ce code décompose le signal en deux parties, réelle et imaginaire. Puis il applique le calcul brut d'une TFD sur des bouts de signal de plus en plus petits jusqu'à ce que la taille du bout de signal soit inférieur ou égale à 1.

Attention, dans ce code, le nombre de points doit être égale à une puissance de 2, car il fonctionne par dichotomie. Je n'ai pas fait le contrôle dans ce bout de code très classique et pédagogique, mais si vous voulez l'utiliser dans vos programmes, je vous recommande vivement d'effectuer le contrôle sous peine de plantage en cas de non-respect de la règle.

Exemple d'application de la FFT en FORTRAN

Le programme de mise en oeuvre de la subroutine FFT n'est pas très compliqué. J'utilise les caractéristiques du signal défini plus. Après le calcul de la FFT et du spectre, je stocke les données dans un fichier pour les tracer plus tard. En effet, F90 comme les autres FORTRAN ne brille pas par ses capacités graphiques. Voici le code :

program FFTTest

  implicit none

  double precision, parameter :: pi = 4.0*atan(1.0d0)

  double complex, allocatable :: s(:)
  double complex, allocatable :: Coeff(:)
  double precision, allocatable :: freq(:)
  double precision, allocatable :: spectre(:)
  double precision :: dt, x, fd
  integer :: i, N

  ! Paramètres du signal
  double precision, parameter :: xmin = 0.0
  double precision, parameter :: xmax = 2.0
  double precision, parameter :: f_0 = 1.0

  N = 32                            ! Nombre de points de la FFT

  allocate(s(0:N-1))
  allocate(Coeff(0:N-1))
  allocate(freq(0:N-1))
  allocate(spectre(0:N-1))

  ! Définition du signal
  dt = (xmax - xmin)/N
  fd = 1./(xmax - xmin)
  do i = 0, N-1
     x = i*dt
     s(i) = cmplx(sin(2*pi*f_0*x))
  enddo

  ! Définition du domaine fréquentiel
    do i = 0, N-1
     freq(i) = i*fd
  enddo

  ! Calcul de la FFT
  call FFT(s, Coeff, N)

  ! Calcul du spectre
    do i = 0, N-1
     spectre(i) = abs(Coeff(i))
  enddo

!  Stockage des donnees pour tracé extérieur
   open(20,FILE = 'Spectre.dat')
   do i=1,n
       write(20,*) freq(i), spectre(i)
   enddo
   close(20)

   deallocate(s)
   deallocate(freq)
   deallocate(Coeff)
   deallocate(spectre)
   print *, 'Fin du calcul'

contains
(....)               ! ici est le code de la subroutine FFT
end program FFTTest

Pour tracer le spectre du signal, vous pouvez utiliser gnuplot. Pour ma part, j'ai écrit un script Python TraceSpectreFFT.py qui fait ça pour moi. En voici le code :

# récupération des données dans le fichier Spectre.dat
nomfile = 'Spectre.dat'
freq, spectre = loadtxt(nomfile,dtype = 'f8,f8', unpack = True)
N = len(freq)
# tracé du spectre
fig1 = figure(figsize=(10,8))
grid(True)
stem(freq[0:N//2], spectre[0:N//2],'r')
title('Spectre du signal - Algo FFT en Fortran')
xlabel(u'Fréquence')
ylabel('Amplitude')

Et le spectre obtenu :

Spectre du signal - FFT FORTRAN

Sans surprise, la raie principale du spectre est à 1 Hz, ce qui correspond à la fréquence de mon signal pur.

L'algorithme FFT et Scipy

Le package scipy de Python dispose des fonctions nécessaires pour calculer des TFD par FFT et même beaucoup plus que cela. Vous trouverez ici la documentation nécessaire. Je vous propose un exemple de mise en oeuvre simple et rapide.

Le script TFDRapide.py montre la mise en oeuvre des fonctions du module scipy.fftpack. Le code commence par l'importation des deux fonctions fft() et fftfreq() :

from scipy.fftpack import fft,fftfreq

La définition du domaine fréquentiel par la fonction fftfreq() :

freq = fftfreq(N,dt)

puis le calcul de la TFD par l'algorithme de FFT :

tfd = fft(s)

En retour du calcul, les coefficents de Fourier calculés sont :

0 (-8.60422844084e-16+0j)
1 (0.0612788888701-0.622174999862j)
2 (3.0495991862-15.3313704239j)
3 (-0.379828054934+1.25212529248j)
4 (-0.272607019011+0.658131562496j)
5 (-0.241152272555+0.451164169154j)
6 (-0.226975604581+0.3396929978j)
7 (-0.219253419244+0.267161064346j)
8 (-0.214566898571+0.214566898571j)
9 (-0.211521769862+0.173591430324j)
10 (-0.209454528981+0.13995304188j)
11 (-0.20801511089+0.111186393217j)
12 (-0.207005169857+0.0857443488361j)
13 (-0.206306972572+0.0625825359348j)
14 (-0.205849192441+0.0409459501916j)
15 (-0.205589449523+0.020248793399j)
16 (-0.205505224099+0j)
17 (-0.205589449523-0.020248793399j)
18 (-0.205849192441-0.0409459501916j)
19 (-0.206306972572-0.0625825359348j)
20 (-0.207005169857-0.0857443488361j)
21 (-0.20801511089-0.111186393217j)
22 (-0.209454528981-0.13995304188j)
23 (-0.211521769862-0.173591430324j)
24 (-0.214566898571-0.214566898571j)
25 (-0.219253419244-0.267161064346j)
26 (-0.226975604581-0.3396929978j)
27 (-0.241152272555-0.451164169154j)
28 (-0.272607019011-0.658131562496j)
29 (-0.379828054934-1.25212529248j)
30 (3.0495991862+15.3313704239j)
31 (0.0612788888701+0.622174999862j)

L'affichage du spectre est obtenu en utilisant la première moitié du vecteur tfd, notre fonction étant réelle :

fig2 = figure(figsize=(10,8))
grid(True)
stem(freq[0:N//2],abs(tfd[0:N//2]),'r')
title('Spectre du signal - Algo FFT')
xlabel(u'Fréquence')
ylabel('Amplitude')

Le spectre obtenu est :

Spectre du signal - FFT scipy

L'algorithme de FFT en Python

Si vous le souhaitez, vous pouvez construire un script de calcul d'une TFD avec un algorithme de FFT, sans faire appel aux fonctions de scipy. L'implémentation est basée sur la récurrence de l'algorithme, exactement comme pour l'implémentation en FORTRAN :

# fonction de calcul de la FFT par récursivité
def fft(x):
    N = len(x)
    if N <= 1: return x    
    even = fft(x[0::2])
    odd =  fft(x[1::2])
    T= [cmath.exp(-2j*cmath.pi*k/N)*odd[k] for k in range(N//2)]
    return [even[k] + T[k] for k in range(N//2)] + \
           [even[k] - T[k] for k in range(N//2)]

Le spectre obtenu avec le script FFTPython.py et le signal défini plus haut  :

Spectre du signal - FFT python pur

Comme attendu, il est quasiment identique à celui obtenu avec la fonction fft() de scipy.

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