Licence CC BY-SA

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

La transformée de Fourier discrète est une transformation qui découle de la transformée de Fourier et est, comme son nom l’indique, adaptée pour des signaux discrets. Dans cette première partie je vous propose de découvrir comment construire la transformée de Fourier discrète puis comprendre pourquoi la transformée de Fourier rapide est utile.

La transformée de Fourier

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 ff, sa transformée de Fourier f^\hat{f} est définie par:

f^(ν)=+f(x)e2iπνxdx\hat{f}(\nu) = \int_{-\infty}^{+\infty}f(x)e^{-2i\pi \nu x}\text{d}x

De la transformée de Fourier à la transformée de Fourier discrète

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 xx, 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,1 ce qui simplifie les représentations graphiques.

Le signal qui nous servira d'exemple
Le signal qui nous servira d’exemple

Plus formellement, on a:

f(x)=ex2,  f^(ν)=πe(πν)2f(x) = e^{-x^2},\;\hat{f}(\nu)=\sqrt{\pi}e^{-(\pi\nu)^2}

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

шT(x)=k=+δ(xkT)ш_T(x) = \sum_{k=-\infty}^{+\infty}\delta(x-kT)

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

Le signal et le signal échantillonné.
Le signal et le signal échantillonné.

La transformée de Fourier de la nouvelle fonction gg s’écrit2 :

g^(ν)=+k=+δ(xkT)f(x)e2iπxνdx=k=++δ(xkT)f(x)e2iπxνdx=k=+f(kT)e2iπkTν\begin{aligned} \hat{g}(\nu) &= \int_{-\infty}^{+\infty} \sum_{k=-\infty}^{+\infty}\delta(x-kT) f(x) e^{-2i\pi x \nu} \text{d}x \\ &= \sum_{k=-\infty}^{+\infty}\int_{-\infty}^{+\infty}\delta(x-kT) f(x) e^{-2i\pi x \nu}\text{d}x \\ &= \sum_{k=-\infty}^{+\infty}f(kT)e^{-2i\pi kT\nu} \end{aligned}

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

g^=k=+f[k]e2iπkννech\hat{g} = \sum_{k=-\infty}^{+\infty}f[k]e^{-2i\pi k\frac{\nu}{\nu_{\text{ech}}}}

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

Transformée de Fourier du signal et de son signal échantillonné.
Transformée de Fourier du signal et de son signal échantillonné.

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.

Transformée de Fourier du signal et de son signal échantillonné, illustration du phénomène de repliement de spectre.
Transformée de Fourier du signal et de son signal échantillonné, illustration du phénomène de repliement de spectre.

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 ff que pour xx compris entre x0-x_0 et +x0+x_0. Cela signifie que l’on multiplie la fonction ff par une fonction porte Πx0\Pi_{x_0} qui vérifie :

