TP TS 1en

TRAVAUX PRATIQUES DE TRAITEMENT DU SIGNAL
Volume 1
Marie Chabert et Corinne Mailhes
Anne dEdition : 2010
Table des matires

Avant-Propos
1 TP Probabilits et Statistiques
1.1
Etude de lhistogramme damplitude . . . . . . . . . . . . . . . . . . . . . . . . .
1.1.1
Variables alatoires discrtes - Estimation de la probabilit . . . . . . . .
1.1.2
Variables alatoires continues - Estimation de la densit de probabilit . .
1.1.3
Estimation de la fonction de rpartition . . . . . . . . . . . . . . . . . . .
12
1.1.4
Travail eectuer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
12
1.2
Loi des grands nombres . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13
1.3
Thorme Central-Limite - Test de Kolmogorov . . . . . . . . . . . . . . . . . . .
13
1.3.1
Thorme central-limite . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13
1.3.2
Test de Kolmogorov . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
14
1.3.3
Travail eectuer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
14
Etude de distributions particulires . . . . . . . . . . . . . . . . . . . . . . . . . .
15
1.4.1
Triangle de Galton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
15
1.5
Test de Neyman-Pearson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
16
1.6
Annexes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
17
1.6.1
Demos de la bote outils stats de Matlab . . . . . . . . . . . . . . . .
17
1.6.2
Rappels de quelques distributions . . . . . . . . . . . . . . . . . . . . . . .
17
1.6.3
Table du Test de Kolmogorov . . . . . . . . . . . . . . . . . . . . . . . . .
18
1.4
2 TP Echantillonnage et Quantification
2.1
2.2
21
Rappels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
21
2.1.1
Echantillonnage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
21
2.1.2
Restitution aprs chantillonnage du signal dorigine . . . . . . . . . . . .
22
2.1.3
Quantification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
25
2.1.4
Annexe : Test de Kolmogorov . . . . . . . . . . . . . . . . . . . . . . . . .
29
Travail raliser . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
31
TABLE DES MATIRES
2.3
2.2.1
31
2.2.2
Quantification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
32
2.2.3
Restitution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
33
2.2.4
Signal de parole
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
33
2.2.5
Dmonstrations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
34
Elements de programmation Matlab . . . . . . . . . . . . . . . . . . . . . . . . .
34
2.3.1
Quelques fonctions Matlab utiles dans le TP . . . . . . . . . . . . . . . .
34
2.3.2
34
2.3.3
Quantification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
35
3 TP Analyse Spectrale - Corrlations et Spectres

3.1
3.2
3.3
3.4
3.5
Rappels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
37
3.1.1
Proprits des fonctions de corrlation . . . . . . . . . . . . . . . . . . . .
37
3.1.2
Algorithmes de calcul des fonctions de corrlation . . . . . . . . . . . . .
38
Exemples dutilisation des fonctions de corrlation . . . . . . . . . . . . . . . . .
39
3.2.1
Dtection dun signal priodique noy dans un bruit . . . . . . . . . . . .
39
3.2.2
Identification dun filtre . . . . . . . . . . . . . . . . . . . . . . . . . . . .
40
Analyse spectrale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
40
3.3.1
Dirents estimateurs . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
40
3.3.2
Intrt du zero-padding . . . . . . . . . . . . . . . . . . . . . . . . . . . .
41
3.3.3
Corrlations thoriques de signaux particuliers . . . . . . . . . . . . . . .
42
3.3.4
Priodogramme dune sinusode bruite et estimation du SNR . . . . . . .
42
Travail eectuer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
43
3.4.1
Autocorrlations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
43
3.4.2
Estimation spectrale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
44
3.4.3
Analyse spectrale de signaux rels . . . . . . . . . . . . . . . . . . . . . .
47
Programmation Matlab (pour ceux qui le souhaitent) . . . . . . . . . . . . . . . .
47
3.5.1
Quelques fonctions Matlab utiles dans le TP . . . . . . . . . . . . . . . .
47
3.5.2
Autocorrlations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
48
3.5.3
Estimations spectrales . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
48
4 Filtrage Numrique
4.1
4.2
37
51
Filtre Rponse Impulsionnelle Finie (RIF) : rappels . . . . . . . . . . . . . . . .
51
4.1.1
Dfinition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
51
4.1.2
Synthse par la mthode de la fentre . . . . . . . . . . . . . . . . . . . .
52
4.1.3
Optimisation de la synthse obtenue . . . . . . . . . . . . . . . . . . . . .
53
Filtre Rponse Impulsionnelle Infinie (RII) : rappels . . . . . . . . . . . . . . .
53
Avant propos
4.2.1
4.3
Dfinition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
53
Travail eectuer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
54
4.3.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
54
4.3.2
Gabarit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
55
4.3.3
Filtres Rponse Impulsionnelle Finie (RIF) . . . . . . . . . . . . . . . .
55
4.3.4
Filtre Rponse Impulsionnelle Infinie (RII) . . . . . . . . . . . . . . . .
57
4.3.5
Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
58
5 Initiation Matlab
59
5.1
Commandes daide . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
59
5.2
Matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
59
5.2.1
Construction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
59
5.2.2
Oprations lmentaires . . . . . . . . . . . . . . . . . . . . . . . . . . . .
61
5.2.3
Oprations lment par lment . . . . . . . . . . . . . . . . . . . . . . . .
62
5.2.4
Fonctions lmentaires . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
62
5.2.5
Fonctions matricielles . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
63
5.3
Gestion de lespace de travail . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
63
5.4
Boucles, tests et relations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
63
5.4.1
Boucles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
64
5.4.2
Test : instruction if
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
64
5.4.3
Relations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
65
5.4.4
Pause et break . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
66
Fichiers .m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
66
5.5.1
Scripts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
66
5.5.2
Fonctions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
66
5.5.3
Textes, entres, messages derreurs . . . . . . . . . . . . . . . . . . . . . .
67
5.5.4
Vectorisation et pr-allocation . . . . . . . . . . . . . . . . . . . . . . . . .
67
5.6
Graphiques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
67
5.7
Application 1 : probabilits et statistiques . . . . . . . . . . . . . . . . . . . . . .
68
5.8
Application 2 : traitement du signal . . . . . . . . . . . . . . . . . . . . . . . . .
69
5.9
Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
69
5.5
Avant propos
Avant-Propos
Un ensemble de logiciels de Travaux Pratiques pour le Traitement du Signal a t labor par les
chercheurs de lquipe Signal et Communication (SC) de lInstitut de Recherche en Informatique
de Toulouse (IRIT), galement membres du TSA (Laboratoire de Tlcommunications Spatiales
et Aronautiques). Les domaines de recherche abords par ce groupe de chercheurs sont varis :
analyse spectrale, classification et reconnaissance des formes, compression de donnes, traitement
statistique du signal, tlcommunications spatiales... Pour plus de renseignements, on peut se
reporter aux sites Web :
http : //www.irit.fr/ Equipe SC
http : //www.tesa.prd.f r
Le Traitement du Signal ne doit pas tre vu seulement de faon thorique, au niveau des cours
et des Travaux Dirigs ; le Traitement du Signal est une matire vivante, actuellement en plein
essor, avec des applications de plus en plus nombreuses, dans des domaines trs varis. Donnons
deux exemples dtudes eectues au SIC dans le cadre de doctorats.
Une premire tude concerne le systme daide latterrissage des avions (ILS : Instrument
Landing System). Ce systme, adopt internationalement par lAviation Civile, utilise pour sa
partie radioalignement de piste une bande de frquences VHF qui est, depuis quelques annes,
trs vulnrable des brouilleurs externes, comme les radios FM, les brouilleurs industriels etc...
Il est trs important davoir une connaissance sur ces brouilleurs, aussi le Service Technique de
la Navigation Arienne a-t-il suscit une tude qui a fait lobjet dune thse. Le problme a t
abord laide des mthodes de Traitement du Signal. Le signal quon dsire conserver est un
signal de faible puissance (le brouilleur) par rapport au signal ILS utilis latterrissage qui doit
tre estim, puis rejet le plus parfaitement possible afin de ne pas nuire lidentification du
brouilleur. Ce travail a mis en jeu des mthodes destimation et de rjection des deux sinusodes
qui constituent le signal ILS, puis des mthodes de segmentation et de classification permettant
de reconnaitre les brouilleurs.
Une deuxime application relve du domaine biomdical. Lusage de llectromyographie
(EMG) en tant quindicateur de ltat du systme neuromusculaire remonte la fin du XIXme
7
Avant propos
sicle. Cette technique est aujourdhui trs importante dans la dtection et le suivi datteintes
nerveuses ou musculaires. Il sagit dobserver lactivit lectrique mise lors dune contraction
naturelle du muscle. Pendant longtemps, ces signaux taient recueillis au moyen dune lectrode quon piquait dans la masse musculaire du patient. Outre laspect traumatisant de la
mthode, le diagnostic savrait parfois dicile. Do lide dutiliser des lectrodes de surface,
nullement erayantes pour le malade. Malheureusement, linspection visuelle du signal recueilli
nest pas susante pour transmettre la caractrisation dune pathologie neuromusculaire. Les
mthodes de Traitement du Signal, en particulier les mthodes danalyse spectrale, de modlisation paramtrique et de classification, ont conduit llaboration dun outil daide au diagnostic
pour des signaux EMG recueillis en surface.
Pour mettre en vidence cet aspect application pratique du Traitement du Signal, les chercheurs
du SIC ont dcid de crer un ensemble de Travaux Pratiques, venant complter lenseignement
du Traitement du Signal de lENSEEIHT.
En deuxime anne, sont prvus quatre TP illustrant les cours de Thorie du Signal (TS)
et de Traitement Numrique du Signal (TNS). Le premier TP TS sintresse aux probabilits lmentaires, le deuxime traite de lchantillonnage et de la quantification. Le
premier TP TNS sintresse au calcul des fonctions de corrlations et des densits spectrales
de puissance et le deuxime au filtrage numrique.
En troisime anne, les TP sont plus spcifiques : deux TP sur lestimation paramtrique et
non paramtrique (cours de Reprsentation et Analyse des Signaux et de Modlisation
Paramtrique) et un TP sur le filtrage adaptatif (cours de Traitement Adaptatif ).
Ces TP ont t organiss autour de logiciels MATLAB. Ces logiciels ont t conus pour tre
utiliss partir de menus droulants, vitant lutilisateur de se plonger dans la lourde tche
de la programmation.
Des indications de programmation Matlab, donnes en fin de chaque TP, associes une initiation Matlab donne en annexe, permettent de reprogrammer facilement certaines fonctionnalits
des logiciels. Ceci peut constituer une aide la ralisation des projets Matlab de deuxime et
troisime anne.
Chapitre 1
TP Probabilits et Statistiques
Lancer prob_stat sous Matlab.
1.1
1.1.1
Etude de lhistogramme damplitude

Variables alatoires discrtes - Estimation de la probabilit
Comment obtenir un estimateur de la loi de probabilit dune variable alatoire discrte partir
de lhistogramme de N ralisations de cette variable? En dduire un estimateur de la fonction
de rpartition.
A laide du logiciel, observer les performances de ces estimateurs en fonction du nombre de
ralisations N pour les variables alatoires :
uniforme sur {1, ..., P 1}
binomiale B(P 1 = n, P 2 = p)
de Poisson de paramtre = P 1.
Le bouton OK permet de gnrer N nouvelles ralisations de la variable alatoire. Noter les
moments estims en fonction de N (moyenne et variance) et comparer la thorie.
1.1.2
Variables alatoires continues - Estimation de la densit

de probabilit
La probabilit pour que la variable (v.a.) X soit comprise entre x et x + dx peut tre approche,
pour dx petit, par p(x)dx, p(x) tant la densit de probabilit de la v.a. X. La probabilit quune
9
10
Chapitre 1. TP Probabilits et Statistiques
ralisation de la v.a. X soit comprise entre mq et (m + 1) q est :

P [m] =
(m+1)q
p(x)dx
mq
Considrons une fonction reprsente graphiquement sur un boulier. Sur ce boulier sont places
des billes qui reprsentent plus ou moins finement la fonction : les tiges peuvent tre assez
rapproches et les billes assez petites pour donner une impression de continuit selon les deux
axes de coordonnes.
Si on place ce boulier dans le plan vertical et quensuite, on lui fait subir une rotation de 90
dans ce plan vertical, les billes vont glisser le long des fils et on aura le nombre de billes sur
chaque fil, image (grossire bien sr) de la densit de probabilit. Cest exactement lopration
que lon ralise lorsquon fait lhistogramme damplitude. Cest ce que nous nous proposons de
faire dans ce T P : utiliser lhistogramme pour estimer la densit de probabilit dune variable
alatoire.
Soit Z(t) la variable indicatrice dfinie par :
x
,x +
Z(t) = 1 si X(t) x
2
x
Z(t) = 0 si X(t)
/ x
,x +
2
x
2
x
2
Calculons la moyenne de cette variable :

E {Z (t)} = 1.P [Z (t) = 1] + 0.P [Z (t) = 0]
x
x
= P [Z (t) = 1] = P X(t) x
,x +
2
2
On obtient finalement :
E {Z (t)}
= p (x) x si x susamment petit
Estimer la densit de probabilit revient construire un estimateur de moyenne de la variable
indicatrice. On dfinit des classes, cest--dire que laxe des y est divis en intervalles de mme
1.1. Etude de lhistogramme damplitude
11
largeur x. On note Zi (t) la variable dfinie de la mme faon que Z(t) sur lintervalle n i.
Lestimateur habituel de la moyenne de cette v.a. Zi (t) est :
N
1 X
Nx
Zi (t) =
N
N
t=1
Nx tant le nombre dobservations de X(t) appartenant lintervalle considr de largeur x

et N tant le nombre total dobservations de X(t). Do lestimation de la densit de probabilit
suivante :
Nx
N x
x
,
x
+
o Nx Nombre dchantillons [x x
2
2 ] reprsente lhistogramme. Les valeurs
pb(x)
minimale et maximale prises par la variable alatoire sur les ralisations considres et le nombre
dintervalles (ou nombre de classes) permettent de dterminer la valeur de x.
Un estimateur est caractris par son biais et sa variance.
Soit b
un estimateur dun paramtre . Le biais de cet estimateur est dfini par :
n o
biais () = E b

