Chapitre 2

Traitement du signal

2.1 Series de Fourier

Avant de rentrer dans le vif sujet et detudier la transformee de Fourier, nous allons nous arreter un instant
sur les series de Fourier. Ces dernieres permettent de comprendre les concepts qui sous-tendent toute lanalyse de
Fourier pour le traitement du signal que nous verrons dans ce cours et sont, selon moi, plus intuitives pour aborder
ces questions.

Definition
Pour une fonction f (x) definie sur [, ], on peut ecrire sa decomposition en serie de Fourier de la facon
suivante :

a0 X
f (x) = + an cos(nx) + bn sin(nx). (2.1)
2 n=1

Nous allons maintenant utiliser Matlab pour obtenir les coefficients an et bn .

Exercice 2.1
R
1. Ecrire f (x) cos(mx) dx a laide de lequation (2.1).
2. Decomposer lexpression precedente et isoler les 3 termes :

a0
Z
P1 = cos(mx) dx,
2

X Z
P2 = an cos(nx) cos(mx) dx,
n=1

X Z
P3 = bn sin(nx) cos(mx) dx.
n=1

3.a Utiliser une representation graphique sous Matlab pour trouver la valeur de P1.
3.b Utiliser une representation graphique sous Matlab pour trouver la valeur de P2.
3.c Utiliser une representation graphique sous Matlab pour trouver la valeur de P3 pour n 6= m.
3.d Utiliser une representation graphique sous Matlab pour trouver la valeur de P3 pour n = m.

1
R
4. Deduire de ce calcul la valeur des coefficients an et monter que an =
f (x) cos(nx) dx, pour tout n.

1
R
5. Effectuer la meme demarche pour obtenir la valeur des coefficients bn et monter que bn =
f (x) sin(nx) dx,
pour tout n.
2 2.2. Application a la synthese de signaux sous Matlab

Notation exponentielle
Dans lexercice precedent nous venons de voir que lon peut decomposer une fonction definie sur [, ] a laide
de fonctions sinus et cosinus. Il est possible de faire cette decomposition de maniere plus compacte a laide de la
notation exponentielle. On ecrira alors f (x) definie sur [, ] sous la forme :

X
f (x) = cn einx . (2.2)
n=

Notons que lon peut etendre a des fonctions f (x) definie sur [L, L] par :

n
X
f (x) = cn ei L x . (2.3)
n=

Exercice 2.2
En utilisant la relation dEuler : ei = cos + i sin , demontrer que :

a. cn = an ib2
n
pour n > 0
an +ibn
b. cn = 2 pour n < 0.
c. c0 = a20 .

2.2 Application a la synthese de signaux sous Matlab

Tres bien, nous avons donc maintenant compris que les series de Fourier correspondent a des decompositions
dun signal sur la base des fonctions trigonometriques. Essayons de faire un exemple sous Matlab.

Exercice 2.3 : Effet la somme non-infinie sur les series de Fourier

1. Creer un script qui va realiser les operation suivantes :

Definir la taille L de lechantillon utilise.

Creer un vecteur time 100L points regulierement espaces entre 0 et L.

Creer un vecteur signal contenant 100L points valant tous 0.

Definir a=[0 0 0 0 0 0 ] ; et b=[0 0 0 0 0 0 ] ;

2. Les vecteurs a et b sont les 6 premieres composantes de la serie de Fourier. Ajouter au script une boucle permet-
tant de generer une serie de Fourier tronquees au 6 premiers coefficients (equation 2.1).

3. Choisir les vecteurs a et b au hasard et afficher le signal en fonction de time.

2
4. Essayer les vecteurs an = 0 et bn = (1)n+1 n et afficher le signal en fonction de time. Quen dites vous ?

2
5. Ajouter maintenant 2 termes dans a et b (toujours an = 0 et bn = (1)n+1 n ) et afficher le signal en fonction de
time. A-t-il change ? Pourquoi ?

6. Realiser une boucle pour ajouter un nombre arbitraire de coefficients dans la serie et faire un test avec an = 0
2
et bn = (1)n+1 n . A partir de combien de coefficients lapproximation dun signal triangulaire vous semble bonne ?

7. Trouver la decomposition de Fourier dune fonction creneau et realiser un script afin de voir leffet du nombre de
coefficients sur la qualite de lapproximation.
 Traitement du signal 3

Exercice 2.4 : Synthe Matlab

Dans cet exercice nous allons realiser un synthetiseur Matlab ! Pour ce faire nous allons proceder comme a
lexercice precedent en definissant les composantes frequentielles a ajouter dans notre signal. Puis nous utiliserons la
fonction sound pour jouer ce magnifique morceau de musique electronique... Cest la premiere fois que nous allons
developper un petit projet Matlab qui utilisera plusieurs fonctions, soyez attentifs a vos noms de fichier.

