Transformée de Fourier rapide

1. Transformée de Fourier discrète

La transformée de Fourier discrète est une transformation mathématique permettant d'obtenir le spectre de fréquence d'un signal échantillonné.

On dispose de N échantillons d'un signal, que l'on note u_k avec k=0,...,N-1. Ces échantillons ont été prélevés avec une période d'échantillonnage T_e. La durée du signal échantillonné est donc T=NT_e.

Pour obtenir le spectre du signal analogique, on cherche à calculer les coefficients de Fourier complexes d'une fonction périodique de période T :

$$C_n=\frac{2}{T}{\int }_0^{T}u\left(t\right)exp(-jn\frac{2\pi }{T}t)dt\text{(1)}$$

La méthode des rectangles permet d'obtenir une approximation de cette intégrale :

$$C_n\approx \frac{2}{N{T}_e}\sum _{k=0}^{N-1}{u}_{k}exp(-jn\frac{2\pi }{N{T}_e}k{T}_e)T_e=\frac{2}{N}\sum _{k=0}^{N-1}{u}_{k}exp(-jn\frac{2\pi }{N}k)\text{(2)}$$

On définit la transformée de Fourier discrète (TFD) par :

Posons :

$$S_{n,k}=exp(-j2\pi \frac{n}{N}k)\text{(4)}$$

La TFD s'écrit :

$$F_n=\sum _{k=0}^{N-1}{u}_k{S}_{n,k}\text{(5)}$$

Le nombre complexe F_n, multiplié par 2/N, est une approximation du coefficient de Fourier complexe C_n de la fonction périodique de période T. La TFD définit le spectre du signal échantillonné. Le nombre complexe F_n, multiplié par 2/N, est la composante de fréquence n/(NT_e) du spectre du signal échantillonné.

La TFD transforme N nombres, que nous supposerons réels, en N nombres complexes. Introduisons une notation faisant apparaître le nombre N :

$$F_n^{\left(N\right)}=\sum _{k=0}^{N-1}{u}_k{S}_{n,k}^{\left(N\right)}\text{(6)}$$

Le nombre N est placé en exposant mais il ne s'agit pas d'un exposant.

La TFD est périodique, de période N :

$$F_{n+N}^{\left(N\right)}=F_n^{\left(N\right)}\text{(7)}$$

En terme de fréquences, cela implique la périodicité du spectre avec une période $\frac{N}{N{T}_e}=\frac{1}{T_e}$, qui est la fréquence d'échantillonnage.

Les échantillons du signal étant des nombres réels, nous avons :

$$F_{N-n}^{\left(N\right)}={\left[F_n^{\left(N\right)}\right]}^{*}\text{(8)}$$

où ${\left[z\right]}^{*}$ désigne le complexe conjugé de z. La représentation du spectre consiste à tracer le module de $F_n^{\left(N\right)}$ en ordonnée et la fréquence n/(NT_e) en abscisse. La propriété :

$$\left|F_n^{\left(N\right)}\right|=\left|F_{N-n}^{\left(N\right)}\right|\text{(9)}$$

implique que la composante du spectre de fréquence n/(NT_e) est égale à celle de fréquence $\frac{N-n}{N{T}_e}=\frac{1}{T_e}-\frac{n}{N{T}_e}$.

2. Transformée de Fourier rapide

La méthode basique de calcul de la TFD consiste à calculer les N nombres $F_n^{\left(N\right)}$ (n=0,...,N-1) tels qu'ils sont écrits dans la TFD. Chacun de ces nombres est la somme de N termes, chacun étant le produit de u_k par une valeur complexe dont l'évaluation nécessite le calcul d'un cosinus et d'un sinus. La complexité temporelle de cet algorithme est donc Θ(N²). Le nombre d'échantillons est de l'ordre de plusieurs milliers à plusieurs millions et cette méthode conduit à des calculs extrêmement longs.

