TP Isj

Télécharger au format pdf ou txt
Télécharger au format pdf ou txt
Vous êtes sur la page 1sur 45

TRAITEMENT du SIGNAL

2022
__________________________________________________________________

TP sous MATLAB
%analyse harmonique d'une tension sinusoïdale redressée
N=2048; %nbre de points
F=50; %fréquence de la sinusoïde
T=1/F; %période
Fe=2000; % fr. ech.
Te=1/Fe; % per ech
t=0:Te:(N-1)*Te; % axe du temps, N points
x=abs(sin(2*pi*F*t)); % tension redressée
fprintf('moyenne de x(t) :\n m = %d\n',mean(x));
fprintf('puissance de x(t):\n P = %d\n',mean(x.^2));
fprintf('valeur efficace de x(t) :\n s = %d\n',sqrt(mean(x.^2)));
% spectre de Fourier de x
y=fft(x)/N; %TFD de x = |Cn|=An/2=sqrt(an^2+bn^2)
z=abs(y); % spectre de module de x
f=0:Fe/N:Fe/2-Fe/N; %axe des fréquences
%affichage
subplot(121);plot(t(1:50),x(1:50));grid on;
xlabel('temps en s');
ylabel('tension redressée');
A=[51.07 9.238 3.469 2.857];
fprintf('taux de distorsion harmonique TO = %d\n',sqrt(sum(A.^2)/(A(1)^2)-1));
subplot(122);plot(f(1:N/2),z(1:N/2));grid on;
xlabel('fréquence en Hz');
ylabel('spectre de module');
TP I INITIATION au logiciel MATLAB
_____________________________________________________________________________________________________

Ce TP a pour but d’apprendre à utiliser le logiciel Matlab afin de pouvoir développer des
applications simples en traitement du signal.
Créer d’abord un répertoire portant votre nom sous la directorie work pour y
développer vos programmes : >>mkdir « votre nom » puis >>cd « votre nom ».

I- Manipulation des variables


On distingue les variables scalaires et les variables vectorielles (matricielles en général).

1- Variables scalaires
Dans un premier temps on génère trois variables scalaires a, b et c de la manière
suivante :
>> a=2 ;
>> b=3;
>> c=4 ;

On peut consulter la valeur d’une variable en entrant son nom :

>> a
La réponse serait :

a=

2
Cela signifie que les valeurs des variables sont mémorisées automatiquement avec leurs
noms.

2- Taille des variables dans la mémoire

La taille d’un scalaire de type « double » est 8 Bytes = 8 Octets = 64 Bits

3- Commandes de base :
 who et whos :
Affiche la taille mémoire et types de toutes les variables utilisées.
 cd :
Affiche le répertoire (directorie) où vous opérez en ce moment.
 what, dir :
Affiche la liste les noms des fichiers contenus dans le répertoire actuel.
 help nom_fonction :
Donne un descriptif de la fonction et ses arguments d’entrée sortie.

4- Opérations sur les variables


>>d=a+b+c
>>e=a+b*c
>>f=(a+b)*c
>>g=(a/b)*c
>>h=a^2

5- Les matrices et les variables vectorielles

