Quelques rappels sur la transformée de Fourier discrète

Ce tutoriel n’a pas vocation à présenter la transformée de Fourier. Cependant, il existe plusieurs définitions de la transformée de Fourier et même au sein d’un unique domaine, il arrive que l’on en utilise plusieurs. Nous utiliserons la suivante : pour une fonction $f$, sa transformée de Fourier $\hat{f}$ est définie par:

Telle que défini dans la section précédente, la transformée de Fourier d’un signal est une fonction continue de la variable $\nu$. Or, pour représenter un signal quelconque, nous ne pouvons utiliser qu’un nombre fini de valeurs. Pour cela on procède en quatre étapes :

1. On échantillonne (ou discrétise) le signal à analyser. Cela signifie qu’au lieu de travailler sur la fonction qui associe la valeur du signal en fonction de la variable $x$, on va travailler sur une suite discrète de valeurs du signal. Dans le cas de la FFT, on échantillonne avec un pas constant. Par exemple si on regarde un signal temporel comme la valeur d’une tension lue sur un voltmètre, on pourrait enregistrer la valeur à chaque tic d’une montre.
2. On fenêtre le signal discrétisé. Cela signifie que l’on garde uniquement un nombre fini de points du signal.
3. On échantillonne la transformée de Fourier du signal pour obtenir la transformée de Fourier discrète.
4. On fenêtre la transformée de Fourier discrète pour le stockage.

Je vous propose de raisonner sur un signal jouet qui aura la forme d’une Gaussienne. Cela rend le raisonnement un peu plus simple car la transformée de Fourier d’une Gaussienne réelle est elle aussi une gaussienne réelle,¹ ce qui simplifie les représentations graphiques.

Plus formellement, on a:

Intéressons-nous tout d’abord à l’échantillonnage. Mathématiquement, on peut représenter le processus par la multiplication du signal $f$ par un peigne de Dirac de péridode $T$, $ш_T$. Le peigne de Dirac est défini ainsi :

Avec $\delta$ la distribution de Dirac. Voici le graphe que l’on peut obtenir si l’on représente $f$ et $g=ш_T\times f$ ensemble :

La transformée de Fourier de la nouvelle fonction $g$ s’écrit² :

Si on pose $f[k]=f(kT)$ le signal échantillonné et $\nu_{\text{ech}} = \frac{1}{T}$ la fréquence d’échantillonnage, on a:

Si on trace la transformée de Fourier du signal de départ $\hat{f}$ et celle du signal échantillonné $\hat{g}$, on obtient le graphe suivant :

On remarque que l’échantillonnage du signal a amené à la périodisation de sa transformée de Fourier. Cela mène à une propriété importante en traitement du signal : le critère de Nyquist-Shanon, et une de ses conséquences, le repliement de spectre, ou aliasing. Je vous laisse consulter l’article Wikipédia à ce propos si cela vous intéresse, mais on peut se faire une rapide idée de ce qu’il se passe si l’on trace le graphe précédent avec un échantillonnage trop large : les cloches de la transformée du signal échantillonné se superposent.

On peut ensuite regarder le processus de fenêtrage. Il existe plusieurs méthodes qui ont chacune leurs avantages, mais nous nous intéresserons ici uniquement à la fenêtre rectangulaire. Le principe est simple : on ne regarde que les valeurs de $f$ que pour $x$ compris entre $-x_0$ et $+x_0$. Cela signifie que l’on multiplie la fonction $f$ par une fonction porte $\Pi_{x_0}$ qui vérifie :

Graphiquement, voici comment on pourrait représenter $h$ et $f$ ensemble.

Concrètement, c’est équivalent à limiter la somme du peigne de Dirac à un nombre fini de termes. On peut alors écrire la transformée de Fourier de $h=\Pi_{x_0} \times ш_T \times f$ :

