TP3 - Tranformée de Fourier¶

I) Théorème d'échantillonage¶

Pour illustrer comment échantillonner un signal continue, on va considérer le signal à valeurs complexes suivant

$$ f_a(t) = e^{-\pi (t-2)^2}e^{2\pi i a(t-2)} = h_a(t-2)\qquad a = 1. $$

Ce signal a pour fréquence centrale $a=1$ et, comme on a un décalage, sa transformée de Fourier est simplement

$$ \mathcal{F}f_a(s) = \mathcal{F}h_a(s) e^{2 \pi i b s}\qquad b=-2, $$

où

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

Comme

$$ \mathcal{F}h_a(s) = e^{-\pi(s-a)^2}, $$

on a au final

$$ \mathcal{F}f_a(s) = e^{-\pi(s-a)^2}e^{2 \pi i b s}. $$

1) Illustrer la fonction $f_a$¶

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

In [1]:
import matplotlib.pyplot as plt
import numpy as np

Définissons la fonction $f_a$.

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

puis une grille de temps.

In [28]:
T = 4.
t = np.linspace(0, 4, 200)

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

In [4]:
plt.plot(t,np.real(f(t)))
plt.show()
No description has been provided for this image

et la partie imaginaire par

In [5]:
plt.plot(t,np.imag(f(t)))
plt.show()
No description has been provided for this image

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

On va prendre $B$ comme étant 2 fois la fréquence centrale pour définir notre grille de temps.

In [210]:
B = 2*a
dt = 1/(2*B)
t_grid = np.arange(0, 4, dt)

Evaluons $f_a$ sur la grille t_grid

In [211]:
f_discrete = f(t_grid)

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

In [212]:
plt.plot(t_grid, np.real(f_discrete), "o")
plt.plot(t, np.real(f(t)))
plt.show()
No description has been provided for this image
In [213]:
plt.plot(t_grid, np.imag(f_discrete), "o")
plt.plot(t, np.imag(f(t)))
plt.show()
No description has been provided for this image

3) Evaluer la fonction suivante sur la grille t_grid et la comparer avec la fonction $f_a$ originale.¶

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

pour différentes valeurs de $B$.¶

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

In [10]:
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$.

In [11]:
f_reconstruceted  = np.array([f_sample(t) for t in t])
In [12]:
plt.plot(t, np.real(f_reconstruceted))
plt.plot(t, np.real(f(t)))
plt.show()
No description has been provided for this image
In [13]:
plt.plot(t, np.imag(f_reconstruceted))
plt.plot(t, np.imag(f(t)))
plt.show()
No description has been provided for this image

On peut voir que la reconstruction est très bonne. Elle devient moins bonne si on prend $B$ plus petit.

4) Tracer la transformée de Fourier de $f_a$ et tracer les droites verticales d'abscisse $B$. Que peut-on remarquer, et est ce en accord avec l'hypothèse principale du théorème d'échantillonnage?¶

Evaluons la transformée de Fourier de $f_a$ sur une grille fine.

In [14]:
b = -2
s = np.linspace(-1, 2.5, 100)
TF_fa = np.exp(-np.pi*(s-a)**2) * np.exp(2*np.pi*1j*b*s)

Sa partie réelle est donnée par

In [15]:
plt.plot(s, np.real(TF_fa), label='TF')
plt.vlines(0, -0.5, 0.5, label=' freq 0', color='green')
plt.vlines(B, -0.5, 0.5, label=' freq B', color='orange')
plt.legend()
plt.show()
No description has been provided for this image

et sa partie imaginaire par

In [16]:
plt.plot(s, np.imag(TF_fa), label='TF')
plt.vlines(0, -0.5, 0.5, label=' freq 0', color='green')
plt.vlines(B, -0.5, 0.5, label=' freq B', color='orange')
plt.legend()
plt.show()
No description has been provided for this image

Tout d'abord on n'est pas dans les conditions standards du théorème d'échantillonnage car $\mathcal{F}f_a$ n'est pas limité en bande de fréquence (band limited). Cependant, la fréquence centrale est donnée par $a$ dans le terme $e^{2\pi i a t}$. Le théorème d'échantillonnage nous dirait d'échantillonner à $2a = 2$ (2 fois la fréquence max). Avec ce choix, on observe une mauvaise reconstruction du signal. En un sens cette observation est normale car dans notre cas $a$ n'est pas la fréquence maximale portée par le signal à cause de la gaussienne. Le choix $B=2a$ correspond quasiment à une fréquence maximale et donne une reconstruction quasi partfaite du signal original.