Matlab est optimisé pour l’usage matriciel : Eviter les formulations non matricielles.
 [ a b c ] est un vecteur ligne.
 [a ;b ;c ;] est un vecteur colonne.
 V’ est le transposé du vecteur V
 u=1:5 est le vecteur [1 2 3 4 5] (de même que [1 :5] et (1 :5).
 t=0 :2 :15 est le vecteur [0 2 4 6 8 10 12 14]
 sin(t) est le vecteur [sin(0) sin(2) …sin(14)]
 zeros(1,N) est le vecteur ligne nul à N éléments.
 Ones(1,N) est le vecteur ligne à N éléments égaux à 1.

6- Opérations sur les vecteurs


 Il faut respecter les dimensions des vecteurs et matrices.

Génération automatique d’un vecteur

V=début:pas:fin;
 On définit une valeur de début : début
 On définit une valeur de fin : fin
 On définit un pas de progression linéaire ou logarithmique
(incrémentation) : pas ;
 si le pas n’est pas spécifié, il est égal à 1 automatiquement.

>>v3=1:10
>>v4=1:-0.5:-1
>>debut=0;fin=256;pas=8;
>>v5=debut:pas:fin

>>debut=0;fin=2*pi;pas=0.1; où pi désigne le nombre 3.14…


>>v5=debut:pas:fin
Une liste des fonctions les plus courantes est disponible dans l’aide en ligne en tapant la
commande (elfun désigne elementary functions : sin, cos, log, exp…)
>>help elfun
7- Les matrices
Saisir la matrice 3x3 suivante :
a=[1 2 3 ;4 5 6 ;7 8 9]
Noter que deux lignes sont séparées par un point virgule.
Si vous entrez la commande a(:) vous obtenez le vecteur colonne [1 4 7 2 5 8 3 6 9]’.
La fonction eig donne les valeurs propres de la matrice a.
zeros(N) est la matrice nulle NxN, eye(N) est la matrice identité NxN.

II- ROGRAMMATION

1-LES SCRIPTS

Plutôt que de taper les commandes au clavier les unes après les autres pour effectuer
une tâche, ce qui vous oblige à refaire la même chose à chaque utilisation de cette tâche,
il est préférable de grouper les commandes dans un fichier à extension .m, ainsi tous les
programmes auront pour nom name.m. Il suffit alors de taper name pour que la tâche
s’exécute.

Exemple 1 :
Créez un fichier qui s’appelle essai1.m, qui génère un signal sinusoïdal x(t) de n points,
puis visualisez le à l’aide de la commande plot(x).
 Etape 1 : édition du fichier par la commande : edit essai1.m
 Etape 2 : taper les commandes suivantes dans la fenêtre d’édition :

%ce fichier génère et affiche une sinusoïde


t=0:0.1:2*pi; % t est le vecteur temps avec un pas d’échantillonnage 0.1

x=sin(t); % le signal x est un vecteur de même taille que le vecteur t


plot(t,x) ;grid on % dessine à l’écran x en fonction de t
Remarques :
1- Le texte débutant par % est un commentaire de votre choix, il est ignoré par le logiciel.
Vous pouvez supprimer ces commentaires, mais ils sont utiles lorsqu’on a plusieurs
programmes.
2- Chaque instruction doit être suivie d’un point virgule.
 Etape 3 : Sauvegardez le fichier dans le répertoire en cours ;
 Etape 4 : Exécutez le programme en entrant la commande suivante :
>>essai1

2- Utilitaires graphiques

 L’instruction title(‘titre de la courbe ou la figure’) ajoute un titre à la figure


visualisée.
 xlabel(‘titre des abscisses’) affiche un titre horizontal suivant x.
 ylabel(‘titre des ordonnés’) affiche un titre vertical suivant y.
 grid on et grid off quadrille ou non le graphique.

Exemple 2 :

Générez,à l’aide d’un programme prog1.m, deux périodes des deux signaux x(t)=cos(t)
et y(t)=sin(t) et visualisez les deux signaux sur une même figure, l’une en rouge et l’autre
en bleu à l’aide de la fonction plot(t,x,’b’,t,y,’r’). Mettre une légende.

Exemple 3 :
Visualisez sur une même figure,à l’aide d’un programme prog2.m, les quatre signaux
cos(t), sin(t),log10(t) et exp(t) en utilisant la fonction subplot(.,.,.) qui divise l’écran en
quatre sous figures : (2,2,1),(2,2,2),(2,2,3) et (2,2,4).
>>subplot(2,2,1),plot(t,x);

Saisie d’une donnée au clavier :

Pour saisir une variable x à partir du clavier, on utilise l’instruction :

x=input(‘x=’) ;

Affichage d’un message à l’écran :

Pour afficher à l’écran un message personnel suivi d’un retour à la ligne :

fprintf(‘message personnel \n’) ;

Affichage de la valeur d’une variable :


Pour afficher à l’écran la valeur d’une variable x :
fprintf(‘x=’%d) ;

III- LES FONCTIONS


Si plusieurs de vos programmes personnels utilisent en commun une liste d’instructions,
il est préférable de regrouper ces instructions sous forme d’un programme indépendant.
A chaque besoin on appelle le dit programme par son nom : c’est une fonction.
Une fonction possède des paramètres d’entrée et des paramètres de sortie, dont la
syntaxe de déclaration est la suivante :
function [sortie1,sortie2,…]=nom_fonction(entrée1,entrée2,…)
Le programme matlab correspondant à la fonction doit porter le même nom que la
fonction : nom_fonction.m
Exemple 1 :
Calcul de la moyenne arithmétique m d’un vecteur v de dimension n :
1 n
m=  v(i)
n i 1
La fonction moyenne.m comportera les instructions suivantes :
function resultat=moyenne(v,n); %déclaration de la fonction

r=0; %initialisation de la moyenne


for
k=1:n r=r+v(k);
end
resultat=r/n;
On note la présence de la boucle itérative :
for indice=debut:pas:fin
instruction ;
end
Qui est équivalente à :
indice=debut ;
(B) instruction;
debut=debut+pas;
si debut≤fin aller à B;
Mise en œuvre de la fonction
>>moyenne(1:100,100)
Comparer avec la fonction matlab mean.
Exemple 2 :
Créez une fonction puissance.m qui donne la valeur efficace s et la puissance moyenne
p d’un signal.
x étant un vecteur de composantes x(i) , i=1..n, la valeur efficace de x est :

1 n 2
s= p =  x (i)
n i 1
Testez cette fonction sur un signal sinusoïdal d’amplitude A et de période 1.
Comparer avec la fonction matlab std.

Exemple 3 : Diagramme de bode.


1- Créez une fonction spectre(f,H) qui affiche à l’écran le spectre d’amplitude
20log|H(f)| et le spectre de phase Arg(H(f)) en fonction de f en Hertz avec
indication du titre de chaque courbe ; Les paramètres d’entrée sont le vecteur
fréquence f=0:fmax et le vecteur complexe H.

Fonctions conseillées :abs(x),angle(x) et log10(x).

2- Utilisez la fonction spectre dans un programme indépendant myfilter.m pour


tracer les diagrammes de Bode des filtres passe bas et passe haut dont les
fonctions de transferts sont définies par :
f
j.
1 fc
H1(f)= , H2(f)=
f f
1  j. 1  j.
fc fc
fc étant la fréquence de coupure demandée par le programme et fournie par l’utilisateur à
l’exécution.
Exemple 4 : vitesse de convergence d’une série
Considérons la série en N suivante (qui converge vers exponentiel de x) :
N
xk
s ( x, N )  
k 0 k!
Réaliser :
 Une fonction factoriel(n) qui calcule le factoriel d’un entier naturel.
 Une fonction somme(x,n) qui calcule s(x,n).
 Un programme prog3.m qui calcule la valeur de l’exponentiel de 2 avec une
précision (incertitude relative) supérieure à 95% (incertitude inférieur à 5%).

PARTIE THEORIQUE
 
1. En utilisant e j .n  e j . , Trouver les expressions de cos(nα) et sin(nα) pour n=2 et n=3.
n

2. De même retrouver les expressions de cos(α+β) et sin(α+β).



1 1
3. Calculer la somme  s(n) , avec s(n)= 2
n 1
n
 j
3n
j
4. Même question pour s(n)= ( ) n
3
5. Trouver les complexes z tels que : z N  z  N
6. trouver les racines cubiques de 1.

7. Calculer le produit infini e 
n 1
j / 2n

 /2  /2
8. Calculer I n   sin ( x)dx et J n   cos ( x)dx
n n

0 0
________________________________________________________________________

Le compte rendu de la séance doit être rédigé sur la feuille qui vous est fournie et
doit contenir en plus de la partie théorique :
Un listing des programmes :
prog1.m
prog2.m
puissance.m
spectre(f,H)
myfilter.m
factoriel(n)
somme(x,n)
prog3.m

___________________________________________________________________________
TP II NOTIONS DE SIGNAL NUMERIQUE
_________________________________________________________________________________________________________________

I- Définition
Un signal numérique s(k) est une suite de N échantillons régulièrement espacés de T e
secondes :
s(0,),s(Te),s(2Te),…,s((N-1)Te)
La grandeur Fe=1/Te est appelée fréquence d’échantillonnage : c’est le nombre
d’échantillons par seconde.
Le nombre d’échantillons N est le plus souvent une puissance de 2 : N= 2 p où p est un
entier naturel.
Signal numérique

100
80
s(kTe)

60
40
20
0
Te 2Te 3Te 4Te kTe

kTe (secondes)

II- Visualisation du signal s


Pour visualiser s(kTe) correctement à l’écran, il faut préciser l’axe du temps t=k Te :
t=0:Te:(N-1)*Te;
plot(t,s);
Remarque : la commande plot(s), dessine s(k) en fonction de 1,2,…,N ; c’est-à-dire Te=1,
ce qui est rarement le cas.

Exemple 1 : Signal rectangulaire


 On considère le signal rectangulaire Rec(k) égal à l’unité pour
k=1…N/4 et nul pour k=N/4+1…N.
Pour N=512 et Te=0.01s visualisez à l’aide d’un programme
Rec(k) en précisant sur la figure l’axe du temps avec l’unité de
mesure.
 Visualisez deux périodes de l’onde carrée de période 2  fournie
par Matlab square(t) pour plusieurs pas d’échantillonnage :
t=0 :pas :4*pi
Exemple 2 : convolution de deux signaux
x(k) et y(k) étant deux signaux numériques de même durée N.Te, la convolution de x et y
est un signal z(k) de durée (2N-1).Te défini par :
k
z(k)=x*y(k)=  x(m) y (k  m  1) , k=1…2N-1
m 1

Remarque :
 Le nombre de points de z(k) est généralement N+M-1 où N et M désignent
le nombre d’échantillons respectifs de x et y.
 L’expression générale de la convolution est un peu différente si la variable
temps débute à zéro au lieu de un (ce qu’exige Matlab).

1- Réalisez un programme prog1.m qui effectue la convolution de x(k)=Rec(k) avec lui-


même et visualisez Rec(k) et z(k) résultant. On prendra N=512.
Observez la forme de z(k) ainsi que sa largeur par rapport à celle de Rec(k). A quel instant
k, z(k) est-il maximal ? Quel est la valeur de ce maximum ? Interpréter.
2- comparer z(k) et le résultat de z2=conv(x,x)
Exemple 3 : Signaux harmoniques
On considère les deux signaux :
4
1
s1(t)=  2k .sin( 2 .2kf t )
k 1
1

4
(1) k
s2(t)=  k 1 k
. cos(2 .k . f 2 .t )

Avec : N=1024, f1=50Hz, Fe=8Khz pour s1 et f2=1/2 Hz , Fe=100Hz pour s2,.


Générez et visualisez en même temps, à l’aide d’un programme prog2.m, les deux
signaux s1 et s2.
 Quelles remarques pouvez vous faire sur la forme des signaux?
 Que constatez vous lorsque les sommations vont plus loin que 4?
 A quel type de signaux connus pouvez vous comparer les signaux
générés (comparer à sawtooth (t) et square (t))?

III- NOTION de BRUIT

Le bruit b(k) est défini comme étant un signal indésirable se mêlant additivement ou
autrement à un signal s(t) qui intéresse l’observateur.
Le model le plus fréquent de bruit de mesure de grandeurs physiques est le bruit blanc
gaussien : c’est un bruit aléatoire b(k) dont les échantillons sont dé-corrélés (entendez
par cela que même si on connaît b(k1), b(k2) reste aussi imprévisible qu’avant).
La fonction b=randn(1,N) permet de générer un vecteur bruit b de distribution pseudo
normale (de Gauss) de taille N de moyenne nulle et d’écart type 1.
Un coefficient multiplié par randn permet d’augmenter à volonté la puissance du bruit.
Ainsi pour générer un vecteur bruit blanc de taille N, de moyenne m et d’écart type  la
commande est :
bruit= m+  *randn(1,N);
Dont la puissance est Pb  m 2   2 (Puissance de la moyenne plus celle des fluctuations
autour de cette moyenne).

IV- Rapport Signal sur Bruit (SNR )


Observons un signal y(k) bruité additivement : y(k) = s(k)+b(k)
Le SNR en dB définit le rapport de la puissance du signal Ps à celle du bruit Pb :

Ps
SNR= 10. log 10
Pb
On distingue trois cas de figure :
 SNR positif : le signal est plus puissant que le bruit.
 SNR nul : il y a autant de bruit que de signal.
 SNR négatif : le signal est dégradé, il ya plus de bruit que d’information.

L’énergie d’un signal x(k) est fournie sous matlab par sum(x.^2). Concernant la puissance
moyenne, il faut diviser l’énergie par le nombre d’éléments de x(k).

Exemple 1 : sinusoïde bruitée

Réalisez un seul programme prog3.m qui effectue les tâches suivantes :


 Génère une sinusoïde s(k) de fréquence 50Hz, d’amplitude 1,
échantillonnée à 2.5KHz (Te=0.4ms) et de taille n=256.
 Calcule l’énergie de s(k).
 Génère un bruit d’amplitude quelconque.
 Additionne le bruit à la sinusoïde
 Visualise les trois signaux simultanément.
 Affiche le SNR à l’écran.

Exécuter trois exemples donnant des SNR négatif, nul et positif.

PARTIE THEORIQUE

1- s(t) étant un signal défini sur [-T/2,T/2], on peut l’analyser sous la forme :

s (t )  C e
n  
n
j 2nf0t

Avec :

T
2
1
 s(t )e
 j 2nf0t
Cn= dt , f0=T 1
T T

2

(Cn)n constituent le spectre de Fourier du signal s(t) aux fréquences (harmoniques) nf 0 multiples
du fondamental f0.

a) Montrer que C0 est la valeur moyenne du signal s(t) sur [-T/2,T/2].