1. Creer une fonction qui prendre en parametres dentree : T la duree du signal genere et f req la frequence de la
composante. Cette fonction donnera en sortie un vecteur qui sera le signal genere.

2. La frequence dechantillonnage de Matlab par defaut pour le son est 8192 points par seconde. On definit donc une
variable time telle que cette variable soit un vecteur regulierement espace entre 0 et T avec une longueur 8192*T.
On cree aussi une variable signal de meme taille et de valeur 0.

3. Creer un vecteur a de dimension egale a la frequence maximale f req et dont tous les coefficients valent 0 sauf le
coefficient f req qui vaut 1.

4. Realiser une boucle qui permet de faire la synthese du signal (equation 2.1). Attention pour avoir les frequences
en Hz et le temps en seconde, il est important davoir ici lexpression cos(n time 2) dans la serie de Fourier.

5. Tester le signal avec la commande sound(...). PAS plus long que 2 secondes sil vous plait !

6. Modifier le programme pour quil prenne en entree non pas une valeur de frequence mais un vecteur contenant
plusieurs frequences.

7. Faire un script qui appelle votre fonction afin de jouer Au clair de la lune.
4 2.3. Transformee de Fourier

2.3 Transformee de Fourier

2.3.1 Definition
Nous venons de voir quil etait possible de decomposer (ou synthetiser cest la meme chose), la plupart des
signaux analytiques a laide des series de Fourier. Il sagit dune somme des differents coefficients an et bn si lon
utilise la forme trigonometrique ou cn si lon utilise la forme exponentielle. Lorsque lon passe du discret (somme)
au continu (integrale), on dit que lon realise une transformee de Fourier. Il faut bien comprendre que lon parle
ici presque du meme animal mathematiques, la principale differences etant que lon ne realise plus une somme
des differentes composantes. Par analogie avec la relation (2.2), on ecrit pour une fonction f sa decomposition de
Fourier : Z
1
f (x) = eikx F (k) dk. (2.4)
2
Les composantes de Fourier se retrouvent donc dans la fonction F (k). Pour obtenir les coefficients de la decompo-
sition, cest a la fonction F (k), on fait loperation :
Z
1
F (k) = eikx f (x) dx. (2.5)
2
Cest cette operation qui sappelle la transformee de Fourier. La transformee de Fourier possede de nombreuses
proprietes tres interessantes, et je vous suggere de passer quelque temps avec un bon livre de Mathematiques pour
les revoir...

2.3.2 Transformee de Fourier discrete

On vient de voir que pour passer de la serie de Fourier a la transformee de Fourier on passait du discret au
continu. Or, vous le savez maintenant, Matlab (et les ordinateurs en general) travaillent avec des donnees discretes !
Il faut donc utiliser utiliser un algorithme pour faire loperation transformee de Fourier discrete.

La bonne nouvelle est que Cooley et Tukey dans le milieu des annees 60 ont trouve un algorithme tres efficace
pour faire des transformee de Fourier qui sappelle FFT (Fast Fourier Transform). Cette algorithme repose sur
une division du probleme en 2 sous-problemes, eux meme divises en 2 sous-sous-problemes, etc. Il faut donc un
nombre de point en 2N pour que lalgorithme FFT sapplique bien, ce qui est une condition importante lorsque vous
travaillerez avec la commande Matlab fft. Voyons maintenant comment utilise-t-on cette commande.

2.3.3 Un exemple pas a pas de FFT

Realisons ensemble le code suivant :
clear all; close all; clc;

L=20;
n=128;
x2=linspace(-L/2,L/2,n+1); x=x2(1:n);

u=exp(-x.^2);
plot(x,u)

ut=fft(u)
plot(ut) % sur certaine version de Matlab il faut faire plot(real(ut)).

Ajouter maintenant la commande fftshift :

uts=fftshift(ut)
plot(abs(uts))

Il faut maintenant redefinir correctement laxe des frequences :

k=(2*pi/L)[0:n/2-1 -n/2:-1];

2pi/L permet de normaliser au domaine L et non pas au domaine 2pi.

 Traitement du signal 5

Exercice 2.6 : FFT dun signal electrique - Rapport signal sur bruit.
On a mesure le signal electrique U (t) (tension en fonction du temps) avec un echantillonnage temporel dt = 0.1 ms
pour des valeurs de t comprises entre 0 et 10 s. La tension U (t) mesuree est donnee dans la matrice U de dimensions
(100001,1). Elle contient une composante periodique et un bruit aleatoire.
Telecharger ce fichier depuis http://quentinglorieux.fr/matlab/files/U.mat

