Signaux Matlab2

TRAITEMENT du SIGNAL
2018/2019
__________________________________________________________________
SOLUTIONS des TP sous MATLAB
Mohamed SABRI- département de physique- FST Béni Mellal- Université sultan Moulay Slimane
%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 :
n
1
∑ v (i )
m= 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 :
n
s= √p = √ 1
∑
n i=1
x 2 (i)
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
f f
1+ j. 1+ j.
H1(f)= , H2(f)= 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
n
1. En utilisant e j .nα = [ e j. α ] , Trouver les expressions de cos(nα) et sin(nα) pour n=2 et
n=3.
2. De même retrouver les expressions de cos(α+β) et sin(α+β).
∞
1 1
∑ s(n) n
+j n
3. Calculer la somme n=1 , avec s(n)= 2 3
j
( )n
4. Même question pour s(n)= 3
N −N
5. Trouver les complexes z tels que : z̄ =z
6. trouver les racines cubiques de 1.
∞ n
∏ e jπ /2
7. Calculer le produit infini n=1
π /2 π /2
n
I n= ∫ sin ( x )dx J n = ∫ cosn ( x )dx
8. Calculer 0 et 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
TPI : INITIATION à MATLAB
___________________________________________________________________________
factoriel.m
%TP1 : fonction factoriel.m
%calcul du factoriel d’un nombre naturel
function p=factoriel(n)
if n>1
p=factoriel(n-1)*n ;% récursivité
else
p=1 ;
end
end
___________________________________________________________________________
myfilter.m
%TP1 : myfilter.m
% tracage des réponses en fréquence de deux filtres
%utilisant la fonction spectre(f,H)
fc=input(‘fréquence de coupure à -3dB (entre 0 et 100 Hz) ?’) ;
f=0 :0.01 :100 ;
H1=(1+j*f./fc).^(-1) ;
H2=(j*f./fc).*H1 ;
spectre(f,H1,1) ;% l’indice 1 pour tracer les deux spectres de H1 en haut
spectre(f,H2,3) ;% l’indice 3 pour tracer les deux spectres de H2 en bas
___________________________________________________________________________
prog1.m
%TP1 : prog1.m, 2008
%tracage des fonctions sin et cos sur la même figure
Te=0.001;
T=2*pi;
t=0:Te:2*T;
x=sin(t);
y=cos(t) ;
plot(t,x,’b’,t,y,’r’) ;grid on ;
legend(‚Sin(t)’,’Cos(t)’);
___________________________________________________________________________
prog2.m
%TP1 : prog2.m, 2008
%tracage des fonctions sin, cos, log10, exp
Te=0.001;T=2*pi;t=0:Te:2*T;
x=sin(t);y=cos(t) ;t1=Te :Te :5 ;z=log10(t1) ;
t2=-2 :Te :2 ;u=exp(t2) ;
subplot(2,2,1) ;plot(t,x) ;grid on ;legend(‘Sin(t)’) ;
subplot(2,2,2) ;plot(t,y) ;grid on ;legend(‘Cos(t)’) ;
subplot(2,2,3);plot(t1,z);grid on;legend(‚log10(t)’);
subplot(2,2,4);plot(t2,u);grid on;legend(‚Exp(t)’);
_________________________________________________________________________
TP1 : prog3.m
% calcule exp(2) avec une précision > 95%
s=somme(2,1) ;
i=1 ;%nombre d’itérations
while ((exp(2)-s)/exp(2)>0.05) i=i+1 ;% itération suivante s=somme(2,i) ; end
fprintf(‘valeur de exp(2) : %d\n’,s) ;
fprintf(‘nombre d’’itérations effectuées : %d\n’,i) ;
fprintf(‚incertitude : %d\n’,(exp(2)-s)/exp(2));
___________________________________________________________________________
puissance.m
% TP1 : fonction puissance.m
%calcul de la puissance moyenne et de la valeur efficace d’un signal

function r=puissance(v) ;
p=mean(v.^2) ;%la puissance de v
s=sqrt(p) ;%la valeur efficace de v, ou son écart type
r=[p,s];
end
___________________________________________________________________________
somme.m
%TP1 : fonction somme.m
% calcule le développement de Taylor de Exp(x)
function s=somme(x,n) s=1;for k=1:n s=s+x^k/factoriel(k);end
end
___________________________________________________________________________
spectre.m
% TP1: fonction spectre.m
% tracage des spectres de phase et d'amplitude d'une réponse fréquentielle
function spectre(f,H,i)
subplot(2,2,i);plot(f,20*log10(abs(H)));title('spectre d''amplitude');
xlabel('fréquence (Hz)');ylabel('20Log10(|H|)');grid on;
subplot(2,2,i+1);plot(f,angle(H));title('spectre de phase en radian');
xlabel('fréquence (Hz)');ylabel('arg(H)');grid on;

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.
p
Le nombre d’échantillons N est le plus souvent une puissance de 2 : N= 2 où p est
un entier naturel.
Signal numérique
100
80
60
s(kTe)
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 T e :
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.T e, la convolution de x et y
est un signal z(k) de durée (2N-1).Te défini par :
k
∑ x (m) y (k−m+1 )
z(k)=x*y(k)= m=1 , k=1…2N-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
∑ 2k1 . sin(2 π .2 kf 1 t )
s1(t)= k =1
4
(−1)k
∑ k . cos(2 π . k . f 2 . t )
s2(t)= k =1
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);
2 2
Dont la puissance est Pb =m +σ (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 P s à celle du bruit Pb :
Ps
10 . log10
SNR= 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 :
+∞
j 2π nf 0 t
s(t )= ∑ C n e
n=−∞
Avec :
T
2
1 − j 2π nf 0 t
∫ s(t )e dt
T T −1
−
Cn= 2 , f0=T
(Cn)n constituent le spectre de Fourier du signal s(t) aux fréquences (harmoniques) nf0 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).
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).
2
d) On définit le spectre de puissance du signal s(t) par P n=Cn , Cn désigne un coefficient
de Fourier.
Calculer Pn.
2
Pour quel rapport θ /T a-t-on P1=A ?
_____________________________________________________________________________________________________
Un listing des programmes : prog1.m, prog2.m, prog3.m
TPII : SIGNAUX NUMERIQUES
___________________________________________________________________________
prog1.m
%TP2 : prog1.m
% calcule la convolution simple d'un signal rectangulaire avec lui même
n=512; % nombre d'échantillons
for i=1:n/4 rec(i)=1;end % on génére une impulsion rectangulaire Rec de durée n/4
for i=n/4+1:n rec(i)=0;end
x=[rec,zeros(1,n-1)]; % on complète avec des zéros car la convolution x*x admet 2n-1 échantillons
for k=1:2*n-1 % calcule l'échantillon z(k)=x*x(k)
z(k)=0; % initialisation de z(k)
for m=1:k
z(k)=z(k)+x(m)*x(k-m+1); %calcule la convolution proprement dite
end
end
subplot(2,1,1);plot(rec);grid on;title('x : impulsion rectangulaire');
subplot(2,1,2);plot(z);grid on;title('convolution x*x');