On l’appelle aussi la composante continue, s(t)-C0 étant la composante oscillante ou variable.

b) Exprimer s(t) en fonction de sin(2  nf0t) et cos(2  nf0t), n=0…   .

c) Quelle relation lie C-n et Cn?

2- Soit le signal périodique de période T=2  :


s(t)=0 si -  < t <0

s(t)=sin(t) si 0≤ t < 

a) Tracer s(t).
b) Calculer la valeur moyenne de s(t).
c) Calculer la série de Fourier de s(t).

3- Soit le signal carré de période T défini par :


s(t)=A si |t|<T/4

s(t)=0 si T/4<|t|<T/2

a) Tracer s(t).
b) Calculer la valeur moyenne de s(t).
c) Calculer la série de Fourier de s(t).

4- Soit le signal s(t) de période T avec 0<  <T/2 :


s(t)=A.sign(t) si T/4-  /2 < |t| < T/4+  /2

s(t)=0 si |t| < T/4-  /2 ou T/4+  /2 < |t|

a) Tracer s(t).
b) Calculer la valeur moyenne de s(t).
c) Calculer la série de Fourier de s(t).
d) On définit le spectre de puissance du signal s(t) par Pn=Cn 2 , Cn désigne un coefficient de
Fourier.
Calculer Pn.

Pour quel rapport  /T a-t-on P1=A 2 ?

_____________________________________________________________________________________________________
Le compte rendu de la séance doit être rédigé sur la feuille qui vous est fournie et
doit contenir en plus de la partie théorique :
Un listing des programmes : prog1.m, prog2.m, prog3.m
TP III TRANSFORMEE de FOURIER DISCRETE (TFD)
_____________________________________________________________________________________________________

I- Echantillonnage des signaux analogiques


Considérons un signal analogique s(t) avec t  R. Pour pouvoir manipuler ce signal par
un logiciel tel que matlab on doit l’échantillonner, c’est-à-dire prélever à une cadence Te
secondes des valeurs s(k) de s(t) :
s[-(N-1)Te]…s[Te],s(0),s[Te]…S[(N-1)Te]…
Après traitement du signal échantillonné et obtention des résultats projetés, on doit
retrouver à partir des seuls échantillons de s le signal analogique en totalité : c’est
possible si on respecte le théorème de Shannon.
Ce théorème stipule que si s(t) possède une transformée de Fourier de largeur de bande
2Fmax et si la fréquence d’échantillonnage Fe est telle que : 2Fmax<Fe alors on peut
retrouver s(t) à partir des échantillons s(k) de la manière suivante :

s(t)=  s(k ) sin c[ Fe (t  k .Te )]


sin(u )
Où : sin c(u )  .
u

Obstacles pratiques :
 Les signaux physiques observés possèdent le plus souvent des spectres
de Fourier occupant des bandes fréquentielles illimitées, donc Fmax est
infini.
 L’observation d’un signal physique se fait durant un intervalle limité [0,T]
ce qui donne le plus souvent un spectre très étalé.

Pour pallier ces difficultés, il faut noter que le spectre d’un signal physique tend vers zéro
avec la fréquence (l’énergie du signal s’évanouit aux hautes fréquences). Ceci nous
permet de limiter la bande significative du spectre et par la suite appliquer le théorème de
Shannon.
II- Transformée de Fourier Discrète :
La TFD d’ordre N d’un signal numérique s(kTe), k=0…N-1 est définie par :
N 1
n
S(f)=S( .Fe )=  s(kTe) Exp[ j 2kn / N ] , n= -N/2…N/2-1
N k 0

La transformation inverse :
N / 2 1
1
s(t)=s(kTe)=
N
 S (nFe / N ) Exp[ j 2kn / N ] , k=0…N-1
n N / 2

Remarques :
 Par abus d’écriture et pour simplifier on note s(k) et S(n) le signal et
sa TFD en omettant les facteurs Te et Fe.
 N est souvent une puissance de 2 et dans ce cas il ya un algorithme
rapide de calcul des N coefficients de la TFD, on parle de FFT (Fast
Fourier Transform).
Ainsi aux coefficients s(1)…s(N) correspondent par TFD les coefficients S(1)…S(N).
Sous matlab la TFD est donnée par la fonction fft(s,N).
La transformée inverse est donnée par ifft(s,N).
Transformée de Fourier Discrète

35
30
|S(nFe/N)|

25
20
15
10
5
0
Fe/N 2Fe/N 3Fe/N 4Fe/N (1-1/N)Fe

nFe/N (Hertz)

On définit les spectres d’amplitude et de phase du signal s(k) par :

20*log10(abs(fft(s,N)))
angle(fft(s,N))
Repliement de spectre (aliaising)

 Générez une sinusoïde de fréquence f 0 échantillonnée à Fe=100Hz avec


f0<fe/2.Observez le signal et sa TFD.
 Augmentez la fréquence f0 de la sinusoïde de manière à dépasser Fe/2.
En observant le spectre du signal donnez sa fréquence réelle fr.
 Que se passe-t-il lorsque f0=Fe/2 ?
 Exprimez fr en fonction de f0 et Fe.
 Expliquez à l’aide d’un schéma ce qui se passe en général lorsque f 0
dépasse Fe/2.

Signal Chirp linéaire

Générer le signal x(t)= sin( .t 2  2 .t  1) de durée une seconde avec Fe=1kHz. Visualiser
x en fonction de t. Visualiser la TFD de x. Commenter.
La fréquence instantanée de x(t) dépend de t, exprimer cette fréquence. Cela implique
que le théorème d’échantillonnage ne pourra jamais être respecté : aliaising inévitable.
Matlab permet de visualiser le spectre par tranches temporelles en fonction du temps :
C’est le spectrogramme. Observer ce spectrogramme à l’aide de la fonction
specgram(x) ;

prog1.m :
 Générer et visualisez les n/8 points de la TFD d’ordre n d’une sinusoïde
sin(t) ainsi que sa phase. Utiliser le zoom et « data cursor » pour lire les
pics.
 Faites de même pour une somme de sinusoïdes sin[(2k+1)t], k=0…4.

prog2.m
 Générer et visualisez Rec(k), sur 2 secondes. Visualisez les n/8 points des
