Table des matières PDFPython

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 uk avec k=0,...,N-1. Ces échantillons ont été prélevés avec une période d'échantillonnage Te. La durée totale du signal échantillonné est donc T=NTe.

Pour simplifier les notations, on considère que Te constitue l'unité de temps, ce qui revient à poser Te=1. Les instants correspondant à ces échantillons sont alors tk=k.

Considérons une sinusoïde de fréquence n/N, où n est un entier (n=0,...,N-1). La représentation complexe de cette sinusoïde s'écrit :

Les N échantillons de S aux instants tk=k sont :

La transformée de Fourier discrète (TFD) consiste à multiplier ces échantillons par ceux du signal et à calculer la somme des produits :

Le nombre complexe Fn est l'amplitude de la composante de fréquence n/N (en réalité n/(NTe) du signal. Plus précisément, Fn est le coefficient de Fourier de la série de Fourier (écrite sous forme complexe) d'un signal de période T=NTe=N.

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

Attention : le nombre N est placé en exposant mais il ne s'agit pas d'un exposant; ce type de notation est souvent utilisé en calcul numérique.

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

En terme de fréquences réelles, cela implique la périodicité du spectre avec une période , qui est la fréquence d'échantillonnage.

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

désigne le complexe conjugé de z. La représentation du spectre consiste à tracer le module de en ordonnée et la fréquence n/(NTe) en abscisse. La propriété :

implique que la composante du spectre de fréquence n/(NTe) est égale à celle de fréquence .

2. Transformée de Fourier rapide

La méthode basique de calcul de la TFD consiste à calculer les N nombres (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 uk par une valeur complexe dont l'évaluation nécessite le calcul d'un cosinus et d'un sinus, deux fonctions dont l'évaluation est très coûteuse en temps de calcul. La complexité temporelle de cet algorithme est donc O(N2). 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=2p. 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 :

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, mutlipliée par le coefficient suivant :

Notons la TFD des termes de rangs pairs et la TFD des termes de rangs impairs. Chacune de ces TFD contient N/2 nombres :

Ces deux TFD sont en fait définis pour grace à la périodicité de la TFD :

qui fait que, pour on a :

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 :

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 est appliquée N fois. L'ordre de grandeur maximal de la complexité est donc O(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(liste):
    N = len(liste)
    if N==1:
        return liste
    pair = []
    impair = []
    for k in range(0,N,2):
        pair.append(liste[k])
    for k in range(1,N,2):
        impair.append(liste[k])
    tfd_pair = fft(pair)
    tfd_impair = fft(impair)
    tfd = [0]*N
    W = numpy.exp(-1j*2*numpy.pi/N)
    N2 = int(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 avec une fonction comportant 3 harmoniques, échantillonnée sur une période :

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)
         
from matplotlib.pyplot import *
spectre = numpy.absolute(tfd)*2/N
figure(figsize=(10,4))
stem(k,spectre,'r')
            
figA.svgFigure pleine page

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 donc en général une implémentation itérative.

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

inversion.svgFigure 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 à 2q éléments. Il y a 2p-q TFD pour le niveau q.

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

avec et . Les opérations à effectuer pour une TFD à 2 éléments sont schématisées sur la figure suivante :

TFD.svgFigure pleine page

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

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 2p-q groupes de 2q éléments à calculer. Considérons le groupe numéro m, dont les indices vont de i=m2q à i=(m+1)2q-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 m2q à m2q+2q-1-1, et une boucle sur la seconde moitié des indices, soit i variant de m2q+2q-1 à (m+1)2q-1. Les opérations à effectuer sont schématisées sur la figure suivante.

TFD-niveau-q.svgFigure 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')
                    
figB.svgFigure pleine page
Creative Commons LicenseTextes et figures sont mis à disposition sous contrat Creative Commons.