Approximation de la transformée de Fourier spatiale grâce à la FFT

Présentation

Pour une introduction à la FFT en Python, vous pouvez consulter la page Transformation de Fourier, FFT et DFT.

Nous allons considérer l’expression suivante pour la transformée de Fourier :

\(g(k) = \frac{1}{\sqrt{2 \pi}}\int\limits_{-\infty}^{+\infty} f(x) e^{- i k x}dx\)

Lorsqu’on désire calculer la transformée de Fourier de la fonction \(f(x)\) à l’aide d’un ordinateur, ce dernier ne travaille que sur des valeurs discrètes, on est amené à :

  • discrétiser la fonction spatiale,
  • tronquer la fonction spatiale,
  • discrétiser la fonction des fréquences spatiales.

En approchant l’intégrale par une somme d’aires de rectangles de largeur \(\Delta x\) et en limitant la domaine d’intégration à un intervalle de longueur \(n \Delta x\), c’est-à-dire en faisant une somme pour \(n\) valeurs \(x_m = m \Delta x\) avec \(m\) allant de \(0\) à \(n-1\), on obtient :

\(g(k) \approx \frac{\Delta x}{\sqrt{2 \pi}}\sum\limits_{m=0}^{n-1}{ f(x_m) e^{-i k x_m}}\)

Or, si on choisit des fréquences spatiales discrètes telles que \(k_j = j \frac{2\pi}{n \Delta x}\), on a \(k_j x_m = j \frac{2\pi}{n \Delta x} \, m \Delta x = 2\pi\frac{j m}{n}\):

\(g(k_j) \approx \frac{\Delta x}{\sqrt{2 \pi}} \sum\limits_{m=0}^{n-1}{ f(x_m) e^{-i k_j x_m} } \approx \frac{\Delta x}{\sqrt{2 \pi}} \sum\limits_{m=0}^{n-1}{ f(x_m) e^{-2\pi i \frac{j m}{n}}}\approx \frac{\Delta x}{\sqrt{2 \pi}} \,\text{fft}(f)\)

Exemple dans le cas d’une gaussienne

import numpy as np
import matplotlib.pyplot as plt

alpha = 20

nc = 40
dx = 0.1
xmax = (nc-1) * dx
xmin = -nc * dx

# definition d'un signal gaussien
x = np.linspace(xmin, xmax, 2*nc)
f = np.exp(-alpha/2 * x**2)

plt.subplot(411)
plt.plot(x,f)

# on effectue un ifftshift pour positionner le temps zero comme premier element
plt.subplot(412)
a = np.fft.ifftshift(f)
plt.plot(a)

A = np.fft.fft(a)
# on effectue un fftshift pour positionner la frequence zero au centre
g = dx/np.sqrt(2*np.pi)*np.fft.fftshift(A)

# calcul des frequences avec fftfreq
n = x.size
freq = 2*np.pi * np.fft.fftfreq(n, d=dx)
k = np.fft.fftshift(freq)

# comparaison avec la solution exacte
plt.subplot(413)
plt.plot(k, np.real(g), label="fft")
plt.plot(k, 1/np.sqrt(alpha) * np.exp( -k**2 / (2*alpha)), label="exact")
plt.legend()

plt.subplot(414)
plt.plot(k, np.imag(g))

plt.show()

(Source code)

_images/spatial_fft_2015a.png

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 \frac{x^2}{2}}\) est donnée par :

\(\frac{1}{\sqrt{\alpha}}e^{-\frac{k^2}{2\alpha}}\)

Exemple avec visualisation en couleur de la transformée de Fourier

import numpy as np
import matplotlib.pyplot as plt

alpha = 20

nc = 40
dx = 0.1
xmax = (nc-1) * dx
xmin = -nc * dx

# definition d'un signal gaussien
x = np.linspace(xmin, xmax, 2*nc)
f = np.exp(-alpha/2 * x**2)

plt.subplot(211)
plt.plot(x,f)

# on effectue un ifftshift pour positionner le temps zero comme premier element
a = np.fft.ifftshift(f)

A = np.fft.fft(a)
# on effectue un fftshift pour positionner la frequence zero au centre
g = dx/np.sqrt(2*np.pi)*np.fft.fftshift(A)

# calcul des frequences avec fftfreq
n = x.size
freq = 2*np.pi * np.fft.fftfreq(n, d=dx)
k = 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(k, k[0]) # calcul d'une valeur supplementaire
z = np.append(g, g[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()

(Source code)

_images/spatial_fft_2015b.png

Exemple pour le produit d’un cosinus par une gaussienne

import numpy as np
import matplotlib.pyplot as plt

alpha = 3
k1 = 10

nc = 40
dx = 0.1
xmax = (nc-1) * dx
xmin = -nc * dx

# definition d'un signal gaussien
x = np.linspace(xmin, xmax, 2*nc)
f = np.cos(k1*x)*np.exp(-alpha/2 * x**2)

plt.subplot(211)
plt.plot(x,f)

# on effectue un ifftshift pour positionner le temps zero comme premier element
a = np.fft.ifftshift(f)

A = np.fft.fft(a)
# on effectue un fftshift pour positionner la frequence zero au centre
g = dx/np.sqrt(2*np.pi)*np.fft.fftshift(A)

# calcul des frequences avec fftfreq
n = x.size
freq = 2*np.pi * np.fft.fftfreq(n, d=dx)
k = 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(k, k[0]) # calcul d'une valeur supplementaire
z = np.append(g, g[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()

(Source code)

_images/spatial_fft_2015c.png