___________________________________________________________________________
prog2.m
%TP 2 : prog2.m
%génération et comparaisons de signaux harmoniques
f0=50;%fréquence du fondamental pour s1
fe=5000;% fréquence d'échantillonnage de s1
n=1024;te=1/fe;%période d'échantillonnage ou pas
t1=0:te:(n-1)*te;% durée d'observation de s1
I=50; %nombre d'itérations
s1=zeros(1,n);%initialisation
for k=1:I s1=s1+sin(2*pi*f0*(2*k)*t1)/(2*k);end
fe=100;%fréquence du fondamental pour s2
te=1/fe;t2=0:te:(n-1)*te;s2=zeros(1,n);
for k=1:I s2=s2+(sin(pi*k/2)/k)*cos(pi*k*t2);end
subplot(2,1,1);plot(t1,s1);grid on;legend('s1');subplot(2,1,2);plot(t2,s2);grid on;legend('s2');
___________________________________________________________________________
TP2 : prog3.m
% génére une sinusoïde bruitée et calcule le SNR
f=50;%fréquence de la sinusoide
a=1;%amplitude de la sinusoide
fe=2500;%fréquence d'échantillonnage
te=1/fe;n=256;t=0:te:(n-1)*te;% durée d'observation
s=a*sin(2*pi*f*t);%génération du signal
Es=sum(s.^2);%énergie du signal
m=0;%moyenne du bruit
sigma=0.1;%ecart type du bruit
randn('state',0);%initialisation du simulateur de la v.a. gaussienne
b=m+sigma*rand(1,n);%génération du bruit
y=s+b;%bruitage
Ps=Es/n;%puissance du signal
Pb=sigma^2;%puissance du bruit
SNR=10*log10(Ps/Pb);%rapport signal sur bruit
subplot(2,2,1);plot(t,s);title('sinusoide');subplot(2,2,2);plot(t,b);title('bruit');subplot(2,2,3);
plot(t,y);title('sinusoide dégradée');subplot(2,2,4);plot(SNR*ones(1,n));title('SNR');grid on;

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 F e 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)sinc[ F e (t−k .Te )]
sin( πu)
sin c (u)=
Où : π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 F max 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
.F e ∑ s(kTe )Exp [− j 2 π kn /N ]
S(f)=S( N )= k =0 , n= -N/2…N/2-1
La transformation inverse :
N /2−1
1
∑ S(nFe / N ) Exp[ j2 π kn/ N ]
s(t)=s(kTe)= N n=−N /2 , k=0…N-1
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
40
30
|S(nFe/N)|
20
10
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 f r.
 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
2
Générer le signal x(t)= sin( π .t +2 π . t+1 ) de durée une seconde avec F e=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 f 0 en
2
lui rajoutant un bruit blanc centré de puissance σ et observez si l’on
distingue le pic correspondant à la fréquence de la sinusoïde. Testez
2
plusieurs valeurs de σ 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 s m(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).
_____________________________________________________________________________________________________
Un listing des programmes : prog1.m, prog2.m, prog3.m
TPIII : TFD
___________________________________________________________________________
prog1.m
%calcul de la TFD de signaux harmoniques
n=1024;%nombre de points
te=4*pi/(n-1);%période d'échantillonnage connaissant n et 2 périodes d'observation
t=0:te:4*pi;%on observe deux périodes
fe=1/te;f=0:fe/n:(n-1)*fe/n;%on observe la TFD sur[0:fe/n:fe-fe/n]
x=sin(t);%sinusoide
y=fft(x);%sa tfd
moduly=abs(y);phasey=angle(y);
L=4;%nombre d'harmoniques additionnés
s=zeros(1,n);%initialisation de la somme
for k=0:L s=s+sin((2*k+1)*t);end%somme d'harmoniques
z=fft(s);%tfd de s
modulz=abs(z);phasez=angle(z);
%On visualise les n/8 points pour bien voir les pics
subplot(3,2,1);plot(t,x);title('sin(t)');grid on;
subplot(3,2,2);plot(t,s);title('somme de sinusoïdes');grid on;

subplot(3,2,3);plot(f(1:n/8),moduly(1:n/8));
title('spectre d''amplitude de sin(t)');grid on;
subplot(3,2,4);plot(f(1:n/8),modulz(1:n/8));title('spectre d''amplitude de la somme des sinusoïdes');grid

on;
subplot(3,2,5);plot(f(1:n/8),phasey(1:n/8));title('spectre de phase de sin(t)');grid on;
subplot(3,2,6);plot(f(1:n/8),phasez(1:n/8));title('spectre de phase de la somme des sinusoïdes');grid on;
___________________________________________________________________________
prog2.m
%TP 3 : prog2.m
%génére et calcule la tfd d'impulsions rectangulaires
clear all;%iniialisation de la mémoire
n=512;rec=10*[ones(1,n/4),zeros(1,3*n/4)];%impulsion rectangulaire de durée n/4
s=fft(rec);moduls=abs(s);%spectre de module
phases=angle(s);%spectre de phase
te=2/n;%pas d'échantillonnage
fe=1/te;%fréquence d'échantillonnage
t=0:te:(n-1)*te;%durée d'observation
f=0:fe/n:fe/2-fe/n;%on observe les spectres sur 0,fe/2
x=conv(rec,rec);%convolution de rec avec lui même
sx=fft(x);modulx=abs(sx);phasex=angle(sx);
t2=0:te:2*(n-1)*te;%la durée d'observation de la convolution est deux fois plus longue
f2=0:fe/(2*n-1):fe/2;%le pas fréquentiel change aussi, 2n-1 points au total au lieu de n
wave=square(2*pi*3*t);%signal périodique carré de fréquence 3Hz
sw=fft(wave);modulw=abs(sw);
%on visionne les n/8 points pour voir les pics sans zoomer
subplot(4,2,1);plot(t,rec);title('Impulsion Rec');grid on;

