TP1 - Tranformée de Fourier¶
Premières illustrations numériques¶
On va s'intéresser à la fonction porte:
$$ \Pi(t) = \begin{cases} 1 & \text{ si } |t| \leq T/2 \\ 0 & \text{ sinon.} \end{cases} $$
1) Illustrer numériquement la fonction porte.¶
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 prote $\Pi$.
Pi = lambda t, T: 1.*(np.abs(t) < T/2)
# On utilise ici lambda pour déclarer la fonction
# par exemple
# f = lambda x: x**2
# f est alors la fonction qui a x associe x**2
Définissons l'intervalle de temps qu'on va considérer.
t = np.linspace(-10, 10, 200)
T = 4
plt.plot(t,Pi(t,T))
plt.show()
2) Calucler sa transformée de Fourier $\mathcal{F}(\Pi)$.¶
On a pour $s\in \mathbb{R}$ \begin{align*} \mathcal{F}\Pi(s) = \int_{-\infty}^{+\infty} \Pi(t) e^{-2 \pi i s t}dt & = \int_{-T/2}^{T/2} \Pi(t) e^{-2\pi i s t} dt \\ & = \int_{-T/2}^{T/2} e^{-2\pi i s t} dt \\ & = \Big[-\frac{e^{-2\pi i s t}}{2\pi i s}\Big]_{-T/2}^{T/2} \\ & = \frac{-(e^{-\pi i s T}-e^{\pi i s T})}{2\pi i s} \\ & = \frac{1}{\pi s}\frac{e^{\pi i s T}-e^{-\pi i s T}}{2 i} \\ & = \frac{\sin(\pi s T)}{\pi s} \\ & = T \frac{\sin(\pi s T)}{\pi s T} \\ & = T\text{sinc}(sT), \end{align*} où $\text{sinc}(x) = sin(\pi x) /(\pi x)$.
3) Comparer la formule théorique avec une évaluation numérique de la transformée de Fourier. On utilisera pour le moment la fonction quad de la librairie scipy.integrate (On verra comment faire autrement plus tard, l'approche qu'on va utliser n'est pas optimale!).¶
Commençons par importer la fonction quad et $\pi$
from scipy.integrate import quad
from numpy import pi as pi
Un exemple avec la fonction quad. Pour calculer l'integrale de $f$ entre $a$ et $b$ on appele quad(f,a,b). Par exemple pour calculer
$$ \int_{-1}^1 x^2 dx $$
on tape
f = lambda x: x**2
quad(f, -1, 1)
(0.6666666666666666, 7.401486830834376e-15)
le premier chiffre correspond à la valeur de l'integrale, et le second à une estimation de l'erreur absolue. C'est donc le premier chiffre qui va nous intéresser. Prour juste avoir ce terme, on peut procéder ainsi
I, _ = quad(f, -1, 1)
I
0.6666666666666666
Dans notre cas $$ \mathcal{F}\Pi(s) = \int_{-\infty}^{+\infty} \Pi(t) e^{2\pi ist} dt = \int_{-\infty}^{+\infty} \Pi(t) \cos(2\pi ist) dt + i \int_{-\infty}^{+\infty} \Pi(t) \sin(2\pi ist) dt $$ et on va calculer les parties réelles et imaginaires de la transformée de Fourier.
Définissons l'intervalle de fréquences.
s = np.linspace(-2, 2, 200)
Définissons maintenant notre fonction de transformée de Fourier qui prend en entrée une certain fonction et renvoit les parties réeelles et imagniare de la transformée de Fourier.
def TF_func(fct, s):
fct_tmp_re = lambda t: fct(t) * np.cos(2*pi*s*t)
re_tf, _ = quad(fct_tmp_re, -10, 10)
fct_tmp_im = lambda t: fct(t) * -np.sin(2*pi*s*t)
im_tf, _ = quad(fct_tmp_im, -10, 10)
return np.array([re_tf, im_tf])
On définit la fonction qu'on souhaite intégrer pour un $T$ donné.
Pi_tmp = lambda t: Pi(t,T)
On calul la transformée de Fourier en tous les points $s$
TF_Pi = np.array([TF_func(Pi_tmp, s) for s in s])
/tmp/ipykernel_15231/2371038607.py:4: IntegrationWarning: The integral is probably divergent, or slowly convergent. re_tf, _ = quad(fct_tmp_re, -10, 10)
Ci-dessous on voit que la partie imaginaire est bien nulle car notre fonction porte est paire.
TF_Pi[::10,:]
array([[-4.75942143e-14, 0.00000000e+00],
[-1.02185217e-01, 0.00000000e+00],
[ 1.87829609e-01, 0.00000000e+00],
[-2.19214515e-01, 0.00000000e+00],
[ 1.67112298e-01, 0.00000000e+00],
[-2.01885966e-02, 0.00000000e+00],
[-2.10418106e-01, 0.00000000e+00],
[ 4.93897423e-01, 0.00000000e+00],
[-7.93724657e-01, 0.00000000e+00],
[ 1.12643685e+00, 0.00000000e+00],
[ 3.98937483e+00, 0.00000000e+00],
[ 7.08982983e-01, 0.00000000e+00],
[-6.90215307e-01, 0.00000000e+00],
[ 5.13382114e-01, 0.00000000e+00],
[-2.81887914e-01, 0.00000000e+00],
[ 5.90512453e-02, 0.00000000e+00],
[ 1.08222350e-01, 0.00000000e+00],
[-1.93935460e-01, 0.00000000e+00],
[ 1.95978387e-01, 0.00000000e+00],
[-1.33549852e-01, 0.00000000e+00]])
plt.plot(s, TF_Pi[:,0])
plt.plot(s, T*np.sinc(s*T))
plt.show()
La formule théorique est notre simulation sont bien en adéquation.
4) Jouer avec le paramètre $T$. Qu'observe-t-on sur le signal porte et sa transformée de Fourier quand on joue avec ce paramètre.¶
T = 1
T1 = 8
plt.plot(t,Pi(t,T))
plt.plot(t,Pi(t,T1))
plt.show()
plt.plot(s, T*np.sinc(s*T))
plt.plot(s, T1*np.sinc(s*T1))
plt.show()
Plus le signal porte est étroit plus sont spectre est étalé, et inversement. Cela reflète le fait qu'un signal ne pas être localisé en espace et en temps.
5) Evaluer numériquement $\mathcal{F}^{-1} \mathcal{F} \Pi$. Obtient-on ce à quoi on s'attend?¶
On va calculer la transformée de Fourier inverse de $\mathcal{F} \Pi(s) = T\text{sinc}(sT)$. Comme $\mathcal{F}^{-1} \mathcal{F} \Pi = \Pi$, on s'attend à retourver la fonction $\Pi$.
Pour simplifier l'implémentation on rappelle que $\mathcal{F}^{-1} h = \mathcal{F} h^-$. Mais si $h$ est une fonction paire, on a $\mathcal{F}^{-1} h = \mathcal{F} h^- = \mathcal{F} h$, ce qui est le cas de $\text{sinc}$.
Ainsi, on a
$$ \mathcal{F}^{-1} \mathcal{F} \Pi = \mathcal{F}^{-1}(T\text{sinc}(\cdot T)) = T\mathcal{F}^{-1}(\text{sinc}(\cdot T)) = T\mathcal{F}(\text{sinc}(\cdot T)).$$
T = 4
t = np.linspace(-10,10,200)
TF_Pi_tmp = lambda s: T*np.sinc(s*T)
TF_TF_Pi = np.array([TF_func(TF_Pi_tmp, t) for t in t])
/tmp/ipykernel_15231/2371038607.py:4: IntegrationWarning: The maximum number of subdivisions (50) has been achieved. If increasing the limit yields no improvement it is advised to analyze the integrand in order to determine the difficulties. If the position of a local difficulty can be determined (singularity, discontinuity) one will probably gain from splitting up the interval and calling the integrator on the subranges. Perhaps a special-purpose integrator should be used. re_tf, _ = quad(fct_tmp_re, -10, 10) /tmp/ipykernel_15231/2371038607.py:4: IntegrationWarning: The integral is probably divergent, or slowly convergent. re_tf, _ = quad(fct_tmp_re, -10, 10)
TF_TF_Pi[::10,:]
array([[ 3.01393912e-02, 0.00000000e+00],
[-1.50820363e-01, 0.00000000e+00],
[ 6.44224136e-03, 0.00000000e+00],
[ 3.22946833e-04, 0.00000000e+00],
[ 1.89680408e-04, 0.00000000e+00],
[-1.51597256e-05, 0.00000000e+00],
[-5.65440109e-04, 0.00000000e+00],
[-2.58929823e-03, 0.00000000e+00],
[ 1.06909743e+00, 0.00000000e+00],
[ 1.00624714e+00, 0.00000000e+00],
[ 1.00506795e+00, 0.00000000e+00],
[ 1.00660882e+00, 0.00000000e+00],
[-7.04235134e-02, 0.00000000e+00],
[-2.19894109e-03, 0.00000000e+00],
[-4.80079364e-04, 0.00000000e+00],
[-3.33267162e-02, 0.00000000e+00],
[ 4.29873762e-02, 0.00000000e+00],
[ 2.65517359e-04, 0.00000000e+00],
[ 9.45843247e-03, 0.00000000e+00],
[-1.43934819e-01, 0.00000000e+00]])
La partie imaginaire de la transformée de Fourier que l'on vient de calculer est encore nulle car $\mathcal{F}\Pi$ est encore une fonction pare.
plt.plot(t, Pi(t,T))
plt.plot(t, TF_TF_Pi[:,0])
plt.show()
On s'attend à obtenir $\Pi$, car on a la relation $\mathcal{F}^{-1} \mathcal{F} \Pi = \Pi$, mais on voit ici que la reconstruction de la porte n'est pas top. Il y a des parasites, et la porte n'est pas bien capturée sur les bords.
6) Pour la fonction triangle du cours, comparer la formule théorique avec une évaluation numérique de la transformée de Fourier. Evaluer numériquement $\mathcal{F}^{-1} \mathcal{F} \Pi$. Obtient-on ce à quoi on s'attend?¶
On rappelle la fonction triangle, $$ \Lambda(t) = \begin{cases} 1 -|t| & \text{ si } |t| \leq 1 \\ 0 & \text{ sinon,} \end{cases} $$ dont la transformée de Fourier est donnée par $$ \mathcal{F}\Lambda(s) = \text{sinc}^2(s). $$
Lambda = lambda t: (1. - np.abs(t))*(np.abs(t) < 1)
s = np.linspace(-3, 3, 200)
TF_Lambda = np.array([TF_func(Lambda, s) for s in s])
/tmp/ipykernel_15231/2371038607.py:4: IntegrationWarning: The integral is probably divergent, or slowly convergent. re_tf, _ = quad(fct_tmp_re, -10, 10)
TF_Lambda[::10,:]
array([[2.99193757e-15, 0.00000000e+00],
[9.16952950e-03, 0.00000000e+00],
[1.58513366e-02, 0.00000000e+00],
[2.01451833e-03, 0.00000000e+00],
[1.14475344e-02, 0.00000000e+00],
[4.54621480e-02, 0.00000000e+00],
[2.27684000e-02, 0.00000000e+00],
[1.48376714e-02, 0.00000000e+00],
[2.71304178e-01, 0.00000000e+00],
[7.57602796e-01, 0.00000000e+00],
[9.99252545e-01, 0.00000000e+00],
[7.10820027e-01, 0.00000000e+00],
[2.30355352e-01, 0.00000000e+00],
[7.48310420e-03, 0.00000000e+00],
[2.78409717e-02, 0.00000000e+00],
[4.34838005e-02, 0.00000000e+00],
[8.38816838e-03, 0.00000000e+00],
[3.31536648e-03, 0.00000000e+00],
[1.63137438e-02, 0.00000000e+00],
[7.71450293e-03, 0.00000000e+00]])
plt.plot(s, TF_Lambda[:,0])
plt.plot(s, np.sinc(s)**2)
plt.show()
L'adéquation entre la formule théorique et l'approximation numérique est très bonne.
t = np.linspace(-10,10,200)
TF_Lambda_tmp = lambda s: np.sinc(s)**2
TF_Lambda = np.array([TF_func(TF_Lambda_tmp, t) for t in t])
/tmp/ipykernel_15231/2371038607.py:4: IntegrationWarning: The integral is probably divergent, or slowly convergent. re_tf, _ = quad(fct_tmp_re, -10, 10) /tmp/ipykernel_15231/2371038607.py:4: IntegrationWarning: The maximum number of subdivisions (50) has been achieved. If increasing the limit yields no improvement it is advised to analyze the integrand in order to determine the difficulties. If the position of a local difficulty can be determined (singularity, discontinuity) one will probably gain from splitting up the interval and calling the integrator on the subranges. Perhaps a special-purpose integrator should be used. re_tf, _ = quad(fct_tmp_re, -10, 10)
plt.plot(t, Lambda(t))
plt.plot(t, TF_Lambda[:,0])
plt.show()
On s'attend à obtenir $\Lambda$, car on a la relation $\mathcal{F}^{-1} \mathcal{F} \Lambda = \Lambda$, mais on voit ici que la reconstruction du signal triangle est bien meilleur que pour le signal porte. Il y a ici quelques parasites, mais le triangle est bien reconstruit.