V.

La convolution

Exemple 1.
Un système de mesure introduit inéluctablement une déformation de l'objet mesuré.
Cette déformation est appelée réponse du système.
Un système de mesure numérique est encore plus drastique où une portion de l'objet à
mesurer doit être représentée par un point de mesure statistiquement défini (pixel, voxel,
intervalle de temps etc...).
L'exemple ci-dessous représente des portions de l'objet par des pixels. Plus la taille
des pixels est grosse, plus l'image perd les détails. L'objet est convolué par la fonction de
réponse du système de mesure (point spread function en anglais).

20 20

40 40

60 60

80 80

100 100

120 120

140 140

160 160

180 180

200 200
50 100 150 200 250 50 100 150 200 250

800 Originale
Modifiée
20
700
40
600
60

80 500

100 400

120
300
140
200
160

180 100

200
0
50 100 150 200 250 0 50 100 150 200 250 300

Cours IMN 359. Page 68/138

 Image2 - Image1

20

40

60

80

100

120

140

160

180

200
50 100 150 200 250

Exemple 2. Faire la moyenne entre 3 nombres successifs du vecteur a = [0 4 1 3 6 5 2 0]

en faisant la moyenne de chaque valeur avec ses deux voisines. (Par exemple dans un
but de réduire les écarts entre les valeurs expérimentales a):
(0+0+4)/3 = 1.33; (0+4+1)/3=1.67; (4+1+3)/3=2.67.....
On trouve: b= [1.33 1.67 2.67 3.33 4.67 4.33 2.33 0.67].

Exemple 3. La moyenne entre les 3 nombres successifs du vecteur précédent peut se

faire avec l'opération de convolution:
N 1
L'opération de convolution: b = a h=h a. b(k ) ¦ a(n)h(k  n) et 0”N”10-1 avec N
n 0

et M sont les longueurs de a et de h respectivement.

Exemple d'aplication de cette formule:

a=[1 2 3], h=[1 3]. Calculer b=a h. Nous pouvons utiliser les indices n=0:N-1. n est l'indice
du vecteur a. Dans cet exemple, les indices de a sont 0, 1 et 2, et les indices de h sont 0
et 1. Les valeurs de a et de h ailleurs qu'à ces indices valent 0. Comme le vecteur a
possède 3 valeurs et h possède 2 valeurs, le produit de la convolution b aura 3+2-1=4
valeurs ce qui implique que b a pour indices k variant de 0 à 3:
n=0 n=1 n=2
k=0: b(0)= a(0).h(0-0) + a(1).h(0-1) + a(2).h(0-2) = 1.1 + 2.0 + 3.0 = 1
k=1: b(1)= a(0).h(1-0) + a(1).h(1-1) + a(2).h(1-2) = 1.3 + 2.1 + 3.0 = 5
k=2: b(2)= a(0).h(2-0) + a(1).h(2-1) + a(2).h(2-2) = 1.0 + 2.3 + 3.1 = 9
k=3: b(3)= a(0).h(3-0) + a(1).h(3-1) + a(2).h(3-2) = 1.0 + 2.0 + 3.3 = 9

Cours IMN 359. Page 69/138

L'opération peut-être faite en intervertissant les arguments de a et de h comme suit:
N 1
b( k ) ¦ a ( k  n) h( n) et 0”N”10-1.
n 0

n=0 n=1 n=2

k=0: b(0)= a(0-0).h(0) + a(0-1).h(1) + a(0-2).h(2) = 1.1 + 0.3 + 0.0 = 1
k=1: b(1)= a(1-0).h(0) + a(1-1).h(1) + a(1-2).h(2) = 2.1 + 1.3 + 0.0 = 5
k=2: b(2)= a(2-0).h(0) + a(2-1).h(1) + a(2-2).h(2) = 3.1 + 2.3 + 3.0 = 9
k=3: b(3)= a(3-0).h(0) + a(3-1).h(1) + a(3-2).h(2) = 0.1 + 3.3 + 2.0 = 9
Dans ce cas, il n'est pas nécessaire d'évaluer pour n=2 car h(2)=0, et le calcul se réduit à
2 colonnes correspondant aux valeurs de h.

