Méthodes quantitatives en neurosciences

Cours NSC-2006, année 2015
Laboratoire d'introduction au traitement des signaux - Réponses
*Pierre Bellec, Yassine Ben Haj Ali*

Objectif:

Ce laboratoire a pour but de vous familiariser avec la notion de convolution et de système linéaire invariant dans le temps. Nous allons implémenter des convolutions avec différents noyaux, et tester la linéarité du système. Nous allons reproduire la plupart des résultats présentés dans le cours avec matlab pour des signaux temporels, et également appliquer les mêmes techniques à des images. Pour réaliser l'exercice 2 de ce laboratoire, il est nécessaire de récupérer l'archive suivante, et copier l'image lena_gray_512.tif dans votre répertoire de travail.

Ne pas tenir compte et ne pas executer cette partie du code:



In [3]:

    
%matplotlib inline
from pymatbridge import Octave
octave = Octave()
octave.start()
%load_ext pymatbridge









    



Starting Octave on ZMQ socket ipc:///tmp/pymatbridge
Send 'exit' command to kill the server
.Octave started and connected!
The pymatbridge extension is already loaded. To reload it, use:
  %reload_ext pymatbridge

Section 1 : Convolution de signaux temporels

Commençons par générer un signal de type `bloc`, auquel on rajoute un bruit gaussien:



In [4]:

    
%%matlab
time_samples = 0.1*(0:1000);
bloc = zeros(size(time_samples));
bloc(300:500) = 1;
bloc = bloc + sqrt(0.01)*randn(size(time_samples));

Quel est la taille du signal?

Le signal bloc est un vecteur ligne de longueur 1001:



In [ ]:

    
%%matlab
size(bloc)

Quelle est la taille du bloc?

A l’aide de la ligne bloc(300:500) = 1 on voit que le bloc contient 500-300+1=201 valeurs non-nulles, avant addition du bruit.

Quelle est la fréquence d’échantillonnage, en supposant que les temps d’acquisition des échantillons sont représentés par time_samples en secondes?

A l’aide de la ligne time_samples = 0.1*(0:1000); on voit que les échantillons temporels sont séparés de 0.1 s, soit une fréquence d’échantillonnage de 1 Hz.

Quelle est la durée d’acquisition?

L’enregistrement dure 100 secondes.

Générez un graphe représentant bloc avec l’axe temporel adéquat (voir la fonction plot).

Pour générer un graphe avec l’axe temporel adéquat, on utilise plot avec deux arguments, ainsi que les commandes xlabel et ylabel pour labéliser les axes.



In [4]:

    
%%matlab
figure
plot(time_samples,bloc)
xlabel('time')
ylabel('a.u.')

On va maintenant construire un noyau, et appliquer la convolution sur le signal `bloc`



In [5]:

    
%%matlab
noyau = ones(1,3);
noyau = noyau/sum(abs(noyau));
bloc_conv = conv(bloc,noyau,'same');

Faire un graphe avec le noyau (note: stem est plus adapté que plot pour les signaux courts, et axis peut être utilisé pour ajuster les limites des axes).



In [6]:

    
%%matlab
figure
stem(noyau)
axis([0 4 0 0.5])

Quelle opération est associée avec la convolution basée sur ce noyau?

Ce noyau implémente une moyenne mobile.

Superposez les signaux bloc et bloc_conv sur un même graphe. Qu'observez-vous?



In [7]:

    
%%matlab
figure
plot(time_samples,bloc,'b')
hold on
plot(time_samples,bloc_conv,'r')

On observe une diminution de la variance du bruit, liée au moyennage.

Répétez la question précédente en utilisant un noyau de longueur 7, puis 100, puis 300.

Qu'observez-vous?



In [8]:

    
%%matlab
for nn = [7 100 300]
    figure
    noyau = ones(1,nn); noyau = noyau/sum(abs(noyau)); bloc_conv = conv(bloc,noyau,'same');
    hold on, plot(time_samples,bloc,'b'), plot(time_samples,bloc_conv,'r')