En dautres termes, le biais reprsente lerreur systmatique que lon commet lorsquon estime .
La variance de b
permet de mesurer la vitesse de convergence de b
vers lorsque N augmente.
Un estimateur asymptotiquement sans biais
n o
=
lim E b
et dont la variance vrifie :
n o
=0
lim var b
est convergent. Etudions les proprits de lestimateur de la densit de probabilit que nous
avons dfini, pb(x). Pour cela, calculons la moyenne de la variable indicatrice Z(t) :
Z x+ x
2
p(u)du
E {Z(t)} =
x x
2
Pour x, largeur dune classe de lhistogramme susamment petite et p(x) rgulire, par un
dveloppement en srie de Taylor, on a :
x2
p(x) x
E {Z(t)} p(x) +
24
Lestimateur de densit de probabilit possde donc un biais additif :

x2
p(x)
24
Si les chantillons sont indpendants, alors la variance est gale :
1
p(x)
V ar {b
p(x)}
var [Z(t)]
N x2
N x
On retrouve le dilemme biais-variance o, pour N fix, la variance et le biais voluent en sens
b {b
p(x)}
inverse de x.
12
1.1.3
Estimation de la fonction de rpartition
La fonction de rpartition F se dfinit de la faon suivante :

Z x
p(u)du
F (x) = P [X < x] =
Pour estimer la fonction de rpartition F, on considre une autre variable indicatrice Z(t) telle
que :
Z(t) = 1 si X(t) x
Z(t) = 0 si X(t) > x
Alors, on a :
E {Z(t)} = F (x)
Lestimateur de la fonction de rpartition est un estimateur de moyenne dans le cas dchantillons
indpendants. Sa variance est telle que :
n
o
1
var Fb(x) F (x) (1 F (x))
N
Cet estimateur est non biais et possde une variance qui est minimale pour :
F (x) 0 et F (x) 1
Remarquons que les estimateurs de la densit de probabilit et de la fonction de rpartition ont
tous les deux une variance en 1/N.
1.1.4
Travail eectuer
Estimation de la densit de probabilit

Gnrer N = 1000 ralisations dune variable alatoire gaussienne centre rduite. En faire
lhistogramme sur Nc = 300 classes.
Comment obtient-on lestimateur de densit de probabilit partir de lhistogramme? Observer lestimation de sa densit de probabilit. Etudier linfluence du nombre de classes,
nombre dchantillons constant, sur la variance de lestimateur (en utilisant la fonction superposition des tracs).
Observer la densit de probabilit estime dune variable alatoire uniformment rpartie.
Noter les valeurs des moments estims et comparer la thorie.
Le programme moy_e_dp, indpendant du logiciel, permet dtudier les proprits de lestimateur
de densit de probabilit en termes de biais et de variance en fonction du nombre de classes de
lhistogramme et du nombre dchantillons. La moyenne de lestimateur ainsi que sa variance
1.2. Loi des grands nombres
13
sont estimes sur un nombre de ralisations quil est galement possible de modifier. Le programme trace la moyenne (en vert), la moyenne plus ou moins lcart type (en magenta) ainsi
que la densit de probabilit thorique (en rouge). Observer les proprits de lestimateur en
fonction du nombre de classes Nc et du nombre de ralisations de la variable alatoires N .
Estimation de la fonction de rpartition
Dterminer la fonction de rpartition de la gaussienne. Observer le biais et la variance de
cet estimateur (utiliser la fonction superposition des tracs pour superposer plusieurs
ralisations.
Observer la fonction de rpartition dune variable alatoire uniformment rpartie sur [0, 1].
1.2
Loi des grands nombres
Le traitement grands nombres permet de tracer en fonction de N lestimateur de la moyenne

de X calcul sur N chantillons :
N
X
cN = 1
x(k)
M
N
k=1
Observer le comportement de cet estimateur en fonction de N pour des variables alatoires

gnres suivant les distributions suivantes :
loi normale centre rduite,
loi uniformment rpartie,
loi exponentielle
loi de Cauchy.
1.3
1.3.1
Thorme Central-Limite - Test de Kolmogorov

Thorme central-limite
Soient Xk , k = 1, , K, des variables alatoires indpendantes et de mme loi, de moyenne m
et de variance 2 . On pose
SK
1 X
=
Xk
K
k=1
alors
lim
K+
SK m
en loi
= N (0, 1)
cest--dire que lon sait construire une loi normale centre rduite partir de variables alatoires de nimporte quelle loi indpendantes et de mme loi. Ce thorme est trs puissant et
14
trs important. Il va permettre de traiter de nombreux problmes en faisant lhypothse de loi

gaussienne pour le phnomne observ.
1.3.2
Test de Kolmogorov
Le test de Kolmogorov permet de dcider (avec une marge derreur donne) si une variable
alatoire X suit une loi de fonction de rpartition F (x). Le droulement du test est le suivant :
1) Estimation de la fonction de rpartition Fb(x) partir de K observations xk k = 1, ..., K
de la variable alatoire X.
2) Recherche de lcart maximum max entre Fb(x) et F (x).
3) Pour un risque de premire espce donn ( 1% ou 5% gnralement) vrifier que la valeur
donne dans la table de Kolmogorov est suprieure max . Si cest le cas, on dcide que X suit
la loi de fonction de rpartition F (x) avec un risque derreur %.
Rappelons la notion de risque lorsquon fait un test dhypothse H0 contre H1 : 2 risques sont
dfinis :
= P [rejeter H0 sachant H0 vraie]
= P [accepter H0 sachant H0 fausse]
risque de non dtection

risque de fausse alarme
Le test de Kolmogorov, trs simple mettre en oeuvre, ne prend malheureusement en compte

que le risque de premire espce .
1.3.3
Travail eectuer
A laide de la fonction Th.
Central-Limite et du bouton OK, il est possible de gnrer la
somme de K variables alatoires uniformment rparties sur lintervalle [0, 1], X1 ... XK . Soit
SK la variable telle que :
SK =
K
1 X
Xj
K
j=1
Gnrer N = 500 ralisations et choisir un nombre de classes pour lhistogramme Nc = 500. A

laide du thorme central-limite, dterminer la loi de SK quand K tend vers linfini. Observer
la moyenne et la variance ainsi que la densit de probabilit. Expliquer.
A laide du test de Kolmogorov, tudier partir de quelle valeur de N on peut considrer
avoir une loi de Gauss. Le logiciel eectue le test suivant les tapes suivantes :
1 - Fonction Thorme Central-Limite du logiciel :
Gnration et sommation de K squences alatoires (de N chantillons) indpendantes de mme
loi, en loccurence uniformment rparties sur [0, 1] .
Estimation de la fonction de rpartition (sur N c classes) de la somme SK de ces K squences
alatoires.
Estimation de la moyenne m
b et de la variance
b2 de YN .
1.4. Etude de distributions particulires
15
2 - Fonction Kolmogorov normal?

Calcul thorique par le logiciel de la fonction de rpartition de la gaussienne ayant pour moyenne
m
b et pour cart type
b.
Calcul par le logiciel de lcart maximum entre cette fonction de rpartition thorique et la
fonction de rpartition estime de YN .
Comparaison avec lcart maximal donn dans les tables.
1.4
1.4.1
Etude de distributions particulires

Triangle de Galton
Il sagit dun plan inclin hriss de clous rgulirement disposs sur n lignes horizontales, la
ki`eme ligne comportant k clous. Une bille, place au sommet du triangle, oblique, chaque ligne,
vers la droite ou la gauche avec une probabilit
numrotes de
n2
n
2 qui
1
2.
Les billes samassent dans les n + 1 cases
constituent la dernire ligne.
n=4
-n/2
n/2
Triangle de Galton pour n=4

Pour observer la trajectoire dune bille, slectionner dans le menu signal bille de Galton.
Le nombre de lignes est rglable (pour cela modifier le nombre de ralisations). On prendra
n = 50.
On note X la variable alatoire correspondant au point de chute de la bille (indice de la case
de la dernire ligne o la bille est tombe). La bille aboutit dans la ki`eme case en partant de la
gauche si elle a obliqu k fois vers la droite et n k fois vers la gauche. On note Xi la variable
alatoire asssocie chaque rebond
Xi = 0 la bille oblique gauche

Xi = 1 la bille oblique droite
Vrifier que X =
Pn
j=1 Xj
n2 . Quelle est la loi de X +
n
2
16
Grce au menu demosGalton pas pas il est possible de lancer successivement l

billes (cliquer sur OK) et dobserver le nombre de billes tombes dans chacune des cases.
En dduire une approximation de la loi de X pour l grand. Vrifier lhypothse grce au traitement Galton n grand (il permet de faire tomber l billes sans avoir cliquer pour chacune
delle).
Chi-Deux
Soient Xk , k = 1, , K, des variables alatoires indpendantes et de loi N (0, 1). La variable
vK
VK =
K
X
Xk2
k=1
suit une loi du Chi-Deux K degrs de libert. La moyenne et la variance de cette variable
alatoire sont :
E[VK ] = K
var[VK ] = 2K
Gnrer des signaux de lois Chi-Deux de degrs croissant (bouton OK). Comparer les moyennes
et les variances estimes aux valeurs thoriques. Que se passe-t-il lorsque le degr du Chi-Deux
devient trs grand. A partir de quel degr le test de Kolmogorov accepte lhypothse gaussienne?
Box-Muller
Dans le menu distribution, le choix Box-Muller permet de gnrer partir de deux variables
alatoires indpendantes de mmes lois uniformes sur ]0, 1] U et V la variable alatoire X de
la faon suivante :
X=
2 ln U cos (2V )
Quelle est la loi de X ?
1.5
Test de Neyman-Pearson
Soient X1 , ..., Xn n variables alatoires indpendantes de loi N (m, 2 ) avec 2 connue. On

souhaite raliser le test dhypothses :
H0 : m = m0
H1 : m = m1 > m0
Daprs le lemme de Neyman-Pearson, la rgion de rejet de H0 est donne par
f (X1 , ..., Xn |H1 )
> s()
f (X1 , ..., Xn |H0 )
1.6. Annexes
17
o s() est le seuil de dcision, fonction de la probabilit de fausse alarme . Ce test permet de
maximiser la puissance du test = 1 . Dans ce cas, on obtient :
n
X
Xi > ()
i=1
Expliquer intuitivement ce rsultat.

Le programme puissance_th(m0,m1,sigma2,n) renvoie la puissance thorique du test pour
n variables pour les valeurs = 0.01, 0.02, ..., 0.99.
Tester ce programme pour m0 = 0, m1 = 0.5, 2 = 1 et n = 10. Recommencer avec :
m0 = 0, m1 = 1, 2 = 1 et n = 10 ;
m0 = 0, m1 = 0.5, 2 = 2 et n = 10 ;
m0 = 0, m1 = 0.5, 2 = 0.5 et n = 10 ;
m0 = 0, m1 = 0.5, 2 = 1 et n = 20 ;
m0 = 0, m1 = 0.5, 2 = 1 et n = 50 ;
Commenter ces rsultats (utiliser la commande hold on pour superposer les courbes).
Le programme puissance_est(m0,m1,sigma2,n,K) renvoie la puissance du test pour les
valeurs = 0.01, 0.02, ..., 0.99 estime laide de K simulations. Reprendre les questions prcdentes avec K = 50, puis K = 500. Commentaires ?
1.6
Annexes
1.6.1
Demos de la bote outils stats de Matlab
La bote outils Matlab dispose de programmes de dmonstration permettant de gnrer des

valeurs alatoires suivant une loi donne et dobserver leurs distributions et fonctions de rpartition thoriques ou estimes.
Taper disttool ou randtool sous la fentre de commande Matlab.
Pour connatre toutes les fonctions disponibles dans la bote outils statistiques taper help
stats.
1.6.2
Rappels de quelques distributions
Distribution binmiale B (n, p) :

P [X = k] = Cnk pk q nk
(k = 0, 1, ..., n; 0 < p < 1; q = 1 p)
18
Distribution uniforme sur [a, b] : f (x) =

E(X) =
b+a
2
1
ba
pour x [a, b] (0 ailleurs)

V ar(X) =
(b a)2
12
Distribution Normale N (m, 2 ) :
1
(x m)2
f (x) = exp
2 2
2
Distribution du Chi-deux n degrs de libert :

n
x
1
x 2 1 e 2 pour x [0, ] (0 ailleurs)
n
2 ( 2 )
E(X) = n
V ar(X) = 2n
f (x) =
1.6.3
n
2
Table du Test de Kolmogorov
La table suivante donne lcart maximal thorique, kolmo , entre fonctions de rpartition empirique et thorique pour accepter lhypothse H0 avec un risque de 5 ou 1%, en fonction de
N c (nombre de points de calcul de la fonction de rpartition, cest--dire nombre de classes de
lhistogramme). Si lcart maximal mesur max est infrieur kolmo , on accepte lhypothse
H0 avec un risque .
Nc
= 5%
= 1%
0.5633
0.6685
10
0.4087
0.4864
15
0.3375
0.4042
20
0.2939
0.3524
25
0.2639
0.3165
30
0.2417
0.2898
40
0.2101
0.2521
50
0.1884
0.2260
60
0.1723
0.2067
70
0.1597
0.1917
80
0.1496
0.1795
90
0.1412
100
0.1340
Cette table nest valable que pour Nc infrieur 100. Pour Nc plus grand, on sait calculer
P (max |Fb F | > kolmo ) sous H0 (Vladimir N.Vapnik, The Nature of Statistical Learning
1.6. Annexes
19
Theory, Springer Ed. p87 ). Or :

= P [rejeter H0 | H0 vraie] = P (max |FbF | > kolmo ) '
k=+
X
2(1)k1 exp(2k 2 Nc 2komo )
k=1
On mesure tout dabord max = max |Fb F |, cart maximal entre la fonction de rpartition
thorique et la fonction de rpartition estime. On calcule ensuite :