Nous pouvons maintenant passer à la dernière étape : l’échantillonnage de la transformée de Fourier. En effet, on ne peut stocker qu’un nombre fini de valeurs sur notre ordinateur or, telle que définit, la fonction $\hat{h}$ est continue. On sait déjà qu’elle est périodique, de période $\nu_{\text{ech}}$, on peut donc ne stocker que les valeurs entre $0$ et $\nu_{\text{ech}}$. Il reste encore à l’échantillonner, et en particulier à trouver le pas d’échantillonnage adéquat. Il est clair que l’on souhaite que l’échantillonnage soit le plus "fin" possible, pour ne manque aucun détail de la transformée de Fourier ! Pour cela on peut s’inspirer de ce qu’il s’est passé lorsque l’on a échantillonné $f$: sa transformée de Fourier est devenue périodique, de période $\nu_{\text{ech}}$. Or la transformée de Fourier inverse (l’opération qui permet de retrouver le signal à partir de sa transformée de Fourier) possède des propriétés similaires à la transformée de Fourier. Cela signifie que si on échantillonne $\hat{h}$ avec un pas d’échantillonnage $\nu_s$, alors sa transformée de Fourier inverse devient périodique de période $1/\nu_s$. Cela donne une limite basse sur les valeurs que peut prendre $\nu_s$ ! En effet, si la transformée inverse a une période inférieure à la largeur de a fenêtre ($1/\nu_s < 2x_0$), alors le signal reconstruit pris entre $-x_0$ et $x_0$ ne correspondra pas au signal initial $f$ !

On choisit donc $\nu_s = \frac{1}{2x_0}$ pour discrétiser $\hat{h}$. On utilise le même processus de multiplication par un peigne de Dirac pour discrétiser. De cette manière on obtient la transformée de Fourier d’une nouvelle fonction $l$ :

Cette notation est un peu compliquée, et on peut s’intéresser plutôt à $\hat{l}[n]=\hat{l}(n\nu_s)$ :

Pour obtenir la dernière ligne, j’ai ré-indicé $f[k]$ de manière à commencer à 0, en notant $N$ le nombre d’échantillons. J’ai ensuite supposé que la taille de la fenêtre correspondait à un nombre entier d’échantillons c’est-à-dire que $2x_0 = N\times T$, ce qui se réécrit $N\times \nu_s = \nu_{\text{ech}}$. Cette expression est la transformée de Fourier discrète du signal.

Ce problème est résolu assez simplement par un fenêtrage de la transformée de Fourier discrète. Puisque la transformée a été périodisée par l’échantillonnage du signal de départ, il suffit de stocker une période de la transformée pour stocker l’intégralité de l’information contenue dans cette dernière. Le choix qui est généralement fait est de garder tous les points entre O et $\nu_{\text{ech}}$. Cela permet de n’utiliser que des $n$ positifs, et on peut reconstruire facilement le graphe de la transformée au besoin en inversant la première et la seconde moitiée de la transformée calculée. En pratique (pour l’implémentation), la transformée de Fourier discrète est donc donnée par :

Pour conclure sur notre fonction d’exemple, on obtient le graphe suivant :

Nous avons donc à notre disposition l’expression de la transformée de Fourier discrète d’un signal $f$:

Cette expression est celle d’un produit matriciel qui ressemblerait à ceci:

avec

Les connaisseurs remarqueront qu’il s’agit ici d’une matrice de Vandermonde sur les racines de l’unité.

On peut donc implémenter relativement simplement ce calcul !

function naive_dft(x)
	N = length(x)
	k = reshape(0:(N-1), 1, :)
	n = 0:(N-1)
	M = @. exp(-2im * π * k * n / N)
	M * x
end

Et pour vérifier qu’il donne effectivement le bon résultat, il suffit de le comparer avec une implémentation de référence :

using FFTW

let
	a = rand(1024)
	b = fft(a)
	c = naive_dft(a)
	b ≈ c
end

Le dernier bloc s’évalue à true, ce qui confirme que nous ne sommes pas totalement à côté de la plaque !

