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 :eiφ=cosφ+isinφe^{i\varphi}=\cos{\varphi}+i\sin{\varphi}

  • Le formule de Moivre : (cosφ+isinφ)n=(eiφ)n=einφ=cos(nφ)+isin(nφ)(\cos{\varphi}+i\sin{\varphi})^n=(e^{i\varphi})^n=e^{in\varphi}=\cos{(n\varphi)}+i\sin{(n\varphi)}

Exercice

Soit j=ei2π3j=e^{i\frac{2\pi}{3}}.

Calculer 1+j+j21+j+j^2.

Racines de unité

On cherche à résoudre dans C\mathbb{C} l'équation zn=1z^n=1, c'est à dire "trouver les racines de l'unité".

  • Pour n=2n=2, l'équation z2=1z^2=1 possède 22 solutions : 11 et 1-1 soit 1=ei2π2×01=e^{i\frac{2\pi}{2}\times 0} et 1=ei2π2×1-1=e^{i\frac{2\pi}{2}\times 1}

  • Pour n=4n=4, l'équation z4=1z^4=1 possède 44 solutions : 11, ii, 1-1 et i-i soit 1=ei2π4×01=e^{i\frac{2\pi}{4}\times 0}, i=ei2π4×1i=e^{i\frac{2\pi}{4}\times 1}, 1=ei2π4×2-1=e^{i\frac{2\pi}{4}\times 2}, i=ei2π4×3-i=e^{i\frac{2\pi}{4}\times 3}.

GeoGebra - Racines de l'unité

Cas n=3n=3:

On cherche zCz\in\mathbb{C} tel que z3=1z^3=1. En écrivant z=reiθz=re^{i\theta}, il vient r3ei3θ=1ei×0r^3e^{i3\theta}=1e^{i\times 0}. Ainsi il faut pour le réel r=1r=1, et pour 3θ=0  (2π)=2kπ3\theta=0\;(2\pi) = 2k\pi avec kZk\in\mathbb{Z}. Ainsi les 33solutions sont z=ei2kπ3z=e^{i\frac{2k\pi}{3}} avec k=0  ;  1  ;  2k=0\;;\;1\;;\;2.

Cas général :

Les solutions complexes de l'équation zn=1z^n=1 sont les nn nombres complexes : zk=ei2kπnz_k=e^{i\frac{2k\pi}{n}} pour kk de 00 à nn. En notant ω=ei2πn\omega=e^{i\frac{2\pi}{n}}, ces solutions sont les puissances de ω\omega, de k=0k=0 à n1n-1: 1  ;  ω  ;  ω2  ;;  ωk  ;  ;  ωn11\;;\;\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+ω+ω2+ω3+ω4+ω5=01+\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 (x0  ;  x1  ;  x2)(x_0\;;\;x_1\;;\;x_2) une séquence de 3 nouveaux nombres complexes (X0  ;  X1  ;  X2)(X_0\;;\;X_1\;;\;X_2). On procède de la façon suivante :

On pose ω=ei2π3\omega = e^{i\frac{2\pi}{3}}, et alors :

  • X0=x0×ω0×0+x1×ω0×1+x2×ω0×2X_0=x_0\times\omega^{-0\times0}+x_1\times\omega^{-0\times1}+x_2\times\omega^{-0\times2}

  • X1=x0×ω1×0+x1×ω1×1+x2×ω1×2X_1=x_0\times\omega^{-1\times0}+x_1\times\omega^{-1\times1}+x_2\times\omega^{-1\times2}

  • X2=x0×ω2×0+x1×ω2×1+x2×ω2×2X_2=x_0\times\omega^{-2\times0}+x_1\times\omega^{-2\times1}+x_2\times\omega^{-2\times2}

Remarque : on a X0=x0+x1+x2X_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=(ωi×j)M=(\omega^{-i\times j}) b) Calculer M×(x0x1x2)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 nNn\in\mathbb{N}, on pose ω=ei2πn\omega=e^{i\frac{2\pi}{n}}.

TFD(x0  ;  x1  ;    ;  xn1)=(X0  ;  X1  ;    ;  Xn1)TFD(x_0\;;\;x_1\;;\;\dots\;;\;x_{n-1})=(X_0\;;\;X_1\;;\;\dots\;;\;X_{n-1})

Avec Xi=j=0n1xj×ωi×j=x0×ωi×0+x1×ωi×1++xn1×ωi×(n1)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 ω=ei2π4\omega=e^{i\frac{2\pi}{4}} .

  2. Calculer la TFD(1,2,1,3)TFD(1, -2, 1, 3)

  3. Vérifier par le calcul matriciel

Travaux pratiques

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 ii 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 [x0,x1,,xn1][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(x0  ;  x1  ;  x2  ;  ;  xn1)=(X0  ;  X1  ;  X2  ;  ;  Xn1)TFD(x_0\;;\;x_1\;;\;x_2\;;\dots\;;\;x_{n-1})=(X_0\;;\;X_1\;;\;X_2\;;\dots\;;\;X_{n-1}) on a k=0n1xk2=1nk=0n1Xk2\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(x0  ;  x1  ;  x2  ;  ;  xn1)=(X0  ;  X1  ;  X2  ;  ;  Xn1)TFD(x_0\;;\;x_1\;;\;x_2\;;\dots\;;\;x_{n-1})=(X_0\;;\;X_1\;;\;X_2\;;\dots\;;\;X_{n-1}), on peut retrouver (x0  ;  x1  ;  x2  ;  ;  xn1)(x_0\;;\;x_1\;;\;x_2\;;\dots\;;\;x_{n-1}) à partir de (X0  ;  X1  ;  X2  ;  ;  Xn1)(X_0\;;\;X_1\;;\;X_2\;;\dots\;;\;X_{n-1}) en utilisant la transformée inverse donnée par la formule (avec ω=ei2πn\omega=e^{i\frac{2\pi}{n}}) :

xi=1nj=0n1Xj×ωi×j=1n(X0×ωi×0+X1×ωi×1++Xn1×ωi×(n1))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 1n\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 ff est défini par s(t)=sin(2π×f×t)s(t)=\sin{(2\pi\times f\times t)}

Pour échantillonner ce signal sur 1 seconde à une fréquence Fe=1000Fe=1000, on prend 1000 abscisses également espacées de l'intervalle [0  ;  1][0\;;\;1], pour lesquelles on récupère les images par le signal ss. 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 ss 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 ss.

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 +π4+\dfrac{\pi}{4} et π2-\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π×40×t)+sin(2π×90×t+π/4)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')

Last updated

Was this helpful?