subplot(4,2,3);plot(f(1:n/8),moduls(1:n/8));title('module du spectre de Rec');grid on;
subplot(4,2,5);plot(f(1:n/8),phases(1:n/8));title('phase du spectre de Rec');grid on;
subplot(4,2,2);plot(t2,x);title('Impulsion Tri');grid on;
subplot(4,2,4);plot(f2(1:n/8),modulx(1:n/8));title('module du spectre de Tri');grid on;
subplot(4,2,6);plot(f2(1:n/8),phasex(1:n/8));title('phase du spectre de Tri');grid on;
subplot(4,2,7);plot(t,wave);title('wave:carré périodique');grid on;
subplot(4,2,8);plot(f(1:n/8),modulw(1:n/8));title('module du spectre de wave');grid on;
___________________________________________________________________________
prog3.m
%TP3 : prog3.m
%calcul de la TFD d'une sinusoïde bruitée
n=1024;%nombre de points
te=4*pi/(n-1);%période d'échantillonnage connaissant n et 2 périodes d'observation
t=0:te:4*pi;%on observe deux périodes
fe=1/te;
f=0:fe/n:(n-1)*fe/n;%on observe la TFD sur[0:fe/n:fe-fe/n]
sigma=1;%ecart type du bruit
f0=1;
x=sin(2*pi*f0*t)+sigma*randn(1,n);%sinusoide de fréquence f0 bruitée
y=fft(x);%sa tfd
moduly=abs(y);
phasey=angle(y);
%On visualise les n/8 points pour bien voir les pics
subplot(3,1,1);plot(t,x);title('sin(t)+randn(t)');grid on;
subplot(3,1,2);plot(f(1:n/8),moduly(1:n/8));title('spectre d''amplitude de sinus bruité');grid on;
subplot(3,1,3);plot(f(1:n/8),phasey(1:n/8));title('spectre de phase de sinus bruité');grid on;

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 F e.
p
On prendra soin de manipuler des signaux de tailles N=2 :
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 f0 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).
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 F e :
Δ 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.
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 :
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 [f 0- Δ f/2, f0+ Δ f/2] d’un signal s(t) à
énergie finie.
5- Quelle est la puissance ΔP localisée dans la bande [f 0- Δ 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
TP IV : ANALYSE SPECTRALE
___________________________________________________________________________
prog1.m
%TP4 : prog1.m
% calcule l'effet de l'ordre de la TFD sur le spectre
clear all;
f0=100;%fréquence de la sinusoïde
te=1/fe;
ordre=[16 32 64 128 256 512];%les ordres m de la TFD à étudier

for k=1:6
m=ordre(k);%ordre de la TFD
t=0:te:(m-1)*te;%durée observée
f=0:fe/m:fe/2-fe/m;%TFD observée sur 0,fe/2
s=sin(2*pi*f0*t);%génération de la sinusoïde
x=fft(s,m);
y=abs(x);%spectre de module
subplot(2,3,k);plot(f,y(1:m/2));title(['spectre d''ordre M=' num2str(m,'%d')]);grid on;
end
___________________________________________________________________________
prog2.m
%TP 4: prog2.m
% analyse spectrale par fenêtrage de hanning et blackman
clear all;
n=512;
han=hanning(n);%fenêtre de hanning:c'est un vecteur colonne
blak=blackman(n);%fenêtre de blackman
f=0:1/n:1/2-1/n;%observation sur [0,1/2Hz]
mhan=abs(fft(han));
mblak=abs(fft(blak));
subplot(4,2,1);plot(han);title('fenêtre de hanning');grid on;
subplot(4,2,2);plot(blak);title('fenêtre de blackman');grid on;
subplot(4,2,3);plot(f,mhan(1:n/2));title('spectre de hanning');grid on;
subplot(4,2,4);plot(f,mblak(1:n/2));title('spectre de blackman');grid on;
f0=100;%fréquence de la sinusoïde
te=1/fe;t=0:te:(n-1)*te;%durée observée
s=sin(2*pi*f0*t);clear f;%initialisation de f
f=0:fe/n:fe/2-fe/n;%observation du spectre sur [0, fe/2]
sh=s.*han';%fenêtrage
sb=s.*blak';%fenêtrage
msh=abs(fft(sh));msb=abs(fft(sb));
subplot(4,2,5);plot(sh);title('sin(t)*hanning');grid on;subplot(4,2,6);plot(sb);title('sin(t)*blackman');grid
on;subplot(4,2,7);plot(f,msh(1:n/2));title('spectre de sin*hanning');grid on;
subplot(4,2,8);plot(f,msb(1:n/2));title('spectre de sin*blackman');grid on;
___________________________________________________________________________
prog3.m
%TP4 : prog3.m
% Etude de la résolution spectrale de la méthode d'analyse
% par TFD (pouvoir de séparation)
clear all;f0=100;%fréquence de la sinusoïde 1
f1=105;%fréquence prôche de la sinusoïde 2
te=1/fe;
for k=5:8% calcul des TFD d'ordres 256, 512, 1024, 2048
n=2^(3+k);%ordre de la TFD
t=0:te:(n-1)*te;%durée observée
f=0:fe/n:fe/4-fe/n;%TFD observée sur 0,fe/4
s=sin(2*pi*f0*t)+sin(2*pi*f1*t);%génération de la somme des deux sinusoïdes
x=fft(s,n);%TFD de la somme
y=abs(x);%spectre de module
%Affichage
subplot(4,2,2*(k-5)+1);plot(t,s);title('somme des deux sinusoides');grid on;
subplot(4,2,2*(k-5)+2);plot(f,y(1:n/4));title(['TFD d''ordre M = ' num2str(n,'%d')]);grid on;
end
___________________________________________________________________________
prog4.m
%TP4 : prog4.m
%calcul du spectre à partir de la fonction d'autocorrélation
clear all;n=2048;rec=[ones(1,n/4),zeros(1,3*n/4)];%génération du signl rectangulaire
R=xcorr(rec,rec,'biased');%autocorrélation normalisée biaisée
X=fft(rec);% TFD de rec
mx=abs(X);
mx2=mx.^2;%module au carré de la TFD de rec:sa densité spectrale d'énergie
sr=fft(R);%TFD de R
mr=abs(sr);randn('state',0);%initialisation du générateur aléatoire
b=randn(1,n);%bruit blan gaussien
B=xcorr(b,b,'biased');%autocorrélation du bruit
sb=fft(B);% dsp du bruit
mb=abs(sb);
%affichage
subplot(4,2,1);plot(rec);title('rec : impulsion rectangulaire');grid on;
subplot(4,2,3);plot(mx(1:n/16));title('TFD de rec');grid on;
subplot(4,2,5);plot(R);title('R : autocorrélation de rec');grid on;
subplot(4,2,7);plot(mr(1:n/16));title('TFD de R');grid on;
subplot(4,2,8);plot(mx2(1:n/16));title('TFD au carré de rec');grid on;
subplot(4,2,2);plot(b);title('b : bruit blanc');grid on;subplot(4,2,4);plot(B);title('B : autocorrélation de

b');grid on;
subplot(4,2,6);plot(mb(1:n/16));title('TFD de B');grid on;
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
PdB=10.log 10
P étant la puissance par cm2, la puissance en dB s’écrit : P0 , avec
P0=10-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=
ln(1+μ ) , 0 ¿ X ¿ 1 (standard US, ex. μ =255)
1+ln( A . X )
Y=
1+ln( A ) , 1/A ¿ X ¿ 1 (standard EU, ex. A=87.6)
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 .
___________________________________________________________________
doit contenir un listing des programmes : prog1.m, prog2.m, prog3.m
TP V : NOTIONS d’AUDIO
%enregistrement d'un son wav
clear all;fs=44100 ; %fréquence d’échantillonnage
D=8 ; %Durée de l’enregistrement en secondes
input('taper un caractere pour enregistrer:');
input('taper un caractere pour ecouter l''enregistrement:');
play(x);
% Ecouter le son enregistré
___________________________________________________________________________
%lecture d'un son
input('ENTRER UN CARACTERE POUR ECOUTER UN SON WAV:');
son='helicoptere.wav';%nom du fichier son
[x,fs,bits]=wavread(son);%lecture des donnees, frequence et codage
wavplay(x,fs);%ecoute
input('ENTRER UN CARACTERE POUR connaitre le nombre de samples et de canaux:');
siz=wavread(son,'size');%
fprintf('nombre d''échantillons:');siz(1)
fprintf('nombre de canaux:');siz(2)
fprintf('durée du son:');siz(1)/fs
input('ENTRER UN CARACTERE POUR connaitre la taille de l''entête:');
fid=fopen(nom,'r');%ouverture du fichier pour lecture
son_lu=fread(fid,inf,'uint8');%lecture de tout les bytes
fclose(fid);
entete=size(son_lu,1)-siz(1)*siz(2)*bits/8%taille en byte de l'entête
___________________________________________________________________________
prog1.m
%synthèse de sons purs et harmoniques
f0=1000;%le son fondamental
fe=44100;%la fréquence d'échantillonnage
te=1/fe;%pas d'échantillonnage
D=100/f0;%durée du son
n=D*fe+1;%nombre de samples
t=0:te:D;%durée du son
x=sin(2*pi*f0*t);%le son pur
input('taper return pour ecouter le fondamental');

wavplay(x);%écouter
y=x;for k=1:10 y=y+sin(2*pi*k*f0*t);end%sommation des harmoniques
input('taper return pour écouter les harmoniques');
wavplay(y);
___________________________________________________________________________
prog2.m
%enregistrement audio
clear all;fs=44100 ; %fréquence d’échantillonnage
D=16; %Durée de l’enregistrement en secondes
input('taper un caractere pour ecouter l''enregistrement:');
play(x);
% Ecouter le son enregistré
n=2048;y1=y(1:n,1);%extrait canal 1
y2=y(1:n,2);%extrait canal 2
subplot(3,2,1);plot(y1(1:n));title('extrait de l''enregistrement');
subplot(3,2,3);plot(abs(fft(y1)));title('spectre d''amplitude');
subplot(3,2,5);plot(angle(fft(y1)));title('spectre de phase');
subplot(3,2,2);plot(y2(1:n));title('extrait de l''enregistrement');
subplot(3,2,4);plot(abs(fft(y2)));title('spectre d''amplitude');
subplot(3,2,6);plot(angle(fft(y2)));title('spectre de phase');
%modification de l'enregistrement
L=size(y,1);
z=[y(1:L/4,:);y(L/2:L,:)];%on supprime L/4 échantillons
input('taper un caractere pour ecouter l''enregistrement modifié:');
wavplay(z,fs);
_________________________________________________________________________
prog3.m
%TP5 prog3.m
%mixage audio
clear all;
fs=44100 ; %fréquence d’échantillonnage
D=8; %Durée de l’enregistrement en secondes
x1=audiorecorder(fs,nbits,mode) ;% x1 est une variable objet
recordblocking(x1,D) ;
y1=getaudiodata(x1) ; % y1 contient les données son de type double
x2=audiorecorder(fs,nbits,mode) ;% x2 est une variable objet
recordblocking(x2,D) ;
y2=getaudiodata(x2) ; % y2 contient les données son de type double
%modification de l'enregistrement
L=size(y1,1);%nombre d'échantillons
z=[y1(1:L/4,:);y2(1:L/4,:);y1(L/4+1:L/2,:);y2(L/4+1:L/2,:)];%mixage
input('taper un caractere pour ecouter le mixage:');
wavplay(z,fs);
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 s m(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).
~x (t )=x (t )+ j. y(t )
On appelle signal analytique de x(t) le signal complexe :
~x (t )
1- quelle est la T.F. du signal analytique ? 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é.
~ j.ϕ(t )
La forme polaire du signal analytique : x (t )=r(t ).e
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
− j . ω0 t
~ ~
par : r (t )= x (t ). e
On peut interpréter l’enveloppe complexe comme le résultat d’une modulation
d’amplitude du signal analytique.
~r (t )
5- Exprimer les parties réelles et imaginaires a(t) et b(t) de 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 f 1 ≤|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.
_____________________________________________________________________________________________________
doit contenir en plus de la partie théorique : Un listing des programmes prog1.m et
prog2.m
TP VI : MODULATION
prog1.m
%modulation analogique
f0=50;%fréquence du signal original
te=1/fe;
n=512;%ombre de points
t=0:te:te*(n-1);%durée d'observation
x=sin(2*pi*f0*t);%signal à moduler
fp=2500;%fréquence porteuse
p=cos(2*pi*fp*t);%la porteuse
m=(x+2).*p;%signal modulé m=0.5
dm=m.*p;%signal démodulé non filtré
%les spectres d'amplitude
sx=abs(fft(x));
sp=abs(fft(p));
sm=abs(fft(m));
sdm=abs(fft(dm));
f=0:fe/n:fe/2-fe/n;
%visualisation
subplot(4,2,1);plot(x);title('signal original');grid on;
subplot(4,2,2);plot(f,sx(1:n/2));title('le spectre original');grid on;
subplot(4,2,3);plot(p);title('la porteuse');grid on;
subplot(4,2,4);plot(f,sp(1:n/2));title('spectre de la porteuse');grid on;
subplot(4,2,5);plot(m);title('signal modulé');grid on;
subplot(4,2,6);plot(f,sm(1:n/2));title('spectre du signal modulé');grid on;
subplot(4,2,7);plot(dm);title('signal démodulé non filtré');grid on;
subplot(4,2,8);plot(f,sdm(1:n/2));title('spectre démodulé non filtré');grid on;
___________________________________________________________________________
prog2.m
%Modulation démodulation analogique d'amplitude
% d'un signal audio
clear all;
nom=input('nom du son wav à moduler ?');
[x,fe]=wavread(nom);
sx=abs(fft(x));%spectre original
te=1/fe;
n=size(x,1);%taille du son
T=te*(n-1);
fp=5000;%fréquence porteuse
t=0:te:T;%durée du son
y=ammod(x,fp,fe);%modulation d'amplitude
sy=abs(fft(y));%spectre du son modulé
wavwrite(y,fe,'son_module.wav');%le son modulé
z=amdemod(y,fp,fe);%démodulation d'amplitude
sz=abs(fft(z));%spectre du son démodulé
wavwrite(z,fe,'son_demodule.wav');%le son démodulé
f=0:fe/n:fe/2-fe/n;%plage spectrale
input('ENTRER UN CARACTERE POUR ECOUTER LE SON ORIGINAL:');
wavplay(x,fe);
input('ENTRER UN CARACTERE POUR ECOUTER LE SON MODULE:');
wavplay(y,fe);
input('ENTRER UN CARACTERE POUR ECOUTER LE SON DEMODULE:');
wavplay(z,fe);
input('ENTRER UN CARACTERE POUR VISUALISER LES SIGNAUX et leurs spectres:');
subplot(3,2,1);plot(x(1024:2*1024));title('signal audio');
subplot(3,2,2);plot(f,sx(1:n/2));title('spectre du signal audio');
subplot(3,2,3);plot(y(1:512));title('signal audio modulé');
subplot(3,2,4);plot(f,sy(1:n/2));title('spectre du signal modulé');
subplot(3,2,5);plot(z(1:512));title('signal audio démodulé');
subplot(3,2,6);plot(f,sz(1:n/2));title('spectre du signal démodulé');

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):
N
1
μ^ ( 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):
N
1
σ^ 2 ( N )= . ∑ [ x (i)− μ^ (N )]2
1 N i=1
N
1
σ^ 2 ( N )= . ∑ [ x (i )− μ^ ( N )]2
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.
PARTIE THEORIQUE
1- Donner le biais et la variance de l’estimateur de la moyenne empirique :

N
1
μ^ ( N )= . ∑ X (i)
N i =1
2- Donner le biais et la variance de l’estimateur de la variance empirique :

N
1
σ^ 2 ( N )= . ∑ [ x (i)− μ^ (N )]2
1 N i=1
3- Estimation empirique de la fonction de corrélation
Soit N échantillons [x1,…xN] extraits d’un signal stationnaire X centré d’ordre deux de fonction de
corrélation inconnue RX(n)=E[x(k).x(k+n)].
On peut généraliser l’estimation empirique pour n=0…N-1, sous la forme :

N −n
1
R^ X (n )= R^ X (−n)= . ∑ x .x
N −n k=1 k k +n
a- Montrer que cet estimateur est non biaisé.

b- Considérons le nouvel estimateur empirique :
N−n
1
R^ X (n )= R^ X (−n)= . ∑ x k . x k +n
N k =1
Calculer le biais de cette estimateur et montrer que la matrice de corrélation est
définie positive.
_____________________________________________________________________________________________________
Un listing du programme : prog1.m, prog2.m, prog3.m
TP VII : ESTIMATION et DETECTION
prog1.m
% Estime la moyenne arithmétique d'un signal
n=1000;%nombre d'échantillon du signal
T=2*pi;%période du signal
fe=(n-1)/(4*T);%fréquence d'chantillonnage pour observer 4 période
te=1/fe;
t=0:te:4*T;%durée observée
s=sin(2*pi*t/T);%le signal
randn('state',0);%initialisation du générateur aléatoire
mu=5;%moyenne du bruit
b=mu+randn(1,n);%le bruit
x=s+b;%signal bruité
for k=1:n %calcul de la moyenne pour k échantillons
moyenne(k)=mean(x(1:k));
end
e=moyenne-5;%le biais
eqd=e.^2;%l'erreur quadratique
subplot(2,2,1);plot(t,x);
title(['signal x de moyenne mu = ' num2str(mu,'%d')]);grid on;
subplot(2,2,2);plot(moyenne);
title('moyenne arithmétique de x');xlabel('N : nombre de points utilisés');grid on;
subplot(2,2,3);plot(e); title(['Biais de l''estimateur de moyenne arithmétique B='

num2str(mean(e),'%d')]);
xlabel('N');grid on;subplot(2,2,4);plot(eqd);
title(['Erreur quadratique de moyenne arithmétique V=' num2str(mean(eqd),'%d')]);
xlabel('N');grid on;
___________________________________________________________________________
prog2.m
% Estime la variance empirique d'un signal
clear all;
T=2*pi;%période du signal
a=1;%amplitude du signal
fe=(n-1)/(4*T);%fréquence d'chantillonnage pour observer 4 période
te=1/fe;
t=0:te:4*T;%durée observée
s=a*sin(2*pi*t/T);%le signal
sigma=sqrt(1.5);%écart type du bruit
b=mu+sigma*randn(1,n);%le bruit
x=s+b;%signal bruité
Px=sigma^2+a^2/2;%puissance du signal bruité (sa variance)
variance1(1)=0;%initialisation
variance2(1)=0;
for k=2:n %calcul de la variance pour k échantillons
variance1(k)=mean((x(1:k)-mean(x(1:k))).^2);
variance2(k)=k/(k-1)*variance1(k);
end
e1=variance1-Px;%le biais
eqd1=e1.^2;%l'erreur quadratique
e2=variance2-Px;%le biais
eqd2=e2.^2;%l'erreur quadratique
%premier estimateur
subplot(2,4,1);plot(t,x);
title(['x de variance = ' num2str(Px,'%d')]);grid on;
subplot(2,4,2);plot(variance1);
title('variance1');xlabel('N');grid on;
subplot(2,4,3);plot(e1);
title(['Biais 1 =' num2str(mean(e1),'%d')]);
subplot(2,4,4);plot(eqd1);
title(['EQ V1=' num2str(mean(eqd1),'%d')]);
%second estimateur
subplot(2,4,5);plot(variance2);
title('variance2');xlabel('N');grid on;
subplot(2,4,6);plot(e2);
title(['Biais 2 =' num2str(mean(e2),'%d')]);
subplot(2,4,7);plot(eqd2);
title(['EQ V2=' num2str(mean(eqd2),'%d')]);
___________________________________________________________________________
prog3.m
% Détection d'une composante continue par moyennage empirique
T=2;%durée du signal observé
fe=(n-1)/T;%fréquence d'chantillonnage pour observer 4 période
te=1/fe;
t=0:te:T;%durée observée
C=1;%la composante continue à détecter
v=C/2;
b=mu+sqrt(v)*randn(1,n);%le bruit
x=C+b;%le signal bruité
Ce=mean(x);
e=Ce-C;%le biais
eqd=e.^2;%l'erreur quadratique
fprintf('Valeur estimée : %d\n',Ce);
fprintf('Biais : %d\n',e);
fprintf('EQ : %d\n',eqd);
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)|


 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 z−m
