Sommaire
Ce tutoriel est en cours de rédaction / approbation.
Ce tutoriel s’appuie sur la sous-bibliothèque fft de la bibliothèque SciPy.
Vous trouverez des rappels sur la transformée de Fourier et sa version discrétisée sur la page suivante : Transformée de Fourier & FFT
Transformée de Fourier
Depuis la sous-bibliothèque fft de la bibliothèque SciPy, il est possible d’accéder à des algorithmes de calcul de la transformée de Fourier discrète (et son inverse).
from scipy.fft import fft, ifft, fftshift
Calcul de la FFT d’un vecteur de données
Le calcul de la transformée de Fourier discrète d’un vecteur de données, par l’intermédiaire de l’algorithme FFT, se fait par la fonction fft de la sous-bibliothèque fft de SciPy.
x = np.array([1.0, 2.0, 1.0, -1.0, 1.5]) y = fft(x) print(y)
Dans l’exemple précédent, le vecteur y contient alors la transformée de Fourier de la suite de nombres du vecteur x.
[ 4.5 -0.j 2.08155948-1.65109876j -1.83155948+1.60822041j -1.83155948-1.60822041j 2.08155948+1.65109876j]
Il s’agit d’un vecteur de nombres complexes.
Calcul de la transformée de Fourier inverse
Le calcul de la transformée de Fourier inverse discrète d’un vecteur de données, par l’intermédiaire de l’algorithme FFT, se fait par la fonction ifft de la sous-bibliothèque fft de SciPy.
yinv = ifft(y) print(yinv)
Cet exemple donne l’affichage suivant :
[ 1.0+0.j, 2.0+0.j, 1.0+0.j, -1.0+0.j, 1.5+0.j]
Affichage de la FFT
Afin de pouvoir afficher de manière graphique les résultats des calculs de transformées de Fourier, nous allons utiliser la sous-bibliothèque pyplot de la bibliothèque Matplotlib (voir aussi tutoriel : Python / Matplotlib / Premiers pas ).
Pour cela, il faut importer ce module par l’intermédiaire de la commande suivante :
import matplotlib.pyplot as plt
Parties réelles et imaginaires
Dans cet exemple, nous nous intéresserons à un signal dont toutes les valeurs sont à 0 sauf une qui vaudra 1. Pour réaliser ce signal, on utilisera le code suivant :
n = 100 # sur 100 points a = np.zeros(n) a[1] = 1
Affichage du signal initial
plt.figure() plt.plot(a, '*') plt.show()
Ce qui donne l’affichage suivant :
Calcul et visualisation de la FFT
On calcule ensuite la transformée de Fourier de ce signal par la fonction suivante :
A = np.fft.fft(a)
Ce calcul fournit alors un vecteur complexe de même taille que le vecteur de données initial. Il est alors possible de visualiser les valeurs réelles et imaginaires de ce vecteur comme suit :
plt.figure()
plt.subplot(211)
plt.plot(np.real(B))
plt.ylabel("partie reelle")
plt.subplot(212)
plt.plot(np.imag(B))
plt.ylabel("partie imaginaire")
plt.show()
Cet exemple donne l’affichage suivant :
Module et phase
Il est courant également d’afficher le résultat d’une transformée de Fourier à travers le module et la phase du nombre complexe. C’est le cas par exemple pour l’affichage d’un diagramme de Bode, réponse en fréquence d’un système.
On prendra dans cet exemple un signal sinusoïdal (cosinus), sur un nombre de période p. avec un nombre d’échantillons n.
n = 100 p = 5 m = np.linspace(0, 2*p*np.pi-2*np.pi/p, n) a = np.cos(m)
Affichage du signal initial
L’affichage peut se faire par l’intermédiaire du code suivant :
plt.figure()
plt.plot(m, a)
plt.xlabel('Angle (rd)')
plt.ylabel('Cosinus')
Calcul et visualisation de la FFT
On calcule ensuite la transformée de Fourier de ce signal par la fonction suivante :
A = fft(a)
Ce calcul fournit alors un vecteur complexe de même taille que le vecteur de données initial. Il est alors possible de visualiser le module et la phase de ce vecteur comme suit :
plt.figure()
plt.subplot(211)
plt.plot(np.abs(A), '.')
plt.ylabel('Module')
plt.subplot(212)
plt.plot(np.angle(A), '.')
plt.ylabel('Phase (rd)')
Ce qui donne l’affichage suivant :
La transformée de Fourier d’un signal sinusoïdal de fréquence \(f_0\) donne deux pics de Dirac en \(-f_0\) et \(f_0\). On remarque bien ces deux pics sur le résultat précédent (partie module).
Un peu de théorie… et fftshift
La transformée de Fourier d’un signal échantillonné, comme c’est le cas du vecteur de données que l’on fournit à la fonction fft, est périodisée autour de tous les multiple entiers de la fréquence d’échantillonnage \(F_e\).
L’algorithme de calcul utilisé par la fonction fft fournit alors des résultats dans l’intervalle fréquentiel \( [0 , F_e [ \). Or, en physique, nous avons plus l’habitude de visualiser les réponses en fréquence centrées autour de la fréquence nulle. Il existe pour cela une fonction fftshift qui permet de recaler la partie droite du vecteur du côté gauche (périodisation) et ainsi obtenir l’affichage habituel.
Ainsi, le code suivant :
A = fft(a)
As = fftshift(A)
plt.figure()
plt.subplot(211)
plt.plot(np.abs(As), '.')
plt.ylabel('Module')
plt.subplot(212)
plt.plot(np.angle(As), '.')
plt.ylabel('Phase (rd)')
Donne l’affichage suivant :
On remarque aussi que ces pics semblent convolués avec un autre signal… Peut-être un sinus cardinal, signifiant que le signal a été multiplié par une porte…
Reconstruction des axes fréquentiels
A venir…
FFTFREQ
numpy.fft.fftfreq renvoie les fréquences du signal calculé dans la transformée de Fourier discrète
import numpy as np import matplotlib.pyplot as plt
definition du signal
dt = 0.1 T1 = 2 T2 = 5 t = np.arange(0, T1T2, dt) signal = 2np.cos(2np.pi/T1t) + np.sin(2np.pi/T2t)
affichage du signal
plt.subplot(211) plt.plot(t,signal)
calcul de la transformee de Fourier et des frequences
fourier = np.fft.fft(signal) n = signal.size freq = np.fft.fftfreq(n, d=dt)
affichage de la transformee de Fourier
plt.subplot(212) plt.plot(freq, fourier.real, label="real") plt.plot(freq, fourier.imag, label="imag") plt.legend() plt.show()
Ce tutoriel a été co-écrit par Aboubakar BAYERO (promo 2024) dans le cadre d’une semaine spécifique.