Pmax = P (max |Fb F | > max ) '
k=+
X
2(1)k1 exp(2k 2 Nc 2max )
k=1
Or cette probabilit est une fonction dcroissante de max . Donc si Pmax > , on peut en
dduire que max < kolmo . Il nest donc pas ncessaire de calculer kolmo pour prendre la
dcision. En eet, on accepte lhypothse H0 si :
k=+
X
2(1)k1 exp(2k 2 Nc 2max ) >
k=1
Le calcul de cette probabilit est eectue par le logiciel.
20
Chapitre 2
TP Echantillonnage et Quantification
2.1
Rappels
Lchantillonnage et la quantification permettent lacquisition et le codage dun signal en vue

dun traitement ou dun stockage. Le but de ce TP est dillustrer les notions dchantillonnage
et de quantification et danalyser leurs eets sur dirents signaux.
2.1.1
Echantillonnage
Lchantillonnage consiste reprsenter un signal temps continu s (t) par ses valeurs s (nTe )
des instant multiples de Te , Te tant la priode dchantillonnage. La condition de
Shannon permet dchantillonner un signal sans perte dinformation aucune perte dinformation
si la frquence dchantillonnage fe =
1
Te
est au moins 2 fois suprieure la plus grande frquence
intervenant dans le spectre (rpartition de la puissance du signal en fonction des frquences) du

signal. Si cette condition nest pas respecte on observe un repliement de spectre (voir figure
suivante). On note ce signal chantillonn se (t) :
se (t) = s (t)
+
X
(t nTe )
Dans le plan frquentiel :
Se (f ) = S (f )
+
+
1 X
1 X
n
n
=
f
S f
Te
Te
Te
Te
21
22
Chapitre 2. TP Echantillonnage et Quantification
Spectre du signal
-Fmax
Frquence
Fmax
Spectre du signal chantillonn
-Fe
-Fmax 0
Fmax Fe
2Fe
zone de repliement
Eet de repliement de spectre
Lchantillonnage la priode Te a donc introduit une priodicit du spectre du signal chantillonn, de priode Fe . Lorsquon observe des spectres de signaux chantillonns, on prfre
remplacer les frquences (F ) en Hertz par des frquences normalises par rapport la frquence
dchantillonnage (Fe ) :
fe = F/Fe
Remarquons quaprs chantillonnage, le spectre entier du signal est disponible entre 0 et 1
en frquence normalise. En pratique, par raison de symtrie, on ne sintresse qu la partie
centrale entre 0 et 0.5 en frquence normalise.
2.1.2
Restitution aprs chantillonnage du signal dorigine
Plusieurs mthodes dinterpolation existent pour restituer le signal entre deux points dchantillonnage
successifs. Entre deux points dchantillonnage situs entre kTe et (k + 1) Te , comment reconstituer le signal sb (kTe + )?
Par bloqueur :
x
b(kTe + ) = x(kTe )
0 Te
2.1. Rappels
23
(k-1)Te
kTe
(k+1)Te (k+2)Te (k+3)Te (k+4)Te
Par extrapolateur linaire :

x
b(kTe + ) = x(kTe ) +
(x(kTe ) x((k 1)Te ))

Te
0 Te
Les marches descalier sont remplaces par des rampes.
(k-1)Te
kTe
(k+1)Te (k+2)Te (k+3)Te (k+4)Te
Par interpolateur linaire :

x
b(kTe + ) = x(kTe ) +
(k-1)Te
(x((k + 1)Te ) x(kTe ))

Te
kTe
0 Te
(k+1)Te (k+2)Te (k+3)Te (k+4)Te
Dans la restitution, il y aura un retard de Te : il est ncessaire dattendre la valeur en kTe

pour restituer lintervalle [(k 1)Te , kTe ] .
24
Par filtrage
Aprs chantillonnage, le spectre du signal est :
Se (f ) = Fe
nZ
S(f nFe ),
Fe =
1
Te
Si x (t) est spectre born (FM , +FM ), et si la condition de Shannon est vrifie, on peut
restituer le signal par filtrage de faon ne rcuprer que le spectre dordre 0 :
S(f ) =
1
[Se (f )]n=0
Fe
Donc
S(f ) = Se (f )Hr (f )
o Hr (f ) est le filtre de restitution :
1
si f (FM , FM )
Fe
Hr (f ) = 0 si |f | Fe FM
Hr (f ) =
Hr(f)
-Fe
-FM
FM
Fe
Dans le cas limite o Fe = 2FM :

Hr (f ) =
1
F (f )
Fe e
La rponse impulsionnelle associe scrit :

hr (t) =
sin (Fe t)
Fe t
soit, pour le signal restitu :
s (t) =
X
k
sin
se (kTe )
Te
Te
(t kTe )
(t kTe )
La courbe restitue passe exactement par les valeurs des chantillons.
2.1. Rappels
2.1.3
25
Quantification
La quantification consiste allouer aux chantillons un nombre fini de valeurs damplitude. Deux
grandeurs caractrisent un quantificateur : le nombre de niveaux et la dynamique de codage.
On distingue des quantificateurs uniformes (lcart entre chaque valeur quantifie est constant)
et non uniformes (lcart entre chaque valeur quantifie est variable) comme lillustre la figure
suivante.
sortie quantifie
Arrondi
Dynamique
entre
Quantificateur
uniforme
pas de quantification
sortie quantifie
Quantificateur
non uniforme
(loi logarithmique)
entre
Quantificateurs uniforme et non uniforme
On note N le nombre de valeurs quantifies ou nombre de niveaux du quantificateur. Le

codage ncessite alors b bits, o 2b est la puissance de 2 immdiatement suprieure ou gale
N.
Leet de la quantification revient, en premire approximation sous certaines conditions ralises en pratique (approximation dite de Sheppard) ajouter au signal s (t) un signal derreur
e (t) appel bruit de quantification, non corrl avec s(t). Le pas de quantification doit tre choisi
de faon minimiser lerreur en sortie du quantificateur.
e(t)
s(t)
sq(t)
Bruit de Quantification
26
Quantification uniforme
Le pas de quantification est alors constant et vaut :
q=
pleine chelle du quantificateur

N
On peut distinguer deux types de quantificateurs uniformes : type A et type B.

type A :
type B :
1
1
q s (t) n +
q alors sQ (t) = nq
si n
2
2
1
q
nq s (t) (n + 1) q alors sQ (t) = n +
2
Leurs caractristiques sont donnes sur la figure suivante. Dans les deux cas, le domaine de
quantification est symtrique par rapport zro.
sortie quantifie
Quantificateur
Uniforme type A
2q
q
-3q/2 -q/2
entre
q/2
3q/2
-q
5q/2
-2q
3q/2
-2q
-q
sortie quantifie
q/2
entre
-q/2
Quantificateur
Uniforme type B
2q
-3q/2
-5q/2
Quantificateurs uniformes de type A et B

Lamplitude du signal derreur est comprise entre q/2 et q/2. Quand les variations du signal
sont grandes par rapport au pas de quantification, cest--dire que la quantification est faite avec
susamment de finesse, et lorsque le signal nest pas crt, le signal derreur peut tre suppos
distribution uniforme, de densit de probabilit :
1
q
pe (x) 0
pe (x)
q q
pour x [ , ]
2 2
ailleurs
Dans ce cas, la puissance du bruit de quantification est donne par :

2e =
q2
12
2.1. Rappels
27
Rapport signal bruit de quantification

Considrons la quantification uniforme de type B. La gamme des amplitudes quil est possible
de coder est soumise une double limitation : vers les faibles valeurs, elle est limite par le pas
de quantification q et vers les fortes valeurs par
Nq
2 .
Toute amplitude qui dpasse cette valeur ne
peut tre reprsente et il y a crtage du signal. Il sen suit une dgradation par distorsion
harmonique, par exemple si le signal est sinusodal.
Pour dfinir et calculer le rapport signal bruit de quantification, considrons un signal
sinusodal non crt et occupant la pleine chelle du quantificateur.
On appelle puissance crte Pc dun codeur la puissance du signal sinusodal ayant lamplitude
maximale Am admissible sans crtage :
Nq
Am =
2
1
Pc =
2
Nq
2
Le rapport signal bruit de quantification est alors le rapport de la puissance crte et de la

puissance du bruit de quantification.
3
Pc
= N2
2
e
2
soit en dcibels et pour N = 2b :

Pc
= 20 log10 N + 1.76dB ' 6b + 1.76dB
10 log10
2e
Le rapport signal bruit de quantification est proportionnel N 2 . Doubler N (ajouter un bit)
amliore le rapport signal bruit de 6dB.
Quantification non uniforme : codage non linaire suivant une loi segmente
Le rapport signal bruit de quantification varie fortement avec le niveau du signal dans le cas de
la quantification uniforme. La quantification non uniforme vise maintenir lerreur relative
de quantification constante, quelle que soit lamplitude du signal alors que la quantification
uniforme introduit une erreur absolue maximum de 12 q.
La solution thorique est de faire varier le pas de quantification q proportionnellement
lamplitude de lentre. Or ceci est incompatible avec la ncessit davoir un nombre fini de
niveaux de quantifications. En pratique, la quantification non-uniforme utilise une quantification
uniforme avec compression-extension : on ralise tout dabord une compression de la dynamique
du signal dans laquelle le signal x est transform en un signal y, puis une quantification uniforme
du signal y et finalement une extension (opration inverse de la compression) du signal quantifi.
La compression du signal x est destine amplifier les faibles amplitudes et minimiser
leet des fortes amplitudes. En pratique, deux caractrisitiques de compression (correspondant
28
deux approximations de la loi logarithmique) sont normalises au niveau international par le

C.C.I.T.T. : la loi A et la loi .
- la loi A :
1
1 + Ln (A |x|)
pour
< |x| 1
1 + Ln (A)
A
1
A |x|
pour 0 |x|
y = sign (x)
1 + Ln (A)
A
y = sign (x)
- la loi :
y = sign (x)
Ln (1 + |x|)
Ln (1 + )
pour 1 x 1
Les paramtres de A et dterminent laugmentation de la dynamique du codeur. La valeur

retenue pour A est 87.6, celle retenue pour est 255. Ces lois sont ensuite approches par des
segments de droite, de faon coder sur un nombre fini de bits le signal ainsi distordu. La loi A
est normalise par le C.C.I.T.T. comme une loi 13 segments alors que la loi est normalise
15 segments. La caractristique de compression de la loi A 13 segments est donne sur la figure
suivante. La caractristique fait apparatre 7 segments dans le quadrant positif, on en imagine
autant dans le cadrant ngatif et, les deux segments entourant lorigine tant colinaires, on
obtient bien 13 segments. Les lois A et sont utilises en tlphonie. Elles constituent deux
approximations de la loi logarithmique proposes respectivement en Europe et aux Etats-Unis.
Ces lois permettent dobtenir un rapport signal bruit de quantification constant partir dun
seuil (1/A ou 1/).
S/B en dB
Loi linaire
(12 bits)
Loi linaire
(8 bits)
Loi A
Puissance (dB) du signal
Rapport signal bruit de quantification avec quantification non uniforme

Lintrt de telles lois rside dans la diminution du nombre dlments binaires pour coder
les chantillons, tout en conservant la dynamique du signal dentre (au dtriment du RSB aux
forts niveaux).
Si on prend lexemple de la loi A, pour coder lamplitude du signal distordu par cette loi, il
faut :
2.1. Rappels
29
- 1 bit de signe
- 3 bits indiquant le numro de segments
- n autres bits pour coder le niveau o on se trouve dans le segment.
La dynamique du signal tlphonique est denviron 65dB, ncessitant 12 lments binaires
dans le cas dun codage linaire. Avec la compression de la loi A, le nombre de bits descend
8, ce qui signifie que n = 4 bits susent pour un codage correct du niveau des segments.
Approximation de la loi A laide de 13 segments :
1
alors
y = 16x
si 0 |x|
64
1
1
1
|x|
alors
y = 8x +
64
32
8
1
1
1
|x|
alors
y = 4x +
32
16
4
1
3
1
|x|
alors
y = 2x +
16
8
8
1
1
1
|x|
alors
y =x+
8
4
2
1
1
5
1
|x|
alors
y = x+
4
2
2
8
1
1
1
|x|
alors
y =x+
8
4
2
1
3
1
|x| 1
alors
y = x+
2
4
4
y
1
111
110
101
100
011
0.5
010
001
000
1
64
2.1.4
0
1
1
16
32
1
8
1
4
1
2
Annexe : Test de Kolmogorov
Le test de Kolmogorov permet de dcider (avec une marge derreur donne) si une variable
alatoire X suit une loi de fonction de rpartition F (x). Dans ce TP, il est utilis pour ltude
de la distribution de lerreur de quantification.
Le droulement du test est le suivant :
1) Estimation de la fonction de rpartition Fb(x) partir de K observations xk k = 1, ..., K
de la variable alatoire X. Remarque importante : pour se rapprocher au maximum des
30
conditions de validit du test de Kolmogorov, le nombre de classes pour lestimation de la

fonction de rpartition doit tre pris gal au nombre de points de signal.
2) Recherche de lcart maximum max entre Fb(x) et F (x).
3) Pour un risque de premire espce donn (1% ou 5% gnralement) vrifier que la valeur
donne dans la table de Kolmogorov est suprieure max . Si cest le cas, on dcide que X suit
la loi de fonction de rpartition F (x) avec un risque derreur %. Rappelons la notion de risque
lorsquon fait un test dhypothse H0 contre H1 : 2 risques sont dfinis :
= P [rejeter H0 sachant H0 vraie]
= P [accepter H0 sachant H0 fausse]
risque de non dtection