Cependant, ce code est-il efficace ? On peut le vérifier en comparant l’empreinte mémoire et la vitesse d’exécution.

using BenchmarkTools

let
	@benchmark fft(a) setup=(a = rand(1024))
end

BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  26.134 μs …  13.261 ms  ┊ GC (min … max): 0.00% … 56.15%
 Time  (median):     32.502 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   37.547 μs ± 184.201 μs  ┊ GC (mean ± σ):  4.03% ±  0.82%

    ▆▆▄▅▅▇▇▆██▅▅▆▅▃▁              ▁▁▁▁▂▂▂▃▂▃▂▂▂▂▁▁  ▁▁         ▃
  ▇▆████████████████████████▇█████████████████████▇████▆▇▅▆▆▇▇ █
  26.1 μs       Histogram: log(frequency) by time      61.5 μs <

 Memory estimate: 34.19 KiB, allocs estimate: 34.

let
	@benchmark naive_dft(a) setup=(a = rand(1024))
end

BenchmarkTools.Trial: 52 samples with 1 evaluation.
 Range (min … max):  54.047 ms … 131.510 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     96.858 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   96.784 ms ±  10.861 ms  ┊ GC (mean ± σ):  0.92% ± 2.47%

                                              ▁ █▂              
  ▃▁▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▃▁▇█▃██▃▃▇▃▁▃▄▃▄▁▃▄ ▁
  54 ms           Histogram: frequency by time          109 ms <

 Memory estimate: 16.03 MiB, allocs estimate: 4.

Notre implémentation est donc vraiment lente (environ 10000 fois plus) et possède une empreinte mémoire très élevée (environ 500 fois) par comparaison avec l’implémentation de référence ! Pour améliorer cela, nous allons implémenter la transformée de Fourier rapide.

Avant de replonger les mains dans le cambouis, posons-nous d’abord la question : est-il bien nécessaire de chercher à améliorer cet algorithme ?

Avant de répondre directement, intéressons-nous à quelques applications de la transformée de Fourier et de la transformée de Fourier discrète.

La transformée de Fourier a d’abord une foultitude d’applications théoriques, que ce soit pour résoudre des équations différentielles, en traitement du signal ou en physique quantique. Elle possède également des applications pratiques en optique et en spectroscopie.

La transformée de Fourier discrète a également de nombreuses applications, en analyse des signaux, pour la compression de données, la multiplication de polynômes ou le calcul de produits de convolutions.

Notre implémentation naïve de la transformée de Fourier discrète a une complexité en temps et en mémoire en $\mathcal{O}(N^2)$ avec $N$ la taille de l’échantillon d’entrée, cela est dû au stockage de la matrice et au temps de calcul du produit matriciel. Concrètement, si on souhaitait analyser un signal sonore de 3 secondes échantilloné à 44kHz avec des données stockées sur avec des flottants simple précision (4 octets), il faudrait donc environ $2\times(44000\times3)^2\times 4\approx100\;000\;000\;000$ octets de mémoire (un nombre complexe est stocké sur 2 flottants). On peut aussi estimer le temps nécessaire à faire ce calcul. Le temps médian pour 1024 points était de 38.367 ms. Pour notre signal de 3 secondes, il faudrait environ $38.867\times\left(\frac{44000\times3}{1024}\right)^2\approx 637\;537$ millisecondes, soit plus de 10 minutes !

On comprend aisément l’intérêt de réduire la complexité du calcul. En particulier l’algorithme de la transformée de Fourier rapide (utilisé par l’implémentation de référence) possède une complexité en $\mathcal{O}(N\log N)$. D’après notre benchmark, l’algorithme traite une entrée de 1024 points en 23.785µs. Il devrait donc traiter le signal sonore en environ $23.785\times\frac{44000\times3\times\log(44000\times3)}{1024\log1024}\approx 5\;215$ micro-secondes, soit environ 120000 fois plus rapidement que notre algorithme. On peut vraiment dire que le fast de Fast Fourier Transform n’est pas volé !