spectres d’amplitude et de phase de la TFD d’ordre n=512 en respectant
les unités temporelle et fréquentielle.
 Générer et visualisez la convolution z(k) de Rec(k) avec lui-même (utilisez
la fonction conv(Rec,Rec)) ainsi que sa TFD et sa phase (Les n/8 points
uniquement).
 Comparez la TFD de z(k) et le carré de la TFD de Rec(k).
 Générez wave(k) une onde carrée à l’aide de la fonction square( pour
l’utiliser correctement tapez help square).
 Visualisez wave et sa TFD.
prog3.m
 Générer et visualisez la TFD d’une sinusoïde bruitée de fréquence f0 en lui
rajoutant un bruit blanc centré de puissance  2 et observez si l’on
distingue le pic correspondant à la fréquence de la sinusoïde. Testez
plusieurs valeurs de  2 et de f0.

PARTIE THEORIQUE

I- Propriétés de la TFD

On considère les signaux causals x(k) et y(k).

1. Déterminer la TFD de x(k-k0).


2. Déterminer la TFD de k.x(k).
3. Déterminer la TFD de x(k)*y(k).
4. Déterminer la TFD du produit xy.

II- Modulation
On considère un signal s(t) dont la transformée de Fourier S(f) occupe la bande de
fréquences [-Fmax,Fmax].

1- Donner l’expression de la transformée de Fourier analogique de sm(t)=s(t).cos(2 


fpt) en fonction de S(f).
2- Donner l’expression de la transformée de Fourier analogique du signal démodulé
mais non filtré : sdm(t)=sm(t).cos(2  fpt), en fonction de S(f).
3- Proposez une méthode pour extraire S(f) de sdm(t).
_____________________________________________________________________________________________________
Le compte rendu de la séance doit être rédigé sur la feuille qui vous est fournie et
doit contenir en plus de la partie théorique :
Un listing des programmes : prog1.m, prog2.m, prog3.m
TP IV CALCUL et ANALYSE SPECTRALE
_________________________________________________________________________________________________________________

I- Analyse du spectre par TFD


On va procéder à une étude de l’effet de la durée M d’observation d’un signal sur sa TFD
ainsi que l’effet de la fréquence d’échantillonnage Fe.
On prendra soin de manipuler des signaux de tailles N=2 p :
2,4,8,16,32,64,128,256,512,1024,2028,4056…
1- Influence de la taille de la fenêtre d’observation (windowing).
Il est possible d’appliquer la TFD sur la totalité ou une partie d’un signal de taille N. La
partie analysée est appelée fenêtre d’observation.
On va mettre en évidence le fait que la précision fréquentielle (résolution) augmente avec
la taille de la fenêtre.

Réalisez prog1.m qui :