L'algorithme de transformée de Fourier rapide (FFT pour Fast Fourier Transform) a été inventé en 1965 par Cooley et Tukey. Il s'agit d'un algorithme de type diviser pour régner, qui fonctionne si le nombre d'échantillons est une puissance de 2, soit N=2^p. En pratique, cette condition n'est pas contraignante car il suffit d'ajouter si nécessaire des zéros pour allonger la liste des échantillons (l'ajout de zéro permet d'ailleurs d'améliorer la précision du spectre).

L'algorithme FFT repose sur le regroupement des termes de rangs pairs d'une part, des termes de rangs impairs d'autre part :

$$\begin{array}{ccc}\hfill F_n^{\left(N\right)}=& \sum _{k=0}^{\frac{N}{2}-1}{u}_{2k}exp(-j2\pi \frac{n}{N}2k)+exp(-j2\pi \frac{n}{N})\sum _{k=0}^{\frac{N}{2}-1}{u}_{2k+1}exp(-j2\pi \frac{n}{N}2k)\hfill & \hfill \text{(10)}\\ \hfill & =\sum _{k=0}^{\frac{N}{2}-1}{u}_{2k}{S}_{n,k}^{(N/2)}+exp(-j2\pi \frac{n}{N})\sum _{k=0}^{\frac{N}{2}-1}{u}_{2k+1}{S}_{n,k}^{(N/2)}\hfill & \hfill \text{(11)}\end{array}$$

Le premier terme est la TFD des échantillons de rangs pairs (TFD de taille N/2). Le second terme est la TFD des échantillons de rangs impairs, multipliée par le coefficient suivant :

$$W_n^{\left(N\right)}=exp(-j2\pi \frac{n}{N})\text{(12)}$$

Notons $P_n^{(N/2)}$ la TFD des termes de rangs pairs et $I_n^{(N/2)}$ la TFD des termes de rangs impairs. Chacune de ces TFD contient N/2 nombres :

$$\begin{array}{ccc}\hfill & P_n^{(N/2)}=\sum _{k=0}^{\frac{N}{2}-1}{u}_{2k}{S}_{n,k}^{(N/2)}pourn=0,\dots \frac{N}{2}-1\hfill & \hfill \text{(13)}\\ \hfill & I_n^{(N/2)}=\sum _{k=0}^{\frac{N}{2}-1}{u}_{2k+1}{S}_{n,k}^{(N/2)}pourn=0,\dots \frac{N}{2}-1\hfill & \hfill \text{(14)}\end{array}$$

Ces deux TFD sont en fait définis pour $n=0,\dots ,N-1$ grace à la périodicité de la TFD :

$$\begin{array}{ccc}\hfill & P_{N/2+n}^{(N/2)}=P_n^{(N/2)}\hfill & \hfill \text{(15)}\\ \hfill & I_{N/2+n}^{(N/2)}=I_n^{(N/2)}\hfill & \hfill \text{(16)}\end{array}$$

qui fait que, pour $n=N/2,\dots ,N-1$ on a :

$$\begin{array}{ccc}\hfill & P_n^{(N/2)}=P_{n-N/2}^{(N/2)}\hfill & \hfill \text{(17)}\\ \hfill & I_n^{(N/2)}=I_{n-N/2}^{(N/2)}\hfill & \hfill \text{(18)}\end{array}$$

La TFD des échantillons de rangs pairs (respectivement de rangs impairs) est calculée récursivement de la même manière puisque N/2 est divisible par 2. Lorsque le calcul de ces deux TFD est terminé, les N valeurs de la TFD sont obtenues par :

$$F_n^{\left(N\right)}=P_n^{(N/2)}+W_n^{\left(N\right)}{I}_n^{(N/2)}pourn=0,\dots ,N-1\text{(19)}$$

La récursion s'arrête lorsqu'on arrive à une TFD de taille 1. La TFD d'un seul nombre est ce nombre.

Voyons la complexité temporelle de cet algorithme. Il y a p=ln(N)/ln(2) niveaux de récursion. À chaque niveau, la relation $\mathrm{(19)}$ est appliquée N fois. L'ordre de grandeur de la complexité est donc Θ(N ln(N)), ce qui en fait un algorithme beaucoup plus efficace que la méthode basique. L'algorithme FFT a permis le développement de l'analyse spectrale numérique, à une époque où les ordinateurs étaient beaucoup moins rapides qu'aujourd'hui.

3. Implémentation récursive

L'implémentation récursive de l'algorithme FFT est immédiate :

def fft(u):
    # u : numpy.ndarray
    N = len(u)
    if N==1:
        return u
    pair = u[0::2]
    impair = u[1::2]
    tfd_pair = fft(pair)
    tfd_impair = fft(impair)
    tfd = numpy.zeros(N,dtype=numpy.complex64)
    W = numpy.exp(-1j*2*numpy.pi/N)
    N2 = N//2
    for n in range(N2):
        tfd[n] = tfd_pair[n]+tfd_impair[n]*W**n
    for n in range(N2,N):
        tfd[n] = tfd_pair[n-N2]+tfd_impair[n-N2]*W**n
    return tfd
            

On teste la fonction fft avec un signal comportant 3 harmoniques, échantillonné sur une période. Si la fréquence d'échantillonnage vérifie la condition de Nyquist-Shannon pour le dernier harmonique, la TFD donne exactement les coefficients de Fourier (aux erreurs d'arrondis près).

import numpy
from matplotlib.pyplot import *

