TP1-Signal
TP1-Signal
TP1-Signal
1. Le chapitre 10 de cet ouvrage (avec les tests et exercices correspondants) me servira très
souvent dans ces TP.
2. Attention : mettre ; derrière une instruction vous évite de voir l’affichage (c’est exactement
le contraire de ce qui se passe avec Maple12 !) ; ne rien mettre derrière vous fait au contraire courir
le risque de voir s’afficher un résultat, ce qui n’est en règle générale pas souhaitable, surtout lorsque
les variables en jeu sont de grands tableaux.
1
2 TP1 : FAMILIARISATION AVEC MATLAB : SIGNAUX ET IMAGES DIGITALES, DFT
déjà effectuée : tapez les premières lettres et vous verrez apparaitre l’ancien
modèle que vous pourrez alors corriger).
– Command History : c’est ici que s’affice l’historique des commandes. C’est
une fenêtre que vous pouvez fermer pour ne pas encombrer votre écran de
console.
– Current Directory : vous voyez dans cet espace le contenu du répertoire dans
lequel vous vous êtes placé (qui devrait être toujours le répertoire TPSignal
que vous venez de créer). C’est en particulier dans ce répertoire que se trou-
veront tous les fichiers .m correspondant à des script (par exemple com-
mençant par [.,., ...,.] = function (.,...,.)) que vous aurez à écrire
ou à utiliser, les fichiers .mat correspondant aux sessions ou aux variables que
vous sauverez éventuellement, ainsi que les fichiers de données numériques
(tableaux de nombres, images jpeg, etc.) préalablement téléchargés.
– Workspace : ici se trouvent recensées les variables actuellement affectées dans
votre présente session de travail.
Les commandes Unix sont valables sous l’environnement MATLAB ; reportez vous au
fascicule de TP 2011-2012 pour l’UE N1MA3003 (page 12), pour un rappel de ces
commandes avec les touches Ctrl + [...]. Ce fascicule est téléchargeable ici :
http://www.math.u-bordeaux1.fr/~ayger/initMS.pdf
Le logiciel MATLAB installé sur les machines du CREMI est augmenté des cinq
toolboxes : Signal Processing, Image Processing, Wavelets, PDE, Statistique.
Les trois premiers pourront être mis à contribution dans ces TP. Est présent aussi
le toolbox libre Wavelab rapatrié depuis ce site
http://www-stat.stanford.edu/~wavelab/Wavelab_850/download.html
Des routines prises dans ce package pourront s’avérer utiles. Pour le rendre opérationnel,
ajouter sous l’onglet Files/Set Path :
/opt/local/toolbox/Wavelab850
Pour vous assurer de la présence de ce package, lancer les commandes suivantes :
>> WavePath
>> S1 = MakeSignal (’Piece-Regular’,1024);
>> S2 = MakeSignal (’Blocks’,1024) ;
>> S3 = MakeSignal (’LinChirp’,1024);
>> S4 = MakeSignal (’TwoChirp’,1024);
>> S5 = MakeSignal (’Chirps’,1024);
>> S6 = MakeGabor (’Gabor’,1024);
>> plot(S1) ; plot (S2) ; .... ; plot(S6);
Comme vous le constatez, la routine MakeSignal nous servira à générer des signaux
1D.
Quelques tests préliminaires. Au contraire de Maple12 (qui est un logiciel in-
terprété de calcul symbolique ou formel), donc entrainé aux calculs sans pertes
lorsque les variables d’entrée sont des entiers, des rationnels, ou des polynômes à
coefficients entiers, rationnels, voire algébriques, MATLAB est un logiciel interprété
de calcul scientifique, où les calculs numériques impliquent en général des nombres
réels ou complexes déclarés en flottants (floating, soit -c’est le cas par défaut- en
TP1 : FAMILIARISATION AVEC MATLAB : SIGNAUX ET IMAGES DIGITALES, DFT 3
double précision double, ou bien en simple précision single) et sont donc tribu-
taires des arrondis inévitables liés à l’erreur machine (en accord toutefois avec la
norme du standard IEEE754). Testez les commandes suivantes :
>> eps
>> realmax
>> realmin
A l’aide de help (exemple : tapez help eps), analysez le sens de ces variables
modifiables ; comment interpréter le résultat de la commande
>> eps(x)
lorsque x est un nombre flottant (faites le test avec eps(pi), avec eps(1), eps(1/2),
eps(10(20) + pi), eps(107 +pi)). Conclusion ? En quoi eps(x) peut il être considéré
comme un indicateur de la précision au niveau du flottant x ? Le logiciel n’affiche
par défaut que quatre chiffres représentatifs pour les nombres réels flottants (en
double précision) ; pour en voir 16 (comme cela est théoriquement possible, aux
arrondis près, cf. la valeur de eps(x)), il faut préalablement taper l’instruction :
>> format long
Testez ceci sur des exemples (par exemple sur les variables non effacables i, j,
sqrt(-1)). Attention : j est le i des physiciens, ce n’est pas ici le exp(2*i*pi/3)
des mathématiciens !
Dans chaque cas où vous n’avez pas reçu un message d’erreur, affichez BB
ou C en faisant
>> BB
>> C
(sans point virgule cette fois) et tentez d’analyser ce qui se passe. Les
bases de calcul matriciel acquises en L1 et en L2 vous seront ici utiles.
(2) Les commandes rand (M,N) et randn (M,N) génèrent des matrices à M
lignes et N colonnes dont les entrées sont des réalisations de variables
aléatoires réelles (mutuellement indépendantes) suivant respectivement :
– en ce qui concerne la fonction rand, la loi uniforme sur [0, 1] (chaque
entrée est un nombre flottant aléatoire entre 0 et 1, tous les choix
étant équiprobables) ;
– en ce qui concerne la fonction randn, la loi normale (moyenne 0 et
écart type égal à 1) sur R.
Générez des matrices suivant les instructions :
>> A = rand(10,8);
>> A
>> Adouble = [A A];
>> Adouble
>> B = randn(12,9);
>> B
>> Bdouble = [B;B];
>> Bdouble
>> AA = A(1:2:10,1:2:8);
>> AA
>> BB = B(1:3:8,2:2:9);
>> ONES = ones (7,5);
>> ONES
>> ZEROS = zeros (10,7);
>> ZEROS
>> EYE = eye (10,5);
>> EYE
>> AAA = flipud (A);
>> AAA
>> BBB = fliplr (B);
>> BBB
Analysez l’opération qu’exécute chacune de ces instructions. Qu’exécute
en particulier l’instruction Mat1=Mat(1:pas1:M,1:pas2:N) si Mat est un
tableau numérique à M lignes et N colonnes, pas1 et pas2 étant des entiers
respectivement entre 1 et M et 1 et N ?
(3) Construisez une matrice Apad à 16 lignes et 16 colonnes en bordant la
matrice A générée à la question (1) de manière symétrique (des deux côtés)
par des blocs de zéros 3. Faites ensuite la même chose avec cette fois des
blocs de 1.
3. Pareille opération, fort utile dans la pratique, soit pour compléter une matrice rectangulaire
en une matrice carrée pour des calculs ultérieurs, soit pour s’accorder une marge de sécurité ,
telle un cadre , autour d’un tableau ou d’une image, s’appelle le zeropadding.
TP1 : FAMILIARISATION AVEC MATLAB : SIGNAUX ET IMAGES DIGITALES, DFT 5
(4) Créez une matrice à 19 lignes et 15 colonnes dont les lignes et les colonnes
sont celles de la matrice A générée à la question (1), mais séparées cette
fois entre elles par des lignes et des colonnes de zéros.
(5) Quelle est la méthode la plus judicieuse (au niveau du temps d’exécution)
pour déclarer sous MATLAB les matrices :
0 0 3 0 1 0 1 1 3 1 1 1
5 0 0 0 0 8 5 1 1 1 1 8
U = 0 2 0 0 1 0 V =
1 2 1 1 1 0 ?
0 0 1 0 8 0 1 1 1 1 8 1
1 0 0 0 0 5 1 1 1 1 1 5
Déclarez ces deux matrices de la manière que vous jugerez la plus efficace
possible (au niveau du temps d’exécution).
Sauvez votre session dans votre répertoire TPSignal (celui dans lequel vous êtes
précisément en train de travailler) avec l’instruction save TP1exo1. Vérifiez que
vous avez bien ainsi créé un fichier TP1exo1.mat dans votre répertoire TPSignal.
Faites ensuite l’instruction clear all pour rendre à nouveau vierge tout votre
Workspace. Vérifiez que c’est bien le cas. Que se passe t’il si vous donnez les ins-
tructions :
>> load TP1exo1
>> who
Faites clear all à nouveau.
Exercice 2 (signaux analogiques discrétisés et transformation de Fourier dis-
crète (DFT)). Les signaux audio constituent des exemples de signaux analogiques
(fonctions d’une variable continue, à savoir le temps) qui ont été discrétisés. Si le pas
temporel d’échantillonnage est de 1/Fs seconde, on dit que le signal analogique 4
est échantillonné à Fs Hertz. La partie entière de Fs constitue donc le nombre
d’échantillons collectés par seconde. Dans un signal audio au format .wav (par
exemple toto.wav) se trouvent encodés la fréquence d’échantillonnage Fs (en Hertz,
c’est-à-dire qu’il y a Fs échantillons par seconde), les échantillons successifs du
signal, ainsi que le nombre de bits B utilisé pour coder chaque échantillon. La
commande
>> [s,Fs,B] = wavread(‘toto.wav’);
génère, outre les valeurs de Fs et B, le signal audio (analogique) discrétisé (c’est-
à-dire les échantillons du signal audio codés initialement sur B bits) et converti au
format numérique, s(k) correspondant à un nombre réel flottant codé en double
précision et normalisé de manière à appartenir au segment [−1, 1] (16 décimales
dans la mantisse suivant la virgule dans le système décimal, vérifiez le en vous
mettant sous format long).
(1) Depuis le répertoire
http://www.math.u-bordeaux1.fr/~yger/MATLABSignal/Audio
téléchargez les 10 fichiers saxo-NOTE-wav, violon-NOTE-wav, où NOTE est
à prendre dans la liste : do, mi, la, ladiapason, si des notes de la
gamme. Convertissez ces 10 fichiers en 10 fichiers digitaux NOTE-saxo,
4. Par exemple la position de la membrane du haut-parleur s’il s’agit d’un signal audio.
6 TP1 : FAMILIARISATION AVEC MATLAB : SIGNAUX ET IMAGES DIGITALES, DFT
2iπ(j − 1)Fs
X
M
t
s : t 7−→ cj+M e N
j=−M +1
de manière à ce que
∀ k = 1, ..., N, S((k − 1)/Fs) = s((k − 1)/Fs).
Vérifiez que le système de N équations à N inconnues que doivent satisfaire
les nombres cj pour qu’il en soit ainsi s’exprime :
X
M
cj+M e2iπ(j−1)(k−1)/N = S((k − 1)/Fs), k = 1, ..., N
j=−M+1
>> CC = fft(s,Fs)/Fs;
Quel est le champ fréquenciel (en Hertz) exploré ? Prenez pour N des va-
leurs diverses (512, 1024, 2048, 4096) et choisissez au hasard une fenêtre
d’observation sN de longueur N pour l’un de vos signaux digitaux NOTE-saxo
ou Note-violon (plutôt au cœur du signal, aux alentours de 20000). Cal-
culez le module et l’argument (phase) des entrées de
>> dft_sN = fftshift(fft(sN,N))/N);
Pourquoi le champ fréquenciel exploré reste-t’il [−Fs/2, Fs/2], mais cette
fois discrétisé avec un pas de Fs/N et non plus de 1 ?
(3) En privilégiant la valeur de N=2048 et en choissant une fenêtre d’observa-
tion judicieusement, précisez (en Hertz) les valeurs des fréquences corres-
pondant aux divers harmoniques de la note enregistrée. Examinez le cas
de diverses notes (en particulier le ladiapason), sur les deux instruments
que sont le saxo et le violon.
(4) Après avoir étudié la fonction de la routine hamming, explicitez ce qu’ef-
fectue le code MATLAB suivant :
function [Omega,PSD] = welch1(s,Fs,fenetre,recouvrement);
NN = length(s);
N = length(fenetre);
Omega = - fix(Fs/2) : Fs/N : fix(Fs/2)-1;
PAS = N-recouvrement;
PSD=zeros(N,1);
q=0;
for j= N+1: PAS : NN,
y= fenetre.* s(j-N:j-1);
PSD=PSD+(abs(fft(y))).^2;
q=q+1;
end
PSD=PSD/(q*(norm(fenetre))^2);
PSD1=PSD(1:fix(N/2));
PSD2=PSD(fix(N/2)+1:N);
PSD=[PSD2;PSD1];
plot(Omega,PSD);
Ce processus d’inspiration statistique correspond à l’implémentation de
ce que l’on appelle la transformation de Welch. Cette transformation cor-
respond à un calcul statistique du carré du module du spectre (sur une
gamme de N fréquences) par simple moyennisation des résultats obtenus
sur des fenêtres d’exploration du signal (chacune de longueur N). Ce cal-
cul correspond à la version discrètisée du calcul de la densité spectrale de
puissance (psd) dont on reparlera.
(5) Générez des signaux se présentant comme des fonctions présentant des
discontinuités, par exemple :
8 TP1 : FAMILIARISATION AVEC MATLAB : SIGNAUX ET IMAGES DIGITALES, DFT
Constatez en faisant des tests que I est un tableau à trois entrées : les deux premières
entrées figurent la position du pixel (ligne, colonne) ; une fois ces entrées fixées (par
exemple i et j, les trois sorties I(i,j,1), I(i,j,2), I(i,j,3) sont des nombres
entiers codés entre 0 et 255 (8 bits car 256 = 28 ), correspondant au codage R,G,B
dans le cube des couleurs. Pour travailler du point de vue numérique sur cette
image, on la convertit d’abord au format Noir et Blanc (gray) :
>> II = rgb2gray(I);
Cette fois II est devenu un tableau à deux entrées, mais toujours au format uint8
(entrées entre 0 et 255, du noir [luminance minimale] au blanc [luminance maxi-
male]. Pour convertir ce tableau en un tableau de nombres flottants entre 0 et 1
(en double précision) entre 0 et 1 (et pouvoir ultérieurement traiter cette image du
point de vue nnumérique), il reste à effectuer la commande :
>> Inumer = double(II);
Testez les commandes :
>> image(I);
>> imagesc(I);
>> image(II);
>> imagesc(II);
>> image(Inumer);
>> imagesc(Inumer);
Téléchargez toujours depuis
http://www.math.u-bordeaux1.fr/∼yger/MATLABSignal/Images
cette fois le fichier (de données) venise. Une fois ceci fait, chargez ce fichier
dans votre Workspace en tapant sous la Command Window les instructions :
>> load venise
>> I=venise;
Grâce à la commande size, vérifiez que la variable ainsi déclarée I est un tableau à
256 lignes et 256 colonnes 5. Testez les valeurs de quelques entrées I(i,j) pour voir
qu’il s’agit en fait d’une matrice dont les entrées sont des nombres entiers positifs
(mais considérés ici comme des nombres flottants en double précision, comme la
réponse 1 à l’instruction isfloat(I) vous le prouvera, vérifiez le).
(1) Visualisez la matrice I de deux manières :
– en utilisant la visualisation en 3D suivant mesh(I) ;
– en utilisant la visualisation en 2D suivant image(I), imagesc(I)
(voire aussi imshow(I,[]) si vous êtes au CREMI, cf. plus loin).
Quelle est, pour ce type de matrice I, de ces deux visualisations, celle qui
est la plus adaptée ? Quel est l’effet de la commande
>> II = imrotate(I,angle);
(où angle désigne la valeur en flottant d’un angle entre 0 et 360 degrés
exprimé en flottant) ? 6. On utilisera l’une des routines du type image
mentionnées plus haut pour visualiser l’image II ainsi qu’une instruction
convenable (faire help axis) pour que l’affichage respecte la taille (ici
5. Notez que 256 est une puissance de 2, ce qui n’est, on le verra plus tard, nullement un
hasard en ce qui concerne les formats d’images jpeg.
6. Cette instruction imrotate bien utile n’est malheureusement pas présente dans la bi-
bliothèque de Scilab.
10 TP1 : FAMILIARISATION AVEC MATLAB : SIGNAUX ET IMAGES DIGITALES, DFT