II) Transformée de Fourier discrète¶

5) Evaluer la tranformée de Fourier discrète de $f_a$ discrétisée sur la grille $(t_n)$ et la comparer graphiquement avec la formule théorique.¶

Définissons une grille de fréquence couvrant $[0,2B]$ de pas de temps $1/T$ et une fonction correspondant à la TFD qu'on va évaluer sur la grille.

In [29]:
s = np.arange(0, 2*B, 1/T)
TFD_f = lambda s: np.sum( f_discrete * np.exp(-2*np.pi*1j*s*t_grid) )
In [18]:
TFD_f_discretise = np.array([TFD_f(s) for s in s])  
TFD_f_discretise
Out[18]:
array([ 1.72852192e-01+1.11022302e-16j, -6.83282922e-01+0.00000000e+00j,
        1.82374885e+00+2.22044605e-16j, -3.28690355e+00-1.38777878e-15j,
        3.99999626e+00+2.24987512e-15j, -3.28690355e+00-1.88737914e-15j,
        1.82374885e+00+1.38777878e-15j, -6.83282922e-01-5.55111512e-16j,
        1.72852192e-01+1.11022302e-16j, -2.95305664e-02+5.55111512e-16j,
        3.40245297e-03-5.55111512e-16j, -2.69009259e-04-2.77555756e-17j,
        2.46532238e-05+1.60383061e-15j, -2.69009259e-04-1.66533454e-16j,
        3.40245297e-03+1.11022302e-16j, -2.95305664e-02-1.11022302e-16j])

Comparons la tranformée de Fourier théorique et ce qu'on obtient avec la TFD. Pour la partie réelle on a

In [19]:
TF_fa = np.exp(-np.pi*(s-a)**2) * np.exp(2*np.pi*1j*b*s)
plt.plot(s, np.real(TFD_f_discretise)*dt, "o")
plt.plot(s, np.real(TF_fa) ,"o")
plt.ylim(-1.5, 1.5)
plt.show()
No description has been provided for this image

et pour la partie imaginaire

In [20]:
plt.plot(s, np.imag(TFD_f_discretise)*dt, "o")
plt.plot(s, np.imag(TF_fa) ,"o")
plt.ylim(-1.5, 1.5)
plt.show()
No description has been provided for this image

On remarquera qu'on a besoin de multiplier par $dt$ pour comparer avec la formulation théorique de la transformée de Fourier.¶

6) Effectuer le même calcul sous forme de produit matriciel (on evitera de faire une boucle for).¶

Le calcul de la DFT fait intervenir les termes

$$ e^{-2 \pi i n m / N} = (e^{-2 \pi i/ N})^{n m} = \omega^{-nm}$$

qui forment la matrice $\underline{\mathcal{F}}$. On va donc créer la matrice des indices $nm$ et l'évaluer par la fonction définie par l'exponentielle complexe.

In [21]:
N = len(t_grid)
m = np.arange(N).reshape((1,N))
n = m.reshape((N, 1))
(n*m).shape
Out[21]:
(16, 16)

On applique cette matrice à l'exponentielle complexe.

In [22]:
M = np.exp(-2 * np.pi * 1j * n * m / N)

La DFT est alors donnée par le produit matrice/vecteur $Mf_{discrete}$

In [23]:
TF_fa_mat = np.dot(M, f_discrete)
TF_fa_mat
Out[23]:
array([ 1.72852192e-01+3.08066025e-17j, -6.83282922e-01-6.88252096e-17j,
        1.82374885e+00+2.82255477e-16j, -3.28690355e+00-1.48981385e-15j,
        3.99999626e+00+2.24987512e-15j, -3.28690355e+00-1.86191077e-15j,
        1.82374885e+00+1.31241985e-15j, -6.83282922e-01-5.19964975e-16j,
        1.72852192e-01+1.38539704e-16j, -2.95305664e-02+5.66374091e-16j,
        3.40245297e-03-5.11594956e-16j, -2.69009259e-04+2.01875289e-17j,
        2.46532238e-05+1.60383061e-15j, -2.69009259e-04-2.68258218e-16j,
        3.40245297e-03+1.85767674e-16j, -2.95305664e-02-1.87189095e-17j])