p = 5
N = 2**p
k = numpy.arange(N)
u = numpy.sin(2*numpy.pi*k/N)+\
    0.5*numpy.sin(4*numpy.pi*k/N)+0.25*numpy.cos(10*numpy.pi*k/N)
tfd = fft(u)
         
from matplotlib.pyplot import *
spectre = numpy.absolute(tfd)*2/N
figure(figsize=(10,4))
stem(k,spectre,'r')
            

Figure pleine page

Comment calcule-t-on la TFD si le nombre d'échantillons N n'est pas une puissance de 2. Si possible, on doit éviter cette situation au moment où l'on choisit ce nombre. Si malgré tout on dispose d'un signal dont le nombre d'échantillons n'est pas une puissance de 2, la meilleure solution dans le cas de l'analyse spectrale d'un signal est d'ajouter des zéros au signal afin d'atteindre un nombre N puissance de 2. Cette méthode est expliquée dans Analyse spectrale d'un signal numérique. Si cela est vraiment nécessaire, il est possible de calculer TFD pour un nombre N quelconque mais ce calcul ne sera pas optimal. Si N est pair, on fait la séparation en termes de rangs pairs et de rangs impairs puis on répète récursivement tant que le nombre d'échantillons obtenus est pair. Lorsqu'on obtient un nombre d'échantillons impair, on calcule la TFD de manière directe. On a donc besoin d'une fonction qui fait un calcul direct de la TFD :

def tfd_directe(u):
    N = len(u)
    F = numpy.zeros(N,dtype=numpy.complex64)
    a = -1j*2*numpy.pi/N
    for n in range(N):
        b = a*n
        for k in range(N):
            F[n] += u[k]*numpy.exp(b*k)
    return F
                

Voici la nouvelle fonction qui calcule par FFT (tant que N est pair) :

def fft2(u):
    # u : numpy.ndarray
    N = len(u)
    if N==1:
        return u
    if N%2==1:
        return tfd_directe(u)
    pair = u[0::2]
    impair = u[1::2]
    tfd_pair = fft2(pair)
    tfd_impair = fft2(impair)
    tfd = numpy.zeros(N,dtype=numpy.complex64)
    W = numpy.exp(-1j*2*numpy.pi/N)
    N2 = N//2
    for n in range(N2):
        tfd[n] = tfd_pair[n]+tfd_impair[n]*W**n
    for n in range(N2,N):
        tfd[n] = tfd_pair[n-N2]+tfd_impair[n-N2]*W**n
    return tfd
            

Voici un exemple :

N = 40
k = numpy.arange(N)
u = numpy.sin(2*numpy.pi*k/N)+\
    0.5*numpy.sin(4*numpy.pi*k/N)+0.25*numpy.cos(10*numpy.pi*k/N)
tfd = fft2(u)
         
spectre = numpy.absolute(tfd)*2/N
figure(figsize=(10,4))
stem(k,spectre,'r')
            

Figure pleine page

Pour cet exemple, le calcul direct de la TFD est effectué sur 4 blocs de 5 échantillons, puis 2 blocs de 10 valeurs sont obtenus puis un bloc de 20.

La fonction numpy.fft.fft effectue le calcul de la TFD de cette manière.

4. Implémentation itérative

Le code précédent n'est pas optimal car il nécessite, à chaque appel de la fonction, la création et le remplissage de deux listes, puis la création d'une nouvelle liste pour contenir la TFD. On préfère en général une implémentation itérative qui évite ces opérations de recopie en mémoire.

Une implémentation itérative consiste à commencer par calculer toutes les TFD à 2 éléments, puis toutes les TFD à 4 éléments, et ainsi de suite jusqu'à la TFD à N éléments. Il faut pour cela déterminer dans quel ordre les échantillons (u_k) doivent être mis pour calculer les TFD à 2 éléments. Voyons le cas particulier p=3. Sur la figure suivante, la liste initiale est représentée en haut, avec les indices en base 10 et en base 2.

Figure pleine page

La première division en une liste d'éléments de rangs pairs et une liste d'éléments de rangs impairs se fait en fonction du bit de poids faible de l'indice initial (c.a.d. le bit situé à droite). La seconde division se fait en fonction du deuxième bit et la dernière division se fait en fonction du troisième bit (toujours de l'indice initial). Pour effectuer la FFT de manière itérative, il faut partir de la dernière ligne, appliquer tout d'abord les 4 TFD à 2 éléments, puis les 2 TFD à 4 éléments et enfin la TFD à 8 éléments. La liste de départ de l'algorithme (en bas sur la figure) contient bien sûr les N échantillons mais dans un ordre différent de l'ordre chronologique. Par exemple, l'échantillon k=4 se trouve à l'indice j=1. On constate que le rangement des échantillons dans la liste de départ se fait en inversant l'ordre des bits de l'indice. Pour savoir où placer l'échantillon d'indice k, il suffit d'inverser l'ordre des bits dans la représentation binaire de k, ce qui donne l'indice j où il faut le placer.

