Exemple simple de DFT

L’objectif de cet exercice est d’étudier numériquement la transformée de Fourier d’une Gaussienne : \(G(x) = \exp(-\frac{x^2}{2\sigma^2})\). Analytiquemet, la T.F. de G vaut \(\tilde G(\omega) = \sqrt{2\pi}\sigma \exp(-\frac{\omega^2\sigma^2}{2})\)

1. Écrire une fonction Python qui retourne une Gaussienne \(\exp(-\frac{x^2}{2\sigma^2})\) sur [-L/2, L/2] et son axe des \texttt{x = L/N*np.arange(N) - L/2}, prenant comme paramètres \(\sigma\), \(L\) et \(N\).

2. Tracer la Gaussienne pour \(\sigma=0.5\) et \(N=1000\) et \(L=5\).

3. Tracer la DFT de la Gausienne pour \(\sigma=0.15\) pour une plage de fréquences inférieures à 20 : (a) Quel taux d’échantillonnage \(\Delta x = L/N\) faut-il utiliser ? (b) Quelle valeur minimum de L faut-il choisir pour obtenir une résolution en fréquence de 0.04 ?

4. Comparer la courbe obtenue avec la formule analytique ? Pourquoi la DFT de la Gaussienne n’est pas possitive ?

5. Que se passe-t-il si \(\sigma=1.5\) ?

6. Appliquer la transformée de Fourier inverse et vérifier que l’on obient bien le signal initial.

import matplotlib.pyplot as plt
import numpy as np
from numpy.fft import fft, ifft, rfft, irfft, rfftfreq, fftfreq

def f(sigma,L,N):
    Delta_x = L/N
    x = Delta_x*np.arange(N)-L/2
    y = np.exp(-x**2/(2*sigma**2))
    return x, y

L=5
N=1000
sigma=0.1
x, y = f(sigma, L, N)
plt.plot(x, y)

[<matplotlib.lines.Line2D at 0x7f7cac6d7d00>]

L = 1/0.04
N = int(2*20*L) # f_max = dx/2
sigma=0.15
x, y = f(sigma, L, N)

f_tilde = fft(y)
f_axis = fftfreq(N, L/N)
plt.plot(f_axis, f_tilde*L/N)

# La TF présente des valeurs négatives car elle n'est pas centrée en j=0 mais en j=N/2
# La phase vaut exp(2 i pi N/2 k/N) donc \pm 1

/home/pierre/.local/lib/python3.8/site-packages/matplotlib/cbook/__init__.py:1335: ComplexWarning: Casting complex values to real discards the imaginary part
  return np.asarray(x, float)

[<matplotlib.lines.Line2D at 0x7f7c78bc6dc0>]

La TF présente des valeurs négatives car elle n’est pas centrée en \(j=0\) mais en \(j=N/2\). La phase vaut \(e^{2 i \pi \frac{N}{2} \frac{k}{N}}\) donc \((-1)^k\)

f_tilde = fft(y)*L/N
f_axis_dft = fftfreq(N, L/N)

plt.plot(f_axis_dft, np.abs(f_tilde), '.')

f_axis = np.linspace(-20, 20, 201)

plt.plot(f_axis, np.sqrt(2*np.pi)*sigma*np.exp(-(2*np.pi*f_axis)**2*sigma**2/2))
plt.xlim(-5, 5)

(-5.0, 5.0)

y_tilde = fft(y)

signal_initial = ifft(y_tilde)
plt.plot(x, signal_initial, 'o')
plt.plot(x, y)

[<matplotlib.lines.Line2D at 0x7f7c78ac3430>]