Sommaire
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]
Il existe des versions “réelles” de ces algorithmes de FFT (rfft et irfft) qui s’appuie sur le fait que les signaux d’entrée sont réels. Certaines propriétés des transformées de Fourier sont alors incluses dans le calcul.
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 = 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(A))
plt.ylabel("partie reelle")
plt.subplot(212)
plt.plot(np.imag(A))
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 de l’axe des fréquences
Il est possible de reconstruire l’axe des fréquences en connaissant le nombre de points utilisés pour réaliser la FFT ainsi que la fréquence d’échantillonnage (ou période d’échantillonnage).
Lorsqu’on réalise l’opération suivante :
A = fft(a, n)
A[0] contient le terme de fréquence nulle.
Selon que \(\)n[\latex] est pair ou impair, les autres termes correspondent à :
Then
From the Numpy.FFT documentationA[1:n/2]contains the positive-frequency terms, andA[n/2+1:]contains the negative-frequency terms, in order of decreasingly negative frequency. For an even number of input points,A[n/2]represents both positive and negative Nyquist frequency, and is also purely real for real input. For an odd number of input points,A[(n-1)/2]contains the largest positive frequency, whileA[(n+1)/2]contains the largest negative frequency. The routinenp.fft.fftfreq(n)returns an array giving the frequencies of corresponding elements in the output. The routinenp.fft.fftshift(A)shifts transforms and their frequencies to put the zero-frequency components in the middle, andnp.fft.ifftshift(A)undoes that shift.
Une méthode possible est alors la suivante (basée sur la fonction fftfreq de scipy.fft) :
def freq_reconstruct(N, Fe=1.0):
if N%2 == 0: # even number (pair)
freq1 = np.linspace(0, (N-2)/(2*N), N//2) * Fe
freq2 = np.linspace(-N/(2*N), -2/(2*N), N//2) * Fe
freq = np.concatenate((freq2, freq1))
print('OK')
else: # odd number (impair)
freq1 = np.linspace(0, (N-1)/(2*N), N//2+1) * Fe
freq2 = np.linspace(-(N-1)/(2*N), -2/(2*N), N//2) * Fe
freq = np.concatenate((freq2, freq1))
return freq
FFTFREQ
scipy.fft.fftfreq renvoie les fréquences du signal calculé dans la transformée de Fourier discrète. Cette fonction reconstruit alors l’axe des fréquences selon les règles précédentes.
import scipy import numpy as np import matplotlib.pyplot as plt dt = 0.1 T1 = 2 T2 = 5 t = np.arange(0, T1*T2, dt) signal = np.cos(2*np.pi/T1*t) + np.sin(2*np.pi/T2*t)
Affichage du signal
plt.subplot(211) plt.plot(t,signal)
Calcul de la transformée de Fourier et des fréquences
fourier = scipy.fft.fft(signal) n = signal.size freq = scipy.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()
Cette méthode affiche alors les fréquences dans le “bon” sens, car le vecteur freq commence par des valeurs positives.
Afin d’optimiser l’affichage, il est possible d’utiliser la fonction fftshift sur les deux vecteurs : celui contenant la TF et celui contenant les fréquences.
freq = scipy.fft.fftshif(freq) fourier = scipy.fft.fftshif(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.