Décomposition en série de Fourier

Ce tutoriel montre comment synthétiser un signal à partir de sa décomposition en série de Fourier.

Home signal_processing .ipynb
In [7]:
%matplotlib inline
from numpy import *
from scipy.signal import *
from matplotlib.pyplot import *

Dans ce tutorial, nous allons apprendre à utiliser Python pour synthétiser un signal périodique à partir de sa décomposition en série de Fourier. Pour synthétiser le signal, nous allons utiliser la librairie Numpy. L'affichage du signal synthétisé sera ensuite réalisé via la librairie Matplotlib.

Décomposition en série de Fourier

Lors du siècle des lumières, J. Fourier s'interessa à la problématique de diffusion de la chaleur. Il publia notamment un ouvrage intitulé "Théorie analytique de la chaleur", rendu aujroud'hui immortel grace à la magie d'internet. Dans son étude, au lieu d'analyser les données experimentales en fonction du temps, notre Frenchy eu la très bonne idée de changer d'espace de répresentation. Cet espace de réprésentation présente l'avantage d'être parcimonieux c-a-d que l'énergie du signal (l'information utile) est concentrée sur un petit nombre de paramètres. Sur le plan philosophique, son approche est résumée dans le discours préliminaire de son ouvrage:

Les causes primordiales ne nous sont point connues; mais elles sont assujetties à des lois simples et constantes, que l'on peut découvrir par l'observation, et dont l'étude est l'objet de la philosophie naturelle.

En gros, même si les phénomènes naturels nous semblent compliquées au premier abord, gardons en tête que la nature fait les choses simplement ! Cette approche est intimmement liée au Rasoir d'Ockham "Pluralitas non est ponenda sine necessitate" que nous pouvons traduire, après de multiples approximations par, "les hypothèses suffisantes les plus simples doivent être préférées".

Sur le plan mathématique, J. Fourier avanca qu'un signal périodique $x(t)$ de période $T_0=\frac{1}{f_0}$ pouvait se décomposer en une somme de signaux simples: des sinusoîdes de fréquence $f_n=kf_0$. En utilisant des exponentielles complexes au lieu de sinusoides (merci Euler), sa proposition se résume par l'équation suivante:

$$x(t)=\sum_{n=-\infty}^{\infty} c_n e^{2j\pi nf_0 t}$$

où $n\in \mathbb{Z}$ et $c_n\in \mathbb{C}$. L'histoire d'aller plus loin, J. Fourier proposa également une technique pour obtenir les coefficients $c_n$ en fonction du signal $x(t)$:

$$c_n=\frac{1}{T_0}\int_{[T_0]} x(t) e^{-2j\pi nf_0 t}dt$$

Damned ! Sa proposition fit couler beaucoup d'encre à l'époque ,notamment de la part de plusieurs de ses contemporains comme Pierre-Simon de Laplace, Joseph-Louis Lagrange, et Siméon Denis Poisson (pas n'importe qui quand même !)...Et oui, Fourier ne peut tout de même pas nous faire croire qu'un signal présentant des discontinuités peut se décomposer en une somme de signaux sinusoidaux.

Cas du Signal Carré

Allure du signal de référence

Dans ce tutorial, nous allons nous focaliser sur la décomposition en série de Fourier d'un signal carré de fréquence $f_0=2$ Hz et d'amplitude crête $A=2$. Ce signal présente des discontinuité en $t=kT_0/2$.

In [8]:
Fe = 1000
f0 = 2
T0 = 1/f0
A = 2

t = arange(Fe)/Fe       # base temps allant de 0 à 1 seconde
x = A*square(2*pi*f0*t)   # génération du signal carré
plot(t,x)
xlabel("temps (s)")
ylabel("$x(t)$");

Décomposition en Série de Fourier

Dans un premier temps, nous devons déterminer la valeur des coefficient $c_n$ à partir du signal $x(t)$. En utilisant la proposition de J. Fourier, le coefficient $c_n$ peut s'obtenir de la manière suivante:

$$\begin{align} c_n&=\frac{1}{T_0}\int_{[T_0]} x(t) e^{-2j\pi n f_0 t} dt\\\\ &=\frac{A}{T_0}\left(\int_{0}^{\frac{T_0}{2}} e^{-2j\pi n f_0 t} dt -\int_{\frac{T_0}{2}}^{T_0} e^{-2j\pi n f_0 t} dt\right)\\\\ &=\frac{A}{2j\pi n}\left(1-e^{-j\pi n} - e^{-j\pi n} +1\right)\\\\ &=\frac{A}{j\pi n}\left(1-e^{-j\pi n}\right) \end{align}$$

En prenant plusieurs valeurs de $n$, nous remarquons rapidement que:

  • $c_n=\frac{2A}{j\pi n}$ lorsque $n$ est impair
  • $c_n=0$ lorsque $n$ est pair.

Nous dirons par la suite que le signal carré ne possède que des harmoniques de rang impair.

Remarquons que pour recontruire le signal $x(t)$ avec la formule de J. Fourier, nous devons sommer une infinité de termes. Toutefois, l'énergie du signal est portée essentiellement par les termes basse-fréquences (c-a-d pour $n$ petit). En limitant la reconstruction aux coefficients $c_n$ à un nombre fini de termes, nous obtenons:

In [9]:
L=9
xr = 0j*zeros(len(t))

for n in range(-L,L+1,2):   #attention, on utilise un pas de 2 pour ne conserver que les harmoniques de rang impair
    cn = 2*A/(1j*pi*n)
    xr += cn*exp(2j*pi*n*f0*t)
 
plot(t,x,label="reference")
plot(t,np.real(xr),label="reconstruction") # on affiche que la partie réelle....la partie imaginaire étant nulle.
xlabel("temps (s)")
legend();

Ca marche plutôt bien ! On pourrait bien sur affiner la reconstruction en ajoutant un plus grand nombre d'harmoniques c-a-d en augmentant la valeur de $L$.

Nous observons donc que le signal peut être synthétisé à partir d'un petit nombre de coefficients $c_n$. Au lieu de réprésentation le signal $x(t)$ en fonction du temps, on préfère alors changer de representation en affichant la valeur des coefficients $c_n$ de sa décomposition en série de Fourier. Comme ces coefficients sont complexes, on affichera à la fois leur module c-a-d $|c_n|$ et leur argument $\arg[c_n]$.

In [10]:
n_vect = arange(-L,L+1,2)
c = 2*A/(1j*pi*n_vect)

stem(n_vect,abs(c))
xlabel("n")
ylabel("$|c_n|$")
figure()
stem(n_vect,angle(c));
xlabel("n")
ylabel("arg[$c_n$]");