coursINTRODUCTION AU FILTRAGE ADAPTATIF ET A L'EGALISATION - Fa - T [PDF]

Zitiervorschau

INTRODUCTION AU FILTRAGE ADAPTATIF ET A L’EGALISATION ENSICAEN ´ Sp´ecialit´e : Electronique et Physique Appliqu´ee Majeure : Signal, Automatique, T´el´ecommunications

[email protected]

1

Table des mati`eres

1 Quelques rappels et objet du cours 1.1 Notations et rappels . . . . . . . . . 1.1.1 La Transform´ee de Fourier . . 1.1.2 L’impulsion . . . . . . . . . . 1.1.3 Le filtrage num´erique lin´eaire 1.1.4 Analyse de corr´elation . . . . 1.2 Objets du cours . . . . . . . . . . . .

I

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

Techniques de filtrage adaptatif

2 Filtrage adaptatif au sens de Wiener 2.1 Principe g´en´eraux . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.1 Contexte de signaux al´eatoires et d´eterministes . . . . . . 2.1.2 Formalisme . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Solution globale . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.1 Recherche des param`etres du filtre . . . . . . . . . . . . . 2.2.2 V´erification du principe d’orthogonalit´e . . . . . . . . . . 2.2.3 Mise en oeuvre . . . . . . . . . . . . . . . . . . . . . . . . 2.2.4 Exemple . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Solution r´ecursive . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.1 Premi`ere approche . . . . . . . . . . . . . . . . . . . . . . 2.3.2 Algorithme du gradient d´eterministe . . . . . . . . . . . . 2.3.3 Propri´et´e de convergence . . . . . . . . . . . . . . . . . . 2.3.4 Exemple . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.5 Algorithme du gradient stochastique . . . . . . . . . . . . 2.3.6 Mise en oeuvre . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Exemples d’application . . . . . . . . . . . . . . . . . . . . . . . . 2.4.1 Application de la solution r´ecursive : annulation d’echo acoustique . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.2 Application de la solution globale : d´ebruitage de la parole

2

5 5 5 6 7 8 11

13 16 16 17 17 18 18 20 20 22 23 23 24 25 27 27 29 30 30 33

3 Filtrage adaptatif au sens des moindres carr´ es 36 3.1 Principe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.2 Exemples d’application . . . . . . . . . . . . . . . . . . . . . . . . 38 3.2.1 Exemple num´erique . . . . . . . . . . . . . . . . . . . . . 38 3.2.2 Identification d’un syst`eme de transmission flexible . . . . 39 3.2.3 Retour sur l’annulation d’echo acoustique . . . . . . . . . 41 3.3 Les moindres carr´es en milieu bruit´e . . . . . . . . . . . . . . . . 41 3.3.1 Le biais . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.3.2 Solutions dans un contexte d’identification . . . . . . . . . 44 3.3.3 Retour sur l’identification d’un syst`eme de transmission flexible . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 4 Filtrage De Kalman 4.1 Principe . . . . . . . . . . . . . . . . . 4.2 Filtre et pr´edicteur de Kalman . . . . 4.2.1 Filtre de Kalman . . . . . . . . 4.2.2 Pr´edicteur de Kalman . . . . . 4.2.3 Mise en oeuvre . . . . . . . . . 4.3 Solution stationnaire . . . . . . . . . . 4.4 Exemples d’application . . . . . . . . . 4.4.1 Filtrage d’un signal d´ecroissant 4.4.2 Position d’un mobile . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

51 51 52 52 54 56 56 59 59 61

II Egalisation dans une chaˆıne de transmission de l’information 64 5 Le probl` eme d’´ egalisation 65 5.1 Mod`ele ´equivalent discret du canal . . . . . . . . . . . . . . . . . 65 5.2 Bruit et interf´erence inter symbole . . . . . . . . . . . . . . . . . 68 6 Egaliseur avec connaissance du canal ou s´ equence d’apprentissage 6.1 Egaliseur lin´eaire . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1.1 Egaliseur zero forcing . . . . . . . . . . . . . . . . . . . . 6.1.2 Egaliseur de Wiener (ou `a erreur quadratique moyenne minimale) . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1.3 Egaliseur par retour de d´ecision . . . . . . . . . . . . . . . 6.2 Egaliseur adaptatif . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2.1 Egaliseur de Wiener adaptatif . . . . . . . . . . . . . . . . 6.2.2 Egaliseur par retour de d´ecision adaptatif . . . . . . . . .

72 72 72 79 85 87 88 89

7 Egaliseurs aveugles (ou autodidactes) 92 7.1 Egaliseurs de type Bussgang . . . . . . . . . . . . . . . . . . . . . 92 7.2 Egaliseurs aveugles par identification . . . . . . . . . . . . . . . . 98 A R´ esultats divers 100 A.1 Lemme d’inversion matricielle . . . . . . . . . . . . . . . . . . . . 100 A.2 Forme quadratique incompl`ete . . . . . . . . . . . . . . . . . . . 100

3

Remerciements `a Guy Binet, Miloud Frikel, Olivier Gehan et Eric Pigeon pour les remarques et critiques.