In [24]:
plt.plot(s, np.real(TF_fa_mat)*dt, "o")
plt.plot(s, np.real(TF_fa), "o")
plt.ylim(-1.5, 1.5)
plt.show()
No description has been provided for this image
In [25]:
plt.plot(s, np.imag(TF_fa_mat)*dt, "o")
plt.plot(s, np.imag(TF_fa), "o")
plt.ylim(-1.5, 1.5)
plt.show()
No description has been provided for this image

7) Effectuer le même calcul à l'aide de la function fft du module numpy.fft.¶

FFT correspond à Fast Fourier Transform et est un algorithme pour calculer le transformée de Fourier discrète.

In [26]:
from numpy.fft import fft
In [27]:
TFD_f_discretise_np = fft(f_discrete)
TFD_f_discretise_np
Out[27]:
array([ 1.72852192e-01+1.70830506e-21j, -6.83282922e-01+1.70830506e-21j,
        1.82374885e+00+1.70830506e-21j, -3.28690355e+00+1.70830506e-21j,
        3.99999626e+00+1.70830506e-21j, -3.28690355e+00+1.70830506e-21j,
        1.82374885e+00+1.70830506e-21j, -6.83282922e-01+1.70830506e-21j,
        1.72852192e-01+1.70830506e-21j, -2.95305664e-02+1.70830506e-21j,
        3.40245297e-03+1.70830506e-21j, -2.69009259e-04+1.70830506e-21j,
        2.46532238e-05+1.70830506e-21j, -2.69009259e-04+1.70830506e-21j,
        3.40245297e-03+1.70830506e-21j, -2.95305664e-02+1.70830506e-21j])
In [28]:
plt.plot(s, np.real(TFD_f_discretise_np)*dt, "o")
plt.plot(s, np.real(TF_fa) ,"o")
plt.ylim(-1.5,1.5)
plt.show()
No description has been provided for this image
In [29]:
plt.plot(s, np.imag(TFD_f_discretise_np)*dt, "o")
plt.plot(s, np.imag(TF_fa) ,"o")
plt.ylim(-1.5,1.5)
plt.show()
No description has been provided for this image

8) Quel est la méthode la plus efficace en terme de temps calcul pour calculer la DFT? On pourra ecrire chacune des méthodes sous forme de fonction, et évaluer le temps de calcul de chacune des méthodes sur des entrées (signaux discrétisés tirés aléatoirement).¶

Les fonctions dont on va évaluer le temps d'execution sont les suivantes:

In [59]:
def TFD_1(x, t_grid):
    return np.array([np.sum( x * np.exp(-2*np.pi*1j*s*t_grid)) for s in s])
In [60]:
def TFD_2(x, M):
    m = np.arange(N).reshape((1,N))
    n = m.reshape((N, 1))
    M = np.exp(-2 * np.pi * 1j * n * m / N)
    return np.dot(M, x)
In [51]:
def TFD_3(x):
    return fft(x)

Pour mesurer le temps d'execution on va utiliser la librairie time.

In [52]:
import time

On va mesurer le temps qu'il faut pour chaque fonction de faire 100000 évaluations, comprenant la génération aléatoire des signaux discérisés (mais qui est la même pour chacune des fonctions).

In [61]:
N_repeat = 100000
time_exec = 0.

for j in range(N_repeat):
    x = np.random.rand(N)
    start_time = time.perf_counter()
    TFD_1(x, t_grid)
    stop_time = time.perf_counter()
    time_exec += stop_time - start_time
In [62]:
time_exec
Out[62]:
10.62068264796335
In [63]:
N_repeat = 100000
time_exec = 0.

for j in range(N_repeat):
    x = np.random.rand(N)
    start_time = time.perf_counter()
    TFD_2(x, M)
    stop_time = time.perf_counter()
    time_exec += stop_time - start_time
In [64]:
time_exec
Out[64]:
1.8500823179479085
In [65]:
N_repeat = 100000
time_exec = 0.

for j in range(N_repeat):
    x = np.random.rand(N)
    start_time = time.perf_counter()
    TFD_3(x)
    stop_time = time.perf_counter()
    time_exec += stop_time - start_time
In [66]:
time_exec
Out[66]:
0.19817016194065218

Sans surprise la fonction fft de numpy est la plus performant pour calculer la DFT. Cependant la formulation matricielle, évitant la boucle for explicite de la première méthode améliore déjà grandement les choses.