risque de fausse alarme
Pour expliquer qualitativement le rle jou par les risques et , imaginons la situation suivante
o on eectue la surveillance militaire dun terrain. Prenons lhypothse H0 : il y a des
ennemis et lhypothse H1 : il ny a pas dennemis. Le risque est le risque de croire
quil ny pas dennemis et donc, de ne pas tirer alors quil y a eectivement des ennemis sur
le terrain... Dangereux!... Le risque est le risque de tirer en croyant quil y a des ennemis
alors quil ny a personne. Le second risque est moins important que le premier qui peut avoir
des consquences dsastreuses... Thoriquement, le rsultat dun test nest valable que si lon
prcise les valeurs des 2 risques. Le test de Kolmogorov, trs simple mettre en oeuvre, ne
prend malheureusement en compte que le risque de premire espce .
La table suivante donne lcart maximal thorique, kolmo , entre fonctions de rpartition
empirique et thorique pour accepter lhypothse H0 avec un risque de 5 ou 1%, en fonction
de N c (nombre de points de calcul de la fonction de rpartition, cest--dire nombre de classes de
lhistogramme). Si lcart maximal mesur max est infrieur kolmo , on accepte lhypothse
H0 avec un risque . Cette table nest valable que pour Nc infrieur 100. Pour Nc plus
grand, on sait calculer P (max |Fb F | > kolmo ) sous H0 (Vladimir N.Vapnik, The Nature of
Statistical Learning Theory, Springer Ed. p87 ). Or :
= P [rejeter H0 | H0 vraie] = P (max |FbF | > kolmo ) '
k=+
X
2(1)k1 exp(2k 2 Nc 2komo )
k=1
On mesure tout dabord max = max |Fb F |, cart maximal entre la fonction de rpartition
thorique et la fonction de rpartition estime. On calcule ensuite :

Pmax = P (max |Fb F | > max ) '
k=+
X
2(1)k1 exp(2k 2 Nc 2max )
k=1
Or cette probabilit est une fonction dcroissante de max . Donc si Pmax > , on peut en
dduire que max < kolmo . Il nest donc pas ncessaire de caluler kolmo pour prendre la
dcision. En eet, on accepte lhypothse H0 si :
k=+
X
k=1
2(1)k1 exp(2k 2 Nc 2max ) >
2.2. Travail raliser
31
Le calcul de cette probabilit sera eectu sous Matlab, la somme infinie pouvant tre approche
par une somme de 1 K trs grand (par exemple K = 1000).
2.2
Nc
= 5%
= 1%
0.5633
0.6685
10
0.4087
0.4864
15
0.3375
0.4042
20
0.2939
0.3524
25
0.2639
0.3165
30
0.2417
0.2898
40
0.2101
0.2521
50
0.1884
0.2260
60
0.1723
0.2067
70
0.1597
0.1917
80
0.1496
0.1795
90
0.1412
100
0.1340
Travail raliser
Lancer ech_stat sous Matlab.
2.2.1
Echantillonnage
Echantillonnage et priodisation du contenu spectral.

Gnrer N=100 points dun signal sinusodal de frquence 5000 Hz chantillonn 100 000 Hz.
Dans une priode du sinus, combien y a-t-il de points ?
Eectuer le spectre du signal (cette fonction sera tudie en dtail dans le TP Corrlations
et Spectres).
Expliquer ce que lon obtient, en particulier la prsence des deux pics et leurs frquences
respectives (utiliser le zoom...).
Changer le nombre de points N=1000. Que se passe-t-il ? (utiliser le zoom)
Passer en frquences normalises (bouton changement chelle frquentielle). Retrouver les
valeurs des frquences intressantes.
32
Eet de la frquence dchantillonnage sur les reprsentations temporelle et

spectrale dun signal : le repliement
Reprendre le signal gnr au dpart : N=100 points, sinus de frquence 5000 Hz, chantillonn
100 000 Hz.
Pour direntes frquences dcantillonnage Fe, observez le signal en temporel et en spectral :
Fe= 50 000 Hz
Fe= 20 000 Hz : combien de points par priode du sinus ?
Fe= Frquence de Shannon : combien de points par priode du sinus ? Que se passe-t-il du
point de vue spectral ?
Fe= 7 500 Hz : sur le spectre, do proviennent les pics observs ? (Faire un dessin). La
frquence du sinus reconstitu est de quelle valeur ? Vrifier sur la reprsentation temporelle
que le nombre de points par priode du sinus correspond bien cette frquence.
Fe= 4 000 Hz : mmes questions.
Pour un signal carr

Gnrer N=1000 points dun signal carr de frquence fondamentale de 800 Hz, chantillonn
100 000 Hz.
Visualiser son spectre. Quobservez-vous ?
Changer la frquence dchantillonnage : Fe=10 000 Hz. Observez le spectre. Que se passet-il ? Un tel signal est-il chantillonnable en thorie ?
Pour un fichier signal

Charger le fichier mf_c1_v6 (fixer la frquence dchantillonnage 40 kHz pour la visualisation spectrale). Ce signal correspond une passe dun navire (Moon Fish) au dessus dun
hydrophone de la base du Brusc (prs de Toulon). Visualiser le spectre. A votre avis, le signal
a-t-il t correctement chantillonn ? Pour rpondre correctement la question, il est ncessaire
de visualiser ce spectre en chelle semilogarithmique.
2.2.2
Quantification
Quantification uniforme grossire

Gnrer N=1000 points dun signal sinusodal de frquence 20 Hz chantillonn 2000 Hz.
Quantifier ce signal laide dun quantificateur uniforme 4 niveaux. Pour voir lerreur de
quantification, utiliser traitement. Que se passe-t-il en temporel ?... Bien caler la dynamique !!!
Visualiser le spectre du signal quantifi. Quel est leet dune quantification grossire sur un
signal sinusodal ?
2.2. Travail raliser
33
Considrer une dynamique de 40 pour une quantification uniforme 4 niveaux. Eectuer la

quantification uniforme de type A et de type B du signal sinusodal prcdent. Conclusions?
Quantification uniforme haute rsolution

Vrifier la formule liant SNR et nombre de niveaux de quantification dans le cas haute rsolution
(nombre de niveaux de quantification 2b avec b = 16, 15, 14...) :
3
Pc
= N2
2
e
2
Montrer en particulier que doubler N amliore le rapport signal bruit de 6dB.
Une des hypothses couramment utilise en quantification est de dire que lerreur de quantification est uniformment rpartie sur [q/2, +q/2].
Gnrer N=1000 points dun signal sinusodal de frquence 2100 Hz chantillonn 100 000
Hz. Eectuer la quantification uniforme de ce sinus sur 14 bits (attention la dynamique !). Observer lerreur de quantification puis sa densit de probabilit. En appuyant sur OK, direntes
ralisations du mme processus alatoires sont gnres. Que penser de cette hypothse ?
Pour sassurer de la validit de cette hypothse, on eectue un test de Kolmogorov. Jusqu
quel nombre de bits de quantification peut-on considrer que cette hypothse est valable ?
2.2.3
Restitution
Gnrer N=100 points dun signal sinusodal de frquence 2100 Hz, chantillonn 100 000
Hz. Comparer les mthodes de restitution par bloqueur, interpolateur linaire, extrapolateur et
filtrage tant sur le plan temporel que sur le plan spectral.
Mmes questions lorsque le signal sinusodal est chantillonn 10 000 Hz.
Mmes questions pour un signal carr de frquence fondamentale 2100 Hz, chantillonn
100 000 Hz.
2.2.4
Signal de parole
Charger le fichier signal allo, fichier de parole. Observer le signal en temporel, remarquer le
caractre non stationnaire du signal.
Eectuer la densit de probabilit du signal et retrouver la validit apparente de lhypothse
de loi double exponentielle souvent nonce pour des signaux de parole.
Comparer la quantification uniforme et non uniforme (logarithmique) de ce signal nombre
de niveaux fix 8. En particulier, comparer la puissance de lerreur de quantification.
34
2.2.5
Dmonstrations
Si le TP se droule en salle TSI (B301), des dmonstrations peuvent tre faite par les enseignants
en traitement du signal et en traitement dimage.
Quantification des images

On utilisera le logiciel de traitement dimages et, en chargeant une image noir et blanc, on mettra
en vidence leet de la quantification sur des images.
Quantification de la parole
Sur le poste quip de haut-parleurs, on observera les eets de la quantification et du souschantillonnage sur les signaux de parole.
2.3
Elements de programmation Matlab
2.3.1
Quelques fonctions Matlab utiles dans le TP
Pour connatre le mode dappel de ces fonctions, penser laide en ligne de Matlab : help nom
de la fonction ou lookfor mot.
fft : transforme de Fourier discrte
ifft : transforme de Fourier inverse
interp1q:
interpolation linaire
quantification uniforme du signal.
quantiz:
lin2mu:
conversion la loi .
mu2lin:
conversion inverse (de la loi la loi linaire).
2.3.2
Echantillonnage
Gnration dune sinusode : Pour gnrer le signal obtenu aprs lchantillonnage dune
sinusoide de frquence 20 KHz chantillonne un rythme de 1 chantillon toute les 10 microsecondes, on peut procder soit en utilisant les frquences physiques :
N = 1000;
Te = 10*10^(-6);
fsin = 20*10^3;
temps = (0:N-1)*Te;
signal = sin(2*pi*fsin*temps);
soit en raisonnant directement en frquences normalises :
2.3. Elements de programmation Matlab
35
N = 1000;
Fe = 1/(10*10^(-6));
fsin_norm = 20*10^3/Fe;
indice = 0:N-1;
signal = sin(2*pi*fsin_norm*indice);
Visualisation de la densit spectrale de puissance : Le calcul et le trac de la densit
spectrale de puissance seront dtaills au cours du TP Corrlations et Spectres.
Linstruction de base pour raliser la transforme de Fourier discrte dun signal x est fft(x).
Les commandes suivantes permettent de tracer la densit spectrale de x soit en frquence normalise :
nfft = 2^nextpow2(N); %calcul de la puissance de 2 immdiatement suprieure N.
plot(linspace(0, 1, nfft), abs(fft(x, nfft)).2./nfft).
soit en frquence physique :
plot(linspace(0, Fe, nfft), abs(fft(x, nfft)).2./nfft).
Restitution : la commande y = kron(x, ones(1,Nrestit)); permet de restituer par bloqueur dordre 0 le signal x chantillonn Te. Le signal y rsultant est une version surchantillonne de x de priode dchantillonnage Te/Nrestit.
La commande interp1q permet de raliser linterpolation linaire (taper help interp1q);
La commande interp permet de gnrer un filtre passe-bas dinterpolation.
2.3.3
Quantification
Quantification uniforme : La fonction quantiz permet de raliser directement la quantification sous Matlab.
Etude statitique de lerreur de quantification : La fonction var donne une estimation de
la variance. Pour lestimation de la densit de probabilit et le test de Kolmogorov se reporter
au TP 1.
Quantification non uniforme selon la loi : voir les fonctions lin2mu, mu2lin et auwrite.
36
Chapitre 3
TP Analyse Spectrale - Corrlations

et Spectres
Le but de ce TP est danalyser des estimateurs de la fonction de corrlation et de la densit
spectrale de puissance (DSP).
La fonction de corrlation constitue une mesure de ressemblance entre deux signaux (intercorrlation) ou, si elle est applique un seul signal, cette fonction, appele alors autocorrlation,
mesure le lien qui existe entre les direntes valeurs du signal des instants dirents. Cette
quantit dpend dun paramtre de dcalage temporel (dphasage entre les deux signaux dans le
cas de lintercorrlation) et possde selon la nature des signaux, plusieurs dfinitions.
La densit spectrale de puissance (transforme de Fourier de la fonction dautocorrlation)
reflte la contribution quapporte chaque frquence la puissance moyenne du signal.
3.1
3.1.1
Rappels
Proprits des fonctions de corrlation
La fonction de corrlation se dfinit de direntes faons suivant la classe de signaux laquelle

on sadresse. Nous donnons ci-aprs les direntes dfinitions de la fonction dintercorrlation.
Pour la fonction dautocorrlation, il sut de faire y = x.
Signaux dterministes
Energie finie
Cxy ( ) =
37
x(t)y (t )dt
38
Chapitre 3. TP Analyse Spectrale - Corrlations et Spectres
Puissance finie
1
T + T
Cxy ( ) = lim
Cas particulier des signaux priodiques :

1
Cxy ( ) =
T
x(t)y (t )dt
x(t)y (t )dt
Signaux alatoires
Cxy ( ) = E[x(t)y (t )]
Lorsque y = x on parle de fonction dautocorrlation (Cxx ( ) = Cx ( )).
Proprits de la fonction dautocorrlation

Parit : Cx ( ) = Cx ( )
Maximum en zro : |Cx ( )| Cx (0)
Puissance moyenne du signal = Cx (0)
Proprits de la fonction dintercorrlation

( ) = C ( )
Symtrie hermitienne : Cxy
xy
Majoration : |Cxy ( )| 12 (Cx (0) + Cy (0))
3.1.2
Algorithmes de calcul des fonctions de corrlation
Pour estimer la fonction dautocorrlation ou dintercorrlation, on utilise deux types destimateurs.

Dans une premire approche, la fonction dautocorrlation est estime comme la valeur moyenne
de x(n)x(n + k) . Si on ne dispose que de N chantillons du signal x(n) (rel), Cxx (k) ne peut
tre estime qu partir de N k valeurs :
bb (k) = 1
C
N
Nk1
X
n=0
x(n)x(n + k) 0 k N 1
(3.1)
Lorsque k tend vers N 1 peu de termes interviennent dans le calcul de la moyenne alors que
le terme de normalisation reste gal
1
N.
Cela a pour consquences dintroduire un biais dans
lestimation : la corrlation est pondre par une fentre triangulaire. Ce premier estimateur
bb (k) est donc biais.
C
Pour liminer ce biais, un second estimateur peut tre dfini de la faon suivante :
bnb (k) =
C
Nk1
X
1
x(n)x(n + k) 0 k N 1
N k n=0
(3.2)
3.2. Exemples dutilisation des fonctions de corrlation
39
Le choix entre ces deux estimateurs se porte normalement sur lestimateur non biais mais,
pour les valeurs k proches de N, la variance de cet estimateur augmente considrablement ce
qui ntait pas le cas du prcdent estimateur. Dautre part, il peut tre intressant de choisir
lestimateur biais ce qui est quivalent prendre la fonction de corrlation multiplie par une
fentre temporelle triangulaire.
Remarque : Le calcul de ces deux estimateursl requiert de lordre de
N2
2
multiplications et
additions. Pour rduire le cot calculatoire, on peut utiliser un algorithme base de transforme
de Fourier rapide (FFT). En eet, les deux estimateurs ci-dessus reprsentent un facteur
multiplicatif prs la mme opration de convolution discrte x(k) x(k). Pour calculer cette
convolution, on peut utiliser le fait que
T F (x(k) x(k)) = X(f )X (f ) = |X(f )|2