1.Tracer U (t).
2.Tracer le spectre S(f ), cest a dire le carre de la norme de la transformee de Fourier, exprime en fonction de la
frequence f pour la fonction U (t).
3. Quelle est la periode de ce signal ?
4. Calculer le rapport signal sur bruit. Ce rapport est defini comme I2II1
1
, avec I1 lintegrale de S pour f entre 39.9
et 40.1 Hz et I2 lintegrale de S pour f entre 0 et 5000 Hz.

Exercice 2.7 : Identification des transformee de Fourier

Associer chaque fonction A, B, C. . . a sa transformee de Fourier 1, 2, 3 . . . en justifiant par les proprietes de la
fonction.
Exemple : G 3 La fonction G est une somme dune sinusode de periode 4 amortie sur un temps caracteristique
de 5 et dune sinusode de periode 0,5. Le spectre 3 : un pic de frequence 0,25 de largeur environ 0,2 et un pic etroit
de frequence 2.

!
6 2.3. Transformee de Fourier

Probleme 2.8 : Filtrage sonore

Nous allons analyser spectralement des enregistrements sonores. On pourra creer, pour gagner du temps, la
fonction TFourier.m suivante :
function [Nu,F] = TFourier(f,dt)
N = max(size(f)); % On determine le nombre de points de la fonction f(t)
Nu = (0:N-1)/(N*dt); % Rappel : Ttotal est donne par N*dt
F = abs(fft(f)); % On prend la valeur absolue de la TF.

On peut ouvrir un fichier .wav dans Matlab en utilisant licone ouvrir de lespace de travail (ou avec la
commande wavplay). Ce fichier consiste en une matrice a une ligne, correspondant a un signal S(t) enregistre a
la frequence dechantillonnage freq cest-a-dire que S(t) est connu pour des valeurs de t separees de 1/f req. On
peut lecouter a la frequence dechantillonnage freq par la commande wavplay(W,freq). Par defaut (si on tape
wavplay(W)) la frequence dechantillonnage est fixee a la valeur standard freq = 11025 Hz.
Telecharger les fichiers flute.wav et violon.wav depuis http://quentinglorieux.fr/matlab/files/

Question 1.
Une onde sonore sinusodale de frequence 784 Hz correspond a la note Sol. Ouvrir les fichiers flute.wav et
violon.wav , ou la note Sol, jouee par un instrument, est enregistree a la frequence dechantillonnage 10000 Hz.
a- Combien de temps dure chaque enregistrement ?
b- Tracer la transformee de Fourier de chaque enregistrement, en faisant attention a laxe des abscisses. Decrire ces
courbes : pic principal a 784 Hz, harmoniques, dedoublement des pics, pics parasites a basse frequence etc.
c- Tracer la courbe de lenregistrement de flute entre 250 et 350 ms. Comment les differentes caracteristiques obser-
vees sur la transformee de Fourier se retrouvent-elles sur le signal ?

Conclusion : cest notamment la presence dharmoniques qui fait la difference, pour une meme note (une meme
frequence), entre deux instruments.

Question 2.
a- Par quelle commande peut-on generer artificiellement un signal sonore Wsol dune duree de 2 secondes, a la
frequence dechantillonnage 11025 Hz, de la note Sol ?

b- Tracer la transformee de Fourier de ce signal.

c- Si une note correspond a une frequence donnee (par exemple, Sol : 784 Hz), la frequence double correspond a la
meme note mais a loctave suivante. Quobtient-on si on ecoute le signal sonore Wsol a la frequence dechantillon-
nage 22050 Hz ? et a la frequence 5512 Hz ?

Question 3.
On souhaite filtrer les hautes frequences du signal flute . Pour cela, on calcule la transformee de Fourier, on
supprime la partie contenant les frequences hautes, puis on effectue la transformee de Fourier inverse pour obtenir
le signal filtre :

F = fft(Wflute); F(300:7200)=0; Wfiltre = real(ifft(F));

Commentez le code ci-dessus puis tracer le signal Wfiltre obtenu et sa transformee de Fourier. Montrer que lon a
reussi a supprimer la composante a 800 Hz et a faire ressortir la composante a 150 Hz.