end

On observe que pour des filtres de taille croissante, le bruit devient très faible, mais que les transitions du bloc deviennent de moins en moins abruptes. Pour un noyau dont la taille dépasse celle du bloc (cas 300), la valeur maximale du bloc après convolution devient plus petite que la valeur moyenne dans le bloc avant convolution.

On va maintenant s'intéresser à un nouveau noyau, que l'on va appliquer sur le même signal `bloc`.



In [9]:

    
%%matlab
noyau = [-ones(1,10) 0 ones(1,10)];
noyau = noyau/sum(abs(noyau));
bloc_conv = conv(bloc,noyau,'same');

Faire un graphe avec le noyau.



In [10]:

    
%%matlab
figure
stem(noyau)
axis([0 22 -0.06 0.06 ])

Quelle opération est associée avec la convolution basée sur ce noyau?

Ce noyau est sensible aux variations de signal entre échantillons temporels voisins.

Superposez les signaux bloc et bloc_conv sur un même graphe. Qu'observez-vous?



In [11]:

    
%%matlab
figure
plot(time_samples,bloc,'b')
hold on
plot(time_samples,bloc_conv,'r')

Le signal convolué comporte deux pics importants, qui coincident avec les bords du bloc. une transition de la ligne de base vers le bloc induit une réponse négative, tandis qu'une transition du bloc vers la ligne de base génère une réponse positive.

On va maintenant construire un signal bloc plus court, et sans bruit, pour illustrer le concept de réponse à une impulsion finie.



In [7]:

    
%%matlab
time_samples = 1*(0:100);
bloc = zeros(size(time_samples));
bloc(30:38) = 1;

On va également construire un noyau qui ressemble à une réponse hémodynamique, simplifiée:



In [8]:

    
%%matlab
temps_pic = 3;
noyau = [linspace(0,1,temps_pic+1) linspace(1,-0.3,temps_pic/2) linspace(-0.3,0,temps_pic/2)];
noyau = [zeros(1,length(noyau)-1) noyau];

Faire un graphe avec le noyau.



In [16]:

    
%%matlab
figure
stem(noyau)

Quelle opération est associée avec la convolution avec le noyau?

Ce noyau est un modèle grossier d'une fonction de réponse hémodynamique cérébrale.

On va maintenant appliquer la convolution avec le noyau au signal `bloc`:



In [17]:

    
%%matlab
bloc_conv = conv(bloc,noyau,'same');

Superposez bloc et bloc_conv sur un même graphe.



In [18]:

    
%%matlab
figure
plot(time_samples,bloc,'b-o')
hold on
plot(time_samples,bloc_conv,'r-o')

Notez l'option d'affichage pour ajouter un marqueur à chaque point de mesure.

Combien de temps met la réponse à atteindre son pic?

On voit que trois points séparent le début du bloc du pic de la réponse. En regardant la définition de time samples, on voit que la fréquence d'échantillonnage est de 1 Hz. La réponse et donc 3 secondes à atteindre son pic.

Construire un signal similaire à bloc, mais avec une amplitude du bloc deux fois plus grande.



In [19]:

    
%%matlab
bloc2 = 2*bloc;

Comparez la réponse au gros bloc, avec le double de la réponse au petit bloc, en faisant un graphe de la différence des deux signaux.

On commence par calculer la réponse au bloc simple, que l'on affiche avec la réponse au bloc double.



In [22]:

    
%%matlab
bloc_conv2 = conv(bloc2,noyau,'same');
plot(time_samples,bloc_conv,'b-o')
hold on
plot(time_samples,bloc_conv2,'r-o')

En examinant les valeurs des pics, il semble bien que la réponse rouge (pic à 4) est le double de la réponse bleue (pic à 2). Vérifions cela en affichant la différence des deux signaux:



In [21]:

    
%%matlab
bloc_conv2 = conv(bloc2,noyau,'same');
plot(time_samples,2*bloc_conv-bloc_conv2,'b-o')