avec X(f ) la transforme de Fourier de x(n). Ainsi, les sommes temporelles coteuses
en cot
h
i
calculatoire contenues dans (3.1) et (3.2) peuvent tre remplaces par F F T 1 |F F T (x(n))|2 ,
ce qui rend le calcul plus rapide.
3.2
Exemples dutilisation des fonctions de corrlation
Nous donnons ci-aprs deux exemples dutilisation des fonctions de corrlation.
3.2.1
Dtection dun signal priodique noy dans un bruit
Par dtection, on entend dtection de prsence. Il ne sagit pas de retrouver la forme du signal
priodique mais de dtecter sa prsence, de savoir si ce signal existe ou non.
Soit x(n) le signal priodique de priode inconnue noy dans un bruit centr b(n) :
y(n) = x(n) + b(n)
La fonction dautocorrlation du signal y(n) scrit :
Cyy (k) = Cxx (k) + Cxb (k) + Cbx (k) + Cbb (k)
Si le bruit b(n) est indpendant du signal x(n), les intercorrlations Cxb (k) et Cbx (k) sont nulles
pour tout k. Si de plus la densit spectrale du bruit est absolument continue, on a :
lim Cbb (k) = 0
k+
Do :
Cyy (k) Cxx (k) lorsque k +
40
En calculant la fonction dautocorrlation de y(n), on va voir apparatre celle du signal priodique

x(n) aux grandes valeurs de k.
3.2.2
Identification dun filtre
Considrons un filtre ou dune faon plus gnrale un systme linaire invariant dans le temps,
de rponse temporelle h(t) et de rponse frquentielle H(f ) (fonction de transfert). La relation
temporelle qui relie la sortie y(t) lentre x(t) est la suivante :
y(t) = x(t) h(t) =
x(u)h(t u)du
On obtient la mme relation sur les fonctions de corrlation :

Cyx ( ) = Cxx ( ) h( )
o Cyx ( ) est lintercorrlation entre la sortie et lentre du filtre.
Si le signal x(t) est un bruit blanc, cest--dire Cxx ( ) = ( ), on a :
Cyx ( ) = ( ) h( ) = h( )
Lintercorrlation entre la sortie et lentre du filtre correspond la rponse temporelle du filtre.
La transforme de Fourier de lintercorrlation permet dobtenir la rponse frquentielle. Une
autre possibilit consiste estimer le spectre de la sortie du filtre :
Y (f ) = X(f )H(f )
La fonction de transfert H(f ) permet une interprtation plus visuelle de la nature du filtre
(passe-haut, passe-bas, passe-bande).
3.3
3.3.1
Analyse spectrale
Dirents estimateurs
Il existe 2 grandes familles destimateurs de la densit spectrale de puissance (DSP) :

les mthodes du PERIODOGRAMME
les mthodes du CORRELOGRAMME
3.3. Analyse spectrale
41
Estimation du spectre de puissance
x(t )
FFT
1
2
S$ x ( f ) =
X( f )
N
X( f )
Signal temporel
Priodogramme
C x ()
FFT
S$ x ( f ) = FFT (C x ( ))
Corrlogramme
Autocorrlation
(paire)
Ces deux catgories de mthodes mettent en oeuvre la transforme de Fourier discrte dfinie
par :
X(f ) =
N1
X
x(n)ei2f n
n=0
Lalgorithme de transforme de Fourier discrte rapide, FFT, requiert un nombre de points

frquentiels Nf gal une puissance de 2 afin dobtenir la meilleure ecacit en termes de
temps de calcul. Le spectre Sx (f ) est calcul aux frquences normalises :
f=
3.3.2
k
Nf
k = 0, , Nf 1
Intrt du zero-padding
A titre dexemple, la transforme de Fourier discrte (TFD) dun sinus

x(n) = Acos(2f0 n + ) n = 0, , N 1
donne :
X(f ) =
sin((f f0 )N )
ei((f f0 )(N1))
sin((f f0 ))
sin((f
+ f0 )N )
i((f +f0 )(N1)+)
+e
sin((f + f0 ))
A
2
Ceci correspond 2 noyaux de Dirichlet (ressemblant 2 sinus cardinaux), lun centr en f0 ,

lautre en f0 . Chaque motif correspondant un noyau de Dirichlet se prsente sous forme
dun lobe principal et de lobes secondaires, avec une pseudo priode de
1
N,
et le rsultat du
calcul est reprsent sur seulement N points. Ceci donne 1 point par pseudo-priode et cest
trop peu pour que, visuellement, on voit clairement apparatre le noyau de Dirichlet centr sur
la frquence dintrt. Afin damliorer la visualisation du rsultat de la TFD, il faut prendre
42
un nombre de points en frquence Nf grand devant le nombre dchantillons du signal et faire

ainsi du "zero-padding".
3.3.3
Corrlations thoriques de signaux particuliers
Sinusode phase alatoire

x(n) = Acos(2f n + ) Cx (k) =
Bruit blanc
3.3.4
A2
cos(2fk)
2
Cx (k) = 2 (k)
(3.3)
(3.4)
Priodogramme dune sinusode bruite et estimation du

SNR
On considre N chantillons dune sinusode de frquence f0 et damplitude A perturbe par un

bruit blanc de puissance 2 . Le priodogramme de ce signal est :
Sx (f ) =
avec sinc(x) =
1
A2
|X(f )|2 2 +
N (sinc2 (N (f f0 )) + sinc2 (N (f + f0 )))
N
4
(3.5)
sin(x)
x .
Lamplitude du pic la frquence f = f0 correspond donc
A2
2
4 N +
ce qui permet destimer
le rapport signal--bruit (SNR) :

SN R = 10log
P uissance de la sinusoide =
P uissance du bruit = 2
A2
2
(3.6)
En eet, lobservation de la DSP obtenue permet de mesurer le niveau de bruit moyen 2

ainsi que lamplitude du sinus A, comme lillustre la figure suivante.
A2
N + b2
4
b2
Exemple de DSP estime dune sinusode bruite en chelle linaire.
3.4. Travail eectuer
3.4
43
Travail eectuer
Lancer analyse_spectrale sous Matlab.

Le but de ce TP est deectuer de lanalyse spectrale par Transforme de Fourier (TF) et de
comprendre les intrts, inconvnients et nature des paramtres des estimateurs existants.
Charger successivement les fichiers Eng8_e5.txt et Eng8_e11.txt (dans le logiciel analyse_spectrale,
faire "choix du signal", "fichier") : ces signaux reprsentent des essais de fatigue de trains picyclodaux. La roue de petite vitesse est menante. Elle possde 56 dents et a une frquence de
rotation de 12,43 Hz. La roue de grande vitesse possde 15 dents et une frquence de rotation
de 46,4 Hz. On dispose de mesures issues de capteurs en position horizontale, la frquence
dchantillonnage est de 6365 Hz. Le premier fichier correspond des mesures eectues le
5me jour et le deuxime fichier des mesures eectues le 11me jour. Le 5me jour, ltat des
dentures est sain tandis que le 11me jour, il est dgrad : la dent 51 est use 100%, les dents
24, 26, 27, 31 et 55 sont uses 75% et la dent 56 est use 15%. La question pose est de
savoir si lanalyse spectrale ou le calcul de corrlation nous permettrait de mettre en vidence
lusure des dents.
Cest pourquoi, dans ce TP, nous allons tudier dirents estimateurs de la fonction dautocorrlation
et de la Densit Spectrale de Puissance (DSP). Mais avant de travailler sur des signaux rels,
nous allons valider les dirents estimateurs sur des signaux tests.
Remarque : toutes les frquences sont donnes en frquences normalises.
Autre remarque : il reste encore des "bugs" dans ce logiciel. En cas de comportement bizarre
(refus dacher un signal qui devrait devoir sacher), utiliser le bouton "reset".
3.4.1
Autocorrlations
Autocorrlation dun sinus

Gnrer 50 chantillons dune sinusode de frquence 0.134 (phase uniformment rpartie sur
[0, 2]). Observer les estimations de son autocorrlation en utilisant le bouton "ok" pour gnrer
chaque fois une nouvelle ralisation de signal.
Pour lestimateur biais, quelle est lallure du biais ?
Pour lestimateur non biais, dans quelle partie de lautocorrlation la variance est-elle la plus
importante ? Expliquer.
Retrouver les caractristiques du signal (puissance et frquence).
Rappel : lexpression thorique de la corrlation est donne par (3.3).
44
Autocorrlation dun bruit blanc

Gnrer un bruit blanc (N = 100) de variance 2 = 4. Utiliser le bouton "ok" pour voir des
ralisations direntes de ce processus alatoire. Observer les estimations de son autocorrlation, utiliser le bouton "ok" pour visualiser direntes estimations correspondant direntes
ralisations du signal.
Des 2 estimations (biaise ou non biaise), laquelle parat la plus satisfaisante ? Expliquer.
Retrouver les caractristiques du bruit (remarquer que, pour plus de facilit, la fentre "Informations et mesures" fournit la valeur et lindice du maximum de la fonction).
Rappel : lexpression thorique de la corrlation est donne par (3.4).
Autocorrlation dun sinus bruit

Gnrer 500 chantillons dune sinusode bruite (bruit blanc additif) de frquence 0.0134
(SN R = 7 dB). Pour cela, sachant que les sinusodes gnres sont damplitude A = 2
(i.e. de puissance unit), rgler le niveau de bruit ncessaire en utilisant la dfinition du SNR
donne par (3.6). Observer le signal et lestimation biaise de son autocorrlation (ne pas oublier dutiliser le bouton "ok" pour gnrer plusieurs ralisations du processus). Montrer que
lestimation de lautocorrlation permet de retrouver les caractristiques du signal : prsence
dun signal priodique, puissances du sinus et du bruit.
3.4.2
Estimation spectrale
Priodogramme
Gnrer une sinusode (N = 128, f = 0.134). Observer le priodogramme en utilisant une taille
de FFT = 128 (en chelle linaire et logarithmique).
Remarque : On rappelle que pour le priodogramme, le seul paramtre fixer est la taille de
la FFT.
Recommencer en faisant du zero-padding avec une taille de FFT de 4096. Expliquer les
dirences observes : lexpression thorique de la DSP est donne par (3.5) et des rappels sur
le zero-padding sont fournis au paragraphe 3.3.2). Retrouver les caractristiques du signal :
frquence et amplitude du sinus.
Charger le fichier quisuisje0.txt et cliquer sur le bouton "signal" pour observer le signal en
temporel. En faire lanalyse spectrale en utilisant le priodogramme : retrouver les caractristiques du signal. Ne pas oublier dobserver lestimateur spectral obtenu en chelle log ! Que
reprsente la DSP la frquence f = 0 ?
45
Priodogramme modifi
Pour une sinusode (N = 100, f = 0.25), analyser et comparer les eets des direntes fentres en
utilisant le priodogramme modifi. Classer-les en fonction de leur pouvoir rduire lamplitude
des lobes secondaires (utiliser la reprsentation en chelle logarithmique et le bouton hold ).
Remarque : Pour le priodogramme modifi, les paramtres choisir sont le type de fentre
dapodisation et la taille de la FFT.
Analyse de signaux sinusodaux : intrt du priodogramme modifi

Charger le fichier quisuisje1.txt et cliquer sur le bouton "signal" pour observer le signal en
temporel. En utilisant le priodogramme modifi et les direntes fentres dapodisation, donner
les caractristiques de ce signal.
Refaire la mme analyse sur le fichier quisuisje2.txt.
Corrlogramme
Gnrer une sinusode (N = 100, f = 0.134). Observer le corrlogramme avec lestimateur de
la corrlation non biaise en utilisant une taille de FFT = 4096. Comparer avec le corrlogramme utilisant la corrlation biaise. Quen concluez-vous ? Comparer corrlogramme biais
et priodogramme (utiliser le bouton "hold").
Remarque : Pour le corrlogramme, les seuls paramtres choisir sont la nature de lestimation
de la corrlation (biaise ou non biaise) et la taille de la FFT.
Blackman-Tukey
Pour le mme signal, utiliser le corrlogramme de Blackman-Tukey.
Remarque : Pour la mthode de Blackman-Tukey, il est prconis de lutiliser sur la corrlation
biaise afin de garantir une DSP positive. Toutefois, il est intressant de voir ce quon obtient
en utilisant lestimation de la corrlation non biaise. Cest pourquoi les paramtres choisir
pour cette mthode sont la nature de lestimation de la corrlation (biaise ou non biaise),
le type de fentre dapodisation, la taille de la fentre dapodisation (attention, applique
sur la corrlation qui est de taille double par rapport au signal) et la taille de la
FFT.
Faire varier tous les paramtres afin den observer leur eet (en particulier sur la positivit de
la DSP estime et sur les lobes secondaires des fentres spectrales quivalentes) : choix du type
de corrlation, choix de la fentre, choix de la taille de la fentre.
46
Gnrer 1000 chantillons dune sinusode bruite (bruit blanc additif) de frquence 0.134
(SN R = 7 dB). Comparer lanalyse spectrale obtenue avec le priodogramme modifi et avec
Blackman-Tukey, utilis par exemple en prenant la corrlation biaise et la fentre de Blackman

(permettant davoir une DSP positive). Que remarquez vous sur la variance de lestimation
spectrale ?
Priodogramme de Bartlett
Gnrer un bruit blanc (N = 4000) de variance 2 = 5. Observer son priodogramme en utilisant
une taille de FFT = 8192. Comparer le rsultat obtenu lexpression thorique de la DSP dun
bruit blanc. En remarquant que la fentre "informations et mesures" fournit la valeur moyenne
de la DSP ainsi que sa variance, expliquer la dirence entre la thorie et lestimation obtenue
de cette DSP. Ne pas hsiter utiliser le bouton "ok" pour gnrer plusieurs ralisations du
processus alatoire.
Pour ce mme type de signal, utiliser le priodogramme de Bartlett en choisissant des tailles
de fentre direntes : 1000 et 100 par exemple. Quobserve-t-on ?
Remarque : Pour le priodogramme de Bartlett, les paramtres choisir sont la taille de la
fentre danalyse et la taille de la FFT.
Lutilisation de la mthode de Bartlett se rvle intressante sur le bruit blanc mais elle
saccompagne dun inconvnient. Afin de le mettre en vidence, recommencer la mme exprience
sur un signal constitu dun sinus bruit (N = 4000, f = 0.2) en conservant la mme variance de
bruit. Faire varier la taille de la fentre danalyse de 4000, 1000 puis 100. Pour ce type de
signal, le logiciel nache pas la valeur moyenne et la variance du spectre mais la variance peut
tre observe en se plaant en chelle logarithmique. En dduire les avantages et inconvnients
de la mthode de Bartlett.
Priodogramme de Welch
Pour amliorer la mthode de Bartlett, on peut utiliser le priodogramme de Welch qui permet
dutiliser une fentre danalyse et un recouvrement entre les tranches.
Remarque : Pour le priodogramme de Welch, les paramtres choisir sont la taille et le type
de fentre danalyse, le recouvrement (indiqu en % ici) des tranches et la taille de la FFT.
Gnrer un bruit blanc de N = 50000 chantillons (beaucoup dchantillons !), de variance
2
= 5. Choisir une fentre danalyse rectangulaire de taille 1000 et une taille de la FFT de
8192.
3.5. Programmation Matlab (pour ceux qui le souhaitent)
47
Pour un recouvrement de 0%, acher plusieurs ralisations successives du priodogramme de