m=0
N
∑ a n z−n
H(z)= n=0 , b=[b0,..,bM], a=[a0,..,aN]
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 :
−1 −1 −5
z z −0 . 5 z
−1 −1 −2 −1 −1
0.5.(1-z ), 0.5.(1+ z ), 0.5.(1- z ), 1−0 .5 z , 1−0. 5 z
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 =−
dω , avec θ(ω)=angleH (ω)
θ(ω)
τ ϕ =−
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 Δ tra Δ Δ at
ns ond t
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 F e=8kHz et N=1024.
Visualisez le signal et son spectre.
 Réalisez le filtre adéquat pour extraire les harmoniques 5.f 0 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 :

−1 −2
H(z)=1+2 z +3 z
−1
1) Montrer que son auto-corrélation G(z)=H(z)H( z ) est à phase linéaire.
2) Donner le group delay.
3) Déterminer la réponse impulsionnelle g(n) de G(z).
_____________________________________________________________________________________________________
Un listing des programmes : prog1.m, prog2.m, prog3.m, prog4.m
TP VIII : FILTRAGE NUMERIQUE
___________________________________________________________________________
prog1.m
% filtrage numérique
%réponses impulsionnelles
clear all
n=32;%nombre d'échantillons
d=zeros(1,n);%impulsion de Dirac
d(1)=1;
%les RI des filtres définis par les bi et ai
h(1,:)=filter([0.5 -0.5],[1 0],d);
h(2,:)=filter([0.5 0.5],[1 0],d);
h(3,:)=filter([0.5 0 -0.5],[1 0 0],d);