Πx0(x)={1  si  x[x0,x0]0  sinon\Pi_{x_0}(x) = \left\{ \begin{aligned} 1 &\;& \text{si}\; x\in[-x_0,x_0] \\ 0 &\;& \text{sinon} \end{aligned} \right.

Graphiquement, voici comment on pourrait représenter hh et ff ensemble.

Signal échantillonné et fenêtré
Signal échantillonné et fenêtré

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=Πx0×шT×fh=\Pi_{x_0} \times ш_T \times f :

h^(ν)=k=k0+k0f[k]e2iπkννech\hat{h}(\nu) = \sum_{k=-k_0}^{+k_0}f[k]e^{-2i\pi k\frac{\nu}{\nu_{\text{ech}}}}

Le choix du fenêtrage n’est absolument pas anodin, et peut mener à des problèmes innatendus si on l’ignore. Là encore je vous conseille de consulter l’article Wikipédia associé au besoin.

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 h^\hat{h} est continue. On sait déjà qu’elle est périodique, de période νech\nu_{\text{ech}}, on peut donc ne stocker que les valeurs entre 00 et νech\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é ff: sa transformée de Fourier est devenue périodique, de période νech\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 h^\hat{h} avec un pas d’échantillonnage νs\nu_s, alors sa transformée de Fourier inverse devient périodique de période 1/νs1/\nu_s. Cela donne une limite basse sur les valeurs que peut prendre νs\nu_s ! En effet, si la transformée inverse a une période inférieure à la largeur de a fenêtre (1/νs<2x01/\nu_s < 2x_0), alors le signal reconstruit pris entre x0-x_0 et x0x_0 ne correspondra pas au signal initial ff !

On choisit donc νs=12x0\nu_s = \frac{1}{2x_0} pour discrétiser h^\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 ll :

l^(ν)=n=+δ(νnνs)k=k0+k0f[k]e2iπknνsνech\begin{aligned} \hat{l}(\nu) = \sum_{n=-\infty}^{+\infty} \delta(\nu-n\nu_s) \sum_{k=-k_0}^{+k_0}f[k]e^{-2i\pi k\frac{n\nu_s}{\nu_{\text{ech}}}} \end{aligned}

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

l^[n]=l^(nνs)=k=k0+k0f[k]e2iπknνsνech=k=0N1f[k]e2iπknN\begin{aligned} \hat{l}[n] = \hat{l}(n\nu_s) &=& \sum_{k=-k_0}^{+k_0}f[k]e^{-2i\pi k\frac{n\nu_s}{\nu_{\text{ech}}}}\\ &=& \sum_{k=0}^{N-1}f[k]e^{-2i\pi k\frac{n}{N}} \end{aligned}

Pour obtenir la dernière ligne, j’ai ré-indicé f[k]f[k] de manière à commencer à 0, en notant NN 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 2x0=N×T2x_0 = N\times T, ce qui se réécrit N×νs=νechN\times \nu_s = \nu_{\text{ech}}. Cette expression est la transformée de Fourier discrète du signal.

Échantillonnage de la transformée de Fourier du signal échantillonné pour obtenir la transformée de Fourier discrète.
Échantillonnage de la transformée de Fourier du signal échantillonné pour obtenir la transformée de Fourier discrète.

On voit que la fréquence d’échantillonnage n’entre pas en compte dans cette équation, et il y a de nombreuses applications où on oublie purement et simplement l’existence de cette fréquence.

Il reste un dernier point à éclaircir: cette transformée discrète est définie pour une infinité (discrète) de valeurs de nn. Comment la stocker sur notre ordinateur ?

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 νech\nu_{\text{ech}}. Cela permet de n’utiliser que des nn 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 :

n=0...(N1),  f^[n]=k=0N1f[k]e2iπknN\boxed{ \forall n=0...(N-1),\; \hat{f}[n] = \sum_{k=0}^{N-1}f[k]e^{-2i\pi k\frac{n}{N}} }

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

Fenêtrage de la transformée de Fourier discrète pour le stockage.
Fenêtrage de la transformée de Fourier discrète pour le stockage.

  1. On dit que les gaussiennes sont des fonctions propres de la transformée de Fourier.
  2. Il faudrait ici justifier que l’on peut inverser les signes somme et intégrale.

Calculer la transformée de Fourier discrète

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

f^[n]=k=0N1f[k]e2iπknN\hat{f}[n] = \sum_{k=0}^{N-1}f[k]e^{-2i\pi k\frac{n}{N}}

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

f^=Mf\hat{f} = \mathbf{M} \cdot f

avec

M=(11111e2iπ1×1/Ne2iπ2×1/Ne2iπ1×(N1)/N1e2iπ1×2×1/Ne2iπ2×2/Nee2iπ(N2)×(N1)/N1e2iπ(N1)×1/Nee2iπ(N1)×(N2)/Ne2iπ(N1)×(N1)/N)\mathbf{M} = \begin{pmatrix} 1 & 1 & 1 & \dots & 1 \\ 1 & e^{-2i\pi 1 \times 1 / N} & e^{-2i\pi 2 \times 1 / N} & \dots & e^{-2i\pi 1\times (N-1)/N} \\ 1 & e^{-2i\pi 1 \times 2 \times 1 / N} & e^{-2i\pi 2 \times 2 / N} & \ddots & \vdots \\ \vdots & \vdots & \ddots & \ddots & e^{e-2i\pi (N-2)\times (N-1) / N} \\ 1 & e^{-2i\pi (N-1) \times 1/N} & \dots & e^{e-2i\pi (N-1) \times (N-2) / N} & e^{-2i\pi (N-1)\times (N-1) / N} \end{pmatrix}

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

La macro @. ligne 5 permet de vectoriser le calcul de l’expression qu’elle englobe (exp(-2im * π * k * n / N)). En effet la fonction exp et les opérateurs de division et multiplication sont définis pour des scalaires. Cette macro permet d’informer Julia qu’il doit appliquer les opérations scalaires terme à terme.

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 !

J’utilise l’opérateur pour comparer plutôt que == afin d’autoriser des petites différences, notamment à cause des erreurs d’arrondies sur les flottants.

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: 
  memory estimate:  34.67 KiB
  allocs estimate:  34
  --------------
  minimum time:     20.943 μs (0.00% GC)
  median time:      23.785 μs (0.00% GC)
  mean time:        28.827 μs (1.14% GC)
  maximum time:     6.736 ms (48.98% GC)
  --------------
  samples:          10000
  evals/sample:     1
let
	@benchmark naive_dft(a) setup=(a = rand(1024))
end
BenchmarkTools.Trial: 
  memory estimate:  16.03 MiB
  allocs estimate:  4
  --------------
  minimum time:     35.746 ms (0.00% GC)
  median time:      38.367 ms (0.00% GC)
  mean time:        39.493 ms (1.33% GC)
  maximum time:     51.803 ms (0.00% GC)
  --------------
  samples:          127
  evals/sample:     1

Comme vous avez pu le constater le temps d’exécution maximal de l’implémentation de référence est plus élevé de deux ordres de grandeur que les temps d’exécution moyen et médian. On doit cela à la compilation Just in time (JIT) de Julia. Si nous étions en train d’écrire une vraie bibliothèque Julia on pourrait envisager d’optimiser notre code pour qu’il compile rapidement. Nous nous contenterons d’ignorer le temps d’exécution maximum dans ce tutoriel, qui correspond uniquement au temps de compilation pour la première exécution du code. Je vous renvoie à la documentation de BenchmarkTools.jl pour plus d’informations.

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.

Pourquoi un algorithme de 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 O(N2)\mathcal{O}(N^2) avec NN 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×(44000×3)2×4100  000  000  0002\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×(44000×31024)2637  53738.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 O(NlogN)\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×44000×3×log(44000×3)1024log10245  21523.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é !


Nous avons vu comment la transformée de Fourier discrète était construite, puis nous avons tenté de manière naïve de l’implémenter. Bien que cette implémentation soit relativement simple à mettre en œuvre (surtout avec un langage tel que Julia qui facilite les manipulations de matrices), nous avons également vu ses limites en termes de temps d’exécution et d’empreinte mémoire.

Il est temps de passer à la FFT proprement dite !