Avec n=3 valeurs pour h, il faut m+n-1=10 opérations de multiplications-additions pour

couvrir les m=8 valeurs de a.
6
a
b
a = [0 4 1 3 6 5 2 0]; h=[1 1 1]; h = h / sum(h); 5
b = conv(a,h); figure;plot(1:length(a),a,'bo-
',(1:length(b))-1,b,'rs-','linewidth',2); 4
set(gca,'xcolor',[0 0 0],'ycolor',[0 0
0],'FontSize',16);
set(gcf,'Color',[1 1 1],'Position',[280 3
200 680 500],'InvertHardCopy','off');
legend('a','b'); print -dmeta convol_ab 2

0
0 1 2 3 4 5 6 7 8 9

Cours IMN 359. Page 70/138

Exemple 4.
Faire le produit de deux polynômes a et b:
a = x4 + 2x3 + 3x2 + x
b = 2x2 + x + 3
c= a * b = (x4 + 2x3 + 3x2 + x)( 2x2 + x + 3)
= 2x6 + 5x5 + 11x4 + 11x3 + 10x2 + 3x.

Faire la convolution de a et b:
a = [1 2 3 1]
b = [2 1 3]
c=a b = [2 5 11 11 10 3]
La convolution est une multiplication polynomiale.

Connaissant c et b, on peut retrouver a par une division polynomiale: a = c/b.

La déconvolution est une division polynomiale. [La division ne donne pas toujours un
quotient sans reste, si bien que le quotient pourrait être très approximatif.]

Exemple 5. L'opération de la convolution vue comme un produit matriciel. Cette approche

est utilisée lors des convolutions non-stationnaires.
a = [1 2 3 1]
b = [2 1 3]
c=a b = [2 5 11 11 10 3]

qui peut s'écrire sous la forme:

Cours IMN 359. Page 71/138

Exemple 6:

Figure. La courbe x(t) a été mesurée à t=W sur un intervalle 'W.

Supposons que les signaux x(t) sont des signaux instantanés delta mesurés pendant 'W.
Ceci équivaut à un signal G de Dirac pondéré par la surface du rectangle x(W)'W :
x(W)'WG(W).
Pour une position quelconque Wk, ceci devient: xk(t)=x(Wk)'WG(t-Wk).

Ce signal xk(t) est transformé par un système de fonction de réponse h pour devenir:
yk(t)=x(Wk)'Wh(t-Wk).

En considérant tous les points de x(t), on trouve:

n
y (t ) ¦ x(W k )h(t  W k )'W
k 0

f
et lorsque 'W tend vers 0: y (t ) ³ x(W )h(t  W )dW
f

Cours IMN 359. Page 72/138

Définition de la convolution
La convolution de deux fonctions f et g définies sur R est donnée par:
f f
f (t ) g (t ) ³ f (W ) g (t  W )dW ³ g (W ) f (t  W )dW
f f

Ici t est une variable quelconque et ne représente pas uniquement le temps.

La convolution existe et est définie si:

1- f(t) et g(t) sont deux fonctions de carré intégrable:
1/ 2 1/ 2
ªf 2 º ªf 2 º
« ³ f (t ) dt »  f et « ³ g (t ) dt » f
«¬ f »¼ «¬ f »¼

ou bien:
1/ 2 1/ 2
ªf 2 º ªf 2 º
f (t ) g (t ) d « ³ f (W ) dW » « ³ g (t  W ) dW » f
«¬ f »¼ «¬ f »¼

2- f(t) est intégrable et g(t) est bornée.

Exemple 1:
Calcul graphique de la convolution.

Cours IMN 359. Page 73/138

Cours IMN 359. Page 74/138
 t2
y (t ) ³ x(t  W )h(W )dW
t1

0
y (t ) ³ x ( 0  W ) 0 dW 0
f

0.5
y (t ) ³1* 2dW 1
0

0.5 1
y (t ) ³ 2 * 2 dW  ³ 1 * 2 dW 3
0 0.5

1 1.5
y (t ) ³ 2 * 2dW  ³1* 0dW 2
0.5 1

1,5 2
y (t ) ³ 2 * 0dW  ³1* 0dW 0
1 1,5

Par calcul numérique:

x = [0 1 2 0];
h = [1 1];
y = conv(x,h)
y = [0 1 3 2 0];

Cours IMN 359. Page 75/138

Exemple 2.
Calcul analytique de la convolution.
­°0 si t  0
f (t ) ® -t œ f (t ) e t u (t ) u(t) = fonction de Heaviside.
°̄e si t t 0

­t si 0 d t d T
g (t ) ® œ g (t ) t >u (t )  u (t  T )@
¯0 ailleurs
f f
f (t ) g (t ) ³ f (W ) g (t  W )dW ³e
W
u (W )(t  W )>u (W )  u (t  W )@dW
f f

­
°
°t  0 f(t) g(t) 0
°
> te @> @
t t t
° W
®0 d t d T f(t) g(t) ³ e (t  W )dW t ³ e W dW  ³We W dW t
 t  e  t (t  1)  1 et  t 1
° 0 0 0
° t
°t t T W
e  (t T ) (T  1)  e  t
° f(t) g(t) ³e (t  W )dW
¯ t T

Soit h = f g. h(t<0) = 0; h(0”W”T) = exp(-t) + t - 1;

h(t•T) = exp(-(t-T))(T-1) + exp(-t).

T=2;ts=0.1;t1=-T:ts:0;h1=0*t1; t2=0: ts:T;h2= exp(-t2) + t2 - 1; t3=T: ts:2*T;h3= exp(-(t3-

T))*(T-1) + exp(-t3); t=[t1 t2 t3]; h=[h1 h2 h3]; figure;plot(t,h,'bo-','linewidth',2);
set(gca,'xcolor',[0 0 0],'ycolor',[0 0 0],'FontSize',16);
set(gcf,'Color',[1 1 1],'Position',[280 200 680
500],'InvertHardCopy','off');
xlabel('t');ylabel('Amplitude'); print -dmeta convol_exp1
1.4

1.2

1
Amplitude

0.8

0.6

0.4

0.2

0
-2 -1 0 1 2 3 4
t

Figure. Graphique de h = f g calculée analytiquement.

Cours IMN 359. Page 76/138

Prenons le même exemple et retrouvons les résultats par la convolution de f(t) et g(t)
numériquement:
T=2; ts=0.1; t1=-T:ts:0; f1=0*t1; g1=0*t1;
t2=0:ts:T; f2=exp(-t2); g2=t2;
t3=T:ts:2*T; f3=exp(-t3);g3=0*t3;
t=[t1 t2 t3];f=[f1 f2 f3]; g=[g1 g2 g3];
figure;plot(t,f,'bo-',t,g,'rs-','linewidth',2); set(gca,'xcolor',[0 0 0],'ycolor',[0
0 0],'FontSize',16);
set(gcf,'Color',[1 1 1],'Position',[280 200 680
500],'InvertHardCopy','off');
xlabel('t');ylabel('Amplitude');legend('f(t)','g(t)'); print -dmeta convol_exp2;

h=conv(f,g); imax=find(h==max(h)); t1=(-imax+1)*ts+2:ts:2; t2=2+ts:0.1:(length(h)-imax-

1)*ts+2+ts; t=[t1 t2];figure;plot(t,h,'bo-','linewidth',2); set(gca,'xcolor',[0 0
0],'ycolor',[0 0 0],'FontSize',16);
set(gcf,'Color',[1 1 1],'Position',[280 200 680
500],'InvertHardCopy','off');
xlabel('t');ylabel('Amplitude');legend('f(t)*g(t)');axis([-4.2 8.1 0 13]); print -dmeta
convol_exp3;
2
f(t) 12 f(t)*g(t)
g(t)

1.5 10

8
Amplitude

Amplitude

1
6

4
0.5
2

0 0
-2 -1 0 1 2 3 4 -4 -2 0 2 4 6 8
t t

Figure. Graphiques de f(t) et g(t). Figure. Graphique de h = f g.

% t de la dernière figure est utilisé ici, et en normalisant h par rapport à ha:
T=2;ts=0.1;t1=min(t):ts:0;h1=0*t1; t2=ts:ts:T;h2= exp(-t2) + t2 - 1; t3=T+ts: ts:max(t);h3=
exp(-(t3-T))*(T-1) + exp(-t3); ta=[t1 t2 t3]; ha=[h1 h2 h3]; figure;plot(ta,ha,'bo-
',ta,h*tot1(ha)/tot1(h),'r*-','linewidth',2); set(gca,'xcolor',[0 0 0],'ycolor',[0 0
0],'FontSize',16);set(gcf,'Color',[1 1 1],'Position',[280 200 680
500],'InvertHardCopy','off');xlabel('t');ylabel('Amplitude'); legend('f(t)*g(t) analytique',
'f(t)*g(t) numérique');axis([-4.2 8.1 0 1.5]); print -dmeta convol_exp4;