4

CHAPITRE

1

Quelques rappels et objet du cours

1.1

Notations et rappels

Dans ce cours, seuls sont trait´es les signaux discrets c’est `a dire ´echantillonn´es. Un signal temporel x(t) sera donc d´ecrit par la valeur xk de ses ´echantillons aux instants kTe . Te d´esigne la p´eriode d’´echantillonnage. On supposera dans la suite que la p´eriode d’´echantillonnage a ´et´e judicieusement choisie de mani`ere `a satisfaire la condition du th´eor`eme de Shannon.

1.1.1

La Transform´ ee de Fourier

Rappelons que la transform´ee d’un signal ´echantillonn´e x(t) est d´efinie par : T F [x(t)] = X(f ) =

+∞ X

xk e−j2πf kTe

k=−∞

P+∞ Cette transform´ee de Fourier existe si k=−∞ |xk |. Dans la suite, lorsqu’on sera sera amen´e `a calculer une transform´ee de Fourier, la condition pr´ec´edente sera suppos´ee v´erifi´ee par d´efaut. De mˆeme, certains calculs n´ecessiteront la permutation d’int´egrales, de sommes et de limites. Par d´efaut, les conditions n´ecessaires `a la validit´e de ces permutations seront suppos´ees satisfaites (voir cours de Math´ematiques). La transform´ee de Fourier d’un signal ´echantillonn´e ´etant p´eriodique de p´eriode 1 ee de Fourier inverse s’exprime par : Te , la transform´ Z xk = Te X(f )ej2πf kTe df 1 Te

Cette propri´et´e de r´eversibilit´e est de premi`ere importance puisqu’elle signifie que x(t) et X(f ) contiennent les mˆemes informations, d´ecrivent les mˆemes comportements mais simplement de mani`ere diff´erente. Ce qui diff`ere c’est simplement la nature des informations directement accessibles. La repr´esentation 5

temporelle d´ecrit les variations au cours du temps du signal, la repr´esentation fr´equentielle d´ecrit la vitesse `a laquelle se font ces variations. Ainsi, suivant l’information recherch´ee, suivant le ph´enom`ene observ´e il peut ˆetre plus int´eressant de traiter l’une ou l’autre des repr´esentations. La repr´esentation graphique de la transform´ee de Fourier se fait la plupart du temps par l’interm´ediaire du module |X(f )| ou encore du module au carr´e. Ce dernier terme |X(f )|2 est appel´e densit´e spectrale et est not´e : Sxx (f ) = X(f )X ∗ (f ) = |X(f )|2

De la mˆeme mani`ere, la densit´e spectrale interspectre est d´efinie par Syx (f ) = Y (f )X ∗ (f )

1.1.2

L’impulsion

Tout comme l’analyse en temps continu, l’analyse en temps ´echantillonn´e utilise couramment la s´equence impulsion. Cette s´equence est not´ee δ(t) et se caract´erise par :  1 pour k = 0 δk = 0 pour k 6= 0 δ(t) 6 1 6

- t

Figure 1.1 – Repr´esentation graphique de l’impulsion δ(t) Cette impulsion permet une formulation du signal x(t) comme suit : x(t) =

+∞ X

k=−∞

xk δ(t − kTe )

La transform´ee de Fourier de l’impulsion est T F [δ(t)] =

+∞ X

δk e−j2πf kTe = 1

k=−∞

Un train d’impulsion de p´eriode T est d´efini par : δT (t) =

+∞ X

i=−∞

δ(t − iT )

La transform´ee de Fourier d’un train d’impulsion est : 1 T F [δT (t)] = δ1/T (f ) T

6

1.1.3

Le filtrage num´ erique lin´ eaire

Rappelons ici la notion de filtrage num´erique lin´eaire. On consid`ere pour cela un filtre H(q) d´efini par : H(q) =

+∞ X

hi q −i

i=0

les termes hi d´efinissent les coefficients de la r´eponse impulsionnelle du filtre et q −1 est l’op´erateur retard : q −1 x(t) = x(t − Te ). Le filtrage de x(t) par H(q) s’exprime comme suit : y(t) = H(q)x(t) y(t) =

∞ X

hi q −i x(t)

i=0

y(t) =

∞ X i=0

hi x(t − iTe )

Ainsi la valeur yk de y(t) `a l’instant t = kTe est : yk =

∞ X

hi xk−i

i=0

Le filtrage lin´eaire s’exprime de mani`ere ´equivalente par l’interm´ediaire du produit de convolution ∗. Si on d´efinit un signal h(t) `a partir des coefficients hi : h(t) =

+∞ X

k=0

hk δ(t − kTe )

le produit de convolution entre h(t) et x(t) est : z(t) = h(t) ∗ x(t) +∞ X

z(t) =

i=−∞

h(iTe )x(t − iTe )

+∞ X +∞ X

z(t) =

i=−∞ k=0

z(t) =

+∞ X

hi x(t − iTe )

+∞ X

hi

i=−∞

z(t) =

hk δ(iTe − kTe )x(t − iTe )

i=−∞

+∞ X

q=−∞

xq δ(t − iTe − qTe )