Welch en notant la variance du spectre (ou en faire une moyenne loeil !).
Refaire de mme pour un recouvrement de 30%, 50% et 80%. Quel est lintrt du recouvrement ? Son inconvnient ?
Pour un recouvrement fix 30% par exemple, faire varier la nature de la fentre danalyse
(pas sa taille). Quobserve-t-on ? Que gagne-t-on utiliser des fentres autre que la rectangulaire
? Ne pas oublier leur inconvnient...
Rsumer linfluence (avantages et inconvnients) des divers paramtres de la mthode de
Welch (la plus utilise) :
taille de la FFT :
type de fentre danalyse :
taille de la fentre danalyse :
pourcentage de recouvrement :
3.4.3
Analyse spectrale de signaux rels
Aprs avoir tudi les estimateurs de corrlation et de DSP, analyser les signaux dengrenage cits
au dbut du TP. Donner une description (sommaire) de ces signaux. A votre avis, pourrait-on
distinguer ces deux signaux et donc signaler lusure des engrenages partir de lanalyse spectrale
?
3.5
Programmation Matlab (pour ceux qui le souhaitent)
3.5.1
Quelques fonctions Matlab utiles dans le TP
randn :
gnration dun bruit blanc normal centr et de variance unit
fft : transforme de Fourier discrte

fftshift :
recentre le spectre
ifft :
transforme de Fourier inverse
xcorr :
autocorrlation
specgram :
psd :
spectrogramme (utile pour le calcul du priodogramme cumul)
densit spectrale de puissance Pour connatre le mode dappel de ces fonctions, penser
laide en ligne de MATLAB : help nom de la fonction ou lookfor mot
48
3.5.2
Autocorrlations
Elle peut tre calcule de la faon suivante :

for k=0:N-1,
autocorb(k+1) = fact*signal(1:N-k)*signal(k+1:N);
end
avec fact =1/N pour lautocorrlation biaise ou 1/(N-k) pour lautocorrlation non biaise.
Il existe sous Matlab la fonction xcorr, avec les options biased et unbiased.
Pour la gnration des dirents signaux (sinus, carr, bruit blanc gaussien) voir les indications
Matlab du TP1.
3.5.3
Estimations spectrales
Transforme de Fourier Discrte

Linstruction de base pour raliser la transforme de Fourier discrte dun signal x est fft(x).
Lalgorithme de transforme de Fourier rapide est utilis par Matlab si et seulement si la longueur
du vecteur x, Nech, est une puissance de 2. Si ce nest pas le cas, on peut utiliser la technique
du zero-padding pour sy ramener :
nfft = 2^nextpow2(Nech); %calcul de la puissance de 2 immdiatement suprieure
Nech
fft(x,nfft);
Exemple : Gnration de 512 points dun sinus de frquence normalise f0 = 0.2
x=sin(2*pi*f0*(0:511));
Calcul de sa FFT et trac du module (fonction abs) de sa FFT :
fft_de_x = fft(x, nfft);
Si laxe des x, nest pas spcifi dans la commande plot :
plot(abs(fft_de_x))
cet axe porte alors par dfaut les indices des lments du vecteur trac (ici 0,1,2,...
nfft-1).
Pour interprter correctement le trac de la Transforme de Fourier Discrte, il est ncessaire de graduer laxe des x, par exemple en frquences normalises par rapport la frquence
dchantillonnage :
axe_des_x = linspace(0, 1, nfft)
plot(axe_des_x, abs(fft_de_x))
Lachage du spectre peut tre galement recentr autour de zro laide de la commande
3.5. Programmation Matlab (pour ceux qui le souhaitent)
49
fftshift, laxe des x doit alors tre modifi :

plot(linspace(0.5, 0.5, nfft), fftshift(abs(fft_de_x)))
Priodogramme et priodogramme de Welch

La densit spectrale de puissance (DSP) peut tre estime laide du priodogramme :
den_puis = abs(fft(x, nfft)).2./Nech;
On peut galement utiliser le priodogramme de Welch. Cela consiste :
couper le signal en tranches de Nt points (Nt doit tre videmment plus faible que la
longueur totale du signal) avec recouvrement possible des tranches.
faire une FFT de chacune des tranches (avec du zero-padding et une fentre dapodisation) et
prendre le module au carr de la FFT,
faire la moyenne de toutes les FFT en module au carr
Ces oprations peuvent tre ralises avec la fonction pwelch (voir le help sous Matlab).
Corrlogramme
Il sut dutiliser successivement xcorr et fft.
50
Chapitre 4
Filtrage Numrique
Ce TP est consacr ltude des filtres Rponse Impulsionnelle Finie (RIF) et des Filtres
Rponse Impulsionnelle Infinie (RII). Dans le cas des filtres RIF, tout chantillon du signal en
sortie est la somme pondre dchantillons du signal en entre. Les filtres RIF sont frquemment
dsigns par le terme de filtres non-rcursifs, car ils ne prsentent pas de boucle de raction de
la sortie vers lentre. Ils peuvent tre synthtiss directement par un dveloppement en srie de
Fourier du gabarit idal. Le rsultat obtenu peut tre ensuite optimis grce la mthode des
moindres carrs ou lalgorithme de Remez. La deuxime partie du TP est consacre ltude
des filtres RII ou filtres rcursifs. La synthse de ces filtres sappuie sur les fonctions modles du
filtrage analogique (Tchebychev, Butterworth,...) par lintermdiaire de la transforme bilinaire,
transformation conforme permettant de passer du plan numrique au plan analogique.
Le but du TP est de comparer pour un gabarit donn les rsultats des direntes synthses.
4.1
Filtre Rponse Impulsionnelle Finie (RIF) :

rappels
4.1.1
Dfinition
Ce sont des systmes rponse impulsionnelle finie, de fonction de transfert H(z), dont les
coecients h(k) sont tels que :
h (k) 6= 0 pour k [0, N 1]
h (k) = 0 sinon
On obtient lexpression de la fonction de transfert dans le plan des z :
H (z) =
N1
X
k=0
51
h (k) z k
52
Chapitre 4. Filtrage Numrique
Le caractre non rcursif apparat clairement sur lquation de rcurrence liant lentre et la
sortie du filtre :
y (n) =
N1
X
k=0
h (k) x (n k)
La mthode de synthse des filtres RIF permet dassurer une phase linaire, caractristique
recherche dans de nombreuses applications. Ceci implique une symtrie de la rponse impulsionnelle. En eet, on veut :
H (f ) = R (f ) ej(f )
avec R (f ) R et (f ) = 2f pg . La rponse impulsionnelle dun tel filtre scrit :
Z +
Z +
j(f ) j2f t
R (f ) e
e
df =
R (f ) ej2f (t pg ) df
h (t) =
On dcompose R (f ) en la somme dune partie paire Rp (f ) et dune partie impaire Ri (f ). La

rponse impulsionnelle h (t) tant relle (dans le cas dun signal rel le spectre damplitude est
pair et la phase est impaire), on a :
h ( pg + t) = 2
Rp (f ) cos (2f t) df = h ( pg t)
Cette relation fait donc apparatre la symtrie de la rponse impulsionnelle par rapport au point
t = pg de laxe des temps.
4.1.2
Synthse par la mthode de la fentre
On se donne un gabarit frquentiel H (f ) respecter. On calcule la rponse impulsionnelle du

filtre recherch par transforme de Fourier inverse de ce gabarit. Il sagit alors de faire une
troncature afin de garder un nombre fini N dlments qui seront les coecients du filtre, puis
deectuer un dcalage afin de rendre le filtre causal cest--dire physiquement ralisable.
Etapes successives de la synthse :
1) Dfinir un gabarit frquentiel en frquences normalises
2) En faire un dveloppement en srie de Fourier
Z 1
2
e
e
H ej2f ej2f k dfe
h (k) =
12
3) Multiplication par une fentre temporelle W (k) de longueur N (ordre du filtre) avec W (k) = 0
pour |k| > N,
4) Raliser un dcalage (translation) afin de satisfaire la condition de causalit. La valeur du

temps de propagation de groupe est :
pg

d
fe
N 1
1
=
=
e
2 df
2
4.2. Filtre Rponse Impulsionnelle Infinie (RII) : rappels
4.1.3
53
Optimisation de la synthse obtenue
Uen fois la synthse obtenue par la mthode classique de synthse RIF par fentre, on peut
optimiser le filtre obtenu.
La premire mthode doptimisation consiste minimiser au sens des moindres carrs la
distance entre le gabarit H (f ) dsir et le gabarit du filtre obtenu par la mthode prcdente.
Lobjectif de la seconde mthode doptimisation est dobtenir la meilleure approximation du
gabarit H (f ) prsentant des ondulations damplitude constante. Elle utilise une technique
itrative : lalgorithme de Remez.
4.2
Filtre Rponse Impulsionnelle Infinie (RII) :

rappels
4.2.1
Dfinition
Ces systmes sont caractriss par des rponses impulsionnelles de dure infinie : les coecients
h(k) sont non nuls sur lintervalle [0, +[ . Ceci est ralis par la prsence de ples dans la
fonction de transfert du filtre :
H (z) =
p1
X
bk z k
k=0
m1
X
ai z i
i=0
Cela se traduit par lquation suivante :

y (n) =
m1
X
i=1
ai y (n i) +
p1
X
k=0
bk x (n k)
La condition de stabilit impose que les ples de H(z) soient lintrieur du cercle unit. La
rponse impulsionnelle infinie permet dobtenir un filtrage plus slectif quun filtre RIF pour
une quantit de calcul infrieure. La linarit de la phase est, en thorie, impossible. Cependant
on peut lobtenir approximativement dans une bande limite.
Synthse
On calcule tout dabord le filtre analogique passe-bas qui lui correspond, puis on synthtise ce
filtre analogique. Enfin, on revient au filtre numrique par la transformation bilinaire. Les
mthodes classiques mettent en oeuvre des transformations partir des quations de filtres
analogiques connus (Butterworth, Tchebychev...). Le passage entre les dirents types de filtre (Passe-Bas, Passe-Haut, Passe-Bande et Coupe-Bande) peut se faire dans le plan des p
54
(analogiques) ou dans le plan des z au cours de la transformation p z. La mthode utilise,
la transformation bilinaire, fait correspondre le plan des z au plan des p suivant la relation :
p=
2 1 z 1
Te 1 + z 1
La synthse dun filtre RII se fait donc selon les tapes suivantes :
1) Synthse du filtre analogique passe-bas de pulsation de coupure c :
|H ()|2 =
Butterworth dordre n :
Tchebychev dordre n :
1+
1
2n
1
|H ()|2 =
2
1 + Tn2 ()
Tn () = cos (nAr cos ()) si || 1

avec
Tn () = ch (nArch ()) si || > 1
2) Transformations frquentielles en vue dobtenir les diverses familles de filtres : passe-haut,

passe-bande, coupe bande.
3) Transforme bilinaire du filtre obtenu en 2).
4) Calcul de la fonction de transfert partir des coecients.
5) Visualisation : DSP, temps de propagation de groupe, rponse impulsionnelle, rponse indicielle. Remarque : la transformation bilinaire introduit une dformation des frquences :
fn
Fe
tan
f =
Fe
avec f frquence analogique, fn frquence numrique et Fe la frquence dchantillonnage.
4.3
Travail eectuer
Lancer filtnum sous Matlab.
4.3.1
Introduction
Le TP permet de synthtiser des filtres numriques RIF et RII. Les filtres RIF peuvent tre
calculs par un dveloppement en srie de Fourier avec direntes fentres (rectangulaire, triangulaire, de Hamming et de Kaiser) et par deux mthodes doptimisation (par les moindres
carrs et par Remez). Les mthodes proposes pour les filtres RII utilisent les quations donnant
des filtres analogiques connus (de Butterworth, de Chebychev I, de Chebyshev II et elliptique).
Lintrt sera port tout dabord sur les caractristiques des filtres RIF, puis sur celles des filtres
RII. Chacune de ces tudes sera suivie par des exemples de filtrage de signaux rels.
4.3.2
55
Gabarit
Le gabarit du filtre est dfini par lutilisateur parmi les quatre catgories suivantes : filtre passebas, filtre passe-haut, filtre coupe-bande et filtre passe-bande. On peut rgler Fe la frquence
dchantillonnage, f la largeur de la bande de transition, 1 lamplitude des ondulations en
bande passante et 2 lamplitude des oscillations en bande attnue. Ces amplitudes sont exprimes en chelle linaire. Il est galement possible dutiliser comme
paramtres
londulation en
bande passante et en bande aaiblie exprimes en dB : dp = 20 log10
H (f )
1+
1+ 1
1 1
et da = 20 log10
1
2 .
GABARIT EN ECHELLE LINEAIRE
1 1
f
2
F1
F2
Fe
f
2
Choisir un gabarit en entrant ces dirents paramtres. Le but du TP va tre de comparer les
rsultats de la synthse des dirents filtres pour ce gabarit.
4.3.3
Filtres Rponse Impulsionnelle Finie (RIF)
Remarque pralable : lutilisateur du logiciel peut choisir lordre du filtre synthtiser. Notons
que pour un ordre 2, la synthse par la fentre triangulaire ne donne rien (problme de dfinition
de fentre) et que la synthse par Remez "plante", lalgorithme itratif associ ncessitant un
ordre minimal de 3.
Synthse par la mthode de la fentre