1. Génère une sinusoïde s(k) de fréquence 100Hz, de taille N=512 échantillonnée à
1kHz.
2. Visualise les modules des TFD d’ordres M=16, 32, 64, 128, 256, 512 et donne à
chaque fois la largeur  f du lobe principal du spectre : Cette largeur est l’écart
entre les fréquences où le max du spectre est divisé par deux (largeur à mi-
hauteur) (utiliser la fonction find() sinon le zoom et « data cursor ».
3. Vérifiez que le produit p=M  f est constant. Conclure.

2- Fenêtre d’analyse :

Une fenêtre d’analyse est une mise en forme que l’on donne au signal avant de l’analyser
par TFD :
On multiplie s(k) terme à terme par une fenêtre w(k) de forme adéquate pour réduire les
lobes secondaires afin que le spectre observé tende vers le résultat théorique qui est pour
une sinusoïde de fréquence f 0 une raie pure à la fréquence f0.
Les analyseurs de spectre (appareils qui donnent le spectre d’un signal) proposent de
nombreuses fenêtres d’analyse, qui ont chacune des effets différents : Blackman,
Hamming, Hanning…
Sous matlab on peut utiliser les fonctions hanning(N) et blackman(N) pour le fenêtrage
(ce sont deux vecteurs colonnes de taille N).

Réalisez prog2.m qui :


4. Visualise sur une même figure les deux fenêtres citées et leurs spectres
d’amplitude sur [0, 0.5Hz]. Conclure sur l’utilité de ces deux fenêtres.
5. Multiplie terme à terme la sinusoïde s(k) du prog1.m par chacune des deux
fenêtres, et visualise les spectres résultants.
6. Observez l’effet du fenêtrage sur le spectre. Comparez l’effet des deux fenêtres
proposées.

II- Résolution fréquentielle de la TFD

La résolution fréquentielle  f nous informe sur la capacité à séparer les valeurs du


spectre pour deux fréquences très proches :
On peut distinguer les spectres aux fréquences f et f+  f mais pas aux fréquences f et f+
 où  <  f. Plus  f est petite plus la résolution est haute.
Cette résolution dans le cas de la TFD d’ordre N est déterminée par le nombre de points
d’analyse N et la fréquence d’échantillonnage Fe :

 f=Fe/N
On va essayer de voir les limites de l’analyse par TFD en essayant de détecter des
sinusoïdes prôches.

Réalisez prog3.m qui :

7. Génère un signal échantillonné à 1kHz de 1024 points, résultat de l’addition de


deux sinusoïdes de même amplitude A et de fréquences respectives 100Hz et
105Hz.
8. Déterminez expérimentalement la durée d’observation en dessous de laquelle on
ne distingue plus les lobes principaux des TFD des deux composantes.
9. Que vaut la taille théorique Nmin qui ne permet pas de séparer les deux raies ?
10. On double simultanément la fréquence d’échantillonnage et la taille de la TFD :
Est-ce qu’on obtient de meilleurs résultats ?
11. Un procédé physique génère un signal comportant deux raies 100Hz et 105Hz
présentes dans une bande [0,500Hz]. Quel serait le choix le plus économique de
la longueur d’observation et de la fréquence d’échantillonnage de ce signal afin de
pouvoir distinguer les deux raies ?

III- Calcul du spectre par auto-corrélation :

Réalisez prog4.m qui :

12. Visualise la TFD X(n) d’un signal rectangulaire x(k), sa fonction d’auto-corrélation
R(n) et la TFD de R. Utiliser xcorr(x,x,’biased’) ;
13. Visualise un bruit blanc b(k), son auto-corrélation B(n) et la TFD de B.
14. En comparant R(n) et X(n), vérifiez le théorème de Wiener-kinchine.
PARTIE THEORIQUE
On désigne par u(t) l’échelon unité, nul pour t négatif et égal à 1 ailleurs.

1- Déterminer l’énergie Es et la densité spectrale de puissance  s du signal s(t)=u(t).exp(-at)


avec a réel positif.
2- Déterminer la puissance moyenne Px et la dsp  x du signal déterministe x(t)=A.sin(2  f0t+ 
).
3- Déterminer la puissance moyenne Px et la dsp  y du signal aléatoire y(t)=A.sin(2  f0t+  ),
où  est une v.a. uniforme sur [0,2  ].
4- Quelle est l’énergie ΔE localisée dans la bande [f0-  f/2, f0+  f/2] d’un signal s(t) à énergie
finie.

5- Quelle est la puissance ΔP localisée dans la bande [f0-  f/2, f0+  f/2] d’un signal s(t) à
puissance moyenne finie. Que se passe-t-il si  f tend vers zéro ?

_________________________________________________________________
Le compte rendu de la séance doit être rédigé sur la feuille qui vous est fournie et
doit contenir en plus de la partie théorique :
Un listing du programme :
prog1.m, prog2.m, prog3.m, prog4.m
15. TP V NOTIONS de TRAITEMENT AUDIO
_____________________________________________________________________________________________________

I- Introduction au système acoustique humain


Le système de réception du son chez l’humain (l’oreille) comporte trois parties
principales :
 Un canal externe en contact direct avec le milieu extérieur (l’air).
 Une membrane (le tympan) sensible aux variations de pression dans le canal
externe.
 Une partie interne hétérogène qui détecte les hautes, moyennes et basses
fréquences.
Un son audible est caractérisé par :
 Intensité : la puissance acoustique du son en Watts par centimètre carré.
Elle dépend de la direction (position relative à la source du son)
 Hauteur : la fréquence de la fondamentale en Hertz. Ne dépend pas de la distance
à la source.
 Timbre : le contenu en harmoniques (les multiples du fondamental). Le même
timbre correspond à plusieurs phases relatives ou formes temporelles.

1- La puissance acoustique

Watt/cm2 Décibels Exemple de son


10-16 0 dB Son bas audible à 3kHz
10-14 20 dB Son bas audible à 10kHz
10-12 40 dB Son bas audible à 100Hz
10-10 60 dB conversation normale
10-7 90 dB Limite du bruit industriel
10-5 110 dB Concert de Rock

P
P étant la puissance par cm2, la puissance en dB s’écrit : PdB  10. log 10 , avec P0=10-
P0
16. La puissance est proportionnelle au carré de la pression subie par le tympan.

L’écart de puissance minimum distinguable par l’oreille est de l’ordre de 1dB soit à peu
près 11% d’écart relatif en puissance (correspond à une pression minimale de 20
micropascales).
L’oreille possède une sensibilité logarithmique : un facteur deux ressenti correspond à un
facteur dix en puissance.

2- La hauteur du son
Le fondamental étant f, les harmoniques sont 2f, 3f, 4f, …liés respectivement aux modes
fondamental (do3), second (do4), troisième (sol4)…
3. Le timbre
Il correspond au contenu spectral du son, peu importe son allure temporelle.
4. Echelle logarithmique
la séparation minimum entre niveaux distinguables étant de 1dB, pour coder la plage
allant de 0dB à 120 dB il faudra coder la plage de puissance allant de 0 à 10 12 soit sur 40
bits.
Du fait que l’oreille ait une sensibilité logarithmique, seuls 120 niveaux sont distinguables,
d’où un codage sur 8 bits uniquement. Cela exige une quantification logarithmique des
échantillons sonores (Echelle de quantification compensée).
X étant la donnée originale, la donnée quantifiée Y s’écrit :
ln(1   . X )
Y , 0  X  1 (standard US, ex.  =255)
ln(1   )
1  ln( A. X )
Y , 1/A  X  1 (standard EU, ex. A=87.6)
1  ln( A)
5. Qualité des données audio

Qualité Bande Taux (kHz) codage Débit Commentaire


Musique HIFI 5Hz…20kHz 44.1 16 bits 706 K Audiophile
Téléphone 200Hz…3.2kHz 8 12 bits 96 K La parole
Téléphone 200Hz…3.2kHz 8 8 bits 64 K La parole : courant
Parole codée LPC 200Hz…3.2kHz 8 12 bits 4K Qualité pauvre

II- Le son sous Matlab


Windows propose un standard des fichiers son : Le format ‘.wav’ (‘.au’ est le standard
chez unix). Ce format est rudimentaire et non compressé. Comme tout fichier formaté, un
fichier .wav comporte une entête : Champs de 44 octets, suivie des données (le son
proprement dit).
L’entête commence dès le premier octet (offset 0) et se compose des éléments suivants :
1 TAG1 4 Octets Constante « RIFF » (Ox52,Ox49,Ox46,Ox46)
2 SIZE1 4 Octets Taille du fichier moins 8 octets
3 FORMAT 4 Octets Format= « WAVE » (Ox57,Ox41,Ox56,Ox45)
4 TAG2 4 Octets Identifiant « fmt » (Ox66,Ox6D,Ox74,Ox20)
5 LGDEF 4 Octets Nombre d’octets utilisés pour définir le contenu
6 FORMAT 2 Octets Format de fichier (1 : PCM,….)
7 NBCANAUX 2 Octets Nombre de canaux : 1 pour mono et 2 pour stéréo
8 FREQ 4 Octets Fréquence d’échantillonnage en Hertz
9 BYTEPERSEC 4 Octets Nombre d’octets par seconde de musique
10 NBRBYTE 2 Octets Nombre d’octets par échantillon
11 NBBITS 2 Octets Nombre de bits par donnée
12 TAG3 4 Octets Constante « data » (Ox64,Ox61,Ox74,Ox61)
13 SIZE2 4 Octets Taille du fichier moins 44 octets

1. Lecture d’un fichier son


siz=wavread(‘son1.wav’,’size’) ;
siz(1) est le nombre d’échantillons et siz(2) le nombre de canaux.
Déduire la durée du son.
y=wavread(‘son1.wav’) ;
y est un vecteur colonne si le son est mono et une matrice à deux colonnes si le
son est stéréo qui contient les données. Les données sont dans la plage
[-1,1].
[x,fs,bits]=wavread(‘son1.wav’) ;
x est identique à y; fs est la fréquence d’échantillonnage en Hertz. bits est le
nombre de bits par sample.
wavplay(x,fs);
pour écouter, fs est la fréquence d’échantillonnage en Hertz.

u=wavread(‘son1.wav’,N) ;
u contient les N premières données ou échantillons.
fid=fopen(‘son1.wav’,’r’) ;
Ouvre le fichier son en tant que flot de données non spécifiés pour sa lecture.
fid est le pointeur du fichier :pointe la donnée actuelle.
son_lu=fread(fid,inf,’int16’);
son_lu est un vecteur colonne constitué de données 16 bits ;
inf demande la lecture jusqu’à la fin du fichier.
fclose(fid); Ferme le fichier.
Testez ces instructions sur un fichier son et en déduire la taille de l’entête.

2. Ecriture d’un fichier son

wavwrite(x,fs,’son2.wav’) ;

Crée un fichier son son2.wav dont les échantillons (compris entre -1 et 1) se


trouvent dans la matrice x avec une fréquence d’échantillonnage fs.

3. Enregistrement d’un fichier son


Brancher un micro à votre PC et saisir les commandes :
fs=44100 ; %fréquence d’échantillonnage
nbits=16 ; %bits par échantillon
mode=2 ; %mode stéréo
D=16 ; %Durée de l’enregistrement en secondes
x=audiorecorder(fs,nbits,mode) ;% x est une variable objet
recordblocking(x,D) ;
% Parler au micro pendant D secondes
play(x) ; % Ecouter le son enregistré
y=getaudiodata(x) ; % y contient les données son de type double
wavplay(y) ; % Ecouter le son enregistré
Ou:
x=audiorecorder (fs,nbits,mode);
record(x) ;
Parler au micro….
stop(x) ;
y=play(x) ;
Ecouter….
stop(y) ;
z=getaudiotata(x,’int16’);

prog1.m : sons purs


 Générer un son pur de fréquence audible x(t)= A. sin( 2f 0 t ) . (A=1).

 Ecouter x. Visualiser x et son spectre.


 Ajouter à x des sons de fréquences multiples de f0 et visualiser le résultat y et son
spectre. Ecouter y.

prog2.m
 Enregistrer un son x stéréo de durée 16 secondes.
 Visualisez sur deux systèmes d’axes différents les N=2048 premiers échantillons
de chaque canal (un extrait y).
 Visualisez les spectres de phase et d’amplitude correspondant à y.
 Modifier le contenu de x, en effaçant quelques passages, puis écouter le résultat.
prog3.m
 Mixage : Enregistrer deux sons différents x1 et x2 et fusionner les données
alternativement (un morceau de x1 suivi d’un morceau de x2…etc) dans un fichier
y.
 Ecouter x1, x2 et y .
___________________________________________________________________
Le compte rendu de la séance doit être rédigé sur la feuille qui vous est fournie et
doit contenir un listing des programmes : prog1.m, prog2.m, prog3.m
TP VI MODULATIONS NUMERIQUES
_____________________________________________________________________________________________________

I- Modulations analogiques

La modulation en général est un procédé technique permettant de transmettre des


signaux occupant une certaine bande spectrale de base [-fmax,fmax] d’un émetteur vers un
récepteur tout en occupant dans le milieu de transmission les séparant une bande
spectrale [fp-fmax,fp+fmax] très loin de la bande de base.
La fréquence fp est dite la fréquence porteuse. Evidemment on doit avoir fmax<fp/2 pour
que les deux bandes ne chevauchent pas, afin que le récepteur puisse restituer le signal
original émis.
Grâce aux propriétés de la TFD, il se trouve qu’en multipliant un signal s(t) par cos(2  fpt)
à l’émission, le résultat est un signal sm(t) ayant la même bande spectrale que s(t), mais
translatée de fp.
A la réception il suffit de multiplier par la même fonction cos(2  fpt) pour isoler la bande
de s(t).
Programme à faire : prog1.m
Le programme réalise les tâches suivantes :
 Génère un sinus, x(t), de fréquence 50 Hz échantillonné à 8kHz sur 512
points.
 Module x(k)+2 par un cosinus, p(t), de fréquence porteuse 2500Hz, ce qui
donne m(t).
 Génère dm(t)=m(t).p(t) et affiche les spectres des trois signaux.
 Comment restituer le signal x(k) à partir de dm(k).

Remarque : La multiplication échantillon par échantillon de deux vecteurs x et y sous


matlab : x.*y

On se propose de moduler un signal audio en utilisant les fonctions Matlab ammod et pmmod.

1- Modulation d’amplitude

prog2.m

Fp = 5000; % fréquence porteuse


[x,Fe] = wavread(‘son1.wav’) ;% signal audio à moduler

y = ammod(x,Fp,Fe); % signal modulé

z=amdemod(y,Fp,Fe) ;%démodulation du signal modulé

Visualiser et écouter les trois signaux x ,y et z ainsi que leurs spectres d’amplitude.

2- Modulation de phase

x = sin(2*pi*t) + sin(4*pi*t);

Fp = 10;

phasedev = pi/2; % deviation de phase pour la modulation de phase

y = pmmod(x,Fc,Fs,phasedev); % Modulation

Visualiser les deux signaux x et y.

Bruit du canal et démodulation

On additionne du bruit blanc gaussien à y (bruit du au canal de transmission), puis on effectue la


démodulation de phase.

y = awgn(y,10,'measured',103); % additive white Gaussian noise à y

z = pmdemod ( y, Fp, Fe, phasedev); % Demodulate.

Visualisez les signaux x, y et z.

II- Modulations numériques

1- Tracé des constellations

Matlab permet d’effectuer le tracé des constellations des modulations les plus courantes (trace
l’enveloppe complexe des symboles de l’alphabet). En voici deux exemples :
PSK-16

M = 16; % Le nombre de symboles de l’alphabet à transmettre

x = int16([0:M-1]);

scatterplot (pskmod(x,M)); % effectue le tracage

QAM-32

M = 32;

x = int16([0:M-1]);

y = qammod(x,M);

scale = modnorm(y,'peakpow',1);

y = scale*y; % Scale the constellation.

scatterplot(y); % Plot the scaled constellation.

2- Modulation QAM-16

M = 16; % taille de l’alphabet

x = randint(5000,1,M); % message aléatoire

y = qammod(x,M); % modulation

Visualisez x et y.

ynoisy = awgn(y,15,'measured'); % bruit additive du canal

scatterplot(ynoisy);

z = qamdemod(ynoisy,M); % démodulation
Visualisez x, y, ynoisy et z

PARTIE THEORIQUE

Enveloppe complexe d’un signal


On considère un signal réel x(t) de transformée de Fourier X(f). Soit y(t) le signal réel
défini par sa transformée de Fourier :
sign( f )
Y( f )  .X ( f )
j
y(t) est la transformée de Hilbert de x(t).
On appelle signal analytique de x(t) le signal complexe : ~ x (t )  x(t )  j. y(t )
1- quelle est la T.F. du signal analytique ~
x (t ) ? Faire un schéma des T.F.
2- Soit x(t )  A. cos(2 . f 0 t   ) .
Déterminer y(t) et ~
x (t ) , et représenter les trois signaux dans le plan complexe. En déduire
l’analogie avec la notion de phaseur en électricité.
La forme polaire du signal analytique : ~ x (t )  r (t ).e j . (t )
Permet de définir l’enveloppe réel r (t ) et la phase instantanées  (t ) d’un signal réel x(t).
3- Donner l’enveloppe et la phase du signal x(t )  A. cos(2 . f 0 t   ) .

4- Donner l’enveloppe et la phase du signal modulé x(t )  a(t ). cos(2 . f 0 t   ) , a(t) étant
un signal réel. Faire un tracé qui illustre l’enveloppe de x(t).
L’enveloppe complexe d’un signal réel x(t) est définie, pour une pulsation arbitraire 0
par : ~ x (t ).e  j .0t
r (t )  ~
On peut interpréter l’enveloppe complexe comme le résultat d’une modulation d’amplitude
du signal analytique.
5- Exprimer les parties réelles et imaginaires a(t) et b(t) de ~
r (t ) en fonction de x(t), y(t)
et les porteuses sin( 0 t) et cos( 0 t).
a(t) et b(t) s’appellent la composante en phase (IN) et en quadrature (Q) du signal réel
x(t).
6- Inversement, exprimer x(t) et y(t) en fonction de a(t), b(t) et les porteuses.
7- On considère un signal x(t) réel, aléatoire ou non, à spectre passe bande symétrique :
 x ( f )  1 , pour f1  f  f 2 Et nul ailleurs.
8- Exprimer les spectres de x(t), a(t) et b(t) en fonction du spectre de l’enveloppe
complexe. Faites un schéma illustratif.
_____________________________________________________________________________________________________
Le compte rendu de la séance doit être rédigé sur la feuille qui vous est fournie et
doit contenir en plus de la partie théorique : Un listing des programmes prog1.m et
prog2.m
TP VII ESTIMATION et DETECTION

I- Estimation de la moyenne

Soit X une variable aléatoire stationnaire d’ordre un et x(1)…x(N) une série de N observations ou
mesures de X décorrélées. On se propose d’étudier le comportement de l’estimateur de la
moyenne  =E[X] suivant (moyenne empirique):

1 N
̂ ( N )  . x(i)
N i 1

prog1.m

 Générez n=1000 points d’un signal aléatoire X en additionnant 4 périodes d’une sinusoïde
sin(t) et un bruit blanc de moyenne  =5.
 Calculez l’estimateur ̂ (N) pour les valeurs de k = 1..n. Utiliser mean().
 Visualisez ̂ (k), l’erreur  (k)=( ̂ -  )(k) et l’erreur quadratique  2(k).
Est ce qu’il y a convergence?

N.B. : Avant chaque génération de bruit il faut initialiser le processus de génération à l’aide de la
commande rand('state',0);

II- Estimation de l’écart type

On considère les deux estimateurs de la variance  2  E ( X   ) 2   suivants (variances


empiriques):

1 N
ˆ 1 2 ( N )  . [ x(i )  ˆ ( N )] 2
N i 1

1 N
ˆ 2 2 ( N )  . [ x(i)  ˆ ( N )] 2
N  1 i 1

prog2.m
 Générez 1000 échantillons décorrélés d’une variable aléatoire X stationnaire au second
ordre en additionnant une sinusoïde et un bruit blanc de moyenne  =5 et de variance  2
=1.5.
 Calculez les deux estimateurs ˆ (N) pour les valeurs de N = 1…1000.
 Visualisez ˆ 1(N), l’erreur  (N)=( ˆ 1-  )(N) et l’erreur quadratique  2(N).
Est ce qu’il y a convergence?

 Faites de même pour le deuxième estimateur.

III- Détection d’une composante continue par moyennage

On considère un signal constant C noyé dans un bruit b(k) blanc centré de variance v, le signal
résultant observé est s(k). On estime la valeur de la constante en prenant la moyenne de s(k).

prog3.m : Réaliser un programme qui estime C, en utilisant la moyenne empirique étudiée au


début du TP, ceci pour v=C/10, C/4, C/2, C.

_____________________________________________________________________________________________________
Le compte rendu de la séance doit être rédigé sur la feuille qui vous est fournie et
doit contenir en plus de la partie théorique :
Un listing du programme : prog1.m, prog2.m, prog3.m
TP VIII FILTRAGE NUMERIQUE

Il existe de nombreuses méthodes de synthèse de filtres numériques. Le passage de


l’analogique au numérique est tributaire de problèmes d’échantillonnage , surtout de
quantification des échantillons, faisant souvent perdre au filtre numérique conçu les
propriétés du modèle analogique.

I- Définition d’un filtre :


Un filtre est défini par son gabarit qui donne les limites de tolérance pour les différents
paramètres du spectre d’amplitude |H(f)| du filtre, à savoir :
 La fréquence de coupure fc
 Le module de ’atténuation de |H(f)| dans la bande coupée :  attdB
 Le module de l’ondulation de |H(f)| dans la bande passante :  ond
 La largeur de la bande de transition, située entre la fréquence de coupure
et la zone atténuée  trans

Remarques :
 il peut y avoir plusieurs fréquences de coupure.
 Les fonctions filter, fir1, fir2, remez…permettent la création de filtres
suivant plusieurs méthodes.
 La réponse fréquentielle (RF) du filtre, égale à la TFD de la RI en dB, doit
être contenue dans le gabarit pour que le filtre soit utilisable correctement ;
sinon on dit que la méthode diverge :il faut chercher une autre procédure.
 Le module de H(f) étant déterminé, le choix de la phase se fait
généralement en prenant tous les pôles et zéros de h(k) à l’intérieur du
cercle unité (exigences de stabilité et de causalité).

Passe haut
|H(f)
tra
|
ond

att
coupée

passante
|
f1 fc Fe/2 f

II- Analyse de filtres numériques

La fonction matlab y=filter(b,a,x) permet de filtrer le signal x(kTe), la sortie du filtre étant y(kTe).
b et a sont les vecteurs dont les éléments sont les coefficients du numérateur et du dénominateur
de la fonction de transfert du filtre :
M

b
m 0
m z m
H(z)= N
, b=[b0,..,bM], a=[a0,..,aN]
a
n 0
n z n

Par conséquent h=filter(b,a,d), où d est l’impulsion unité de Dirac.


prog1.m
Pour chacune des fonctions de transfert H(z) suivantes :
z 1 z 1  0.5 z 5
0.5.(1-z 1 ), 0.5.(1+ z 1 ), 0.5.(1- z 2 ), ,
1  0.5 z 1 1  0.5 z 1
Déterminez et tracez les caractéristiques suivantes du filtre associé (entre parenthèses
les fonctions matlab à utiliser):
 réponse impulsionnelle h (N=32) (filter(b,a,d))
 réponse fréquentielle H (TFD d’ordre N=32) (semilogy,unwrap)
 zéros et pôles (zplane(b,a))
 le retard de groupe en fonction de la fréquence (grpdelay(h))
 réponse à une sinusoïde pure de période 100.
 réponse indicielle (réponse à un échelon unité)
 Pour chaque filtre, indiquez de quel type de filtre s’agit-il (d’après la forme de la
réponse en fréquences) et donner les valeurs des fréquences de coupure ainsi
que les paramètres de tolérance.
Remarques :

1- La fréquence d’échantillonnage est égale à un par défaut.


2- Le retard de groupe  g d’un filtre est le retard moyen introduit par le filtre, c'est-
à-dire la durée moyenne qui sépare l’apparition du premier échantillon en sortie et
l’apparition du premier échantillon à l’entrée. Ainsi le filtre z-1 a un groupe delay
égal à une seconde : C’est un retard pure. Il est défini par :
d ( )
g   , avec  ( )  angleH ( )
d
 ( )
Le retard de phase est défini par :     .

Conséquence : Lorsqu’on veut comparer les données (l’entrée) et les données filtrées (la
sortie), il faut retarder l’entrée de  g secondes (compensation).

II- Synthèse de filtres numériques :


prog2.m : Créez par TFD un passe-bas de caractéristiques :

Fe Fc  trans  ond  att


8kHz 1kHz 200Hz 1dB -40dB
 Dessinez le gabarit théorique sur [0, Fe].
 Visualisez la réponse impulsionnelle h de taille n=64 déduite de
l’échantillonnage du gabarit (prendre pour h la partie réelle de la TFD
inverse des échantillons) et vérifiez si la RF finale est bien dans le gabarit.
 Filtrez deux sinusoïdes de fréquences respectives F2 supérieure et F1
inférieure à Fc. Quel effet a le filtre ?
 Mesurez sur l’allure temporelle des signaux filtrés le décalage induit par le
filtre, le retard de groupe, (fonction grpdelay (h)).
 Comparez ce décalage à la longueur de la RI du filtre.
Est-il identique pour les deux fréquences F1 et F2?
 Quelle partie utile du signal filtré de même taille que le signal original peut
on garder ?

III- Filtrage d’harmoniques :


L’opération de filtrage est réalisée par convolution du signal à filtrer x(k) avec la RI du filtre
h(k).
En numérique la convolution y(k)=x(k)*h(k) introduit un décalage de l’ordre du filtre qui
correspond au temps de réponse du filtre en analogique.
En considérant les spectres en dB, il est facile de vérifier que celui de la sortie du filtre est
égal à la somme de la réponse fréquentielle du filtre et du spectre du signal d’entrée.
La distorsion harmonique est une forme de distorsion courante en musique. Elle consiste
à ajouter des harmoniques (fréquences multiples de la fréquence fondamentale) au signal
de départ.
prog3.m
 Ajoutez à une sinusoïde pure (fondamental f 0=100Hz) ses 5 premières
harmoniques impaires de même amplitude avec Fe=8kHz et N=1024.
Visualisez le signal et son spectre.
 Réalisez le filtre adéquat pour extraire les harmoniques 5.f0 et 11.f0 en plus
du fondamental f0.
 Effectuez le filtrage par convolution.
 Visualisez le signal filtré ainsi que son spectre.
 Notez l’amplitude du signal filtré et interprétez sa valeur.

IV- Réduction du bruit :

On se propose de rehausser une sinusoïde noyée dans le bruit. Un bruit blanc possède
un spectre étalé sur toute la bande fréquentielle, alors que celui du sinus est concentré
autour de sa fréquence.
prog4.m

 Générez un sinus de fréquence 1kHz échantillonné à 10kHz


 Ajoutez un bruit centré au sinus avec un RSB=10dB.
 Quel gabarit pour atténuer le bruit ?

Réalisez le filtrage par moyenne mobile sur 5 points et comparez les spectres et
allures temporelles des signaux bruité et débruité.
PARTIE THEORIQUE

I- Filtre à phase linéaire

On considère le filtre FIR défini par :

H(z)=1+2 z 1 +3 z 2

1) Montrer que son auto-corrélation G(z)=H(z)H( z 1 ) est à phase linéaire.


