9. Transformation de Fourier Discrète

Des mathématiques pour "débruiter" un signal ...

É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$.

Exercices de 1 à 3 page 300

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}$.

GeoGebra - Racines de l'unité

https://www.geogebra.org/classic/qagjxphn

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 ?

ATTENTION : ERREUR DE RECOPIE POUR X1 - LE RESULTAT FINAL EST JUSTE

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

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)$

Travaux pratiques

Pour coder du Python en ligne, on utilisera https://repl.it/ Pour effectuer des calculs formels en ligne, on utilisera https://www.xcasenligne.fr/ sous Firefox

359KB

DébuterEnPYTHON.pdf

pdf

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 !

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}]$.

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) :

Obtenu avec Python

Obtenu avec Xcas

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)+3fft(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.

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 ?

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.

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 ?

Débruitage d'un signal

Dans les programmes PYTHON qui vont suivre, nous importeront les deux packages : from matplotlib.pyplot import * from numpy import *

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)

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$.

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.

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 :

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é. :

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 ?

''' ? '''
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')