pour t = kTe z(kTe ) =

+∞ X

i=−∞

hi

+∞ X

q=−∞

xq δ(kTe − iTe − qTe ) 7

zk =

+∞ X

hi xk−i

i=−∞

et puisque les hi sont nuls pour i < 0, on a : zk =

+∞ X

hi xk−i = yk

i=0

Ainsi on retrouve bien l’expression pr´ec´edente sur yk . Remarquons que la transform´ee de Fourier de y(t) s’´ecrit : ! ∞ +∞ X X T F [y(t)] = hi xk−i e−j2πf kTe k=−∞

T F [y(t)] =

i=0

∞ X +∞ X

hi xk−i e−j2πf iTe e−j2πf (k−i)Te

k=−∞ i=0

T F [y(t)] = T F [h(t)]T F [x(t)] La transform´ee de Fourier du filtre est not´ee T F [h(t)] : T F [h(t)] =

+∞ X

hi e−j2πf iTe = H(ej2πf Te )

i=−∞

En terme d’analyse spectrale on a Syy (f ) = |H(ej2πf Te )|2 Sxx (f ) et Syx (f ) = H(ej2πf Te )Sxx (f )

1.1.4

Analyse de corr´ elation

Dans ce cours on distingue deux types de signaux : les signaux d´eterministes et les signaux al´eatoires : – Les signaux d´eterministes sont les signaux caract´eris´es par leur ´evolution temporelle, par exemple l’´echelon, la rampe, une sinuso¨ıde, etc. Ces signaux sont parfaitement connus et pourvu qu’on connaisse la loi d’´evolution et un ´echantillon `a un instant donn´e, il est possible de retrouver l’int´egralit´e du signal. Cette classe de signaux tend `a mod´eliser les signaux type excitation, entr´ee de syst`eme, point de fonctionnement, etc. – Les signaux al´eatoire correspondent aux signaux qu’on va dire ”incertain”. Ils sont caract´eris´es par une densit´e de probabilit´e et la connaissance d’un ou plusieurs ´echantillons `a un instant donn´e ne suffira jamais `a retrouver l’int´egralit´e du signal. Cette classe de signaux tend `a mod´eliser les bruits tels que les bruits de mesure ou encore les dynamiques n´egligeables d’un syst`eme.

8

Les Signaux d´ eterministes Les signaux d´eterministes se r´epartissent en deux classes : les signaux `a ´energie finie et les signaux `a puissance finie. Pour un signal x(t), l’´energie Exx et la puissance Pxx sont d´efinies par Exx =

∞ X

xk x∗k

k=−∞ N X 1 xk x∗k N →∞ 2N + 1

Pxx = lim

k=−N

Une extension de la notion d’´energie (puissance) est la fonction d’autocorr´elation not´ee Rxx . Cette fonction est d´efinie de diff´erentes mani`eres suivant qu’on consid`ere un signal `a ´energie finie ou `a puissance finie : Signaux a `´ energie finie Rxx (i) =

∞ X

xk+i x∗k

k=−∞

On a ainsi Rxx (0) = Exx . Signaux a ` puissance finie N X 1 xk+i x∗k N →∞ 2N + 1

Rxx (i) = lim

k=−N

On a ainsi Rxx (0) = Pxx . Si on consid`ere un second signal y(t) on d´efinit la fonction d’intercorr´elation comme suit : Signaux a `´ energie finie Ryx (i) =

∞ X

yk+i x∗k

k=−∞

Signaux a ` puissance finie N X 1 yk+i x∗k N →∞ 2N + 1

Ryx (i) = lim

k=−N

On peut montrer que les fonctions de corr´elation et les spectres sont li´es par les relations suivantes : T F [Rxx (t)] = Sxx (f ) T F [Ryx (t)] = Syx (f ) d’o` u les figures 1.2 et 1.3.

9

tranform´ee de Fourier  -

x(t)

X(f )

autocorr´elation

module2

? Rxx (t)

? Sxx (f )

tranform´ee de Fourier  -

Figure 1.2 – relations entre signal, corr´elation, transform´ee de Fourier et spectre

x(t) y(t)

tranform´ee de Fourier  -

X(f ) Y (f )

intercorr´elation

produit et conjugu´e

? Ryx (t)

? Syx (f )

tranform´ee de Fourier  -

Figure 1.3 – relations entre signal, corr´elation, transform´ee de Fourier et spectre Les Signaux al´ eatoires Pour les signaux al´eatoires les fonctions d’autocorr´elation et d’intercorr´elation sont d´efinies par : Γxx (t1 , t2 ) = E{x(t1 )x(t2 )} Γyx (t1 , t2 ) = E{y(t1 )x(t2 )} Un signal al´eatoire ergodique stationnaire est n´ecessairement `a puissance finie. On montre alors que les fonctions de corr´elations d´efinies ci-dessus sont telles que N X 1 xk+i x∗k N →∞ 2N + 1

Γxx (t1 , t2 ) = Rxx (i) = lim

k=−N

10

N X 1 yk+i x∗k N →∞ 2N + 1

Γyx (t1 , t2 ) = Rxx (i) = lim