2) Donner le group delay.
3) Déterminer la réponse impulsionnelle g(n) de G(z).
_____________________________________________________________________________________________________
Le compte rendu de la séance doit être rédigé sur la feuille qui vous est fournie et
doit contenir en plus de la partie théorique :
Un listing des programmes : prog1.m, prog2.m, prog3.m, prog4.m
___________________________________________________________________________
TP IX NOTIONS de TRAITEMENT d’IMAGES
_____________________________________________________________________________________________________

I- DEFINITIONS

1- Une image à deux dimensions est définie comme une matrice I(N,N) de N ligne et N colonnes.

2- Le point (i,j) est appelé pixel (picture cell), sa valeur I(i,j) réelle traduit le degré de luminosité
de l’image au point (i,j) : On l’appelle niveau de gris de l’image au point (i,j).

3- Généralement N est une puissance de 2 (32,64,…,1024..) et le niveau de gris des pixels varie
entre 0 et 255=256-1=2 8 -1.

0 correspond au noir (luminosité minimum) et 255 au blanc (luminosité maximum) ; Tous les
niveaux sont codables sur 8 bits (octet ou byte). Ceci pour une image noir et blanc, pour une
image couleur, il faut prévoir trois matrices qui correspondent aux niveaux du rouge, vert et bleu
(RVB ou RGB).