Question 4.
On souhaite generer artificiellement la note de musique jouee par un violon. Pour cela, reperer sur le spectre du
Sol du violon la hauteur des harmoniques 1 a 6, que lon notera A1 a A6 , et creer pour f = 784 Hz le signal
A1 cos(2f t) + A2 cos(4f t) + A3 cos(6f t) + A4 cos(8f t) + ...
Faire de meme pour la flute. Ecouter ces deux signaux : a-t-on correctement reconstruit les sons des instruments ?
 Traitement du signal 7

Probleme 2.9 : Spectroscopie de Fourier (exercice dexamen 2010)

On cherche le spectre dun signal lumineux. Ce spectre est note S(k) avec k le nombre donde defini comme
k = 1 et exprime en m1 .
On envoie le signal lumineux dans un interferometre de Michelson (on rappelle son schema ci-dessous a titre indicatif
mais sa comprehension nest pas necessaire pour lexercice). Lexperience consiste en la mesure de la puissance
lumineuse P en sortie de linterferometre a l aide dun detecteur (typiquement une photodiode) pour differentes
valeurs de la difference de marche . On donne le resultat P = P () P ( = +) de cette mesure (en unite
arbitraire) dans la matrice deltaP en fonction de (en m) dans la matrice delta. Les matrices sont a charger sur
http://quentinglorieux.fr/matlab/files/.

diffrence
de marche

Dtecteur

Question 1.
Tracer P () et placer les legendes.

Question 2.
La technique de spectroscopie par transformee de Fourier repose sur le fait que le spectre recherche S(k) nest autre
que la norme de la transformee de Fourier de la fonction P () mesuree. Tracer donc |S(k)|.

Question 3.
Combien cette courbe comporte-t-elle de pics ? En quelles longueurs donde sont-ils situes ? Lesquels vous paraissent
significatifs, lesquels attribuez-vous a des artefacts ?

Question 4.
Quelle est la resolution sur k (cest-a-dire lechantillonnage dk entre deux valeurs de k successives) que nous pou-
vons atteindre dans ces conditions experimentales ? Quaurait-il fallu changer dans lexperience pour ameliorer la
resolution ?

Question 5.
On sait que :
(k1 + k2 )x (k1 k2 )x
cos(k1 x) + cos(k2 x) = 2 cos cos( ).
2 2
Quel lien peut-on faire entre la presence dun doublet (deux pics tres proches) dans la transformee de Fourier
et la presence de battements (modulation des oscillations par une enveloppe sinusodale) de P () ?
8 2.4. Corriges

2.4 Corriges
Exercice 2.3 - Correction
On cree un signal en fonction de la variable time definit a laide de ses composantes de Fourier. On affiche ici
les 10 premieres composantes.

L=20;
time=linspace(0,L,100*L);
signal=zeros(100*L,1);

a=[1 1 0 0 0 0 0 0 0 0 ];
b=[0 0 0 0 0 0 0 0 0 0 ];

for n=1:length(a)
signal=signal+a(n)*cos(n*time)+b(n)*sin(n*time);
end

plot(time,signal)

2
On affiche ensuite la serie de Fourier avec an = 0 et bn = (1)n+1 n , avec les 6 premiers coefficients.

L=20;
time=linspace(0,L,100*L);
signal=zeros(100*L,1);

for n=1:6
signal=signal+(-1)^(n+1)*2/(pi*n)*sin(n*time);
end

plot(time,signal)

Puis on trace la meme courbes avec les 20 premiers coefficients.

Voila ce que lon obtient :

1.5

0.5

0.5

1.5
0 2 4 6 8 10 12 14 16 18 20

1.5

0.5

0.5

1.5
0 2 4 6 8 10 12 14 16 18 20

2
Figure 2.1: Series de Fourier avec an = 0 et bn = (1)n+1 n pour 6 (haut) et 20 (bas) coefficients.

Enfin, voila la decomposition dune fonction creneau :

for n=0:10
signal=signal+(-1)^(n)*2/((2*n+1)*pi)*cos((2*n+1)*time);
end
 Traitement du signal 9

Exercice 2.4 - Correction

Dans un premier temps, on cherche a synthetiser un son a 440 Hz. Voila le code permettant de faire cette
operation sous Matlab.
% SYNTHESE SONORE
clear all; clc;
L=1;
time=linspace(0,L,8192*L);
signal=time; signal=0;

a=zeros(1000,1);
b=zeros(1000,1);
a(440)=1; %choix de la frequence ici

for n=1:length(a)
signal=signal+a(n)*cos(n*time*2*pi)+b(n)*sin(n*time*2*pi);
end
signal=signal/max(abs(signal));
plot(time,signal)
sound(signal)
On peut transformer ce code relativement simplement en une fonction qui prenne en entree les parametres T et
freq :
function [signal]=note(T,freq)
L=T;
time=linspace(0,L,8192*L);
signal=time; signal=0;