k=−N

avec i = t1 − t2 . La transform´ee de Fourier d’un signal al´eatoire n’´etant pas d´efinie, les figures 1.2 et 1.3 sont remplac´ees par les figures 1.4 et 1.5. x(t)

autocorr´elation

? Γxx (t)

transform´ee de Fourier  -

Sxx (f )

Figure 1.4 – relations entre signal, corr´elation et spectre

x(t) y(t)

intercorr´elation

? Γyx (t)

transform´ee de Fourier  -

Syx (f )

Figure 1.5 – relations entre signal, corr´elation et spectre

1.2

Objets du cours

Le filtrage adaptatif regroupe un ensemble de techniques d´edi´ees `a la r´esolution d’un probl`eme important en sciences de l’ing´enieur, `a savoir :

11

La reconstruction d’un signal en milieu bruit´e.

Nous entendons ici par reconstruction le fait de lisser, filtrer o` u pr´edire un signal. Les techniques pr´esent´ees ici sont bas´ees sur une certaine connaissance de propri´et´es du signal et du processus le g´en´erant. Le traitement sera r´ealis´e par le biais d’un filtrage du signal lui-mˆeme ou d’un second signal. Ce document est scind´e en deux parties : Partie 1 : cette partie est d´edi´ee `a la pr´esentation de plusieurs types de m´ethodes de filtrage adaptatif : le filtrage au sens de Wiener, le filtrage au sens des moindres carr´es et le filtrage de Kalman. La mise en oeuvre de l’une ou l’autre technique d´ependra de l’objectif `a atteindre, des informations disponibles et de la mani`ere dont sera mis en oeuvre le filtre. Partie 2 : cette partie s’int´eresse plus particuli`erement au probl`eme d’´egalisation, probl`eme pos´e en t´el´ecommunication. Diff´erentes variantes, s’exprimant en terme de filtrage adaptatif, sont propos´ees et discut´ees.

12

Premi` ere partie

Techniques de filtrage adaptatif

13

Les techniques de filtrage adaptatif trouvent tout leur sens dans les probl`emes pour lesquels la composante de bruit ou le processus ont un comportement spectral inconnu. Consid´erons par exemple le cas d’un signal perturb´e par un parasite sinuso¨ıdal `a la fr´equence de 50Hz. Ce signal peut ˆetre filtr´e efficacement par un filtre classique coupe bande centr´e sur 50Hz. En revanche consid´erons le cas de la mesure sur ´electrocardiogramme du rythme cardiaque d’un b´eb´e encore dans le ventre de sa m`ere. Le signal va ˆetre parasit´e par le rythme cardiaque de la m`ere. Ce signal parasite est a priori de contenu spectral inconnu et il risque mˆeme de se superposer en partie au signal correspondant au b´eb´e. Le filtrage classique est donc ici inefficace alors que le filtrage adaptatif va se r´ev´eler performant. On distingue quatres classes d’applications : – l’identification : La figure 1.6 illustre le contexte du probl`eme d’identification. Celui-ci consiste en la d´etermination d’un filtre mod´elisant au mieux le comportement d’un processus inconnu. Seuls sont connus les signaux d’entr´ee/sortie de ce processus. Le filtre repr´esentant le mod`ele sera estim´e `a partir de l’observation de la diff´erence entre la sortie du processus et son estimation `a la sortie du filtre. sortie processus ?

-

entr´ ee

7

+

?

−

6

filtre adaptatif

Figure 1.6 – Principe de l’identification – la pr´ediction : La figure 1.7 illustre le contexte du probl`eme de pr´ediction. Ce probl`eme consiste en l’estimation de la valeur future d’un signal `a partir des informations pass´es. Ce peut ˆetre par exemple pour pr´evoir la position futur d’un objet, ou pour anticiper l’´evolution future d’un grandeur afin de prendre au plus vite une d´ecision.

signal de r´ ef´ erence ?

-

-

retard

7

+

?

−

6

filtre adaptatif pr´ ediction

Figure 1.7 – Principe de la pr´ediction – l’annulation d’interf´erence : La figure 1.8 illustre le contexte du probl`eme d’annulation d’interf´erence. Le probl`eme de l’´electrocardiogramme cit´e auparavant est un probl`eme typique d’annulation d’interf´erence. On dispose 14

d’un signal primaire (´electrocardiogramme du b´eb´e) parasit´e par un signal de r´ef´erence d´eform´e. Ce signal de r´ef´erence est l’´electrocardiogramme de la m`ere. Le filtrage adaptatif va permettre une compensation de l’influence de l’´electrocardiogramme de la m`ere sur l’´electrocardiogramme du b´eb´e. signal de r´ ef´ erence d´ eform´ e + signal primaire ?

signal de r´ ef´ erence

-

7

+

?

−

6

filtre adaptatif

Figure 1.8 – Principe de l’annulation d’interf´erence – la mod´elisation inverse : La figure 1.9 illustre le contexte du probl`eme de mod´elisation inverse. Il s’agit ici de reconstruire au mieux un signal de r´ef´erence qui a ´et´e ”d´eform´e” par un processus inconnu. Le filtre adaptatif doit permettre une compensation des d´eformations induites par le processus. En T´el´ecom ce probl`eme est d´esign´e sous le nom de probl`eme d’´egalisation.