h(4,:)=filter([0 1],[1 -0.5],d);
h(5,:)=filter([0 1 0 0 0 -0.5],[1 -0.5 0 0 0 0],d);
%visualisation des RI
c=input('Taper un caractère pour voir la réponse impulsionnelle');
for i=1:5
subplot(2,3,i);stem(h(i,:));title(['h' num2str(i,'%d')]);
end
%visualisation des RF
c=input('Taper un caractère pour voir la réponse fréquentielle');
for i=1:5
H(i,:)=abs(fft(h(i,:)));
subplot(2,3,i);plot([-1/2+1/n:1/n:1/2],H(i,:));title(['H' num2str(i,'%d')]);grid on
end
c=input('Taper un caractère pour voir les zéros et pôles');
%Calcul des pôles et zéros
subplot(2,3,1);zplane([0.5 -0.5],[1 0]);title('Pôles et zéros de h1');
subplot(2,3,2);zplane([0.5 0.5],[1 0]);title('Pôles et zéros de h2');
subplot(2,3,3);zplane([0.5 0 -0.5],[1 0 0]);title('Pôles et zéros de h3');
subplot(2,3,4);zplane([0 1],[1 -0.5]);title('Pôles et zéros de h4');
subplot(2,3,5);zplane([0 1 0 0 0 -0.5],[1 -0.5 0 0 0 0]);title('Pôles et zéros de h5');
%visualisation du retard de groupe
c=input('Taper un caractère pour afficher le retard de groupe');
for i=1:5
subplot(2,3,i);plot(grpdelay(h(i,:)));title(['Retard de groupe de H' num2str(i,'%d')]);grid on
end
%calcul de la réponse à un harmonique

