Transformation de Fourier, FFT et DFT#
Introduction à la FFT et à la DFT#
La Transformée de Fourier Rapide, appelée FFT Fast Fourier Transform en anglais, est un algorithme qui permet de calculer des Transformées de Fourier Discrètes DFT Discrete Fourier Transform en anglais.
Parce que la DFT permet de déterminer la pondération entre différentes fréquences discrètes, elle a un grand nombre d’applications en traitement du signal, par exemple pour du filtrage. Par conséquent, les données discrètes qu’elle prend en entrée sont souvent appelées signal et dans ce cas on considère qu’elles sont définies dans le domaine temporel. Les valeurs de sortie sont alors appelées le spectre et sont définies dans le domaine des fréquences. Toutefois, ce n’est pas toujours le cas et cela dépend des données à traiter.
Il existe plusieurs façons de définir la DFT, en particulier au niveau du signe que l’on met dans l’exponentielle et dans la façon de normaliser. Dans le cas de NumPy, l’implémentation de la DFT est la suivante :
\(A_k=\sum\limits_{m=0}^{n-1}{a_m\exp\left\{ -2\pi i\frac{mk}{n} \right\}}\text{ avec }k=0,\ldots ,n-1\)
La DFT inverse est donnée par :
\(a_m=\frac{1}{n}\sum\limits_{k=0}^{n-1}{A_k\exp\left\{ 2\pi i\frac{mk}{n} \right\}}\text{ avec }m=0,\ldots ,n-1\)
Elle diffère de la transformée directe par le signe de l’argument de l’exponentielle et par la normalisation à 1/n par défaut.
Exemples simples#
Visualisation de la partie réelle et imaginaire de la transformée#
import numpy as np
import matplotlib.pyplot as plt
n = 20
# definition de a
a = np.zeros(n)
a[1] = 1
# visualisation de a
# on ajoute a droite la valeur de gauche pour la periodicite
plt.subplot(311)
plt.plot( np.append(a, a[0]) )
# calcul de A
A = np.fft.fft(a)
# visualisation de A
# on ajoute a droite la valeur de gauche pour la periodicite
B = np.append(A, A[0])
plt.subplot(312)
plt.plot(np.real(B))
plt.ylabel("partie reelle")
plt.subplot(313)
plt.plot(np.imag(B))
plt.ylabel("partie imaginaire")
plt.show()
Visualisation des valeurs complexes avec une échelle colorée#
Pour plus d’informations sur cette technique de visualisation, voir Visualisation d’une fonction à valeurs complexes avec PyLab.
import numpy as np
import matplotlib.pyplot as plt
n = 20
# definition de a
a = np.zeros(n)
a[1] = 1
# visualisation de a
# on ajoute a droite la valeur de gauche pour la periodicite
plt.subplot(211)
plt.plot( np.append(a, a[0]) )
# calcul de k
k = np.arange(n)
# calcul de A
A = np.fft.fft(a)
# visualisation de A - Attention au changement de variable
# on ajoute a droite la valeur de gauche pour la periodicite
plt.subplot(212)
x = np.append(k, k[-1]+k[1]-k[0]) # calcul d'une valeur supplementaire
z = np.append(A, A[0])
X = np.array([x,x])
y0 = np.zeros(len(x))
y = np.abs(z)
Y = np.array([y0,y])
Z = np.array([z,z])
C = np.angle(Z)
plt.plot(x,y,'k')
plt.pcolormesh(X, Y, C, shading="gouraud", cmap=plt.cm.hsv, vmin=-np.pi, vmax=np.pi)
plt.colorbar()
plt.show()
Exemple avec a[2]=1#
Exemple avec a[0]=1#
Exemple avec cosinus#
import numpy as np
import matplotlib.pyplot as plt
n = 20
# definition de a
m = np.arange(n)
a = np.cos(m * 2*np.pi/n)
# visualisation de a
# on ajoute a droite la valeur de gauche pour la periodicite
plt.subplot(311)
plt.plot( np.append(a, a[0]) )
# calcul de A
A = np.fft.fft(a)
# visualisation de A
# on ajoute a droite la valeur de gauche pour la periodicite
B = np.append(A, A[0])
plt.subplot(312)
plt.plot(np.real(B))
plt.ylabel("partie reelle")
plt.subplot(313)
plt.plot(np.imag(B))
plt.ylabel("partie imaginaire")
plt.show()
Exemple avec sinus#
Exemple avec cosinus sans prise en compte de la période dans l’affichage
import numpy as np
import matplotlib.pyplot as plt
n = 20
# definition de a
m = np.arange(n)
a = np.cos(m * 2*np.pi/n)
# visualisation de a
plt.subplot(211)
plt.plot(a)
# calcul de A
A = np.fft.fft(a)
# visualisation de A
plt.subplot(212)
plt.plot(np.real(A))
plt.ylabel("partie reelle")
plt.show()
Fonction fftfreq#
numpy.fft.fftfreq renvoie les fréquences du signal calculé dans la DFT.
Le tableau freq renvoyé contient les fréquences discrètes en nombre de cycles par pas de temps. Par exemple si le pas de temps est en secondes, alors les fréquences seront données en cycles/seconde.
Si le signal contient n pas de temps et que le pas de temps vaut d :
freq = [0, 1, …, n/2-1, -n/2, …, -1] / (d*n) si n est pair
freq = [0, 1, …, (n-1)/2, -(n-1)/2, …, -1] / (d*n) si n est impair
import numpy as np
import matplotlib.pyplot as plt
# definition du signal
dt = 0.1
T1 = 2
T2 = 5
t = np.arange(0, T1*T2, dt)
signal = 2*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 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()
Fonction fftshift#
>>> n = 8
>>> dt = 0.1
>>> freq = np.fft.fftfreq(n, d=dt)
>>> freq
array([ 0. , 1.25, 2.5 , 3.75, -5. , -3.75, -2.5 , -1.25])
>>> f = np.fft.fftshift(freq)
>>> f
array([-5. , -3.75, -2.5 , -1.25, 0. , 1.25, 2.5 , 3.75])
>>> inv_f = np.fft.ifftshift(f)
>>> inv_f
array([ 0. , 1.25, 2.5 , 3.75, -5. , -3.75, -2.5 , -1.25])
Approximation de la transformée de Fourier grâce à la FFT#
Lorsqu’on désire calculer la transformée de Fourier d’une fonction \(x(t)\) à l’aide d’un ordinateur, ce dernier ne travaille que sur des valeurs discrètes, on est amené à :
discrétiser la fonction temporelle,
tronquer la fonction temporelle,
discrétiser la fonction fréquentielle.
\(X(f) = \int\limits_{-\infty}^{+\infty} x(t) e^{-2\pi i f t}dt\)
En approchant l’intégrale par une somme d’aires de rectangles de durée \(\Delta t\) et en limitant la durée d’intégration à un intervalle de longueur \(n \Delta t\), c’est-à-dire en faisant une somme pour \(n\) valeurs \(t_m = m \Delta t\) avec \(m\) allant de \(0\) à \(n-1\), on obtient :
\(X(f) \approx \Delta t \sum\limits_{m=0}^{n-1}{ x(t_m) e^{-2\pi i f \, t_m}}\)
Or, si on choisit des fréquences discrètes telles que \(f_k = k \frac{1}{n \Delta t}\), on a \(f_k t_m = k \frac{1}{n \Delta t} m \Delta t = \frac{mk}{n}\):
\(X(f_k) \approx \Delta t \sum\limits_{m=0}^{n-1}{ x(t_m) e^{-2\pi i f_k t_m} } \approx \Delta t \sum\limits_{m=0}^{n-1}{ x(t_m) e^{-2\pi i \frac{m k}{n}}}\approx \Delta t \,\text{fft}(x)\)
Exemple dans le cas d’une gaussienne#
import numpy as np
import matplotlib.pyplot as plt
alpha = 10
nc = 40
dt = 0.1
tmax = (nc-1) * dt
tmin = -nc * dt
# definition d'un signal gaussien
t = np.linspace(tmin, tmax, 2*nc)
x = np.exp(-alpha * t**2)
plt.subplot(411)
plt.plot(t,x)
# on effectue un ifftshift pour positionner le temps zero comme premier element
plt.subplot(412)
a = np.fft.ifftshift(x)
plt.plot(a)
A = np.fft.fft(a)
# on effectue un fftshift pour positionner la frequence zero au centre
X = dt*np.fft.fftshift(A)
# calcul des frequences avec fftfreq
n = t.size
freq = np.fft.fftfreq(n, d=dt)
f = np.fft.fftshift(freq)
# comparaison avec la solution exacte
plt.subplot(413)
plt.plot(f, np.real(X), label="fft")
plt.plot(f, np.sqrt(np.pi/alpha) * np.exp( -(np.pi*f)**2 / alpha), label="exact")
plt.legend()
plt.subplot(414)
plt.plot(f, np.imag(X))
plt.show()
Pour vérifier notre calcul, nous avons utilisé une transformée de Fourier connue. En effet, pour la définition utilisée, la transformée de Fourier d’une gaussienne \(e^{-\alpha t^2}\) est donnée par :
\(\sqrt{\frac{\pi}{\alpha}}e^{-\frac{(\pi f)^2}{\alpha}}\)
Exemple avec visualisation en couleur de la transformée de Fourier#
import numpy as np
import matplotlib.pyplot as plt
alpha = 10
nc = 40
dt = 0.1
tmax = (nc-1) * dt
tmin = -nc * dt
# definition d'un signal gaussien
t = np.linspace(tmin, tmax, 2*nc)
x = np.exp(-alpha * t**2)
plt.subplot(211)
plt.plot(t,x)
# on effectue un ifftshift pour positionner le temps zero comme premier element
a = np.fft.ifftshift(x)
A = np.fft.fft(a)
# on effectue un fftshift pour positionner la frequence zero au centre
X = dt*np.fft.fftshift(A)
# calcul des frequences avec fftfreq
n = t.size
freq = np.fft.fftfreq(n, d=dt)
f = np.fft.fftshift(freq)
# visualisation de X - Attention au changement de variable
# on ajoute a droite la valeur de gauche pour la periodicite
plt.subplot(212)
x = np.append(f, f[0]) # calcul d'une valeur supplementaire
z = np.append(X, X[0])
X = np.array([x,x])
y0 = np.zeros(len(x))
y = np.abs(z)
Y = np.array([y0,y])
Z = np.array([z,z])
C = np.angle(Z)
plt.plot(x,y,'k')
plt.pcolormesh(X, Y, C, shading="gouraud", cmap=plt.cm.hsv, vmin=-np.pi, vmax=np.pi)
plt.colorbar()
plt.show()
Exemple avec translation#
import numpy as np
import matplotlib.pyplot as plt
alpha = 10
nc = 40
dt = 0.1
tmax = (nc-1) * dt
tmin = -nc * dt
# definition d'un signal gaussien
t = np.linspace(tmin, tmax, 2*nc)
x = np.exp(-alpha * (t-1)**2)
plt.subplot(211)
plt.plot(t,x)
# on effectue un ifftshift pour positionner le temps zero comme premier element
a = np.fft.ifftshift(x)
A = np.fft.fft(a)
# on effectue un fftshift pour positionner la frequence zero au centre
X = dt*np.fft.fftshift(A)
# calcul des frequences avec fftfreq
n = t.size
freq = np.fft.fftfreq(n, d=dt)
f = np.fft.fftshift(freq)
# visualisation de X - Attention au changement de variable
# on ajoute a droite la valeur de gauche pour la periodicite
plt.subplot(212)
x = np.append(f, f[0]) # calcul d'une valeur supplementaire
z = np.append(X, X[0])
X = np.array([x,x])
y0 = np.zeros(len(x))
y = np.abs(z)
Y = np.array([y0,y])
Z = np.array([z,z])
C = np.angle(Z)
plt.plot(x,y,'k')
plt.pcolormesh(X, Y, C, shading="gouraud", cmap=plt.cm.hsv, vmin=-np.pi, vmax=np.pi)
plt.colorbar()
plt.show()
