TP - Echantillonnage et Tranformée de Fourier Discrète

I) Théorème d'échantillonage

On va considérer le signal à valeurs complexes $$ f(t) = e^{-2 i \pi t}, $$ dont la transformée de Fourier est $$ \mathcal{F}f(s) = \delta(s-1). $$ On a donc 1 fréquence qui est $1$.

1) Illustrer la fonction les parties réelle et imaginaire de $f$

Importons les librairies dont nous aurons besoin pour faire des graphiques et du calcul numérique.

import matplotlib.pyplot as plt
import numpy as np

Définissons la fonction $f$.

f = lambda t: np.exp(-2*1j*np.pi*t)

puis une grille de temps.

T = 2.
t = np.linspace(-T, T, 200)

La partie réelle du signal est donnée par

plt.plot(t,np.real(f(t)))
plt.show()

plt.plot(t,np.imag(f(t)))
plt.show()

2) Créer une discrétisation $(t_n)$ de l'intervalle de temps $[-2,2]$ de pas de temps $1/(2B)$ (noté t_grid), où $B$ est un paramètre avec lequel on jouera. Illustrer l'échantillonnage de $f$ sur cette grille.

On pourra tester des valeurs de $B$ autour de $1$.

B = 1.5
dt = 1/(2*B)
t_grid = np.arange(-T, T+dt, dt)

Evaluons $f_a$ sur la grille t_grid

f_discrete = f(t_grid)

Illustrons le signal $f_a$ et ses points d'échantillonnage

plt.plot(t_grid, np.real(f_discrete), "o")
plt.plot(t, np.real(f(t)))
plt.show()

plt.plot(t_grid, np.imag(f_discrete), "o")
plt.plot(t, np.imag(f(t)))
plt.show()

3) Evaluer la fonction suivante sur la grille t_grid et la comparer avec la fonction $f$ originale:

$$ f_{sample}(t) = \sum_{n=0}^{N-1} f(t_n) \, \text{sinc}\Big(\frac{(t-t_n)}{dt}\Big)$$

pour différentes valeurs de $B$.

Définissons la fonction $f_{sample}$.

f_sample = lambda t: np.sum( f_discrete * np.sinc( (t-t_grid)/dt ) )

et évaluons la sur la grille fine t définie plus haut pour la comparer à $f_a$.

f_reconstruceted  = np.array([f_sample(t) for t in t])

plt.plot(t, np.real(f_reconstruceted))
plt.plot(t, np.real(f(t)))
plt.show()

plt.plot(t, np.imag(f_reconstruceted))
plt.plot(t, np.imag(f(t)))
plt.show()

On peut voir que la reconstruction est correcte pour $B$ plus grand que $1$, et si on prend $B$ suffisamment grand, on recontruit le signal très bien.

4) Recommencer l'analyse précédente mais avec le signal

$$ f(t) = e^{-2\pi i t} + e^{-6\pi i t}. $$

Que peut-on observer sur $B$ pour avoir une bonne reconstruction à partir de $f_{sample}$?

f = lambda t: np.exp(-2*1j*np.pi*t)+np.exp(-6*1j*np.pi*t)

B = 3.5
dt = 1/(2*B)
t_grid = np.arange(-T, T+dt, dt)

f_discrete = f(t_grid)

plt.plot(t_grid, np.real(f_discrete), "o")
plt.plot(t, np.real(f(t)))
plt.show()

plt.plot(t_grid, np.imag(f_discrete), "o")
plt.plot(t, np.imag(f(t)))
plt.show()

f_sample = lambda t: np.sum( f_discrete * np.sinc( (t-t_grid)/dt ) )

f_reconstruceted  = np.array([f_sample(t) for t in t])

plt.plot(t, np.real(f_reconstruceted))
plt.plot(t, np.real(f(t)))
plt.show()

plt.plot(t, np.imag(f_reconstruceted))
plt.plot(t, np.imag(f(t)))
plt.show()

II) Transformée de Fourier discrète

On considère le signal

$$ f(t) = e^{-\pi(t-2)^2}e^{2i\pi t}. $$

La fréquence maximal supportée par ce signal est autour de $2$.

5) Illustrer ce signal graphiquement sur l'intervalle $[0,4]$.

f = lambda t: np.exp(- np.pi*(t-2)**2) * np.exp(2*1j*np.pi*t)

t = np.linspace(0, 4, 200)

plt.plot(t,np.real(f(t)))
plt.show()

plt.plot(t,np.imag(f(t)))
plt.show()

6) Evaluer la transformée de Fourier de $f$ sur $[0,2B]$ avec une grille de pas de temps $1/T$ avec $T=4$ qui est en gros la durée du signal.

from scipy.integrate import quad # on importe la fonction permettant le calcul d'intégrale

T = 4
B = 2
freq = np.arange(0, 2*B, 1/T)

def T_Fourier(s):
    
    f_re_tmp = lambda t: np.exp(-np.pi * (t-2)**2) * np.cos(-2*np.pi*(s-1)*t)   
    tf_re, _ = quad(f_re_tmp, -10, 10)
    
    f_im_tmp = lambda t: np.exp(- np.pi* (t-2)**2) * np.sin(-2*np.pi*(s-1)*t)
    tf_im, _ = quad(f_im_tmp, -10, 10)
    
    return np.array([tf_re, tf_im]) 

TF = np.array([T_Fourier(s) for s in freq])

plt.plot(TF[:,0], "o")

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

plt.plot(TF[:,1], "o")
plt.ylim(-1.5,1.5)

(-1.5, 1.5)

7) Evaluer $f$ sur l'intervalle $[0,T]$ à l'aide d'une grille de pas de temps $1/(2B)$ avec $B=2$. On obtiendra un vecteur f_grid. Comparer ce qu'on vient d'obtenir avec ce que retourne la fonction fft du module numpy.fft. On appliquera fft(f_grid).

from numpy.fft import fft

dt = 1/(2*B)
t_grid = np.arange(0, T, dt)
f_grid = f(t_grid)

TFD_f_grid_np = fft(f_grid)

plt.plot(np.real(TFD_f_grid_np)*dt, "o")
plt.plot(TF[:,0] ,"o")
plt.ylim(-1.5,1.5)
plt.show()

plt.plot(np.imag(TFD_f_grid_np)*dt, "o")
plt.plot(TF[:,1] ,"o")
plt.ylim(-1.5,1.5)
plt.show()

On a besoin de rajouter un facteur $dt$ pour que le résultat de notre transformée de Fourier soit en adéquation avec ce que nous renvoie la function fft.