pixel
(1, N )
(1,1)

( N ,1) (N , N )

II- Mémoire occupée par une image :

Les images noir et blanc ou couleur occupent beaucoup d’espace mémoire, ceci est ressenti lors
de la conception d’un programme manipulant plusieurs images.

Ainsi une image telle que celle définie ci-dessus, va occuper NxNx8 bits ou encore NxN octets (un
Kilo octets=1Kb=1024 octets=2 10 octets):
N Mémoire en bytes Mémoire en Kby

128 16.384 16

256 65.536 64

512 262.144 256

1024 1.048.576 1024=1Mby

III- Codage d’une image :

On parle aussi de compression d’une image. Vue la taille occupée par une image, surtout en
couleur et lorsque les images sont animées, on a inventé des méthodes qui permettent de réduire
considérablement la taille mémoire d’une image après compression : On distingue les images
bitmap, tif, gif, epg, jpeg…suivant la méthode de compression utilisée.

MATLAB Dénomination

‘jpg' or 'jpeg' Joint Photographic Experts Group (JPEG)

'tif' or 'tiff' Tagged Image File Format (TIFF)

'gif' Graphics Interchange Format (GIF)

'bmp' Windows Bitmap (BMP)

'png' Portable Network Graphics (PNG)

'hdf' Hierarchical Data Format (HDF)