Cours IMN 359. Page 77/138

 1.5
f(t)*g(t) analytique
f(t)*g(t) numérique

Amplitude

0.5

0
-4 -2 0 2 4 6 8
t

Figure. h(t) analytique et numérique.

Exemple 3:
Calcul analytique de la convolution.

1 / 2
y (t ) ³ x(t  W )0dW 0
f

Cours IMN 359. Page 78/138

 t 1
y (t ) ³ (1  t  W ) *1dW
1 / 2

4t 2  12t  9
8

t
y (t ) ³ (1  t  W ) *1dW 
1 / 2
1/ 2
 4t 2  3
³ (1  t  W ) *1dW 4
t

1/ 2
y (t ) ³ (1  t  W ) *1dW
t 1

4t 2  12t  9
8

1/ 2
y (t ) ³ (1  t  W ) * 0dW 0
t 1

0.7 y(t)

0.6

0.5
Amplitude

0.4

0.3

0.2

0.1

0
-3 -2 -1 0 1 2 3
t

% Tracer y(t) analytique comme définie ci-dessus:

ts=0.1; t1=-3:ts:-3/2; y1=0*t1; t2=-3/2+ts:ts:-1/2; y2=(9+12*t2+4*t2.^2)/8; t3=-
1/2+ts:ts:1/2; y3=(3-4*t3.^2)/4; t4=1/2+ts:ts:3/2; y4=(9-12*t4+4*t4.^2)/8; t5=3/2+ts:ts:3;
y5=0*t5; t=[t1 t2 t3 t4 t5]; y=[y1 y2 y3 y4 y5]; figure;plot(t,y,'bo-','linewidth',2);
set(gca,'xcolor',[0 0 0],'ycolor',[0 0 0],'FontSize',16);set(gcf,'Color',[1 1 1],'Position',[280 200
680 500],'InvertHardCopy','off');xlabel('t');ylabel('Amplitude'); legend('y(t)');axis([-4.2 8.1 0
1.5]); print -dmeta convol_triang1;