Une fois les échantillons placés dans le bon ordre, on peut commencer la boucle des p niveaux de TFD. Le premier niveau consiste à calculer les TFD à 2 éléments, le deuxième les TFD à 4 éléments et ainsi de suite jusqu'au p-ième niveau, qui calcule la TFD à N éléments. Les TFD de niveau q (allant de 1 à p) sont des TFD à 2^q éléments. Il y a 2^p-q TFD pour le niveau q.

Voyons comment se fait le premier niveau. Une TFD à 2 éléments est définie par la relation :

$$F_n^{\left(2\right)}=P_n^{\left(1\right)}+W_n^{\left(2\right)}{I}_n^{\left(1\right)}pourn=0,1\text{(20)}$$

avec $P_1^{\left(1\right)}=P_0^{\left(1\right)}$ et $I_1^{\left(1\right)}=I_0^{\left(1\right)}$. Les opérations à effectuer pour une TFD à 2 éléments sont schématisées sur la figure suivante :

Figure pleine page

Une TFD à 4 éléments est définie par la relation :

$$F_n^{\left(4\right)}=P_n^{\left(2\right)}+W_n^{\left(4\right)}{I}_n^{\left(2\right)}pourn=0,1,2,3\text{(21)}$$

Les opérations à effectuer sont schématisées ci-dessus.

Supposons que les TFD de niveau q-1 soient stockées dans un tableau à N éléments. Nous avons besoin d'un second tableau à N éléments pour générer les TFD de niveau q. Au total, l'algorithme peut donc fonctionner avec seulement deux tableaux à N éléments (on procédera à un échange des références à chaque niveau). Soit A le tableau contenant les TFD de niveau q-1 et B celui devant recevoir les TFD de niveau q. Il y a 2^p-q groupes de 2^q éléments à calculer. Considérons le groupe numéro m, dont les indices vont de i=m2^q à i=(m+1)2^q-1. Il se calcule avec les éléments de mêmes indices du tableau A. Pour ce calcul, il faut écrire une boucle sur la première moitié des indices, soit i variant de m2^q à m2^q+2^q-1-1, et une boucle sur la seconde moitié des indices, soit i variant de m2^q+2^q-1 à (m+1)2^q-1. Les opérations à effectuer sont schématisées sur la figure suivante.

Figure pleine page

Pour effectuer l'inversion des bits des indices, nous avons besoin d'une fonction qui convertit un nombre en représentation binaire à p bits, effectue l'inversion des bits puis renvoie la valeur en base 10 du résultat. Voici cette fonction :

def inversion_bits(i,p): # p = nombre de bits
    c = "{0:0%db}"%p
    binaire = c.format(i)
    inv = binaire[::-1]
    return int(inv,2)
                    

Voici la fonction fft(u,p), qui effectue le calcul itératif de la TFD. Les tableaux u, A et B sont des tableaux de numpy (ndarray) d'éléments de type numpy.complex64.

import numpy
                    
def fft(u,p):
    N=len(u)
    if N!=2**p:
        print("Erreur de taille")
        return
    A = numpy.zeros(N,dtype=numpy.complex64)
    B = numpy.zeros(N,dtype=numpy.complex64)
    for k in range(N):
        j = inversion_bits(k,p)
        A[j] = u[k]
    for q in range(1,p+1): # FFT à 2**q éléments
        taille = 2**q
        taille_precedente = 2**(q-1)
        nombre_tfd = 2**(p-q) # nombre de TFD à calculer
        for m in range(nombre_tfd):
            position = m*taille
            phi = -1j*2*numpy.pi/taille
            for i in range(taille_precedente):
                W = numpy.exp(phi*i)
                B[position+i] = A[position+i] + W*A[position+taille_precedente+i]
            for i in range(taille_precedente,taille):
                W = numpy.exp(phi*i)
                B[position+i] = A[position+i-taille_precedente] + W*A[position+i]
        (A,B)=(B,A) # échange des références des tableaux
    return A

import numpy
from matplotlib.pyplot import *

p = 5
N = 2**p
u = numpy.zeros(N,dtype=complex)
k = numpy.arange(N)
u = numpy.sin(2*numpy.pi*k/N)\
    +0.5*numpy.sin(4*numpy.pi*k/N)+0.25*numpy.cos(10*numpy.pi*k/N)
tfd = fft(u,p)
         

spectre = numpy.absolute(tfd)*2/N
figure(figsize=(10,4))
stem(k,spectre,'r')
                    

Figure pleine page