a=zeros(freq,1);
b=zeros(freq,1);

a(freq)=1;

for n=1:length(a)
signal=signal+a(n)*cos(n*time*2*pi)+b(n)*sin(n*time*2*pi);
end
signal=signal/max(abs(signal));
sound(signal)
Et voila comment modifier le code pour en faire une fonction acceptant un vecteur pour T et pour freq. Par
exemple note([0.5 0.5 0.5],[400 800 400]);.
function [signal]=note(T,freq)

for i=1:length(T)
L=T(i);
time=linspace(0,L,8192*L);
signal=time; signal=0;
a=zeros(freq(i),1);
b=zeros(freq(i),1);
a(freq(i))=1;
for n=1:length(a)
signal=signal+a(n)*cos(n*time*2*pi)+b(n)*sin(n*time*2*pi);
end
signal=signal/max(abs(signal));
sound(signal)
pause(T(i))
end
10 2.4. Corriges

Exercice 2.6 - Correction

Apres avoir charger le fichier on affiche la fonction U , puis la norme carree de sa transformee de Fourier (son
spectre).

plot (U)
Ut=abs(fft(U)).^2;
k=(1/10)*[0:100000];
plot(k,Ut)
axis([1 100 0 5e7])

On observe alors un pic autour de 40 Hz, cest a dire 25 ms de periode.

On calcule ensuite le rapport signal sur bruit (RSB) :

signal = Ut(401);
bruit = sum(Ut([1:50001])) - signal;
RSB=signal/bruit

On applique ensuite un filtre gaussien pour eliminer le bruit :

center=40;
sigma=5;
gauss=exp(-(k-center).^2/sigma^2);
subplot (3,1,2); plot(k,gauss)
axis([1 100 0 1])

subplot (3,1,3); plot(k,gauss.*Ut)

axis([1 100 0 1e7])

Voila ce que lon obtient a laide du filtre gaussien :

6
x 10
10

0
10 20 30 40 50 60 70 80 90 100

0.8

0.6

0.4

0.2

0
10 20 30 40 50 60 70 80 90 100

6
x 10
10

0
10 20 30 40 50 60 70 80 90 100

Figure 2.2: De haut en bas, on obtient le spectre brut, le filtre gaussien, et le spectre filtre.

On peut alors effectuer la TF inverse.

 Traitement du signal 11

signalf=ifft(gauss.*Ut);
figure(2)
subplot (2,1,2);plot (real(signalf))
axis([1 1000 -1e5 1e5])
subplot (2,1,1);plot (U)
axis([1 1000 0 15])

Voila le signal filtre :

15

10

0
100 200 300 400 500 600 700 800 900 1000

5
x 10
1

0.5

0.5

1
100 200 300 400 500 600 700 800 900 1000

Figure 2.3: De haut en bas, on obtient le signal brut, et le signal filtre.

12 2.4. Corriges

Exercice 2.7 - Correction

Rsumons les proprits des diffrentes fonctions dans un tableau (on distingue une
composante priodique 1 et une ventuelle composante priodique 2) :

Fonction A B C D E F G H I

priode 1 1 1 0,5 1 0,5 1 0,5 4 4

amortissement 1 infini 5 5 5 infini infini infini infini infini

1 sinusodale ? Non Non Non Oui Non Oui Non Oui Non

priode 2 - - 4 - 4 - 4 - -

amortissement 2 - - infini - infini - 5 - -

2 sinusodale ? - - Non - Non - Non - -

Par ailleurs on a les transformes de Fourier suivantes (on distingue les pics 1 et 2) :

T. Fourier 1 2 3 4 5 6 7 8 9

frquence 1 0,25 1 2 1 1 2 1 0,25 2

largeur 1 0 0 0 0 0,2 0 0,2 0 0,2

Harmoniques de 1 ? Non Oui Non Non Non Non Oui Oui Non

frquence 2 - - 0,25 - - 0,25 - - 0,25

largeur 2 - - 0,2 - - 0 - - 0

Harmoniques de 2 ? - - Non - - Non - - Non

Pour les temps damortissement et les largeurs des pics, on sait que lun est de lordre
de lautre, mais on ne donne ici que des ordres de grandeur ; pour donner des valeurs
plus prcises, il faudrait avoir dfini prcisment la notion de largeur : par exemple
largeur mi-hauteur etc

Sachant quune composante de priode T, amortie sur un temps T2, correspond un pic
la frquence 1/T de largeur 1/T2 on conclut :

Fonction A B C D E F G H I

T. Fourier 4 5 9 7 6 2 3 8 1