retard

signal de r´ ef´ erence ?

-

-

processus

7

+

?

−

6

filtre adaptatif estimation

Figure 1.9 – Principe de la mod´elisation inverse

15

CHAPITRE

2

Filtrage adaptatif au sens de Wiener

2.1

Principe g´ en´ eraux

Le principe de synth`ese du filtre est pr´esent´e sur la figure 2.1. Il est suppos´e ici que le signal qu’on souhaite estimer, y(t), est corr´el´e `a un second signal u(t). Cette estimation de y(t) sera not´ee yb(t) et prendra la forme suivante : avec

yb(t) = H(q)u(t)

(2.1)

H(q) = h0 + h1 q −1 + h2 q −2 + . . . + hn q −n =

n X

hi q −i

i=0

−1

o` u q est l’op´erateur retard. H(z) est la fonction de transfert du filtre qui est suppos´e ici ˆetre un filtre FIR d’ordre n. Cette fonction de transfert s’´ecrit : H(z) = h0 + h1 z −1 + h2 z −2 + . . . + hn z −n =

n X

hi z −i

i=0

Remarque : Ce type de filtre a l’avantage d’ˆetre toujours stable. D’un point de vue mise en oeuvre, c’est une garantie appr´eciable. yb(t) s’exprime donc comme suit : yb(t) =

n X i=0

hi u(t − iTe )

o` u Te est la p´eriode d’´echantillonnage.

16

y(t)

-

u(t)

7

y(t) b

H(z)

+ -? −

y(t) e

Figure 2.1 – Principe du filtrage de Wiener

2.1.1

Contexte de signaux al´ eatoires et d´ eterministes

– Dans le cas des signaux al´eatoires, le filtrage de Wiener consiste en une d´etermination des coefficients hi du filtre via la minimisation d’une fonction de coˆ ut JW faisant intervenir l’erreur quadratique moyenne entre y(t) et son estimation : JW = E{(y(t) − yb(t))2 } (2.2) Remarquons que dans le cas des signaux al´eatoires stationnaires ergodiques, ce crit`ere peut ˆetre reformul´e de la mani`ere suivante en faisant intervenir la moyenne temporelle : t 1X JW = lim (yk − ybk )2 (2.3) t→∞ t k=1 La minimisation de l’une ou l’autre des fonctions de coˆ ut importe peu en terme de r´esultat final, n´eanmoins les ´ecritures seront plus l´eg`eres pour le crit`ere de nature probabiliste, `a savoir (2.2)

– Dans le cas des signaux d´eterministes, le crit`ere (2.2) n’a pas de sens, n’est alors consid´er´e que le crit`ere (2.3) : t 1X JW = lim (yk − ybk )2 t→∞ t k=1

Dans la suite, afin de ne pas alourdir les ´ecritures, nous d´eterminerons uniquement le vecteur des param`etres minimisant le crit`ere (2.2), on se placera par cons´equent dans un contexte al´eatoire. La r´esolution du second crit`ere ´etant analogue `a celle du premier, le r´esultat sera simplement ´enonc´e et non d´emontr´e.

2.1.2

Formalisme

Dans la suite on notera ye(t) l’erreur d’estimation : ye(t) = y(t) − yb(t) ye(t) = y(t) −

n X i=0

hi u(t − iTe )

Ceci peut aussi se r´e´ecrire sous le formalisme usuel des moindres carr´es comme suit : ye(t) = y(t) − φ(t)T θ 17

avec 

  φ(t) =  

u(t) u(t − Te ) .. . u(t − nTe )

    



  θ= 

et

h0 h1 .. . hn

    

Remarque : Dans certaines situations il peut ˆetre n´ecessaire de faire de la pr´ediction de signaux (pr´esence de retard, anticipation sur une d´ecision), ceci revient alors `a r´ealiser l’estimation d’une valeur future y(t + jTe ) `a partir des donn´ees disponibles jusqu’`a l’instant t. Cette estimation prend la forme : n X yb(t + jTe ) = hi u(t − iTe ) i=0

o` u j d´esigne le pas de pr´ediction. Ceci peut se r´e´ecrire comme suit n X yb(t) = hi u(t − iTe − jTe ) = H(q)u(t) i=0

o` u H(z) est le filtre d´efini par :

H(z) = z −j (h0 + h1 z −1 + h2 z −2 + . . . + hn z −n ) = z −j

n X

hi z −i

i=0

Le crit`ere de Wiener reste le mˆeme, `a savoir (2.2) ou (2.3) suivant le type de signaux trait´es.

2.2

Solution globale

2.2.1