c=input('Taper un caractère pour voir la réponse à une sinusoïde de période 100');
T=100;
t=0:2*T;
x=sin(2*pi*t/T);
y(1,:)=filter([0.5 -0.5],[1 0],x);
y(2,:)=filter([0.5 0.5],[1 0],x);
y(3,:)=filter([0.5 0 -0.5],[1 0 0],x);
y(4,:)=filter([0 1],[1 -0.5],x);
y(5,:)=filter([0 1 0 0 0 -0.5],[1 -0.5 0 0 0 0],x);
%visualisation de la réponse à un harmonique
subplot(2,3,1);plot(t,x);title('entrée : sinus de période 100');grid on;
for i=1:5
subplot(2,3,i+1);plot(y(i,:));title(['réponse de h' num2str(i,'%d')]);grid on;
end
c=input('Taper un caractère pour voir la réponse indicielle');
%calcul de la réponse indicielle
ech=ones(1,n);%échelon unité
d(1,:)=filter([0.5 -0.5],[1 0],ech);
d(2,:)=filter([0.5 0.5],[1 0],ech);
d(3,:)=filter([0.5 0 -0.5],[1 0 0],ech);
d(4,:)=filter([0 1],[1 -0.5],ech);
d(5,:)=filter([0 1 0 0 0 -0.5],[1 -0.5 0 0 0 0],ech);
%visualisation
for i=1:5
subplot(2,3,i);plot(d(i,:));title(['réponse indicielle de h' num2str(i,'%d')]);
end
__________________________________________________________________________
prog2.m
% Synthèse par TFD d'un filtre passe bas
clear all;
fe=8000;%fréquence d'échantillonnage donnée
fc=1000;%fréquence de coupure imposée
dtrns=200;%largeur de la bande de transition
ft=fc+dtrns;%la première fréquence de la bande d'atténuation
dond=1;%amplitude de l'ondulation dans la bande passante en dB
min=10^(-dond/20);%c'est la valeur à la fréquence fc de |H(f)|
datt=-40;%atténuation relative minimale dans la bande atténuée en dB
max=10^(datt/20);
inter=(max+min)/2;%une valeur intermédiaire dans la bande de transition pour adoucir
n=(9-1)*fe/fc;% 64 échantillons : fc correspond au 9ieme échantillons
f=[0:fe/n:fe/2-fe/n];%axe des fréquences 0,fe/2
%échantillonnage du gabarit : prélévement de n échantillons
Hg=[ones(1,8) min inter max max/2 max/4 zeros(1,n-26) max/4 max/2 max inter min ones(1,8)];
%FFT inverse du gabarit échantillonné
h=real(ifft(Hg));%la partie réelle sera la RI estimée
H=fft(h);%FFT de la partie réelle tronquée puisque périodique sera la RF estimée
mH=abs(H);
pH=angle(H);
%visualisations
subplot(2,2,1);stem(f,Hg(1:n/2));title('le gabarit de |H(f)|');
subplot(2,2,2);stem(h);title('RI estimée : partie réelle décalée de la fft inverse');