Evaluation de lordre : Possdant les caractristiques du gabarit respecter, il est possible
de calculer quel devra tre lordre N (trs approximatif et en gnral sous-estim) :
1
Fe
2
b
N = log10
3
10 1 2 f
o Fe est la frquence dchantillonnage, f la largeur de la bande de transition, 1 lamplitude

des ondulations en bande passante et 2 lamplitude des oscillations en bande attnue.
56
Calculer lordre donn par lapproximation pour le gabarit choisi.

b en utilisant la mthode de
Influence de la fentre : Synthtiser le filtre RIF cet ordre N
synthse RIF par fentre. Observer les fonctions de transfert obtenues par les quatre types de
fentres disponibles : valeur de la pente, position et amplitude du premier lobe doscillation.
Cette comparaison sera faite uniquement dun point de vue qualitatif, en utilisant le bouton
"hold" (et ne pas hsiter aussi utiliser le bouton "Ylog"). Quels sont les avantages et les
inconvnients de chacune des fentres ? Observer le temps de propagation de groupe.
Influence de lordre : On appelle ordre optimal Nopt dun filtre numrique, lordre minimal tel
que le gabarit frquentiel soit respect. Synthtiser le filtre RIF laide de la fentre rectangulaire
b (N
b < Nopt !) et en augmentant
en faisant varier lordre du filtre depuis la valeur calcule N
progressivement jusqu dpasser Nopt . Mesurer linfluence de lordre sur les paramtres suivants
:
la raideur de la pente
la position de la frquence de coupure 3dB.(on visualise |H(f )|2 (f ) donc la frquence de

coupure correspond un module carr gal 0.5)
la pseudo-priode des lobes doscillation
lallure du temps de propagation
Comparer lamplitude des ondulations en bande passante et en bande aaiblie. Quelle remarque
peut-on faire sur lamplitude des oscillations ?
Ralisation dun filtre particulier

Pourrait-on raliser un filtre passe-tout? Comment faudrait-il modifier le gabarit pour synthtiser un tel filtre ? Quel serait son intrt?
Optimisation du filtre RIF

La mthode des moindres carrs permet dagir diremment sur les ondulations en bande passante et en bande aaiblie. Comment? Quelle proprit spcifique possde lamplitude des
oscillations lorsquon applique la mthode de Remez ? Observer le temps de propagation de
groupe.
Rponse impulsionnelle et rponse indicielle

Observer la rponse impulsionnelle dun filtre RIF. Comment retrouve-t-on les coecients du
filtre partir de la rponse impulsionnelle ? Justifier son aspect symtrique. Retrouver le
temps de propagation. Quelles remarques peut-on faire sur la rponse indicielle (allure gnrale,
nombre doscillations en fonction de lordre, ...) ?
57
Filtrage de signaux
Vrifier la nature du filtre en lui imposant en entre un bruit blanc (utiliser le bouton "DSP"
qui calcule les DSP de lentre et de la sortie du filtre).
Filtrer un signal sinusoidal. Quobserve-t-on sur la sortie ?
Filtrer une somme de 2 sinus, compose dun sinus dans la bande passante du filtre et dun
autre sinus dans la bande attnue. Quobserve-t-on sur la sortie ?
4.3.4
Filtre Rponse Impulsionnelle Infinie (RII)
A la dirence des filtres RIF, on peut calculer de manire assez prcise lordre minimal ncessaire
pour un type de filtre analogique donn permettant de respecter le gabarit. Ces calculs sont
faits pour chaque type de filtre dans le logiciel par le bouton "Auto".
Synthse avec un Butterworth

Faire varier lordre du filtre synthtis en dmarrant en dessous de lordre donn par la formule
(bouton "auto") et en augmentant progressivement, jusqu dpasser lordre "auto". Observer
linfluence de lordre sur la pente et sur le temps de propagation de groupe. Comment voluent,
pour le temps de propagation, la valeur des maxima et les frquences correspondant ces
maxima? Pourquoi peut-il tre trs gnant davoir des temps de propagation trs dirents dans
la bande passante?
Comparaison entre les dirents filtres RII

Observer les fonctions de transfert obtenues par les quatre mthodes proposes (Butterworth,
Chebyshev I et de Chebychev II et elliptique) : valeur de la pente, position du premier lobe
doscillation, amplitude du lobe secondaire.
Quelle est la dirence entre les mthodes de Chebyshev I et de Chebychev II ?
Quel est lintrt de la synthse dun filtre elliptique ?
Comparer les temps de propagation de groupe des dirents filtres en ne tenant pas compte
des "glitchs" obtenus par la fonction matlab (grpdelay calcule le temps de propagation de groupe
partir des coecients des filtres mais prsente quelques bugs). Quelles consquences cela va-t-il
entraner ?
Rponse impulsionnelle
Quelle est la dirence par rapport la rponse impulsionnelle dun RIF ?
58
Filtrage dun signal

Vrifier la nature du filtre en lui imposant en entre un bruit blanc.
4.3.5
Conclusion
Donner les avantages et les inconvnients respectifs des filtres RIF et des filtres RII.
Chapitre 5
Initiation Matlab
Matlab, abrviation de MATrix LABoratory, est un langage trs performant pour le calcul
numrique. Les utilisations classiques sont : le calcul mathmatique, le dveloppement dalgorithmes,
la modlisation et la simulation, lanalyse de donnes et la visualisation, le dveloppement
dapplications telles que les interfaces utilisateur.
5.1
Commandes daide
Deux instructions seront trs utiles par la suite :

lookfor mot_cl qui permet de trouver une commande laide de mots cls. Exemple :
lookfor filter. Matlab donnera en rponse la liste des fonctions intgres dont laide associe
comporte ce mot cl.
help commande donne une aide sur la commande en question. Contrairement lookfor,
il est ncessaire de connatre le nom exact de la commande. Exemple : help mean.
help toolbox donne toutes les commandes de la toolbox par catgories. Exemple : help
signal; help stats; help optim; help nnet;
5.2
5.2.1
Matrices
Construction
Matlab travaille essentiellement sur des matrices, un vecteur ou un scalaire tant considr
comme une matrice avec une ligne et/ou une colonne. La saisie dune matrice peut seectuer
de direntes faons : soit en la rentrant comme une liste dlments, soit en la gnrant dynamiquement laide de fonctions, soit encore en la chargeant depuis un fichier. La construc59
60
Chapitre 5. Initiation Matlab
tion lment par lment peut se faire de la faon suivante :

A = [1 2 3; 4 5 6; 7 8 9]
Cette instruction cre une matrice 3 3 et laecte la variable A. Les lments dune ligne
sont spars par des espaces ou des virgules et les lignes sont spares par des point-virgules.
Matlab traite aussi les matrices complexes :
A
[1 2; 3 4] + i [5 6; 7 8]
[1 + 5 i
2 + 6 i; 3 + 7 i
4 + 8 i]
i ou j peuvent tre utiliss indiremment pour la notation imaginaire. Dans le cas o ils seraient
utiliss comme variables, on peut les recrer laide de :
ii = sqrt(1)
Les lments dune matrice sont rfrencs de la faon suivante :
A(2,3) est llment de la deuxime ligne et troisime colonne de la matrice A.
si lon dfinit un vecteur B, alors B(4) est le quatrime lment de ce vecteur, indiremment
quil soit ligne ou colonne.
les indices des vecteurs et matrices dbutent 1 et non 0. On peut aecter une valeur
scalaire un lment dune matrice comme suit :
A(2, 4) = 8
qui va aecter la valeur 8 llment de la deuxime ligne et quatrime colonne de A. Certaines
fonctions permettent de crer des matrices. Les fonctions matricielles les plus couramment
utilises sont :
zeros
matrice de zros
ones
matrice de uns
eye
matrice identit
linspace
vecteur en progression arithmtique
logspace
vecteur en progression logarithmique
diag
cration ou extraction de la diagonale dune matrice
Par exemple, zeros(m,n) cre une matrice mn remplie de zros, et zeros(n) cre une matrice
remplie de zros et de taille n n. Si x est un vecteur, diag(x) est la matrice diagonale avec
x sur la diagonale, les autres lments tant nuls. Si A est une matrice carre, diag(A) est un
vecteur compos de la diagonale de A. Les matrices peuvent tre construites par blocs. Soit A
une matrice 3 3 :
B = [A, zeros(3, 2); zeros(2, 3), eye(2)]
5.2. Matrices
61
va construire une matrice 5 5. Les gnrateurs de nombres alatoires (taper help stats
pour en avoir la liste) permettent de gnrer des matrices dont les lments suivent une loi
donne. La commande rand(n), par exemple, cre une matrice n n dont les lments sont
uniformment distribus entre 0 et 1 ; rand(m,n) cre une matrice de taille m n suivant la

mme distribution. La notation : permet dexploiter la structure matricielle de Matlab et
de raliser lconomie de boucles par rapport la plupart des langages de programmation. Par
exemple :
1:5 est le vecteur ligne [1 2 3 4 5].
Bien sr, les incrments peuvent ne pas tre entiers, la syntaxe est la suivante :
dbut :
pas :
fin (le pas tant par dfaut gal 1).
Exemple :
0.2:0.2:1.2 est le vecteur [0.2 0.4 0.6 0.8 1.0 1.2],
et 5:-1:1 est [5 4 3 2 1].
Lexemple suivant cre un signal sinusodal :
temps = [0:0.1:2];
y = sin(2*pi*f0*temps);
La notation : permet aussi daccder aux sous-matrices dune matrice. Par exemple, A(1:4,3)
est le vecteur colonne compos des 4 premiers lments de la 3me colonne de la matrice A.
A(:,3) est la troisime colonne de la matrice A et A(1:4,:) est compose des quatre premires
lignes de la matrice A. Cette notation permet aussi de remplacer des morceaux de matrices :
A(:,[2 4 5]) = B(:,1:3) remplace les colonnes 2, 4 et 5 de A par les trois premires colonnes
de B. La fonction reshape(A,m,n) reconstruit la matrice A avec m lignes et n colonnes.
5.2.2
Oprations lmentaires
Les oprations suivantes sont disponibles dans Matlab :

+
addition
conj
conjugu uniquement
soustraction
.0
transpos uniquement
multiplication
division gauche
puissance
division droite
conjugu transpos
Ces commandes sappliquent aussi aux scalaires. Si les tailles des matrices sont incompatibles,
alors un message derreur apparat, sauf dans le cas dune opration entre un scalaire et une
matrice o loprateur va sappliquer au scalaire et aux lments de la matrice. Attention : La
division matricielle est particulire. Soit A une matrice et b un vecteur respectivement colonne
ou ligne, alors :
x=A\b est la solution de A x = b, et
62
x=A/b est la solution de x A = b.
Les divisions droite et gauche sont lies par b/A = (A\b). Pour rsoudre le systme, Matlab
calcule linverse de A quand A est inversible et sa pseudo-inverse lorsque A est singulire sans
pour autant vous en avertir.
5.2.3
Oprations lment par lment
Les oprations prcdentes sont les oprations matricielles classiques. Les oprations *, ^, \,
/ peuvent aussi sappliquer lment par lment sur une matrice ou un vecteur en les faisant
prcder dun point. Par exemple : [1 2 3 4].*[1 2 3 4] ou [1 2 3 4].^2 donnent [1 4 9
16]
5.2.4
Fonctions lmentaires
Les fonctions lmentaires disponibles sous Matlab, sappliquent aussi bien des matrices qu
des scalaires. Les plus courantes sont :
sin, cos, tan, asin, acos, atan
sinh, cosh, tanh, asinh, acosh, atanh
exp, log, log10
rem, abs, angle, real, imag, conj, sqrt, sign
round, floor, ceil, fix
Dautres fonctions renvoient un scalaire lorsquelles sappliquent sur des vecteurs. Dans le
cas o on les appliquent des matrices, elles agissent gnralement sur chacune des colonnes
indpendamment, le rsultat est alors un vecteur. Ce sont :
max, min,
sort,
sum, prod, median, mean, std
any, all
Par exemple, llment maximal dune matrice A est donn par max(max(A)), et non pas par
max(A) qui renvoie un vecteur (consulter laide).
5.3. Gestion de lespace de travail
5.2.5
63
Fonctions matricielles
Les plus courantes sont :

eig
valeurs propres
poly
polynome caractristique
det
dterminant
inv
inverse
expm
matrice exponentielle
sqrtm
matrice racine carre
size
taille dune matrice
5.3
Gestion de lespace de travail
Les fonctions suivantes permettent davoir une vue globale de lespace de travail :
who
liste les variables
whos
liste les variables, leur taille et leur type
what
liste les fichiers dextension .m et .mat
Les variables peuvent tre supprimes laide de la commande clear nom_variable.
Attention la commande clear qui utilise seule eace toutes les variables de lespace de travail.
Lorsque lon quitte Matlab, toutes les variables sont perdues. On peut les sauver avec la
commande save nom_fichier nom_des_variables_a_sauver.
Le format du fichier est alors
en *.mat. Lorsque lon relance Matlab, on peut les recharger avec load nom_fichier.
Certaines commandes DOS ou Unix classiques sont valables dans Matlab :
pwd
indique le chemin du rpertoire courant
cd
change de rpertoire
dir
liste tous les lments dun rpertoire
what
liste les fichiers .m, .mat et .mex dun rpertoire
path
obtention ou initialisation du chemin
Dautres commandes DOS ou Unix peuvent tre excutes si elles sont prcdes de !. Exemple : !del nom_fichier.
5.4
Boucles, tests et relations
Dans Matlab, ce type de fonctions opre quasiment de la mme faon que pour les langages de
type Pascal ou C. Mais attention : Remplacer une boucle par un traitement matriciel, lorsque
cela est possible, conduit un gain de temps dexcution trs important.
64
5.4.1
Boucles
Lexcution rptitive dune suite dinstruction peut-tre ralise laide de for ou while.
La syntaxe de la boucle for est la suivante :
for condition
instructions
end
Par exemple, soit n donn :
x = [
];
%cre une matrice vide
for (i = 1 : n)
x = [x, i2]
end
La syntaxe de la boucle while est la suivante :
while expression
traitement
end
Le traitement va tre rpt tant que lexpression est vraie. Les groupes dinstructions qui vont
tre excutes ou non selon la condition, peuvent tre dlimites par les instruction case et
otherwise.
5.4.2
Test : instruction if
Cette instruction permet dexcuter une suite dinstructions conditionnelles. Sa syntaxe est la
suivante :
if condition
instructions
end
Le traitement va tre excut uniquement si la condition est vraie. Des branchements multiples
sont aussi possibles :
if condition1
instructions1
elseif condition2
instructions
else
instructions
end
5.4. Boucles, tests et relations
5.4.3
65
Relations
Les principaux oprateurs relationnels sont les suivants :