Recherche des param` etres du filtre

Les coefficients du filtre sont obtenus par minimisation de la fonction de coˆ ut JW , par cons´equent ils se caract´erisent par l’annulation de la d´eriv´ee de cette fonction de coˆ ut : dJW =0 dθ Dans un contexte al´eatoire, JW est de la forme :  JW = E (y(t) − φ(t)T θ)2 Il vient alors

 dJW = 2E φ(t)(y(t) − φ(t)T θ) = 0 dθ  E {φ(t)y(t)} = E φ(t)φ(t)T θ

18

(2.4) (2.5)

D’apr`es la forme de φ(t) on a :  u(t)y(t)     u(t − Te )y(t)  E {φ(t)y(t)} = E  ..   .    u(t − nT e )y(t)  E {u(t)y(t)}  E {u(t − Te )y(t)}  = ..  .

nTe )y(t)}  E {u(t − Γyu (0)  Γyu (1)    =  ..   . Γyu (n)

               

o` u Γyu (i) est la fonction de corr´elation entre y(t) et u(t − iTe ) dont l’expression est rappel´ee `a la page 10. De la mˆeme mani`ere on a :   Γuu (0) Γuu (1) · · · Γuu (n)   .. .    Γuu (−1) Γuu (0) . . . T   = E φ(t)φ(t)  .. . . . . . .   . . . . Γuu (−n) ··· · · · Γuu (0)

Dans le cas de signaux r´eels, la fonction d’autocorr´elation est paire, par cons´equent Γuu (−i) = Γuu (i) et il vient alors :   Γuu (0) Γuu (1) · · · Γuu (n)   .. .    Γuu (1) Γuu (0) . . . T   E φ(t)φ(t) =  .. . . . . . .   . . . . Γuu (n) ··· · · · Γuu (0) On notera dans la suite les matrices :   E φ(t)φ(t)T = Γuu E {φ(t)y(t)} = Γyu

L’´equation (2.5), appel´ee ´equation de Wiener-Hopf, peut donc ˆetre exprim´ee de la mani`ere suivante : ´equation de Wiener-Hopf : Γyu = Γuu θ

(2.6)

Le vecteur des param`etres solutions est alors d´efini par : θ = (Γuu )−1 Γyu

(2.7)

19

Remarque : si on souhaite faire de la pr´ediction lin´eaire, c’est `a dire un filtrage de la forme yˆ(t) = H(q)y(t − 1), on constate que la solution pr´ec´edente correspond au syst`eme d’´equations de Yule-Walker vu pour la mod´elisation AR d’un signal (voir cours de traitement num´erique du signal).

2.2.2

V´ erification du principe d’orthogonalit´ e

D’apr`es (2.4), on constate que la solution de l’´equation de Wiener-Hopf v´erifie :  E φ(t)(y(t) − φ(t)T θ) = 0

Ainsi l’erreur d’estimation ye(t) = y(t) − φ(t)T θ est orthogonale au vecteur des entr´ees du filtre φ(t). Ceci correspond au principe d’orthogonalit´e, principe g´en´eral `a toute estimation quadratique. L’estimation yb(t) apparaˆıt ˆetre ainsi la projection orthogonale de y(t) sur le plan engendr´e par φ(t) comme illustr´e sur la figure 2.2. 6 y(t) e

.... .... .... .... .... y(t) .. .. .. .. .. .. .. . y(t) b

6

1

-

j

plan φ(t)

Figure 2.2 – Principe d’orthogonalit´e

2.2.3

Mise en oeuvre

L’algorithme d’application de la solution globale est pr´esent´e sur la figure 2.3.

d´etermination de Γyu et Γuu ? θ = (Γuu )−1 Γyu ? Figure 2.3 – Solution de Wiener-Hopf Cet algorithme relativement simple souffre de quelques d´efauts. Tout d’abord il n´ecessite la connaissance des termes de corr´elation Γuu (i) et Γyu (i) pour i ∈ [0; n], ceci est dˆ u au fait que la solution (2.9) est obtenu par minimisation du crit`ere (2.2) de nature probabiliste. En g´en´eral on ne dispose que de r´ealisations des signaux al´eatoires u(t) et y(t) et non des termes Γuu (i) et Γyu (i).

20

Remarquons que l’hypoth`ese d’ergodicit´e ´etant v´erifi´ee il est possible de substituer au crit`ere (2.2) le crit`ere (2.3). On peut ais´ement d´emontrer que le vecteur des param`etres minimisant (2.3) est d´etermin´e par : θ = (Ruu )−1 Ryu avec

(

(2.8)

PN Ruu (i) = limN →∞ N1 k=1 uk+i uk P N Ryu (i) = limN →∞ N1 k=1 yk+i uk

Puisque les signaux sont al´eatoires stationnaires ergodiques on a  Ruu (i) = Γuu (i) Ryu (i) = Γyu (i) La solution (2.8) est donc bien ´equivalente `a la solution de l’´equation de WienerHopf (2.7). D’un point de vue pratique, ne disposant des donn´ees que jusqu’`a l’instant t = N Te on notera ( PN buu/t (i) = 1 R N P k=1 uk+i uk byu/t (i) = 1 N yk+i uk R k=1 N b uu/t et En substituant respectivement Γuu (ou Ruu ) et Γyu (ou Ryu ) par R b Ryu/t dans l’´equation (2.7) on obtient : b uu/t )−1 R b yu/t θbt = (R

(2.9)

o` u θbt d´esigne une estimation de θ `a partir des donn´ees jusqu’`a l’instant t = N Te . L’algorithme de mise en oeuvre est pr´esent´e sur la figure 2.4. b yu/t et R b uu/t calcul de R ? b uu/t )−1 R b yu/t θbt = (R ?

Figure 2.4 – Algorithme de mise en oeuvre de la solution de Wiener-Hopf Ensuite, cet algorithme n´ecessite l’inversion d’une matrice Γuu (ou Ruu ), probl`eme plus d´elicat puisqu’il s’agit d’une op´eration coˆ uteuse en temps de calcul et risqu´ee en terme de fiabilit´e. Enfin, si apr`es une premi`ere estimation du vecteur des param`etres, on dispose de nouvelles donn´ees que l’on veut exploiter pour am´eliorer la premi`ere estimation, on est oblig´e de reprendre tous les calculs depuis le d´ebut. 21

Ces deux derniers probl`emes peuvent ˆetre contraignant dans le cas de l’utilisation d’un calculateur ne permettant pas l’inversion de la matrice Γuu (ou Ruu ) ou encore dans le cas d’un probl`eme n´ecessitant une application en temps r´eel pour laquelle il est n´ecessaire d’actualiser `a chaque instant la valeur des param`etres du filtre. Un algorithme r´ecursif sera propos´e dans la suite afin de permettre un contournement intuitive des deux probl`emes pr´ec´edents.

2.2.4

Exemple

Dans le cas d’un filtre d’ordre 0, l’estimation yb(t) prend la forme yb(t) = h0 u(t) et la fonction de coˆ ut devient : JW = E{(y(t) − h0 u(t))2 }

JW = Γy (0) − 2h0 Γyu (0) + h20 Γuu (0) avec

  Γyy (0) = E{y(t)2 } Γuu (0) = E{u(t)2 }  Γyu (0) = E{y(t)u(t)}

JW correspond `a une fonction du second degr´e en h0 . Sa d´eriv´ee par rapport `a h0 s’exprime somme suit : dJW = −2Γyu (0) + 2h0 Γuu (0) dh0 d’o` u dJW Γyu (0) = 0 ⇒ h0 = dh0 Γuu (0) Pour y(t) et u(t) tels que :   Γyy (0) = 10 Γuu (0) = 1  Γyu (0) = 3

(2.10)

on a h0 = 3. Les informations sur Γuu (0) et Γyu (0) sont g´en´eralement estim´ees `a partir de s´eries de donn´ees. Afin d’illustrer ceci, nous avons utilis´e des s´eries de donn´ees correspondant aux r´ealisations de signaux al´eatoires satisfaisants les relations (2.10). L’hypoth`ese d’ergodicit´e ´etant par nature v´erifi´ee ici, Γuu (0) et Γyu (0) peuvent ˆetre estim´es, pour t ≫ 1, comme suit : buu/t (0) Γuu (0) ≃ R byu/t (0) Γyu (0) ≃ R

La figure 2.5 pr´esente l’estimation de h0 `a partir de ces r´ealisations et en fonction du nombre de donn´ees utilis´ees. Il apparaˆıt que cette estimation converge vers la valeur vraie de h0 .

22

3.04

3.03

estimation de h

0

3.02

3.01

3

2.99

2.98

2.97

2.96

2.95

0

0.5

1

1.5

2

2.5 3 nombre de points

3.5

4

4.5

5 4

x 10

Figure 2.5 – Estimation de h0

2.3

Solution r´ ecursive

2.3.1

Premi` ere approche

Consid´erons le cas envisag´e dans une des remarques pr´ec´edentes `a savoir : on suppose disposer d’une premi`ere estimation du vecteur des param`etres (not´ee θbt−Te ) et vouloir am´eliorer cette estimation en exploitant de nouvelles informations disponible `a l’instant t. Ces nouvelles informations sont y(t) et u(t). La solution r´ecursive propose une actualisation de θbt−Te sans passer par l’ensemble des calculs n´ecessaires dans la solution globale et notamment sans passer par l’inversion de la matrice Γuu .

Remarque : Signalons tout de suite que la solution sur θ propos´ee par l’algorithme r´ecursif n’est qu’une approximation de la solution propos´ee par l’´equation de WienerHopf. Il sera vu plus loin les conditions n´ecessaires pour que cette approximation ne soit pas trop grossi`ere. Intuitivement la nouvelle estimation du vecteur des param`etres va prendre la forme : θbt = θbt−Te + terme correctif

Le terme correctif va prendre en compte la d´eriv´ee de la fonction de coˆ ut par rapport `a l’ancien jeu de param`etres θbt−Te , ceci afin de se d´eplacer dans le sens d’une erreur la plus faible possible. Ceci est illustr´e sur la figure 2.6. Notons JW (t/t − Te ) la valeur de la fonction de coˆ ut `a l’instant t, connaissant les param`etres estim´es `a l’instant t − Te : JW (t/t − Te ) = E{(y(t) − φ(t)T θbt−Te )2 }

– Si la d´eriv´ee

dJW (t/t−Te ) dθbt−Te

est positive (c’est le cas sur la figure de gauche) alors

on va diminuer θbt−Te de mani`ere `a se d´eplacer vers le minimum de la fonction de coˆ ut. Le terme correctif sera donc choisi n´egatif. – Si la d´eriv´ee dJW b(t/t−Te ) est n´egative (c’est le cas sur la figure de droite) dθt−Te

alors on va augmenter θbt−Te de mani`ere `a se d´eplacer vers le minimum de la fonction de coˆ ut. Le terme correctif sera donc choisi positif. 23

JW (t/t − Te )

JW (t/t − Te )

6

6 dJW (t/t−Te ) b dθ t−Te

. .. . .. .. .. .

 dJW (t/t−Te ) b dθ t−Te

-

b θ t−Te

.. .. .. w .. ..

b θ

b θ t−Te

-

b θ

Figure 2.6 – Une premi`ere version possible de l’algorithme est donc : ! µ dJW (t/t − Te ) b b θt = θt−Te − signe 2 dθbt−T e

avec µ > 0.

Le terme signe



dJW (t/t−Te ) dθbt−T e



va permettre une ´evolution des param`etres de

mani`ere `a minimiser la fonction de coˆ ut, le terme µ2 permet l’ajustement de l’amplitude du saut de l’ancienne `a la nouvelle estimation. Une difficult´e sur l’algorithme pr´ec´edent est le fait que cette amplitude est constante, d’o` u le risque de voir l’algorithme osciller au bout d’un moment autour de la solution sans parvenir `a se stabiliser.

2.3.2

Algorithme du gradient d´ eterministe

Afin de prendre en compte l’amplitude de la d´eriv´ee et ´eviter dans la mesure du possible ce ph´enom`ene d’oscillation, une autre forme possible de l’algorithme est la suivante :

On a

µ dJW (t/t − Te ) θbt = θbt−Te − 2 dθbt−Te dJW (t/t − Te ) = −2E{φ(t)(y(t) − φ(t)T θbt−Te )} dθbt−T e

dJW (t/t − Te ) = −2(Γyu − Γuu θbt−Te ) dθbt−T

(2.11)

e

Ceci permet l’expression suivante sur θbt : algorithme du gradient : θbt = θbt−Te + µ(Γyu − Γuu θbt−Te )

24

(2.12)

Cet algorithme est aussi appel´e algorithme du gradient d´eterministe puisque les quantit´es Γyu et Γuu sont parfaitement connues et d´etermin´ees. Dans le cas de signaux al´eatoires stationnaires ergodiques ou dans le cas de signaux d´eterministes, cet algorithme s’´enonce de la mani`ere suivante : algorithme du gradient : θbt = θbt−Te + µ(Ryu − Ruu θbt−Te )

2.3.3

(2.13)

Propri´ et´ e de convergence

Pour l’instant, rien ne garanti la convergence des param`etres de la solution r´ecursive θbt vers les param`etres de la solution globale. Il s’agit d’´etablir ici les conditions sur µ garantissant une telle convergence.

Notons θet = θ − θbt l’erreur entre le vecteur des param`etres de la solution globale θ (fournit par l’´equation de Wiener-Hopf) et le vecteur des param`etres de la solution r´ecursive θbt (fournit par l’algorithme du gradient). D’apr`es ce qui a ´et´e vu auparavant on a : θet = θ − θbt−Te − µ(Γyu − Γuu θbt−Te ) et Γyu = Γuu θ

Il vient :

d’o` u

θet = (1 − µΓuu )(θ − θbt−Te ) θet = (1 − µΓuu )θet−Te

ou encore

θet = (1 − µΓuu )N θe0

o` u t = N Te et θe0 d´esigne l’erreur au temps initial t = 0.

Γuu est une matrice sym´etrique d´efinie positive, on peut donc la d´ecomposer comme suit Γuu = P DP T avec D de la forme :  λ1 0 · · · 0  .. . .  0 λ2 . .  D= . . . .. ..  .. 0 0 · · · 0 λn+1

     

o` u les λi d´esigne les valeurs propres de Γuu . θet peut donc ˆetre mis sous la forme suivante : θet = P (1 − µD)N P T θe0

25

avec

θet = P D′N P T θe0 

  D =   ′

1 − µλ1

0

0 .. .

1 − µλ2 .. .

0

···

··· 0 .. .. . . .. . 0 0 1 − µλn+1

     

Il apparaˆıt deux cas figure suivant la nature des valeurs propres 1 − µλi : – Si elles sont toutes de valeur absolue inf´erieure `a 1 alors : lim P D′N P −1 = 0 t→∞ d’o` u lim θet = 0 t→∞ On a donc convergence du vecteur des param`etres de la solution r´ecursive vers la solution globale ; – Si au moins une valeur propre est de valeur absolue sup´erieure `a 1 alors lim P D′N P −1 = ∞ t→∞ d’o` u lim θet = ∞ t→∞

On a donc divergence du vecteur des param`etres θbt .

La condition de convergence de l’algorithme apparaˆıt donc ˆetre : −1 < 1 − µλi < 1 , ∀ i 0