Est ce que la propriété de mise à l'échelle est respectée?

On voit que la différence entre `bloc` et `bloc2` est nulle, les deux signaux sont donc identiques. La propriété de mise à l'échelle semble être ici respectée.

Construisez un nouveau stimulus `impulsion`, contenant une seul impulsion unitaire à l'échantillon numéro 50. Calculez la réponse du système à ce nouveau stimulus, appelée `impulsion_conv`, et représentez-la sur un graphe.



In [19]:

    
%%matlab
impulsion = zeros(size(bloc));
impulsion(50) = 1;
impulsion_conv = conv(impulsion,noyau,'same');
figure
plot(time_samples,impulsion,'b-o')
hold on
plot(time_samples,impulsion_conv,'r-o')

Calculer la réponse à la somme de `bloc` et `impulsion`, et représentez la sur un graphe.



In [20]:

    
%%matlab
somme = bloc+impulsion;
reponse_somme = conv(somme,noyau,'same');
figure
plot(time_samples,somme,'b-o')
hold on
plot(time_samples,reponse_somme,'r-o')

Comparez la réponse à `bloc+impulsion` avec la somme de `bloc_conv` et `impulsion_conv`, en faisant un graphe de la différence des deux signaux. La propriété d'additivité est-elle respectée?



In [21]:

    
%%matlab
somme_reponse = impulsion_conv + bloc_conv;
figure
plot(time_samples,reponse_somme - somme_reponse,'b-o')

La réponse à la somme des signaux est identique à la somme des réponses. La propriété d'additivité est ici respectée.

Répétez la question 8 avec une impulsion localisé à l'échantillon 70. Est ce que la propriété d'invariance dans le temps est respectée?



In [22]:

    
%%matlab
impulsion2 = zeros(size(bloc));
impulsion2(70) = 1;
impulsion2_conv = conv(impulsion2,noyau,'same');
figure
plot(time_samples,impulsion_conv,'r-o')
hold on
plot(time_samples,impulsion2_conv,'b-o')

On observe que la réponse à `impulsion2` est identique à `impulsion`, décalée dans le temps. La propriété d'invariance temporelle semble respectée.

Section 2 : Convolution d'images

En premier lieu, il faut charger l'image lena_gray_512.tif à l'aide de la fonction imread():



In [23]:

    
%%matlab
img = imread ('lena_gray_512.tif');
img = double(img);









    





warning: your version of GraphicsMagick limits images to 8 bits per pixel

La fonction imagesc() nous permet de visualiser une image binaire comme suit:



In [24]:

    
%%matlab
figure
imagesc(img,[0 255]) % le 2ieme argument ajuste le minimum/maximum 
colorbar % Ajoutez une echelle de couleur
colormap gray % Utilisez une echelle en niveau de gris

Quelle est la taille de cette image? son type? Combien y-a-t-il de niveaux de gris dans cette image?



In [25]:

    
%%matlab
whos img









    





Variables in the current scope:

   Attr Name        Size                     Bytes  Class
   ==== ====        ====                     =====  ===== 
        img       512x512                  2097152  double

Total is 262144 elements using 2097152 bytes

l'image et de dimention 250x250 et de type binaire, elle a 250 niveau de gris

Nous allons maintenant créer un noyau pour la convolution, qui est une autre petite image définissant les poids d'une moyenne pondérée dans chaque petit carré (ici, 5x5) de l'image:



In [26]:

    
%%matlab
noyau1 = ones(5,5);
noyau1 = noyau1/sum(sum(abs(noyau1)));

Visualisez le noyau avec la commande imagesc.



In [27]:

    
%%matlab
imagesc(noyau1)
colorbar
colormap gray

Appliquons la convolution entre img et le noyau:



In [28]:

    
%%matlab
img1 = conv2(img,noyau1,'same');