<
infrieur
>=
suprieur ou gal
>
suprieur
==
gal ( = est une aectation)
<=
infrieurou gal
dirent
Ils peuvent tre combins avec des oprateurs logiques suivants :

& et,
| ou,
non
Lorsque une relation est applique des scalaires, le rsultat est alors gal 1 ou 0 selon quelle
soit vraie ou fausse et lorsquelle est applique des matrices ou des vecteurs de mme taille, le
rsultat est une matrice dont les lments prendront les valeurs 1 ou 0.
Une relation entre matrices est interprte par while et if comme vraie si elle est vrifie par
tous les lments de la matrice. Par consquent, si lon veut excuter un traitement quand des
matrices A et B sont gales, on peut faire :
if A = = B
traitement
end
tandis que si on veut lexcuter lorsquelles ne sont pas gales, on peut faire :
if any(any(A = B))
traitement
end
ou tout simplement
if A = = B else
traitement
end
Il faut noter que :
if A = B
traitement
end
ne va pas donner le rsultat escompt, car le traitement va tre excut seulement si tous les
lments de A et B sont dirents.
Les fonctions any et all permettent de rduire les relations entre matrices des vecteurs ou
des scalaires. Consulter laide leur sujet.
Voir aussi la fonction find qui renvoie lindice des lments obssant une condition. Par
exemple, soit y un vecteur, alors find(y>2) va renvoyer le rang de tous les lments de y qui
sont strictement suprieurs 2.
66
5.4.4
Pause et break
Lors de lexcution dun programme, on peut, laide de linstruction pause, suspendre son
excution et la relancer si lutilisateur appuie sur une touche.
Linstruction break arrte lexcution dune boucle while ou for.
5.5
Fichiers .m
Matlab peut excuter un ensemble dinstructions contenues dans un programme. Ces programmes ont une extension .m. Il y a deux types de programmes .m : les scripts et les fonctions.
% permet dinsrer des commentaires.
Si le dernier caractre dun traitement est un point virgule, alors lachage sur la fentre de
commande Matlab est supprim.
5.5.1
Scripts
Un fichier script est compos dun ensemble dinstructions Matlab. Si le fichier a le nom prog.m,
alors la commande prog va lexcuter. Les variables dun programme script sont globales, et
son excution va changer leurs valeurs dans lespace de travail.
5.5.2
Fonctions
Les variables dune fonction sont locales, mais peuvent aussi tre dclares comme globales. Soit
la fonction :
function y
moyenne(x)
%calcul de la valeur moyenne.
sum(x)/length(x);
Cette fonction va calculer la moyenne des lments dun vecteur. Dans la ligne de commande,
on peut alors crire :
vect
[1 2 3 4 5 6];
moy
moyenne(vect)
qui va renvoyer la valeur moyenne des lments du vecteur vect.

La premire ligne du script dune fonction doit obligatoirement tre de la forme :
function [resultat1, resultat2, ...] = nom_fonction(argument1, argument2,...)
5.6. Graphiques
5.5.3
67
Textes, entres, messages derreurs
Les textes doivent tre encadrs par des quotes : s = texte.

Ils peuvent tre achs sur la fentre de commande Matlab : disp(0 affichage0 ).
Les messages derreurs peuvent tre achs laide de la commande error qui arrte aussitt
lexcution du programme : error(0 desole!0 ).
On peut aussi saisir des paramtres lors de lexcution : iter = input(0 nombre d iterations? 0 )
5.5.4
Vectorisation et pr-allocation
Dans le but dacclrer lexcution de programmes, il est important de vectoriser les algorithmes
dans les fichiers .m. En eet, l o dautres langages utilisent des boucles, Matlab peut utiliser
des oprations vectorielles ou matricielles. Considrons lexemple suivant :
for t = 1 : N
y(t) = sin(2 pi f0 t)
end
Si on vectorise ce programme, il sut dcrire:
t = 1 : N;
y = sin(2 pi f0 t);
Dans le cas o il nest pas possible de vectoriser, il est possible de rendre les boucles for plus
rapides en pr-allouant des vecteurs ou des matrices dans lesquels les rsultats seront stocks.
Par exemple :
y = zeros(1, 100);
for (i = 1 : 100)
y(i) = det(Xi);
end
En eet, si on ne pr-alloue pas un vecteur, Matlab va ractualiser la taille du vecteur y chaque
boucle ditration, ce qui prend normment de temps.
5.6
Graphiques
Matlab permet de raliser des graphiques 2-D et 3-D. Pour avoir un aperu des capacits de
Matlab en matire de graphismes, utiliser la commande demo. Dans les travaux pratiques, nous
utiliserons essentiellement des tracs 2-D. Soit deux vecteurs x et y de mme taille. La commande
68
plot(x,y) va tracer y en fonction de x. Soit tracer la fonction sinus de 0 4 :

temps = 0 : 0.01 : 4;
sinus = sin(2 pi f0 temps); plot(temps,sinus);
Si lon souhaite tracer un deuxime graphique sur une figure dirente, taper figure ou figure(i + 1)
si la figure courante porte le numro i.
Pour tracer un deuxime graphique sur la mme figure, taper hold on. Leet de cette commande
sannule avec hold off.
Lcriture de texte sur les axes dun graphe ou sur le graphe seectue laide des commandes :
title
titre du graphe
xlabel
titre de laxe des abcisses
ylabel
titre de laxe des ordonnes
gtext
placement dun texte avec la souris
text
positionnement dun texte spcifi par des coordonnes
La commande grid place une grille sur un graphe.

Les commandes prcdentes doivent tre saisies aprs la commande plot.
Par dfaut les courbes sont traces par des lignes continues, mais il y dautres choix : - -, :,
- ., ., +, *, o et x.
De mme, on peut prciser les couleurs avec : y, m, c, r, g, b, w et k
Par exemple, plot(x,y,r*) va tracer en rouge et avec des * la courbe de y en fonction de

x.
La commande subplot permet de partitionner un graphe en plusieurs graphes : subplot(m, n, p)
divise la fentre de la figure en une matrice m n de petits sous-graphes et slectionne le pe`me
sous-graphe pour la figure courante.
Les commandes semilogx et semilogy la place de plot permettent davoir respectivement

laxe des abcisses ou des ordonnes en chelle logarithmique.
5.7
Application 1 : probabilits et statistiques
Les commandes normrnd ou randn permettent de gnrer des valeurs suivant une loi gaussienne :
Exemple :
ou
bruit=sqrt(sigma2)*randn(1,N);
bruit=normrnd(0,sqrt(sigma2),1,N);
Remarque : Il est possible de gnrer un bruit non gaussien (voir help stats pour dautres
gnrateurs de nombres alatoires. Si on souhaite ajouter un bruit color il faudra dabord gnrer
un bruit blanc (avec randn si il est gaussien) et, ensuite, le filtrer en utilisant la fonction filter
Exemple : y = filter(B,A,z);
5.8. Application 2 : traitement du signal
69
Si on utilise un filtre de type RIF on prendra A=1.

Question 1 : En utilisant la commande randn, gnrer 1000 ralisations indpendantes
x1 , ..., x1000 , dune variable alatoire gaussienne centre de variance 4. Visualiser lhistogramme
et dterminer la moyenne et lcart type de la squence gnre laide des fonctions hist, mean
et std.
Calculer la moyenne et lcart type de la squence x21 , ..., x21000 . Comparer avec les rsultats
thoriques.
Question 2 : A laide de la commande rand, gnrer 1000 ralisations indpendantes de la
variable alatoire discrte X dfinie par :
P [X = 1] = 1/4
P [X = 0] = 1/2
P [X = 1] = 1/4
Calculer la moyenne de la squence gnre.
5.8
Application 2 : traitement du signal
Question 1 : Reprsenter dans le plan polaire (commande polar) les ples (commande roots)
du filtre suivant :
1
1
=
1
A(z)
1 1.45z + 0.81z 2
Question 2 : Gnrer le signal :
x(t) = 10 sin(0.4t + ) + b(t),
t = 0, ..., 127
o est une variable alatoire uniformment distribue sur [0, 2) et b(t) est un bruit blanc
gaussien de variance unit. Etudier la moyenne, la fonction dautocorrlation (fonction xcorr)
et le spectre de x(t) (commande fft). Comparer les rsultats obtenus aux rsultats thoriques.
Question 3 : En crant une fonction, gnrer un signal N.R.Z. soumis une modulation
damplitude (multiplication par cos(2f0 t)). Etudier le spectre du signal avant et aprs modulation. Comparer les rsultats obtenus aux rsultats thoriques.
5.9
Conclusion
Ceci ntait quun aperu des capacits de Matlab. Ce logiciel possde de nombreuses botes
outils lies au traitement du signal et aux probabilits lmentaires : signal (traitement du
signal), stats (statitiques), wavelets (ondelettes), nnet (rseaux de neurones), optim (optimisation), hosa (moments dordres suprieurs)... Pour avoir un aperu des possibilits de Matlab
et de ces direntes botes outils, taper demo.

TP TS 1en

Transféré par

Droits d'auteur :

Formats disponibles

TP TS 1en

Transféré par

Informations du document

Description originale:

Titre original

Copyright

Formats disponibles

Partager ce document

Partager ou intégrer le document

Options de partage

Avez-vous trouvé ce document utile ?

Ce contenu est-il inapproprié ?

Droits d'auteur :

Formats disponibles

TP TS 1en

Transféré par

Droits d'auteur :

Formats disponibles

TRAVAUX PRATIQUES DE TRAITEMENT DU SIGNAL

Table des matires

Etude de lhistogramme damplitude . . . . . . . . . . . . . . . . . . . . . . . . .

Variables alatoires discrtes - Estimation de la probabilit . . . . . . . .

Variables alatoires continues - Estimation de la densit de probabilit . .

Estimation de la fonction de rpartition . . . . . . . . . . . . . . . . . . .

Loi des grands nombres . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Thorme Central-Limite - Test de Kolmogorov . . . . . . . . . . . . . . . . . . .

Etude de distributions particulires . . . . . . . . . . . . . . . . . . . . . . . . . .

Demos de la bote outils stats de Matlab . . . . . . . . . . . . . . . .

Rappels de quelques distributions . . . . . . . . . . . . . . . . . . . . . . .

Table du Test de Kolmogorov . . . . . . . . . . . . . . . . . . . . . . . . .

Restitution aprs chantillonnage du signal dorigine . . . . . . . . . . . .

Annexe : Test de Kolmogorov . . . . . . . . . . . . . . . . . . . . . . . . .

TABLE DES MATIRES

Elements de programmation Matlab . . . . . . . . . . . . . . . . . . . . . . . . .

Quelques fonctions Matlab utiles dans le TP . . . . . . . . . . . . . . . .

3 TP Analyse Spectrale - Corrlations et Spectres

Proprits des fonctions de corrlation . . . . . . . . . . . . . . . . . . . .

Algorithmes de calcul des fonctions de corrlation . . . . . . . . . . . . .

Exemples dutilisation des fonctions de corrlation . . . . . . . . . . . . . . . . .

Dtection dun signal priodique noy dans un bruit . . . . . . . . . . . .

Identification dun filtre . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Corrlations thoriques de signaux particuliers . . . . . . . . . . . . . . .

Priodogramme dune sinusode bruite et estimation du SNR . . . . . . .

Analyse spectrale de signaux rels . . . . . . . . . . . . . . . . . . . . . .

Programmation Matlab (pour ceux qui le souhaitent) . . . . . . . . . . . . . . . .

Quelques fonctions Matlab utiles dans le TP . . . . . . . . . . . . . . . .

Filtre Rponse Impulsionnelle Finie (RIF) : rappels . . . . . . . . . . . . . . . .

Synthse par la mthode de la fentre . . . . . . . . . . . . . . . . . . . .

Optimisation de la synthse obtenue . . . . . . . . . . . . . . . . . . . . .

Filtre Rponse Impulsionnelle Infinie (RII) : rappels . . . . . . . . . . . . . . .

Filtres Rponse Impulsionnelle Finie (RIF) . . . . . . . . . . . . . . . .

Filtre Rponse Impulsionnelle Infinie (RII) . . . . . . . . . . . . . . . .

Oprations lment par lment . . . . . . . . . . . . . . . . . . . . . . . .

Gestion de lespace de travail . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Boucles, tests et relations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Textes, entres, messages derreurs . . . . . . . . . . . . . . . . . . . . . .

Application 1 : probabilits et statistiques . . . . . . . . . . . . . . . . . . . . . .

Application 2 : traitement du signal . . . . . . . . . . . . . . . . . . . . . . . . .

Etude de lhistogramme damplitude

Variables alatoires continues - Estimation de la densit

Chapitre 1. TP Probabilits et Statistiques

ralisation de la v.a. X soit comprise entre mq et (m + 1) q est :

Calculons la moyenne de cette variable :

1.1. Etude de lhistogramme damplitude

Nx tant le nombre dobservations de X(t) appartenant lintervalle considr de largeur x

et dont la variance vrifie :

Lestimateur de densit de probabilit possde donc un biais additif :

Chapitre 1. TP Probabilits et Statistiques

Estimation de la fonction de rpartition

La fonction de rpartition F se dfinit de la faon suivante :

Estimation de la densit de probabilit

1.2. Loi des grands nombres

Loi des grands nombres

Le traitement grands nombres permet de tracer en fonction de N lestimateur de la moyenne

Observer le comportement de cet estimateur en fonction de N pour des variables alatoires

Thorme Central-Limite - Test de Kolmogorov