subplot(2,2,3);plot(f,mH(1:n/2));title('RF estimée : module de la fft de la partie réelle');grid on;
subplot(2,2,4);plot(f,pH(1:n/2));title('phase de la fft de la partie réelle');grid on;
%pôles et zéros
c=input('taper un caractère pour voir les pôles et zéros');
subplot(1,1,1);zplane(h);title('pôles et zéros de la RI estimée');
%Filtrage d'harmoniques
c=input('taper un caractère pour voir les réponses à deux harmoniques');
te=1/fe;%pas d'échantillonnage
f1=400;%fréquence dans la bande passante
e1=sin(2*pi*f1*t);
f2=2500;%fréquence dans la bande atténuée
e2=sin(2*pi*f2*t);
s1=real(ifft(H.*fft(e1)));%réponse par tfd
s2=real(ifft(H.*fft(e2)));
%affichage
subplot(2,2,1);plot(t,e1);title('sinus de fréquence 400<fc');grid on;
subplot(2,2,2);plot(t,s1);title('réponse du filtre');grid on;
subplot(2,2,3);plot(t,e2);title('sinus de fréquence 2000>fc');grid on;
subplot(2,2,4);plot(t,s2);title('réponse du filtre');grid on;
________________________________________________________________________
prog3.m
%Filtrage d'harmoniques
te=1/fe;
f0=100;%le fondamental
n=1024;
x=zeros(1,n);%initialisation
for k=0:5 x=x+sin(2*pi*(2*k+1)*f0*t);end;%somme des harmoniques impaires
f=0:fe/n:fe/2-fe/n;%bande 0,fe/2
sx=abs(fft(x));%module du spectre
%visualisation
subplot(3,2,1);plot(t,x);title('somme d''harmoniques impaires');grid on;
subplot(3,2,2);plot(f,sx(1:n/2));
title('spectre de la somme des harmoniques impaires');grid on;
%gabarit pour extraire le fondamental et deux harmoniques 5 et 11
Hg=[zeros(1,12) 1 1 1 zeros(1,48) 1 1 1 zeros(1,73) 1 1 1 zeros(1,n-142)];
h=real(ifft(Hg));%réponse impulsionnelle estimée
H=fft(h);%RF estimée
mH=abs(H);
xf=conv(h,x);%signal filtré
Xf=abs(H.*fft(x));%spectre du signal filtré
%visualisation
subplot(3,2,3);plot(f,mH(1:n/2));title('RF du filtre extracteur');grid on;
subplot(3,2,4);plot(f,Xf(1:n/2));title('harmoniques extraits');grid on;
subplot(3,2,5);plot(t,xf(n:2*n-1));title('signal extrait');grid on;
subplot(3,2,6);stem(t,h);title('RI du filtre');grid on;
___________________________________________________________________________
prog4.m
%Filtrage par moyenne mobile d'une sinusoïde bruitée

clear all;
te=1/fe;
f0=1000;%le fondamental
n=256;
x=sin(2*pi*f0*t);%sinusoide
rsb=10;%rapport signal sur bruit
sigma=sqrt(var(x)*(10^-(rsb/10)));%ecart type du bruit déduit du snr
b=sigma*randn(1,n);%bruit
y=x+b;%signal bruité
sx=abs(fft(x));%spectre original
f=0:fe/n:fe/2-fe/n;
subplot(3,2,1);plot(t,x);title('sinusoide');
subplot(3,2,2);plot(t,b);title('bruit');
subplot(3,2,3);plot(f,sx(1:n/2));title('spectre de la sinusoide');
subplot(3,2,4);plot(t,y);title('sinusoide bruitée');
m=5;%taille de la fenêtre de moyennage
z=filter(ones(1,m)/m,1,y);%moyennage mobile
subplot(3,2,5);plot(t,z);title('sinusoide débruitée');
sz=abs(fft(z));
subplot(3,2,6);plot(f,sz(1:n/2));title('spectre débruité');
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
8
entre 0 et 255=256-1=2 -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
10
octets (un Kilo octets=1Kb=1024 octets=2 octets):
N Mémoire en bytes Mémoire en Kby
128 16.384 16
256 65.536 64
512 262.144 256
102 1.048.576 1024=1Mby