Cours IMN 359. Page 79/138

Les propriétés de la convolution:
La commutativité: g1(t) ¤ g2(t) = g2(t) ¤ g1(t).
f
g1 (t ) g 2 (t ) ³ g1 (W ) g 2 (t  W )dW
f

En faisant le changement de la variable de t-W = x:

f f

³ g1 (t  x) g 2 ( x)dx ³ g 2 ( x) g1 (t  x)dx g 2 (t ) g1 (t )
f f

L'associativité: [g1(t) ¤ g2(t)] ¤ g3(t) = g1(t) ¤ [g2(t) ¤ g3(t)].

f
f (t ) g1 (t ) g 2 (t ) ³ g1 ( y) g 2 (t  y)dy
f

f f ªf º
f (t ) g 3 (t ) ³ f ( x) g 3 (t  x)dx ³ «³ 1 2
« g ( y ) g ( x  y ) dy »g 3 (t  x)dx
f f ¬ f ¼»
En substituant z=x-y:
f ªf º
f (t ) g 3 (t ) ³ 1 «« ³ g 2 ( z ) g 3 (t  y  z )dz »»dy
g ( y )
f ¬ f ¼
f
Soit: h(t ) ³ g 2 ( z ) g 3 (t  z )dz
f

qui donne avec t = t-y:

f
h(t  y ) ³ g 2 ( z ) g 3 (t  y  z )dz
f

f
f (t ) g 3 (t ) >g1 (t ) g 2 (t )@ g 3 (t ) ³ g1 ( y)h(t  y)dy g1 (t ) h(t )
f
g1 (t ) >g 2 (t ) g 3 (t )@
La distributivité: g1(t) ¤ [g2(t) + g3(t)] = g1(t) ¤ g2(t) + g1(t) ¤ g3(t).

La convolution et les transformées de Fourier:

Z = 2Sf.
TF[g1(t) ¤ g2(t)] = G1(Z)G2(Z) et TFI[G1(Z)G2(Z)] = g1(t) ¤ g2(t).
fªf º  iwt f ªf º
TF[g1(t) ¤ g2(t)] ³ «« ³ g1(x)g2(t - x)dx»»e dt ³ g1( x ) « ³ g2(t - x)e  iwt dt »dx
- f¬- f ¼ -f «¬- f »¼

Cours IMN 359. Page 80/138

à l'aide de la propriété du décalage dans le temps (Time shifting):
f
 iwt
³ g2(t - x)e dt G 2( w)e  iwx
-f

qui est substituée dans l'équation précédente:

f ªf º
 iwx
TF[g1(t) ¤ g2(t)] ³ g1( x)G 2(w)e dx « ³ g1( x)e  iwx dx »G 2( w)
-f «¬- f »¼
ªf º
« ³ g1(t )e  iwt dt »G 2( w) G1( w)G 2( w)
«¬- f »¼

La convolution et les transformées de Fourier inverses:

TFI[G1(w) ¤ G2(w)] = 2Sg1(t)g2(t)
et TF[g1(t)g2(t)] = [G1(w) ¤ G2(w)]/2S.

Le théorème de Parseval:
Soit TF[g(t)] = G(Z) et g*(t) est le conjugué de g(t) ==>
* *

³ >g (t )e @ > @
f f ªf º
TF >g * (t )@ ³ g * (t )e
 iZt
dt iZt
dt « ³ g (t )e  i ( Z )t dt » G * (Z )
f f «¬ f »¼

De la propriété de la convolution dans l'espace des fréquences qui s'écrit:

f
1
TF[g1(t)g2(t)] [G1(Z ) ¤ G2(Z )]/2S
2S ³ G1( y)G 2(Z  y)dy
f

qui s'écrit à son tour:

f f
- iZt 1
³ [g1(t)g2(t)]e dt 2S ³ G1( y)G 2(Z  y)dy
f f

en faisant Z = 0,
f f
1
³ g1(t)g2(t)dt 2S ³ G1( y)G 2( y)dy
f f

puis en changeant y en Z dans le membre de droite:

f f
1
³ g1(t)g2(t)dt 2S ³ G1(Z )G 2(Z )dZ
f f

En prenant le conjugué g2*(t) de g2(t), nous obtenons selon ce qui a été démontré ci-haut
TF[g2*(t)] = G2*(-Z):

Cours IMN 359. Page 81/138

f f
1
³ g1(t)g2 * (t)dt 2S ³ G1(Z )G 2 *[(Z )]dZ
f f

qui devient:
f f
1
³ g1(t)g2 * (t)dt 2S ³ G1(Z )G 2 * (Z )dZ
f f

Qui est l'identité de Parseval.

Si on considère g(t)=g1(t) et g*(t)=g2*(t), l'équation précédente devient:

f f
1
³ g(t)g * (t)dt 2S ³ G(Z )G * (Z )dZ
f f

qui peut s'écrire:

f f
1
³ g(t) dt ³ G(Z )
2 2
dZ
f
2S f

Qui est l'identité de Plancherel.

Dans le monde discrèt, la normalisation par ( 2S ) est remplacé par la taille du signal de
période N; c’est-à-dire
N 1 N 1
1
¦ g[n] ¦ | G[k ] |
2 2

n 0 N k 0

Cours IMN 359. Page 82/138

 1- La corrélation
