TangenteX.com

L'analyse de Fourier

C'est lorsqu'il résolut l'équation de diffusion thermique que Fourier introduisit, en 1807, les séries trigonométriques qui portent son nom, et qu'il construisit ainsi l'analyse harmonique, que l'on nomme analyse de Fourier en son honneur.

L'équation de diffusion thermique unidimensionnelle, appelée improprement aussi "équation de la chaleur" est de la forme :

\( \dfrac{\partial T}{\partial t} = D \dfrac{\partial^2 T}{\partial x^2} \)

Où T est une fonction de x et de t et D le coefficient de diffusivité thermique.

C'est une équation aux dérivées partielles linéaire, dont la solution générale s'obtient par combinaison linéaire de solutions particulières. Ces solutions particulières sont de la forme \( e^{u(x,t)} \), u étant une forme linéaire à déterminer de x et de t. la forme linéaire la plus simple est \( u(x,t) = \omega t + ikx \), forme qui vous est sans doute familière, avec la relation de dispersion \( \omega  = -Dk^2 \).

Dans son ouvrage "Théorie analytique de la chaleur", en 1822, Fourier montra que la fonction solution de l'équation de diffusion thermique pouvait s'écrire sous la forme d'une somme de cosinus :

\( \displaystyle T(x,t) = \sum_{n \in \mathbb{N}} p_n e^{-Dn^2q^2t} \cos(nqx) \)

avec \( \displaystyle \left( 0 \lt x \lt a, q = \dfrac{\pi}{a} \right) \) et  [0, a] la période de la fonction.

Les coefficients \( p_n \) sont les coefficients de Fourier, définis par :

\( \displaystyle  p_n = \dfrac{2}{a} \int_0^a \cos(nqx) \: T(x,0) \: dx  \)

Nous avons là une série de Fourier, objet mathématique que je vais utiliser dans la suite.

Une remarque avant d'aller plus loin : une hypothèse fondamentale de l'analyse de Fourier est que les équations différentielles soient linéaires. Et comme l'analyse de Fourier est un outil très puissant d'analyse des systèmes, vous comprendrez mieux que chaque fois que cela est possible, on fait les approximations nécessaires pour linéariser le système. Toutefois, il est possible d'utiliser l'analyse de Fourier sur des systèmes non-linéaires autour de leurs points fixes, là où il est possible de linéariser localement le système. Il peut aussi être intéressant d'appliquer l'analyse de Fourier sur des solutions connues d'équations différentielles non linéaires, mais cela dépasse de très loin notre scope.

Les séries de Fourier

Soit f une fonction réelle définie partout sur \( \mathbb{R} \), sauf éventuellement en un nombre fini de points. On dit que la fonction f est périodique s'il existe un nombre L tel que \( f(x + nL) = f(x) \)  avec \( \forall \: n \in \mathbb{Z} \). On appelle L la période de la fonction f et l'on dit que f est une fonction L-périodique.

Une remarque de physicien : dans notre monde, une fonction périodique est une vue de l'esprit mathématicien. Il n'existe pas de phénomène infini qui soit modélisable strictement par une fonction périodique. Cependant, car sinon aucune physique ne serait possible, on admet qu'une fonction est périodique sur un intervalle de temps ou d'espace donné lorsque sa période est très petite par rapport à l'intervalle de temps. Dans la pratique, nous ferons en sorte de toujours nous trouver dans cette situation.

Soit donc notre fonction f(x) L-périodique, que nous voulons décomposer en une somme de fonctions élémentaires, qui formeraient une base de la fonction f(x).

Nous pourrions décomposer la fonction f(x) en utilisant un développement limité, de Taylor par exemple. Nous obtiendrions une somme de polynômes, qui ne mettrait pas en avant le caractère périodique de f(x), ce qui serait sans grand intérêt.

la forme trigonométrique d'une série de Fourier

Une autre voie est plus riche : utiliser les fonctions périodiques simples que nous connaissons comme les sinus et les cosinus. Le problème revient donc à exprimer la fonction f(x) par une combinaison linéaire de sinus et de cosinus périodiques sur L. Leur argument sera donc de la forme \( \dfrac{2\pi nx}{L}\). Le terme \( \dfrac{1}{L} \) représente la fréquence de la fonction et le produit \(\dfrac{2\pi nx}{L} \) peut donc s'écrire \( n\omega x\) avec la pulsation \( \omega = 2 \pi \dfrac{1}{L} \). Nous pouvons écrire la combinaison linéaire suivante, qui constitue la définition d'une série de Fourier :