4
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 a 1,3
[ a2,1 a2,2 a 2,3

a3,1 a3,2 a 3,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 )= ∑ ∑ am , 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
____________________________________________________________________________________________________
doit contenir :
Un listing des programmes : prog1.m, prog2.m, prog3.m, prog4.m, prog5.m
TP IX : NOTIONS de TRAITEMENT d’IMAGES
___________________________________________________________________________
prog1.m
%calcul des deux premiers moments statistiques
i=imread('cameraman.tif');%lecture de l'image
x=double(i);%conversion du type pour manipuler
sup=max(x(:));%niveau max
inf=min(x(:));%niveau min
m=mean(x(:));%moyenne
mc=mean(x);%moyenne des colonnes
ml=mean(x');%moyenne des lignes
d=(x-m).^2;%matrice d'écart quadratique
v=mean(d(:));%variance
s=sqrt(v);%ecart type
fprintf('moyenne arithmétique de l''image : %d\n',m);
fprintf('Ecart-type de l''image : %d\n',s);
fprintf('Niveau de gris maximal de l''image : %d\n',sup);
fprintf('Niveau de gris minimal de l''image : %d\n',inf);
___________________________________________________________________________
prog2.m
%TP 9 prog2.m
%calcul de la moyenne locale 3x3 d'une image
clear all;
isource=imread('cameraman.tif');%lecture de l'image
x=double(isource);%conversion du type pour manipuler
n=256;%taille de l'image
for i=2:n-1%balayage de la ligne i
for j=2:n-1%balayage de la colonne j
w=[x(i-1:i+1,j-1:j+1)];%le masque de taille 3x3

mx(i-1,j-1)=mean(w(:));%le filtrage (une moyenne arithmétique)
end
end
%affichage
subplot(1,2,1);imshow(isource);title('image originale');
subplot(1,2,2);imshow(uint8(mx));title('image filtrée par moyennage');
___________________________________________________________________________
prog3.m
%TP 9 prog3.m
%Filtrage median
clear all;
n=256;%taille de l'image
for i=2:n-1%balayage de la ligne i
for j=2:n-1%balayage de la colonne j
w=[x(i-1:i+1,j-1:j+1)];%le masque de taille 3x3
mx(i-1,j-1)=median(w(:));%le filtrage (une moyenne arithmétique)
end
end
%affichage
subplot(1,2,2);imshow(uint8(mx));title('image filtrée par MEDIAN');
___________________________________________________________________________
prog4.m
%Seuillage d'une image

clear all;
s=input('seuil : ');%valeur du seuil
m=mean(x(:));%moyenne de l'image
sx=255*(x>s);%seuillage par s
sx2=255*(x>m);%seuillage par m
%affichage
subplot(2,2,2);imshow(uint8(sx));title(['image seuillée s= ' num2str(s,'%d')]);
subplot(2,2,3);imshow(uint8(sx2));title(['image seuillée par la moyenne= ' num2str(m,'%d')]);
___________________________________________________________________________
prog5.m
%calcule le contour d'une image par le gradient
clear all;
n=256;%taille originale
%calcul du gradient horizontal et vertical et leur moyenne
for i=1:n-1
for j=1:n-1
ch(i,j)=abs(x(i,j+1)-x(i,j));%gradient horizontal
cv(i,j)=abs(x(i+1,j)-x(i,j));%gradient vertical
end
end
mch=mean(ch(:));%moyenne pour seuiller et binariser

mcv=mean(cv(:));
ym=(ch+cv)/2;%le gradient moyen
%binarisation par seuillage avec la moyenne
ch1=255*(ch>mch);%binarisation du gradient h
cv1=255*(cv>mcv);%binarisation du gradient v
ym1=255*(ym>(mch+mcv)/2);%binarisation du gradient moyenne
%affichage
subplot(2,4,2);imshow(uint8(ch));title('gradient horizontal');
subplot(2,4,3);imshow(uint8(cv));title('gradient vertical');
subplot(2,4,4);imshow(uint8(ym));title('moyenne des deux gradients');
subplot(2,4,5);imshow(uint8(ch1));title('gradient horizontal binarisé');
subplot(2,4,6);imshow(uint8(cv1));title('gradient vertical binarisé');
subplot(2,4,7);imshow(uint8(ym1));title('moyenne des deux gradients binarisée');

Signaux Matlab2

Transféré par

Droits d'auteur :

Formats disponibles

Signaux Matlab2

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

Signaux Matlab2

Transféré par

Droits d'auteur :

Formats disponibles

TRAITEMENT du SIGNAL

SOLUTIONS des TP sous MATLAB

I- Manipulation des variables

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

2- Taille des variables dans la mémoire

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

5- Les matrices et les variables vectorielles

6- Opérations sur les vecteurs

Génération automatique d’un vecteur

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

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

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

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

Saisie d’une donnée au clavier :

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

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 :

r=0; %initialisation de la moyenne

Testez cette fonction sur un signal sinusoïdal d’amplitude A et de période 1.

Exemple 3 : Diagramme de bode.

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

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

TPI : INITIATION à MATLAB

%TP1 : fonction factoriel.m

%calcul du factoriel d’un nombre naturel

% tracage des réponses en fréquence de deux filtres

%utilisant la fonction spectre(f,H)

fc=input(‘fréquence de coupure à -3dB (entre 0 et 100 Hz) ?’) ;

spectre(f,H1,1) ;% l’indice 1 pour tracer les deux spectres de H1 en haut

spectre(f,H2,3) ;% l’indice 3 pour tracer les deux spectres de H2 en bas

%TP1 : prog1.m, 2008

%tracage des fonctions sin et cos sur la même figure

%TP1 : prog2.m, 2008

%tracage des fonctions sin, cos, log10, exp

% calcule exp(2) avec une précision > 95%

while ((exp(2)-s)/exp(2)>0.05) i=i+1 ;% itération suivante s=somme(2,i) ; end

fprintf(‘valeur de exp(2) : %d\n’,s) ;

fprintf(‘nombre d’’itérations effectuées : %d\n’,i) ;

% TP1 : fonction puissance.m

%calcul de la puissance moyenne et de la valeur efficace d’un signal

s=sqrt(p) ;%la valeur efficace de v, ou son écart type

%TP1 : fonction somme.m

% calcule le développement de Taylor de Exp(x)

function s=somme(x,n) s=1;for k=1:n s=s+x^k/factoriel(k);end

% TP1: fonction spectre.m

% tracage des spectres de phase et d'amplitude d'une réponse fréquentielle

xlabel('fréquence (Hz)');ylabel('20Log10(|H|)');grid on;

subplot(2,2,i+1);plot(f,angle(H));title('spectre de phase en radian');

xlabel('fréquence (Hz)');ylabel('arg(H)');grid on;

II- Visualisation du signal s

Exemple 1 : Signal rectangulaire

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

III- NOTION de BRUIT

IV- Rapport Signal sur Bruit (SNR )

Exemple 1 : sinusoïde bruitée

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

for k=1:2n-1 % calcule l'échantillon z(k)=xx(k)

for k=1:I s1=s1+sin(2pif0(2k)t1)/(2k);end

for k=1:I s2=s2+(sin(pik/2)/k)cos(pikt2);end

for k=0:L s=s+sin((2k+1)t);end%somme d'harmoniques

n=512;rec=10[ones(1,n/4),zeros(1,3n/4)];%impulsion rectangulaire de durée n/4

t2=0:te:2(n-1)te;%la durée d'observation de la convolution est deux fois plus longue

wave=square(2pi3*t);%signal périodique carré de fréquence 3Hz