Définition: La corrélation croisée entre g1(t) et g2(t) permet de mesurer le degré de
similitude entre les deux fonctions. Si la corrélation croisée vaut 0, les fonctions sont dites
non-corrélées.
La corrélation croisée est donnée par:
f f
r12( s ) g1 $ g 2 ³ g1(t ) g 2(t  s)dt et r 21( s ) g 2 $ g1 ³ g 2(t ) g1(t  s)dt
f f

L'autocorrélation de la fonction g1(t) est donnée par:

f
r11( s ) ³ g1(t ) g1(t  s)dt
f

L'autocorrélation peut être normalisée comme suit:

f

³ g1(t ) g1(t  s)dt

J ( s) f et dans ce cas J(0) = 1.
f

³ >g1(t )@ dt
2

f

Les définitions ne changent pas en changeant t en t-s dans les corrélations croisées et les
autocorrélations.
f f
r12( s ) ³ g1(t ) g 2(t  s)dt ³ g1(t  s) g 2(t )dt
f f

f f
r 21( s ) ³ g 2(t ) g1(t  s)dt ³ g 2(t  s) g1(t )dt
f f

f f
r11( s ) ³ g1(t ) g1(t  s)dt ³ g1(t  s) g1(t )dt
f f

Propriétés des corrélations:

r12(s) = r21(-s)
r11(s) = r11(-s)

Fonctions discrètes:
De façon similaire à la convolution, on utilise un filtre de corrélation h, de préférence avec
un nombre impaire d'éléments, soit 2N+1. La corrélation s'écrit:
N
(h $ g )(n) ¦ h(i) g (n  i)
i N

Cours IMN 359. Page 83/138

Exemple 1.
g1=[5 4 2 3 7 4 6 5 3 6].
On désire faire la moyenne de 3 éléments consécutifs de g1 et obtenir un nouveau vecteur
g2.
Le noyau serait alors h=[1/3 1/3 1/3], et l'opération de corrélation est similaire à celle de
la convolution, mais sans inversion du noyau. Aussi on doit centrer le noyau sur chaque
élément de g1:

g2(1)=h(1)*g1(0)+ h(2)*g1(1)+ h(3)*g1(2) = [0*1/3 + 5*1/3 + 4*1/3]

..
..
g2(4)= h(1)*g1(3)+ h(2)*g1(4)+ h(3)*g1(5) = [2*1/3 + 3*1/3 + 7*1/3]
..
..

Puisque la moyenne sur le 1er élément de g1 se fait seulement sur 2 éléments, on doit
soit ajouter un 0 à gauche du 1er élément et à droite du dernier élément de g1 (zero
padding); ou bien répéter une fois le 1er élément de g1 à gauche et le dernier élément de
g1 à droite. Ces arrangements dépendent de l'application sous étude.

Exemple 2.
Calcul des dérivées dans un vecteur: h=[-1/2 0 1/2] qui permet de calculer [g1(i+1) - g1(i-
1)]/2.

Exemple 3.
Retrouver un signal dans un vecteur (ou dans une image). Soit un modèle h du signal
recherché. On le fait parcourir sur le vecteur g1 pour localiser la partie de g1 qui
ressemble (qui est corrélée) à h. Pour quantifier la ressemblance, on calcule la somme
des carrés des différences entre h et la partie correspondante de g1:
2

¦ >h(i) 2  g1(n  i) 2  2h(i) g1(n  i)@

N N
¦ >h(i)  g1(n  i)@
i N i N
N N N
¦ h(i) 2  ¦ g1(n  i) 2  2 ¦ h(i) g1(n  i)
i N i N i N

Cours IMN 359. Page 84/138

Cela revient à calculer la somme du carré du noyau, la somme du carré de la partie de g1
correspondante et 2 fois la corrélation. Plus ce nombre calculé est faible, plus h est corrélé
avec cette partie de g1.
g1=[5 4 2 3 7 4 6 5 3 6]; h=[3 7 4];
N
¦ h(i) 2 = 32+72+42 = 74
i N

N
¦ g1(n  i) 2 = [..............................]
i N

N
2 ¦ h(i) g1(n  i) = [......................................]
i N

Quelle est la partie de g1 la plus corrélée à h?

Cours IMN 359. Page 85/138