\( \displaystyle f(x) =  a_0 + \sum_{n=1}^{\infty} \left( a_n \cos(n\omega x) + b_n \sin(n\omega x) \right) \)       (1)

Pourquoi le terme \( a_0 \) est-il sorti de la somme ? Parce que cos(0) = 1 et que sin(0) = 0 et donc que pour n = 0, la somme se réduit à \( a_0 \).

Les termes négatifs n'ont pas été pris en compte. C'est dû aux propriétés de parité des sinus et cosinus que vous connaissez bien : sin(-x) = -sin(x) et cos(-x) = cos(x). On retrouve donc les termes négatifs dans l'expression des termes positifs.

Les coefficients \( a_0 \),  \( a_n \) et  \( b_n \) sont les coefficients de Fourier de la série, définis par :

\( \displaystyle a_0 = \dfrac{1}{L} \int_{x_0}^{x_0 + L} f(x) \: dx \)

\( \displaystyle a_n = \dfrac{2}{L} \int_{x_0}^{x_0 + L} \cos(n\omega x) \: f(x) \: dx \)

\( \displaystyle b_n = \dfrac{2}{L} \int_{x_0}^{x_0 + L} \sin(n\omega x) \: f(x) \: dx \)

Notons que si la fonction f(x) est paire, alors les coefficients  \( b_n \) seront nuls et que si la fonction f(x) est impaire, alors ce sont les coefficients \( a_0 \) et  \( a_n \) qui seront nuls.

Les fonctions sinus et cosinus employées dans la combinaison linéaire qui forme la série de Fourier sont orthogonales. Leur produit scalaire est nul, comme vous pourrez vous en convaincre en faisant quelques petits calculs trigo. Elles forment une base complète.

Sa forme exponentielle

L'écriture de la série de Fourier sous sa forme trigonométrique en (1) est lourde, sans compter que les calculs trigonométriques sont eux aussi assez indigestes. Heureusement, il y a la formule d'Euler, qui nous permet une écriture sous forme exponentielle beaucoup plus simple :

\( \displaystyle f(x) = \sum_{n=-\infty}^{\infty} c_n e^{i n\omega x} \)       (2)

avec les coefficients \( c_n \) :

\( \displaystyle c_n = \dfrac{1}{L} \int_{x_0}^{x_0 + L} e^{-i n\omega x} \: f(x) \: dx  \)

ce qui nous simplifie quand même pas mal la vie !

L'ensemble des coefficients \( c_n \) forment le spectre de Fourier de la fonction f(x).

La transformée de Fourier

Dans la définition d'un série de Fourier, la période L d'intégration est une grandeur finie. En faisant tendre L vers l'infini, nous obtenons l'expression d'une intégrale de Fourier :

\( \displaystyle f(x) = \dfrac{1}{2\pi} \int_{-\infty}^{\infty} C_k\: e^{ikx} \: dk \)      

avec :

\( \displaystyle C_k = \int_{-\infty}^{\infty} f(x) \: e^{-ikx} \: dx \)                          

k est un variable conjuguée, qui vaut \( k =  n \omega \).

Cette intégrale va permettre la définition de la transformée de Fourier, c'est à dire d'une opération qui transforme une fonction intégrable f en une fonction \( \widehat{f} \) qui décrit le spectre fréquentiel de f. La fonction f n'est pas nécessairement périodique. En effet, la transformée de Fourier permet d'étendre l'analyse harmonique aux fonctions de période infinie, c'est à dire aux fonctions non périodiques.

Il existe plusieurs définitions de la transformée de Fourier, qui différent par la présence ou non de coefficients constants et par les signes dans les exponentielles. Il s'agit de conventions, et il suffit d'en choisir une et de s'y tenir. Ici, je choisirai par habitude la définition des physiciens spécialistes des systèmes :

\( \displaystyle  \widehat{f}(\omega) = \dfrac{1}{\sqrt{2\pi}}  \int_{-\infty}^{\infty} f(t) \: e^{-i\omega t} \: dt \)        (3)