Visualisez les deux images, avant et après la convolution, en prenant soin d'utiliser les mêmes valeurs min et max dans les deux cas. Quelle opération implémente ce noyau? Que remarquez vous sur l'image convoluée?



In [29]:

    
%%matlab
figure
imagesc(img,[0 255]); colormap gray; colorbar
figure
imagesc(img1,[0 255]); colormap gray; colorbar

Ce noyau implémente un moyenne mobile, soit une opération de lissage.



In [30]:

    
%%matlab
noyau2 = ones(20,20);
noyau2 = noyau2/sum(sum(abs(noyau2)));
img2 = conv2(img,noyau2,'same');
figure
imagesc(img,[0 255]); colormap gray; colorbar
figure
imagesc(img2,[0 255]); colormap gray; colorbar

Le lissage est plus agressif, il y a une dégragation du bord des objets, ainsi que du contraste dans l'image. On voit aussi des artefacts (bandes sombres) sur les bords de l'image.

Construisons maintenant un nouveau noyau:



In [31]:

    
%%matlab
noyau3 = [ -ones(5,2) zeros(5,1) ones(5,2) ]
imagesc (noyau3); colorbar; colormap gray









    





noyau3 =

  -1  -1   0   1   1
  -1  -1   0   1   1
  -1  -1   0   1   1
  -1  -1   0   1   1
  -1  -1   0   1   1

Appliquez la convolution sur l'image img avec le nouveau noyau noyau3 et enregistrez la nouvelle image convoluée dans la variable img3. Visualisez les images avant/après convolution, en utilisant l'échelle de couleur gray. Que remarquez vous?



In [32]:

    
%%matlab
noyau3 = noyau3/sum(abs(noyau3(:)));
img3 = conv2(img,noyau3,'same');
figure
imagesc(img,[0 250]); colormap gray; colorbar
figure
imagesc(img3); colormap gray; colorbar

Ce noyau permet une detection des bords verticaux dans l'image

Répétez les deux questions précédentes avec ces autres noyaux:

1- noyau4:



In [33]:

    
%%matlab
noyau4 = noyau3'
figure
imagesc (noyau4); colorbar; colormap gray









    





noyau4 =

  -0.050000  -0.050000  -0.050000  -0.050000  -0.050000
  -0.050000  -0.050000  -0.050000  -0.050000  -0.050000
   0.000000   0.000000   0.000000   0.000000   0.000000
   0.050000   0.050000   0.050000   0.050000   0.050000
   0.050000   0.050000   0.050000   0.050000   0.050000



In [34]:

    
%%matlab
noyau4 = noyau4/sum(abs(noyau4(:)));
img4 = conv2(img,noyau4,'same');
figure
imagesc(img,[0 250]); colormap gray; colorbar
figure
imagesc(img4); colormap gray; colorbar

C'est un noyau qui detecte le bords horizontaux de l'image

2- noyau5:



In [35]:

    
%%matlab
noyau5 = tril(-ones(size(noyau4))) + triu(ones(size(noyau4)))
figure
imagesc (noyau5)
colorbar
colormap gray









    





noyau5 =

   0   1   1   1   1
  -1   0   1   1   1
  -1  -1   0   1   1
  -1  -1  -1   0   1
  -1  -1  -1  -1   0



In [36]:

    
%%matlab
noyau5 = noyau5/sum(abs(noyau5(:)));
img5 = conv2(img,noyau5,'same');
figure
imagesc(img,[0 250]); colormap gray; colorbar
figure
imagesc(img5); colormap gray; colorbar

Ce noyau détecte les bords correspondant à un angle de -45 degrés

Soit img4 (respectivement img5) le résultat de la convolution avec noyau4 (respectivement noyau5). Visualisez et interprétez l'image abs(img3)+abs(img4)+abs(img5).



In [37]:

    
%%matlab
img_sum = abs(img3)+abs(img4)+abs(img5);
figure
imagesc(img_sum); colormap gray; colorbar

Cette image détecte l'ensemble des bords, sans distinction de direction.