'pcx' Windows Paintbrush (PCX)

'xwd' X Window Dump (XWD)

'cur' Windows Cursor resources (CUR)

'ico' Windows Icon resources (ICO)

IV- Format d’une image :


Comme tout fichier formaté, une image en mémoire est un fichier qui possède deux champs :

 Une entête : c’est une suite d’octets spécifiant le nombre de lignes et de colonnes
de l’image, le nombre d’octets par pixel, le genre de codage utilisé pour compresser
l’image…
C’est justement le format induit par la méthode de compression : tif, gif,

mpeg…

 Les données : c’est la suite des niveaux de gris de tout les pixels de l’image tenant
compte de la méthode de compression, ce qui exige pour lire l’image un programme
spécifique connaissant le code de l’image.

V- Lecture et écriture d’une image :

La fonction imread (‘nom de l’image’) permet de lire une image :

Si R est une matrice NxN, la commande imwrite (R,’nom.tif’) crée une nouvelle image

NxN nom.tif dont les pixels ont les niveaux de gris figurant dans la matrice r.

Exemple 1

i=imread (‘cameraman.tif’) ; %lecture de l’image cameraman de format tif

Cette instruction stock l’image cameraman.tif dans une matrice i de type uint8.

Visualisation d’une image

imshow (i) %affichage de l’image

On peut visualiser une matrice R de type uint8 avec la même commande.

Si R est de type double, il faut la convertir uint8(R) pour la visualiser.


Exécutez la commande whos pour voir la taille de l’image, le nombre de lignes et de colonnes.

Pour connaître les niveaux de gris de l’image, il suffit de taper la commande i.

Les niveaux de gris de l’image sont disponibles maintenant dans la matrice i, ce qui nous permet
d’effectuer toutes les transformations possibles sur cette image et de sauver les résultats du
traitement dans une autre image (il faut la convertir au type double avant de la manipuler
x=double(i)).

VI- Opérations sur les images

Les opérations mathématiques nécessaires s’effectuent sur des variables de type double ; il serait
indispensable de convertir les images sources avant leur traitement de la manière suivante :

isource = fread (‘nom_image’) ;

id=double (isource) ;

C’est sur id que seront appliqués les traitements.

1- calcul de la moyenne et de l’écart type

prog1.m :

Calculer la moyenne arithmétique et l’écart type des niveaux de gris de l’image


cameraman.tif ainsi que les caractéristiques suivantes :

m=mean (i( :) ); % donne la moyenne de tous les niveaux de gris :un scalaire

mc=mean (i) ; % donne la moyenne de chaque colonne de l’image : un vecteur

ml=mean (i) ; % donne la moyenne de chaque ligne de l’image : un vecteur

d=(i-m).^2 ; % écart quadratique : une matrice

v=mean(d(:)) ; % variance des niveaux de gris : un scalaire


s=sqrt(v) ; % écart type des niveaux de gris : un scalaire

inf=min(x( :)) ; % Le niveau de gris le plus faible: un scalaire

sup=max(x (:)) ; % Le niveau de gris le plus élevé: un scalaire

2- Filtrage d’une image

On considère le pixel x(i,j) d’une image donnée à laquelle on a additionné un bruit gaussien
(utiliser rand(N,N)), et on s’intéresse à ses 8 pixels proches voisins : x(i-1,j-1), x(i,,j-1), x(i+1,j-1),
x(i-1,j), x(i+1,,j), x(i-1,j+1), x(i,j+1) et x(i+1,j+1). Ces neuf pixels forment une fenêtre locale :

 x(i  1, j  1) x(i  1, j ) x(i  1, j  1) 


 x(i, j  1) x(i, j ) x(i, j  1) 

 x(i  1, j  1) x(i  1, j ) x(i  1, j  1)

On considère d’autre part « un masque » de taille 3x3 à coefficients ai,j donnés :

 a1,1 a1, 2 a1,3 


 
a 2,1 a 2, 2 a 2,3 
 a3,1 a 3, 2 a3,3 

Un filtrage linéaire de l’image x(.,.) consiste à créer une nouvelle image y(.,.) en glissant le masque
point par point et en effectuant une combinaison pondérée linéaire des niveaux de gris :
m 3 n 3
y (i, j )   a m,n .x(i  m  2, j  n  2)
m 1 n 1

On voit bien que y est une convolution bidimensionnelle de x et a.

Filtrage par la moyenne

Les coefficients du masque sont tous égaux à 1/9 :


Réalisez un programme prog2.m qui effectue le filtrage local 3x3 et visualisez le résultat pour
plusieurs images. Quel est l’effet du filtre ?

Filtrage médian

Le filtrage médian est une opération non linéaire, qui consiste en deux étapes :

 Ranger les niveaux de gris de la fenêtre locale de x dans l’ordre croissant.


 Affecter à y(i,j) le niveau de gris médian (qui se classe au cinquième rang parmi les neuf
niveaux).

prog3.m

Réalisez le filtrage médian et expliquez son effet en le comparant au moyennage.

VII- Seuillage d’une image (Binarisation)

Le seuillage d’une image consiste à choisir un nombre réel s (threshold) et de tester les pixels
de l’image point par point et ligne par ligne : tout niveau de gris inférieur au seuil s est mis à zéro
(noir) sinon il est mis à 255 (blanc) ou 1 : Binarisation.

Matlab permet en une seule instruction de comparer les éléments d’une matrice M à un seuil
donné s : a=(M>s), a est une matrice composée de 1 et de 0.

prog4.m

 Réalisez un programme qui effectue le seuillage d’une image.


 Choisissez plusieurs seuils, y compris la moyenne de l’image.
 Effectuez le seuillage sur plusieurs images.
 Afficher l’image originale et l’image seuillée en même temps

VIII- Détection du contour d’une image

Le contour d’une image Im(N,N) est lui-même une image C de même taille NxN :
L’image C contient les limites ou frontières ou périmètres des objets contenus dans Im. Par
exemple si Im est la photo d’un disque plein alors C sera une image contenant un cercle de même
centre et rayon que le disque.

La détection du contour est l’opération mathématique et informatique qui permet de déduire C à


partir de Im.

Méthode du gradient :

Le calcul des pixels de C s’effectue en faisant la soustraction des pixels voisins de Im :

 Horizontalement (gradient horizontal): C(i,j)=Im(i, j+1)-Im(i,j)


 Verticalement (gradient vertical): C(i,j)=Im(i+1, j)-Im(i,j)
Il faut prendre la valeur absolue des différences !!

Le parcours d’une image ligne par ligne (ordre lexicographique) dans un programme s’effectue
de la manière suivante :

for i=1 :N %parcours des lignes i

for j=1 :N %parcours des colonnes j

traitement1 ;

traitement2 ;…

end

end

prog5.m

 Réalisez un programme qui donne le contour horizontal d’une image et visualisez


le résultat.
 Réalisez un programme qui donne le contour vertical d’une image et visualisez le
résultat.
 Générez et visualisez l’image moyenne des deux contours précédents.
 Testez le programme sur plusieurs images.

Original Saturn Image Edge Map

Original Image Corrupted Image Filtered Image

____________________________________________________________________________________________________
Le compte rendu de la séance doit être rédigé sur la feuille qui vous est fournie et
doit contenir :
Un listing des programmes : prog1.m, prog2.m, prog3.m, prog4.m, prog5.m
___________________________________________________________________________

Vous aimerez peut-être aussi