où \(\omega \) est la pulsation du signal et t le temps. Chez les physiciens, f(t) est appelée le "signal" et \( \widehat{f}(\omega) \) le "spectre" du signal. Autre point de terminologie : on dit que f est définie dans le domaine "temporel" et \( \widehat{f} \) dans le domaine fréquentiel, ou pour simplifier temps et fréquence.

Dans la littérature, on trouve aussi fréquement l'expression de la transformée de Fourier dont la variable est la fréquence \( \nu \) :

\( \displaystyle \widehat{f}(\nu) = \int_{-\infty}^{\infty} f(t) \: e^{-i2\pi \nu t} \: dt \)

Je n'expose pas ici les principales propriétés de la transformée de Fourier, vous les trouverez dans tout bon bouquin d'analyse. A ce sujet, je vous suggère la lecture du cours d'analyse de Jean-Michel Bony et l'excellent manuel "Des mathématiques pour les sciences" de Claude Aslangul (un pavé de 1250 pages pleines de trésors).

La transformée de Fourier discrète

La transformée de Fourier opère sur des fonctions intégrables, qui sont au moins continues par morceaux sur le domaine d'intégration. C'est bien en mathématique, mais cela n'arrange pas vraiment nos affaires de physiciens numériciens. Nous ne travaillons pratiquement jamais sur des fonctions intégrables mais plutôt sur des suites de nombres qui résultent de l'échantillonnage par un instrument de mesure d'un signal continu.

Autre remarque sur l'utilisation pratique d'une transformée de Fourier : ses bornes d'intégration sont l'infini. C'est à dire qu'il faudrait considérer que le signal que nous voulons analyser, la fonction f(t), soit défini et intégrable sur le domaine \([-\infty,+\infty ]\). Gênant s'agissant d'un signal issu d'une acquisition par un instrument de mesure !

Nous avons deux limitations dans l'utilisation de la transformée de Fourier en physique. La transformée de Fourier est inopérante et il faut se tourner vers sa version discrète, la célèbre Transformée de Fourier Discrète ou TFD.

Sa définition

Supposons donc que nous disposions d'un signal, défini par la fonction f(t), qui est acquis entre 0 et L secondes. Implicitement, je considère donc que la fonction f(t) est L-périodique, sa périodicité étant construite par recopie du morceau de signal entre 0 et L.

Lors de son acquisition, ce signal a été échantillonné en N+1 points séparés par un intervalle de temps \( \Delta t \). Le signal est donc découpé en N intervalles de temps. Bien sur, j'ai pris la précaution de vérifier que la fréquence d'échantillonnage \( \dfrac{1}{\Delta t} \) est bien supérieure d'au moins deux fois à la fréquence de mon signal, qui par construction est égal à \( \dfrac{1}{L} \), pour respecter la condition de Shannon-Nyquist. J'obtiens donc une suite de nombres représentant la valeur acquise du signal \( f(t_n) \) pour n variant de 0 à N. Dans la suite, je poserai \( y_n = f(t_n) \).

Pour calculer les coefficients de Fourier de notre signal sur l'intervalle [0,L], je vais intégrer numériquement l'équation (3) en chacun des points \( t_n \), en fixant les bornes d'intégration à 0 à L. Classiquement, on utilise une intégration par la méthode des trapèzes.

Tous calculs faits, nous obtenons la TFD du signal, qui est l'expression des coefficients de Fourier discrets :

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

J'ai choisi de garder le facteur \(\dfrac{1}{\sqrt{2\pi}} \) dans l'expression le TFD, ce qui relève de la convention comme je l'ai dit plus haut. Vous trouverez souvent la définition de la TFD sans ce facteur (dans Wikipedia par exemple) mais cela ne gêne en rien la nature de la TFD.

Les coefficients de Fourier obtenus sont aussi les coefficients du polynôme d'interpolation trigonométriques, interpolation réalisée sur les points \(( t_n, y_n) \).

Sa signification physique

La transformée de Fourier, qu'elle soit continue ou discrète, décompose la fonction f(t), notre signal, en une somme de fonctions trigonométriques. Plus les termes de la somme sont nombreux, et plus la somme de ces fonctions est proche de la fonction initiale f(t). Nous avons vu ce principe dans une autre page de TangenteX.

