> For the complete documentation index, see [llms.txt](https://arnaud-lierville.gitbook.io/bts-sn/llms.txt). Markdown versions of documentation pages are available by appending `.md` to page URLs; this page is available as [Markdown](https://arnaud-lierville.gitbook.io/bts-sn/cours/transformation-de-fourier-discrete-tdf.md).
# 9. Transformation de Fourier Discrète
## Échauffement
On rappelle
* Quelques **valeurs remarquables** de trigonométrie :\
* La **formule d'Euler** :$$e^{i\varphi}=\cos{\varphi}+i\sin{\varphi}$$\
* Le **formule de Moivre** : $$(\cos{\varphi}+i\sin{\varphi})^n=(e^{i\varphi})^n=e^{in\varphi}=\cos{(n\varphi)}+i\sin{(n\varphi)}$$
#### Exercice
Soit $$j=e^{i\frac{2\pi}{3}}$$.
Calculer $$1+j+j^2$$.
{% hint style="success" %}
Exercices de 1 à 3 page 300
{% endhint %}
## Racines de unité
On cherche à résoudre dans $$\mathbb{C}$$ l'équation $$z^n=1$$, c'est à dire **"trouver les racines de l'unité"**.
* Pour $$n=2$$, l'équation $$z^2=1$$ possède $$2$$ solutions : $$1$$ et $$-1$$\
soit $$1=e^{i\frac{2\pi}{2}\times 0}$$ et $$-1=e^{i\frac{2\pi}{2}\times 1}$$
* Pour $$n=4$$, l'équation $$z^4=1$$ possède $$4$$ solutions : $$1$$, $$i$$, $$-1$$ et $$-i$$ \
soit $$1=e^{i\frac{2\pi}{4}\times 0}$$, $$i=e^{i\frac{2\pi}{4}\times 1}$$, $$-1=e^{i\frac{2\pi}{4}\times 2}$$, $$-i=e^{i\frac{2\pi}{4}\times 3}$$.
{% embed url="" %}
GeoGebra - Racines de l'unité
{% endembed %}
{% hint style="success" %}
{% endhint %}
**Cas** $$n=3$$:
On cherche $$z\in\mathbb{C}$$ tel que $$z^3=1$$. En écrivant $$z=re^{i\theta}$$, il vient $$r^3e^{i3\theta}=1e^{i\times 0}$$.\
Ainsi il faut pour le réel $$r=1$$, et pour $$3\theta=0;(2\pi) = 2k\pi$$ avec $$k\in\mathbb{Z}$$.\
Ainsi les $$3$$solutions sont $$z=e^{i\frac{2k\pi}{3}}$$ avec $$k=0;;;1;;;2$$.
**Cas général :**
Les solutions complexes de l'équation $$z^n=1$$ sont les $$n$$ nombres complexes :\
$$z\_k=e^{i\frac{2k\pi}{n}}$$ pour $$k$$ de $$0$$ à $$n$$.\
En notant $$\omega=e^{i\frac{2\pi}{n}}$$, ces solutions sont les puissances de $$\omega$$, de $$k=0$$ à $$n-1$$:\
$$1;;;\omega;;;\omega^2;;\dots;;\omega^{k};;\dots;;;\omega^{n-1}$$.\
On les appelles les **racines nième de l'unité.**
#### Exercice
1. Déterminer les racines sixièmes de l'unité et les représenter dans le plan complexe.
2. En remarquant la somme de termes consécutifs d'une suite géométrique, démontrer que si $$\omega$$ est une racine sixième de l'unité, alors $$1+\omega+\omega^2+\omega^3+\omega^4+\omega^5=0$$.
## Transformée de Fourier Discrète
### Définition
#### Exemple : TFD pour n = 3
La transformée de Fourier Discrète consiste à obtenir à partir d'une séquence de 3 nombres complexes $$(x\_0;;;x\_1;;;x\_2)$$ une séquence de 3 nouveaux nombres complexes $$(X\_0;;;X\_1;;;X\_2)$$.\
On procède de la façon suivante :
On pose $$\omega = e^{i\frac{2\pi}{3}}$$, et alors :
* $$X\_0=x\_0\times\omega^{-0\times0}+x\_1\times\omega^{-0\times1}+x\_2\times\omega^{-0\times2}$$
* $$X\_1=x\_0\times\omega^{-1\times0}+x\_1\times\omega^{-1\times1}+x\_2\times\omega^{-1\times2}$$
* $$X\_2=x\_0\times\omega^{-2\times0}+x\_1\times\omega^{-2\times1}+x\_2\times\omega^{-2\times2}$$
*Remarque* : on a $$X\_0=x\_0+x\_1+x\_2$$
#### Exemple corrigé :
1. Calculer la transformée de Fourier Discrète : TFD(2, 0, 1)
2. a) Écrire la matrice $$M=(\omega^{-i\times j})$$\
b) Calculer $$M\times \left(\begin{array}{c}x\_0\x\_1\x\_2\\\end{array}\right)$$\
c) Que retrouve-ton ?
#### Cas général
Pour $$n\in\mathbb{N}$$, on pose $$\omega=e^{i\frac{2\pi}{n}}$$.
$$TFD(x\_0;;;x\_1;;;\dots;;;x\_{n-1})=(X\_0;;;X\_1;;;\dots;;;X\_{n-1})$$
Avec $$X\_i=\sum\_{j=0}^{n-1}x\_j\times\omega^{-i\times j}=x\_0\times\omega^{-i\times 0}+x\_1\times\omega^{-i\times 1}+\cdots+x\_{n-1}\times\omega^{-i\times(n-1)}$$
#### Exercice
1. Donner la forme algébrique de $$\omega=e^{i\frac{2\pi}{4}}$$ .
2. Calculer la $$TFD(1, -2, 1, 3)$$
3. Vérifier par le calcul matriciel
{% hint style="success" %}
**Exercices :**\
À la main : $$TFD(1;;;2)$$, $$TFD(0;;;1)$$, $$TDF(4;;;2)$$, $$TFD(1;;;0;;;0)$$, $$TFD(1;;;0;;;1)$$\
Avec une calculatrice : $$TFD\left(0;;;\frac{1}{2};;;\frac{\sqrt{3}}{2}\right)$$
{% endhint %}
## Travaux pratiques
{% hint style="success" %}
Pour coder du **Python** en ligne, on utilisera \
Pour effectuer des calculs formels en ligne, on utilisera sous Firefox
{% endhint %}
{% file src="/files/u6RQlW4CZfRwaHa8dwYE" %}
### TP 1 : Xcas, Python et propriétés de la TFD
#### La Transformée de Fourier Discrète en PYTHON
1. Le **package cmath** permet d'effectuer des calculs avec les nombres complexes.\
Le nombre complexe $$i$$ est codé par `1j` . **Essayez !**
```python
from cmath import *
print(1j)
print((1j)**2)
print((2-3j)*(-1+2j))
print(e**(1j*pi/2))
print((sqrt(2)/2+sqrt(2)*1j/2)**2)
```
2\. Recopier et compléter le programme suivant. La **procédure TFD** doit calculer la Transformée de Fourier Discrète de la séquence $$\[x\_0, x\_1, \dots, x\_{n-1}]$$.
```python
from cmath import *
def tfd(x): # x est un tableau [a, b, ..., c]
n = len(x) # nombre d'éléments du tableau
X = []
w = exp(2*pi*1j/n)
for i in range(n): # i = 0 à n-1
Xi = 0
for j in range(n):
Xi = Xi +x[j]*...
X.append(...) # Ajoute un élément au tableau X
return X
print(tfd([1,2,3]))
```
Vérifier en comparant les résultats du programme avec ceux obtenus par Xcas pour différentes séquences comme ici pour **TFD(1,2,3)** :
#### Les propriétés de la TDF
Avec Xcas comparer :
* 4\*fft(1,-0.5) **et** fft(1,-0.5)
* fft(1,-1,2)+fft(2,1,-1) **et** fft(3,0,1)
Comment écrire fft(1,0,0)+2fft*(0,1,0)+3*fft(0,0,1) en **une seule commande fft** ?\
Cette propriété s'appelle **linéarité de la TFD**.
**La formule de Bessel**
Pour $$TFD(x\_0;;;x\_1;;;x\_2;;\dots;;;x\_{n-1})=(X\_0;;;X\_1;;;X\_2;;\dots;;;X\_{n-1})$$\
on a $$\sum\_{k=0}^{n-1}|x\_k|^2=\frac{1}{n}\sum\_{k=0}^{n-1}|X\_k|^2$$
Nous allons vérifier cette propriété à l'aide de programmes PYTHON.\
\
Recopiez ce programme, analysez-le et utilisez-le avec la procédure TFD pour vérifier cette formule pour la séquence \[1 ; 2 ; 3 ]. Essayez avec d'autres séquences.
```python
def somme_des_carrés(x):
S = 0
for xi in x:
S = S + abs(xi)**2
return S
print(somme_des_carrés([1,2,3]))
```
#### La transformée de Fourier Discrète Inverse
Pour $$TFD(x\_0;;;x\_1;;;x\_2;;\dots;;;x\_{n-1})=(X\_0;;;X\_1;;;X\_2;;\dots;;;X\_{n-1})$$, on peut retrouver $$(x\_0;;;x\_1;;;x\_2;;\dots;;;x\_{n-1})$$ à partir de $$(X\_0;;;X\_1;;;X\_2;;\dots;;;X\_{n-1})$$ en utilisant la **transformée inverse** donnée par la formule (avec $$\omega=e^{i\frac{2\pi}{n}}$$) :
$$x\_i=\dfrac{1}{n}\sum\_{j=0}^{n-1}Xj\times\omega^{i\times j}=\dfrac{1}{n}\left(X\_0\times\omega^{i\times 0}+X\_1\times\omega^{i\times 1}+\cdots+X{n-1}\times\omega^{i\times(n-1)}\right)$$
Les **moins en exposant** ont *disparu* et il est *apparu* **un facteur** $$\dfrac{1}{n}$$.
Étudiez ce petit programme. Que permet-il de faire ?
```python
liste1 = [1, 2, 3, 4]
liste2 = [element/2 for element in liste1]
print(liste1)
print(liste2)
```
En copiant et modifiant le programme **tfd,** codez le programme **tfd\_inv** et vérifiez qu'il fonctionne bien avec la séquence \[1, 2, 3]. Essayez avec d'autres séquences.
```python
def tfd_inv(X):
n = len(X)
x = []
w = ...
for i in range(n):
xi = 0;
for j in range(n):
xi = xi + ...
x.append(xi)
x = [xi/... for xi in x]
return x
print(tfd_inv(tfd([1,2,3])))
```
### TP2 : Débruiter un signal
Comment passer de ça à ça ?
{% hint style="info" %}
Dans les programmes PYTHON qui vont suivre, nous importeront les deux packages :\
\
from matplotlib.pyplot import *\**\
from numpy import \*
{% endhint %}
#### A - Définir un signal, l'échantillonner, le représenter
Un signal sinusoïdal de fréquence $$f$$ est défini par $$s(t)=\sin{(2\pi\times f\times t)}$$
Pour **échantillonner** ce signal sur 1 seconde à une fréquence $$Fe=1000$$, on prend 1000 abscisses également espacées de l'intervalle $$\[0;;;1]$$, pour lesquelles on récupère les images par le signal $$s$$.\
\
1\. Taper dans la console repl.it la commande `linspace(0, 0.999, 1000)`
2\. La commande\
`def s(t):`\
`return sin(2*pi*t)` \
permet de définir un signal de fréquence 1.\
\
**Copier-coller** le programme ci-dessous et **le compléter** pour qu'il représente un signal sinusoïdal de fréquence 15 hertz, sur 1 seconde, échantillonné à 1000 Hertz :\
*(le fichier signal.png contiendra la représentation attendue : l'ouvrir après l'exécution du programme)*
```python
import os
os.environ['MPLCONFIGDIR'] = os.getcwd() + "/configs/"
from matplotlib.pyplot import *
from numpy import *
l = ... # Longueur du signal (secondes)
Fe = ... # Freq échantillonage
Te = 1/Fe # Période échantillonage
N = int(l*Fe) # Nombre d'échantillons
# signal à échantillonner
def s(t):
return ...
# variable t sur [0, l-Te] au pas l/N
t = linspace(0, l-Te, N)
''' Signal '''
fig = figure()
plot(t, s(t), label='Signal')
legend(loc='upper center', shadow=True, fontsize='x-small')
title("Signal")
xlabel("Temps")
fig.savefig('signal.png')
```
#### B - Transformée de Fourier Discrète d'un signal
En PYTHON, à l'aide du package **numpy**, on peut calculer la transformée de Fourier discrète d'une séquence : `fft.fft([1,2,3])`
Essayez dans repl.it en ayant chargé le package numpy :\
`from numpy import *`
Obtenir **le spectre des fréquences** d'un signal $$s$$ consiste à
* échantillonner ce signal
* calculer la transformée de Fourier discrète de cet échantillon
* représenter la séquence obtenue
Utiliser le programme suivant pour obtenir le spectre des fréquences du signal $$s$$.
```python
import os
os.environ['MPLCONFIGDIR'] = os.getcwd() + "/configs/"
from matplotlib.pyplot import *
from numpy import *
l = 0.6 # Longueur du signal (secondes)
Fe = 1000.0 # Freq échantillonnage
Te = 1/Fe # Période échantillonnage
N = int(l*Fe) # Nombre d'échantillons
# signal à échantillonner
def s(t):
return sin(2*pi*40*t) + .5*sin(2*pi*90*t+pi/4) + .7*sin(2*pi*120*t-pi/2)
# variable t sur [0, l-Te] au pas l/N
t = linspace(0, l-Te, N)
# On ne garde que la moitié des points (symétrie)
N1 = floor(N/2);
# Fréquences en abscisse
fr = linspace(0, int(N1 -1), int(N1))*Fe/N
''' Spectre du signal '''
y = fft.fft(s(t)) # TFD du signal
fig = figure()
plot(fr[0:int(N1 -1)], abs(y[0:int(N1 -1)]))
title("Spectre des fréquences du signal")
xlabel("Fréquences")
fig.savefig('spectre-signal.png')
```
**Expliquer le spectre des fréquences obtenu :**
1. quelles fréquences contient ce spectre ?
2. où retrouve-t-on ces fréquences dans le signal ?
3. comment expliquer les différences d'amplitude d'une fréquence à une autre dans le spectre ?
4. quel est l'effet sur le spectre des déphasages $$+\dfrac{\pi}{4}$$ et $$-\dfrac{\pi}{2}$$ dans les composantes du signal ?
#### C - Bruitage d'un signal
On souhaite maintenant **"bruiter"** un signal et obtenir un spectre de fréquences comme celui-ci :
Notre signal de départ est le suivant :
$$s(t)=\sin(2\pi\times40\times t) + \sin(2\pi\times90\times t+\pi/4)$$
La commande `random.normal(0, 1, N)` fournit une séquence aléatoire de N nombres distribués selon la loi normale.\
**Essayez-la dans la console sans oublier le package numpy**.
**Copier-coller** et **compléter** le programme suivant pour qu'il affiche le spectre du signal bruité. Pour bruiter le signal, on ajoutera cette séquence aléatoire au signal initial.
```python
import os
os.environ['MPLCONFIGDIR'] = os.getcwd() + "/configs/"
from matplotlib.pyplot import *
from numpy import *
l = 0.6 # Longueur du signal (secondes)
Fe = 1000.0 # Freq échantillonnage
Te = 1/Fe # Période échantillonnage
N = int(l*Fe) # Nombre d'échantillons
# signal à échantillonner
def s(t):
return ...
# variable t sur [0, l-Te] au pas l/N
t = linspace(0, l-Te, N)
# On ne garde que la moitié des points (symétrie)
N1 = floor(N/2);
# Fréquences en abscisse
fr = linspace(0, int(N1 -1), int(N1))*Fe/N
''' Spectre du signal bruité '''
fig = figure()
sb = s(t) + ... # bruitage du signal
yb = ... # TFD du signal bruité
plot(fr[0:int(N1 -1)], abs(yb[0:int(N1 -1)]))
title("Spectre des fréquences du signal bruité")
xlabel("Fréquences")
fig.savefig('spectre-signal-bruité.png')
```
#### D - Débruitage du signal
**Travail préparatoire :**\
copier-coller le code suivant et répondre aux questions :
```python
from numpy import *
# Que fait zeros ?
l1 = zeros(shape= 5)
print('l1')
print(l1)
print()
# random.normal(0, 1, 5) fournit une séquence aléatoire
l2 = random.normal(0, 1, 5)
print('l2')
print(l2)
print()
# Que fait abs ?
l3 = zeros(shape= 5)
for k in range(l3.size):
l3[k] = abs(l2[k])
print('l3')
print(l3)
print()
# Expliquez le résultat
l4 = zeros(shape= 5)
for k in range(l4.size):
l4[k] = (abs(l2[k]) > 0.7)
print('l4')
print(l4)
print()
# Expliquez le résultat
l5 = zeros(shape= 5)
for k in range(l5.size):
l5[k] = l2[k]*(abs(l2[k]) > 0.7)
print('l5')
print(l5)
print()
```
Le principe du **"débruitage"** d'un signal, consiste à ne garder que les fréquences du spectre supérieures à un certain **seuil** (fixé ici à 55).\
À l'aide des instructions découvertes précédemment, copiez et complétez le programme suivant. Il doit afficher le spectre des fréquences du signal **débruité**. :
```python
import os
os.environ['MPLCONFIGDIR'] = os.getcwd() + "/configs/"
from ... import *
from ... import *
l = 0.6 # Longueur du signal (secondes)
Fe = 1000.0 # Freq échantillonnage
Te = 1/Fe # Période échantillonnage
N = int(l*Fe) # Nombre d'échantillons
# signal à échantillonner
def s(t):
return sin(2*pi*40*t) + sin(2*pi*90*t+pi/4)
# variable t sur [0, l-Te] au pas l/N
t = linspace(0, l-Te, N)
# On ne garde que la moitié des points (symétrie)
N1 = floor(N/2);
# Fréquences en abscisse
fr = linspace(0, int(N1 -1), int(N1))*Fe/N
# bruitage du signal
sb = s(t) + random.normal(0, 1, N)
# TDF du signal bruité
yb = fft.fft(sb)
''' Spectre du signal débruité '''
fig = figure()
# débruitage
seuil = 55
y_debruite = zeros(shape=(yb.size), dtype=clongdouble)
for k in range(...):
y_debruite[k]=(yb[k])*...
plot(fr[0:int(N1 -1)], abs(y_debruite[0:int(N1 -1)]))
title("Spectre des fréquences du signal débruité")
xlabel("Fréquences")
fig.savefig('spectre-signal-débruité.png')
```
La commande `ifft` permet de calculer la Transformée de Fourier Discrète Inverse. Ajouter ces lignes de code en fin de programme. Exécuter-le plusieurs fois. Que fait-il ?
```python
''' ? '''
s_debruite = real(fft.ifft(y_debruite))
fig2 = figure()
plot(t, s(t), label='Signal pur')
plot(t, s_debruite, ':', label='Signal débruité')
legend(loc='upper center', shadow=True, fontsize='x-small')
title("?")
xlabel("Temps")
fig2.savefig('what-is-it.png')
```