Dans cette décomposition, les coefficients de Fourier \( c_k \) mesurent, par leur valeur relative, l'importance de la fonction de fréquence correspondante au coefficient. On obtient ainsi le spectre fréquentiel de notre signal. Pour le dire autrement, notre signal est une somme d'ondes de fréquence donnée, et la TFD nous permet d'identifier ces ondes et indique leur point relatif dans la somme.

Algorithmes de TFD bruts

Un algorithme "naïf"

Le code de la TFD

Le calcul numérique de la TFD est assez trivial. Le code de la fonction TFD est une implémentation directe de l'équation :

def TFD(u):
    N = len(u)
    tfd = zeros(N, complex)
    for k in range(N):
        for n in range(N-1):
            tfd[k] += u[n]*cmat.exp(-2j*cmat.pi*k*n/N)
    return tfd

L'algorithme comporte deux boucles imbriquées, l'une sur l'indice k fréquentiel et l'autre sur l'indice n temporel. La boucle interne d'indice n calcule la valeur du coefficient de Fourier pour une fréquence (ou pulsation) donnée et la boucle externe d'indice k provoque le calcul du coefficient de Fourier pour chaque pas temporel.

J'ai laissé tombé le facteur \( \dfrac{1}{\sqrt{2\pi}} \), pour ne pas alourdir le calcul. Comme je le remarquais plus haut, cela ne nuit en rien au principe du calcul de la TFD, mais les amplitudes du spectre ne seront pas normalisées, ce qui n'est pas important ici.

Une remarque concernant la complexité de cet algorithme "brut". Pour un signal échantillonné sur N points, le calcul de la TFD nécessite \( (N - 1)^2 \) multiplications complexes et \( N(N - 1) \) additions complexes, sans parler des \( N(N - 1) \) calculs d'exponentielles. On dit que sa complexité est de l'ordre de N2. Pour un calcul d'une TFD sur 1024 points, très courant, cela représente plus d'un million d'opérations, à supposer que l'on ait tabulé les valeurs des exponentielles ! Heureusement, nous verrons dans une autre page un algorithme, le célèbre FFT, qui permet de résoudre ce très sévère inconvénient lorsque l'on veut calculer une TFD en temps réel, au sens propre du terme, comme c'est souvent le cas.

Calculer la TFD d'un signal simple

Soit le signal \( s(t) = a_0 \sin(2 \pi f_0 t) \), d'une durée L et échantillonné sur N points. Soit le code suivant pour définir ce signal :

# 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 = 20              # nombre de points d'échantillonnage
dt = L/N            # intervalle temporel entre deux points
fe = 1./dt          # fréquence d'échantillonnage

puis le construire :

t = linspace(0.0,L,N)
s = a0*sin(2*pi*f0*t)

Je calcule la TFD de ce signal en appelant la fonction TFD :

tfd = TFD(s)

Les coefficients de Fourier discrets obtenus sont :

0 (1.44328993201e-15+0j)
1 (0.0928311859668-0.586113041007j)
2 (2.95082833029-9.08171577316j)
3 (-0.628708537442+1.23390998048j)
4 (-0.441980587317+0.608334089582j)
5 (-0.389165620936+0.389165620936j)
6 (-0.366074482172+0.265968679715j)
7 (-0.354135972365+0.180441290501j)
8 (-0.347625697966+0.112950436186j)
9 (-0.344318401248+0.0545346772751j)
10 (-0.343300433624-1.79210106224e-15j)
11 (-0.344318401248-0.0545346772751j)
12 (-0.347625697966-0.112950436186j)
13 (-0.354135972365-0.180441290501j)
14 (-0.366074482172-0.265968679715j)
15 (-0.389165620936-0.389165620936j)
16 (-0.441980587317-0.608334089582j)
17 (-0.628708537442-1.23390998048j)
18 (2.95082833029+9.08171577316j)
19 (0.0928311859668+0.586113041007j)

Pour chaque ligne, le premier nombre indique le rang k du coefficient et le second nombre indique la valeur complexe du coefficient. Voici une représentation schématique de leur répartition, en tracant le module de chaque coefficient :

Coefficients de Fourier

L'examen des valeurs de ces coefficients appellent deux remarques :

  • Le coefficient d'ordre 0 (k=0) est quasi-nul. Il représente la valeur moyenne du signal, qui dans le cas de la fonction sinus est nulle. Bien souvent, on élimine ce terme.
  • La liste des coefficients complexes présente une symétrie intéressant, qui est dûe au fait que le signal traité est une fonction réelle. Vous noterez que tfd[k] = tfd*[N-k], son complexe conjugué. C'est une observation importante. Lorsque le signal est une fonction réelle, ce qui est souvent le cas, il n'est pas utile de faire le calcul sur N points, un calcul sur N/2 points suffit.

Le code suivant me permet de tracer le spectre fréquentiel du signal :

# domaine fréquentiel pour affichage du spectre
freq = arange(N)*fd
# tracé du spectre
fig2 = figure(figsize=(10,8))
grid(True)
plot(freq[0:N//2],abs(tfd[0:N//2]),'or')
title('Spectre du signal')
xlabel(u'Fréquence')
ylabel('Amplitude')

ce qui me donne :

Spectre TFD

Comme attendu, le pic du spectre indique la fréquence de base de notre signal, 1 Hz. Quelques remarques à propos du code :

  • J'ai défini le domaine fréquentiel à partir de la durée d'acquisition du signal, la période L, plus exactement à partir de son inverse, c'est à dire de sa fréquence.
  • Pour tracer le spectre, j'ai utilisé la première moitié des données, suite à la remarque faite ci-dessus à propos de la symétrie constatée pour une fonction réelle.
  • Théoriquement, les coefficients de Fourier différents du coefficient correspondant à la fréquence de 1 Hz devraient être nuls, car mon signal est un sinus pur. Ils ne le sont pas parce que notre intégration par la méthode des trapèzes est imprécise. Mathématiquement pourtant, ces coefficients de TFD devraient être exactement égaux aux coefficients de la transformée de Fourier.

Le script TFDBrut.py est téléchargeable dans la bibliothèque de codes de TangenteX.

Un algorithme matriciel

Le script TFDCalculBrutVectorise.py propose une autre méthode de calcul d'une FFT, une méthode matricielle basée sur la définition de la TFD.

Pour cet exemple, j'ai choisi un signal composé de la somme de deux sinus, d'amplitude et de fréquence différentes, tel que :

s = a0*sin(2.0*pi*f0*t) + a1*sin(2.0*pi*f1*t)

En voici le code :

# définition du signal à traiter
Fe = 40.0 # fréquence d'échantillonnage en Hz
dF = 0.1 # résolution fréquentielle en Hz

# le nombre d'échantillons est déterminé par la résolution fréquentielle
# hypothèse N = L
N = int(Fe/dF) # nombre d'échantillons

# la durée du signal à analyser (la période d'échantillonnage Te) dépend du
# nombre d'échantillons à acquérir et de la fréquence d'échantillonnage
Te = 1.0/dF # durée du signal à échantillonner

# construction du signal échantillonné s
t = arange(0.0,Te,Te/N) # vecteur temps discrétisé
a0 = 5.0; f0 = 10.0 # amplitude et fréquence de la sinusoïde
a1 = 2.0; f1 = 5.0
s = a0*sin(2.0*pi*f0*t) + a1*sin(2.0*pi*f1*t)

# construction du vecteur fréquence
f = arange(N)*dF

# définition de la matrice W
n = arange(N)
k = n.reshape((N,1))
W = exp(-2j*pi*k*n/N)

# calcul de la transformée de Fourier discrète par simple produit
# matriciel X = M.s
X = dot(W,s)

# calcul du spectre SP normalisé
SP = (2.0/N)*abs(X)
# pour affichage en dB
SP1 = 20*log10(SP/SP.max())

# affichage du spectre
fig = plt.figure(figsize=(14,8))
plt.grid(True)
plt.title(u"Spectre de fréquence d'un signal composite")
plt.xlabel(u'Fréquence', fontsize = 16)
plt.ylabel('Amplitude (en dB)', fontsize = 16)
plt.plot(f[1:N/2],SP1[1:N/2])

Le tracé du spectre obtenu :

Spectre du signal - FFT python matriciel

On observe très clairement, l'affichage est en dB, les pics des fréquences f0= 5 Hz et f1 = 10 Hz dans le spectre. On observe aussi bien le bruit de "calcul" dans le reste du spectre.

Le script TFDCalculBrutVectorise.py est téléchargeable dans la bibliothèque de codes de TangenteX.

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