Filtrage Optimal [PDF]

Zitiervorschau

Filtrage optimal par

Mohamed NAJIM Professeur à l’École nationale supérieure d’électronique et de radioélectricité de Bordeaux (ENSERB)

1. 1.1 1.2 1.3

Filtre adapté ............................................................................................... Définitions ..................................................................................................... Filtre adapté dans le cas de bruit blanc ...................................................... Cas d’un bruit coloré ....................................................................................

— — — —

2 2 3 4

2. 2.1 2.2 2.3 2.4 2.5 2.6

Filtre de Wiener ......................................................................................... Position du problème................................................................................... Équation de Wiener-Hopf............................................................................. Calcul de l’erreur dans le cas du filtre non physiquement réalisable ...... Équation de Wiener dans le cas de spectres rationnels............................ Filtre de Wiener discret ................................................................................ Application du filtre de Wiener non causal discret pour le débruitage de la parole ...................................................................................................

— — — — — —

6 6 7 8 8 9

—

10

Filtre de Kalman ........................................................................................ Propriétés du filtre de Kalman..................................................................... Algorithme du filtre de Kalman................................................................... Initialisation et mise en œuvre du filtre de Kalman................................... Filtrage de Kalman d’un signal de parole noyé dans un bruit blanc gaussien ........................................................................................................

— — — —

11 11 12 14

—

16

4. 4.1 4.2 4.3

Introduction au filtrage adaptatif ........................................................ Différents algorithmes ................................................................................. Algorithme du type gradient stochastique ou (LMS) ................................ Application à l’annulation du bruit .............................................................

— — — —

18 18 18 20

5.

Annexe. Quelques propriétés de la matrice R ..................................

—

22

3. 3.1 3.2 3.3 3.4

ingénieur doit souvent considérer le cas courant où l’on souhaite, à partir d’un message brut ou signal observé m (t), contenant un signal utile − signal désiré − et un bruit, à déterminer le meilleur récepteur − optimal − permettant de discriminer le signal du bruit. Par récepteur ou filtre optimal, nous entendons un filtre qui satisfait à un certain critère d’optimalité sous des hypothèses que nous préciserons. Par filtre, nous entendons une description mathématique des opérations de traitement que subit le signal mélangé au bruit. Auparavant, nous devons préciser : 1°) que les entrées de ces filtres seront soit des processus aléatoires, soit une combinaison de signaux déterministes et aléatoires. Nous disposerons en général d’un nombre minimal d’informations caractérisant ces entrées ; 2°) que nous ne considérons uniquement les systèmes stationnaires linéaires. Dans les cas où une réalisation matérielle est recherchée, il y aura lieu de considérer la réalisabilité du filtre. Il sera souvent utile de connaître le système optimal, même s’il n’est pas physiquement réalisable. Sa connaissance permettra de mesurer et d’apprécier les performances des systèmes réalisables mais non optimaux.

L’

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite. © Techniques de l’Ingénieur, traité Mesures et Contrôle

R 7 228 − 1

FILTRAGE OPTIMAL _____________________________________________________________________________________________________________________

Nous traiterons trois types de filtres : 1 - le filtre adapté ; 2 - le filtre de Wiener ; 3 - le filtre de Kalman. Ces différents filtres correspondent respectivement à une solution dans les cas où : 1°) le signal désiré est de forme connue. Il est mélangé soit à un bruit blanc, soit à un bruit coloré ; 2°) le signal est, à l’instar du bruit, un processus aléatoire. Le filtre développé par Norbert Wiener constitue une solution non récursive, difficile à implanter sur ordinateur ; 3°) le signal et le bruit sont aléatoires. Le filtre de Kalman est une solution récursive du problème du filtrage qui généralise le filtrage de Wiener. Le développement de ces filtres suppose que l’on dispose d’informations a priori à la fois sur les signaux et sur les bruits. Il s’agit, en particulier, de la connaissance des fonctions ou des matrices d’autocorrélation. Dans le cas où leur connaissance nous fait défaut, on aura comme alternative l’utilisation des filtres adaptatifs. Ces derniers « apprennent » les caractéristiques des signaux au fur et à mesure que ceux-ci se déroulent. On a néanmoins montré récemment qu’une famille de filtres adaptatifs couramment utilisés dans les applications peut, elle aussi, être considérée comme optimale.

1. Filtre adapté

Le filtre étant linéaire, nous aurons : t1

s2 ( t1 ) =

∫

s1 ( τ ) h ( t1 Ð τ ) d τ

(2)

b1 ( τ ) h ( t1 Ð τ ) d τ

(3)

Ð∞

1.1 Définitions

t1

Soit le message observé m (t ) = s1 (t ) + b1 (t ). Le signal s1 (t ), mélangé à un bruit b1 (t ), est supposé être de forme connue. On se propose de réaliser un filtre de réponse impulsionnelle h (t ) qui permette de détecter la présence du signal s1 (t ). Le filtre délivre un signal s2 (t ) mélangé à un bruit b2 (t ) (figure 1). Le filtre cherché doit maximaliser le rapport Ro (t ) :

∫

Ð∞

La borne supérieure de l’intégrale est prise égale à t1 au lieu de l’infini pour que le filtre h (t ) soit physiquement réalisable. En effet, le filtre est réalisable si :

puissance de s 2 ( t ) R o ( t ) = ---------------------------------------------------------puissance de b 2 ( t )

h (t ) = 0 pour t < 0 En prenant t1 comme borne supérieure de l’intégrale de convolut 1 Ð τ > 0 et, de ce fait, le tion au lieu de l’infini, on est sûr que ∀τ filtre obtenu sera réalisable.

à un instant t = t1 c’est-à-dire :

s 22 ( t 1 ) R o ( t 1 ) = -----------------------------E { b 22 ( t 1 ) }

b2 ( t1 ) =

(1)

On supposera dans la suite de l’exposé que s (t ) et b (t ) sont des processus stationnaires. La stationnarité signifie que les paramètres statistiques (moyenne, variance et l’ensemble des moments) ne dépendent pas du temps.

En d’autres termes, on se propose de déterminer quel est le filtre de réponse impulsionnelle h (t ) qui maximalise le rapport signal à bruit Ro (t1 ).

On se suffit souvent, dans la pratique, de l’hypothèse des processus stationnaires du deuxième ordre pour lesquels seuls les deux premiers moments sont indépendants du temps. En outre, pour les processus stationnaires, si l’on note K b (τ , σ ) ∆ = E { b 1 ( τ ) b 1 ( σ ) } , la fonction d’autocorrélation des bruits, celle-ci ne dépend que de l’écart entre les instants τ et σ . Elle s’écrira :

où

E{ }

désigne l’espérance mathématique

m(t ) = s1(t ) + b1(t )

h (t )

Figure 1 – Signal délivré par le filtre h (t )

R 7 228 − 2

s 2(t ) + b2(t )

kb ( τ , σ ) = kb ( τ Ð σ )

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite. © Techniques de l’Ingénieur, traité Mesures et Contrôle

(4)

____________________________________________________________________________________________________________________

Si nous l’appliquons au numérateur de la relation (8) en prenant f ( ω ) = h ( ω ) et g ( ω ) = s 1 ( ω ) exp j ωt 1 il vient :

La relation (1) devient, compte tenu de (2) et (3) : t1

2

∫

s 1 ( τ ) h ( t 1 Ð τ ) dτ

+∞

Ð∞

R o ( t 1 ) = ----------------------------------------------------------------------------------------------------------------------t t

∫ ∫ 1

1

Ð∞

Ð∞

(5)

k b ( τ , σ ) h ( t 1 Ð τ ) h ( t 1 Ð σ ) dσ d τ

∫

■ Calcul de Ro (t1) dans le cas de bruit blanc Si b (t ) est un bruit blanc de densité spectrale de puissance No N ϕ ( ω ) = ------o- et de fonction d’autocorrélation k b ( τ ) = ------- δ ( τ ) , il 2 2 vient :

Ð∞

Ð∞

∫ ∫ N = ------o2

t1

t1

Ð∞

Ð∞

∫ ∫

s1 ( ω ) 2 d ω

Ð∞

∫

s1 ( ω ) 2 d ω

(11)

Ð∞

Ce majorant de Ro (t1) ne dépend pas de h (t ). Il ne dépend que de l’énergie du signal et du bruit, par conséquent Ro (t1) sera maximum quand l’égalité aurla lieu. Le filtre obtenu dans ce cas sera le filtre optimal hop que nous cherchions.

h op ( ω ) s 1 ( ω ) exp { j ωt 1 } = h op ( ω ) h op ( Ð ω ) = s 1 ( ω ) s 1 ( Ð ω ) d′où h op ( ω ) = s 1 ( Ð ω ) exp { Ð j ωt 1 } = s 1* ( ω ) exp { Ð j ωt 1 } car s 1* ( ω ) = s 1 ( Ð ω ) +∞

∫

s 1 ( Ð ω ) exp { Ð j ω ( t 1 Ð t ) } d ω

Ð∞

δ ( τ Ð σ ) h ( t 1 Ð τ ) h ( t 1 Ð σ ) dτ d σ

soit h op ( t ) = s 1 ( t 1 Ð t )

h 2 ( t1 Ð τ ) d τ

(6)

Ð∞

(12)

Il apparaît ainsi que le filtre est adapté à la forme du signal. Ceci justifie le nom de filtre adapté. ■ Réalisabilité du filtre

Supposons que h (t ) remplit effectivement la condition de réalisabilité, dès lors :

∫

Revenons à la réalisabilité de h op ( t ) :  0 pour t < 0 h op ( t ) =   s 1 ( t 1 Ð t ) pour t > 0

+∞

t1

∫

devient

Ð∞

et R o ( t 1 ) devient :

Ð∞

+∞

∫

s 1 ( τ ) h ( t 1 Ð τ ) dτ (7)

∫

Cette dernière relation signifie qu’à l’instant t1 tout le signal a traversé le filtre. Nous avons le signal à la sortie : t

s2 ( t ) =

Cette écriture nous permet de passer du domaine temporel au domaine fréquentiel où la résolution est plus aisée.

+∞

t

s2 ( t ) =

∫

s1 ( τ ) s1 ( τ Ð t + t1 ) d τ

Ð∞

2

Cette dernière relation présente une grande similitude avec la fonction d’autocorrélation de s (t ) :

h ( ω ) s 1 ( ω ) exp ( j ωt 1 ) dω

Ð∞

R o ( t 1 ) = --------------------------------------------------------------------------------------------------+∞ No ⁄ 2

∫

(8) 1 k s ( τ ) = lim ------T → ∞ 2T

h ( ω ) 2 dω

Ð∞

or h ( ω ) 2 = h ( ω ) h ( Ð ω ) = h ( ω ) h∗ ( ω )

(9)

ω2

2


t 1 , s 1 ( t ) = 0 à la sortie du filtre.

2

Ð∞ R o ( t 1 ) = ------------------------------------------------------------------------No + ∞ 2 ------h ( t 1 Ð τ ) dτ 2 Ð∞

∫

∫

+∞

Ro ( t1 ) < ( 2 ⁄ No )

où h op ( t ) =

t1

∫

∫

+∞

h (ω) 2dω

Ð∞

Ð∞

k b ( τ , σ ) h ( t 1 Ð τ ) h ( t 1 Ð σ ) dτ d σ

= ( No ⁄ 2 )


0  h op ( t ) =  =  a exp { Ð α ( t 1 Ð t ) } t t1

Pour ce faire, nous allons remplacer h (t ) par :

0 < t < t1

t 0

1.3 Cas d’un bruit coloré

la condition nécessaire donnant le minimum de Q est :

1.3.1 Position du problème

∂Q (γ ) -----------------∂γ

Nous allons considérer ici le cas d’un signal noyé dans un bruit coloré de fonction d’autocorrélation kb.

(17)

= 0 γ=0

Quel que soit l’incrément δ h , cette relation conduit à A = 0, soit :

Reprenons l’expression (5) : t1

∫

Q  ∂------- =  ∂γ 

2

s1 ( τ ) h ( t1 Ð τ ) d τ

∫ ∫ Ð∞

Ð∞

Ð∞

[ h op ( t 1 Ð τ ) δh ( t 1 Ð σ ) + h op ( t 1 Ð σ ) δh ( t 1 Ð τ )

+ 2 γδh ( t 1 Ð τ ) δh ( t 1 Ð σ ) ] k b ( τ Ð σ ) d σ d τ

1

kb ( τ , σ ) h ( t1 Ð τ ) h ( t1 Ð σ ) d σ d τ

t1

Ð∞

s 22 ( t 1 ) = -----------------------------E { b 22 ( t 1 ) }

е (14)

Quand R o ( t 1 ) aura atteint sa valeur maximale, cela signifiera que la puissance du bruit E { b 22 ( t 1 ) } sera minimale. Ainsi, si l’on considère la relation :

Q = E { b 22 ( t 1 ) } Ð µs 22 ( t 1 )

R 7 228 − 4

t1

∫ ∫

Ð∞ R o ( t 1 ) = --------------------------------------------------------------------------------------------------------------------------t t 1

t1

(15)

∫

δh ( t 1 Ð τ ) s 1 ( τ ) d τ

Ð∞

or t1

t1

Ð∞

Ð∞

t1

t1

Ð∞

Ð∞

∫ ∫

h op ( t 1 Ð σ ) δh ( t 1 Ð τ ) k b ( σ Ð τ ) d τ d σ =

∫ ∫

h op ( t 1 Ð τ ) δh ( t 1 Ð σ ) k b ( τ Ð σ ) d τ d σ

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite. © Techniques de l’Ingénieur, traité Mesures et Contrôle

(18)

____________________________________________________________________________________________________________________

Il vient, en tenant compte de (17) et de la parité de kb : kb ( τ Ð σ ) = kb ( σ Ð τ ) : t1

Ð∞

∫

2 δh ( t 1 Ð τ ) h op ( t 1 Ð σ ) k b ( τ Ð σ ) d τ d σ

Ð∞

h op ( t 1 Ð σ ) exp { Ð j ωσ } ϕ b ( ω ) dσ = s ( ω )

t1

∫

δh ( t 1 Ð τ ) s 1 ( τ ) d τ = 0

Ð∞

En opérant le changement de variable défini par x = t 1 Ð σ , l’équation (23) devient : +∞

ou encore :

∫

exp { Ð j ωt 1 } ϕ b ( ω )

t1

∫

δh ( t 1 Ð τ )d τ

Ð∞

2 h op ( t 1 Ð τ ) k b ( τ Ð σ ) d σ Ð µ s 1 ( τ ) = 0 (19)

Ð∞

Cette relation doit être satisfaite quel que soit δh , donc (19) se réduit à : t1

Ð∞

µ h op ( t 1 Ð σ ) k b ( τ Ð σ ) d σ =  --- s 1 ( τ ) 2

Ð ∞ < τ < t1

µ Le facteur  --- modifie uniquement le gain du filtre, il aura donc la  2 même incidence sur le bruit et sur le signal. Si nous le prenons égal à 1, h op sera alors connue à un facteur constant près, et l’équation générale du filtre adapté pour un bruit coloré devient :

h op ( x ) exp { j ωx } dx = s ( ω )

Ð∞

t1

∫

(23)

Ð∞

е

∫

Si ϕ b ( ω ) est la densité spectrale du bruit b (t ), l’équation (22) devient : ∞

t1

∫ ∫

FILTRAGE OPTIMAL

ou encore :

s (Ð ω) H op ( ω ) = ------------------- exp { Ð j ωt 1 } ϕb ( ω ) Ceci peut être une solution approchée dans certains cas où l’on ne se pose pas le problème de la réalisabilité du filtre, ou bien dans le cas où l’on dispose d’enregistrements. En effet, si l’on opère sur des enregistrements, il est évident qu’à chaque instant on connaît le passé et le futur du signal et il n’y a plus de difficulté à écrire l’intégrale de { Ð ∞ à + ∞ } .

1.3.3 Solution du filtre adapté par les techniques du blanchissement ou « blanchiement »

t1

∫

h op ( t 1 Ð σ ) k b ( τ Ð σ ) d σ = s 1 ( τ )

(20)

Ð∞

Pour nous assurer que ∂Q (γ ) -----------------∂γ

= 0 γ=0

Nous introduirons successivement la factorisation spectrale et la mise en forme spectrale avant d’aborder la résolution de l’équation (23) proprement dite.

est aussi une condition suffisante, il faut calculer la quantité : ∂2 Q ---------- = ∂γ 2

t1

t1

Ð∞

Ð∞

∫ ∫

En réalité, le problème délicat reste la résolution de l’équation (20) et naturellement la difficulté majeure provient du fait que la borne supérieure de l’intégrale soit t1 au lieu de + ∞ . Nous exposerons la technique du blanchissement couramment utilisée dans la résolution des équations de Fredholm qui jouent un rôle important dans divers domaines de la théorie de la détection et de l’estimation.

2 δh ( t 1 Ð τ ) δ h op ( t 1 Ð σ ) k b ( τ Ð σ ) d τ d σ

1.3.3.1 Factorisation spectrale

et s’assurer qu’elle est positive. Ceci est en effet le cas, car k b est définie positive. ∂2 Q ---------- sera positive ou nulle pour δ h = 0 . ∂γ 2

Si l’on considère un processus aléatoire de densité spectrale de puissance (dsp) ϕ ( ω ) et de fonction d’autocorrélation k (t ), le théorème de Wiener-Kinchin donne : +∞

ϕ (ω) =

∫

k ( τ ) exp { Ð j ωτ } d τ

Ð∞

1.3.2 Filtre adapté non physiquement réalisable ■ Rappelons que : Si l’on ne s’impose pas de contrainte de réalisabilité, l’équation du filtre adapté devient :

+∞

+∞

∫

1°) ϕ ( ω ) est une fonction réelle ; en effet :

k b ( τ Ð σ ) h op ( t 1 Ð σ ) d σ = s ( τ )

Ð ∞ 0 pour s’assurer également que c’est une condition suffisante. Nota : nous n’avons fait aucune hypothèse sur la réalisabilité du filtre. Dans le cas d’un traitement en temps réel, où il y a nécessité d’utiliser un filtre réalisable, il faut modifier les limites de l’équation intégrale, la borne Ð ∞ devant être remplacée par 0. Cela ne simplifie naturellement pas la résolution de l’équation de Wiener-Hopf, car on exclut l’utilisation de la transformée de Fourier.

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite. © Techniques de l’Ingénieur, traité Mesures et Contrôle

R 7 228 − 7

FILTRAGE OPTIMAL _____________________________________________________________________________________________________________________

Dans le cas où l’on ne se pose pas le problème de la réalisabilité, l’équation de WienerHopf devient, après application de la transformée de Fourier :

ϕ ms ( ω ) ( ω ) = ---------------------ϕ mm ( ω )

H+

ϕ ϕss

(34)

ϕbb

et dans le cas de signaux non corrélés, ϕ bs ( ω ) = 0 :

ϕ ss ( ω ) H op ( ω ) = --------------------------------------------ϕ ss ( ω ) + ϕ bb ( ω )

(35)

0

f

Figure 4 – Signaux s (t ) et b (t ) dont les densités spectrales n’ont pas de support commun

2.3 Calcul de l’erreur dans le cas du filtre non physiquement réalisable ϕ ϕss

Si l’on reprend l’équation (27) donnant l’erreur :

ϕbb

+∞

E { e 2 ( t ) } = k ss ( 0 ) Ð 2

∫

k ms ( τ ) h ( τ ) d τ

Ð∞

0

+∞

∫∫

+

k mm ( σ Ð τ ) h ( σ ) h ( τ ) d σ d τ

Ð∞

on peut faire apparaître sous l’intégrale double k ms en utilisant l’équation de Wiener-Hopf (33), on obtient : +∞

E { e 2 ( t ) } = k ss ( 0 ) Ð

∫

k ms ( τ ) h ( τ ) d τ

Ð∞

∫

E { e2 ( t ) } =

+∞

ϕ ss (f ) d f Ð

Ð∞

∫

ϕ ms (f ) H ( Ð f ) d f

Ð∞

+∞

=

∫

Ð∞

2.4 Équation de Wiener dans le cas de spectres rationnels

ϕ ms ( ν ) ξ ( ν ) = ---------------------- est une fraction rationnelle, qui peut être ϕ mm ( ν ) décomposée en éléments simples. Elle peut alors se mettre sous la forme :

ϕ sö s (f ) Ð ϕ ms (f ) H ( Ð f ) d f

+∞

min E { e 2 ( t ) } =

ξ (ν) = ξ + (ν) + ξ − (ν) ξ+

Si l’on tient compte de l’expression de la fonction de transfert optimale dont l’expression est donnée par la relation (35) il vient :

ϕ ss (f ) ϕ mm (f ) Ð ϕ ms (f ) 2 -------------------------------------------------------------------------- df ϕ mm ( f ) Ð∞

∫

ou, si s (t ) et b (t ) ne sont pas corrélés : +∞

min E { e 2 ( t ) } =

ϕ ss (f ) ϕ bb (f ) -------------------------------------------- df ϕ ( f ) + ϕ bb ( f ) Ð ∞ ss

∫

Cette relation est intéressante pour commenter les résultats obtenus sur le filtre de Wiener. En effet, supposons que les signaux s (t ) et b (t ) soient tels que leurs densités spectrales n’aient pas de support commun comme cela est représenté dans la figure 4. Le produit ϕ ss . ϕ bb sera quasi nul et conduira à une erreur d’estimation égale à zéro. Dans ce cas, on n’aura pas besoin d’utiliser le filtre de Wiener pour séparer le signal du bruit. Il suffira d’utiliser un simple filtre passe-bas pour éliminer le bruit. Par contre, dans le cas de signaux réels, nous aurons un recouvrement des spectres du signal et du bruit, tel que cela est schématisé dans la figure 5. Dans ce cas, seul le filtre de Wiener permettra une estimation du signal. On constate que le produit ϕ ss . ϕ bb dans la bande de fréquence qui nous intéresse ne sera pas nul et, de ce fait, l’erreur d’estimation ne sera jamais nulle.

R 7 228 − 8

Figure 5 – Signaux s (t ) et b (t ) réels, avec recouvrement des spectres

Dans le cas de spectres rationnels, la quantité :

ou, en remplaçant k ms par sa transformée de Fourier notée ϕ ( ω ) : +∞

f

−

( ν ) et ξ ( ν ) ont respectivement leurs singularités, c’est-àoù dire les pôles et les zéros, dans le demi-plan supérieur et le demiplan inférieur du plan complexe. L’expression (34) de la fonction de transfert du filtre optimal devient : 1 H ( ω ) = -----------------------ϕ mm ( ω )

∞

∫

+∞

exp { Ð 2πj ft } d t

∫

ξ ( ν ) exp ( 2 πj νt ) dν (36)

Ð∞

0

où les bornes d’intégration 0 et l’infini correspondent à la partie réalisable. On peut faire abstraction de la première intégrale, si l’on ne considère, dans ξ ( ν ) , que la partie qui donne bien un filtre physiquement réalisable, c’est-à-dire ξ + ( ν ) . On adopte, traditionnellement, les notations suivantes : [ ξ ] + = partie réalisable de ξ ( ν ) = ξ + ( ν ) L’équation (36) devient :

ϕ ms ( ω ) 1 H ( ω ) = ----------------------- -----------------------+ ϕ mm ( ω ) ϕ mm ( ω )

+

Nota : nous avons porté notre attention essentiellement sur les spectres rationnels. Les spectres non rationnels peuvent être approchés par des spectres rationnels. On peut montrer qu’une condition nécessaire et suffisante de « factorabilité » est que l’intégrale : +∞

∫

Ð∞

log ϕ mm ( ω ) ----------------------------------- d f 1+f2

soit convergente.

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite. © Techniques de l’Ingénieur, traité Mesures et Contrôle

____________________________________________________________________________________________________________________

Exemple : soit à déterminer le filtre optimal, sachant que ϕ b ( ω ) = k 2 (bruit blanc) et : 1 ϕ ss ( ω ) = ------------------ avec ( p = jω ) 1 Ð p2

La figure 7 montre la sortie yö ( n ) du filtre comparée à une sortie désirée d (n ) ; cette dernière joue le rôle de la sortie réelle y ( n ) que l’on cherche à approximer. Ainsi, on peut noter indifféremment l’erreur :

e ( n ) = [ y ( n ) Ð yö ( n ) ]

b2 Ð p2 1 + k2 avec b 2 = -----------------On a : ϕ mm = ϕ ss + ϕ bb = k 2 --------------------k2 1 Ð p2

et :

b + jω b Ð jω ϕ mm = k  ------------------- k  ------------------ = ϕ ++ .ϕ Ð Ð  1 + jω  1 Ð jω 1 1 ϕ ms = -----------------2- = -------------------21+ω 1Ðp

ϕ ms ----------Ð ϕ mm

+

ou encore : e ( n ) = [ d ( n ) Ð yö ( n ) ]

T (n) X (n) HN N

l’erreur sera : T (n) X (n) = d (n) Ð XT (n) H (n) e ( n ) = d ( n ) Ð HN N N N

J (n) = ε (n) = E

1 1 d’où : H ( jω ) = -------------------------- -----------------k2 ( 1 + b ) ( b + p )

{ e2 ( n ) }

= E

(39)

{ d 2 ( n ) } Ð 2 H NT ( n ) E { d ( n ) X N ( n ) } T (n) E [X (n) HT (n)] H (n) + HN N N N

2.5 Filtre de Wiener discret

T (n) r T J = ε ( n ) = E [ d 2 ( n ) ] Ð 2 HN dx + H N ( n ) R . H N ( n )

2.5.1 Filtre de Wiener à réponse finie où ■ Nous allons ici rechercher un filtre numérique de réponse impulsionnelle H nT (figure 6) : T = [ h ( 1 ), ....., h ( N ) ] HN

(40)

r dx = E [ d ( n ) X N ( n ) ]

vecteur 1 x N d’intercorrélation

R = E

vecteur N x N d’autocorrélation

{ X N ( n ). X NT ( n ) }

J est une forme quadratique définie de H N . ■ Le minimum de J sera obtenu en calculant : ∂J ∂J ∂J ------- = --------- , ..., ---------∂ h1 ∂ hN ∂H

où les h (i ) représentent les composantes de la réponse impulsionnelle du filtre. Ces composantes seront ajustées par le filtre de Wiener. Soit par ailleurs N échantillons du signal d’entrée notés :

x (n), x (n − 1), ..... x (n − (N − 1)) qui sont concaténés dans le vecT (n) : teur X N x ( n Ð 1 )..... x ( n Ð N + 1 ) ]

Le filtre de Wiener nous donnera une estimation du signal notée yö ( n ) qui a pour expression : T (n) X (n) ■ yö ( n ) = H N N

H Tn

{ ( d ( n ) Ð H NT ( n ) X N ( n ) ) 2 }

qui donne après développement :

ε (n) = E

x (n )

(38)

■ L’erreur quadratique moyenne, EQM, notée ε ( n ) ou J ( n ) est :

1 1 = ------------------------ -----------------k (1 + b) (1 + p)

T (n) = [x (n) XN

(37)

La sortie du filtre étant :

ϕ ms 1 1 1  1  1Ðp d’où : -------------Ð - =  ------------------- --------------------------- = ---------------------------- ---------------- + --------------ϕ mm 1 Ð ω2 k ( b Ð p ) k ( 1 + b ) 1 + p b Ð p et :

FILTRAGE OPTIMAL

où

T

= 0

h1, ....., hN sont les composantes du filtre H ,

et : 2 RH∗ Ð 2 r dx = 0

H∗ est le filtre optimal : H∗ = R Ð1 r dx (solution optimale au sens de Wiener)

(41)

Cette équation, qui donne les coefficients du filtre qui minimisent l’EQM, est l’équation de Wiener-Hopf ou équation normale. Elle est l’équivalent, dans le cas discret, de l’équation (33) de Wiener-Hopf que nous avons déjà établie dans le cas continu. On peut calculer l’erreur minimale en remplaçant H∗ par son expression (41) dans l’équation (40) soit :

y (n )

Figure 6 – Filtre numérique de réponse impulsionnelle H nT

J min- = E J min = E

T r { d 2 ( k ) } + ( R Ð1 r dx ) T R . R Ð1 r dx Ð 2 R Ð1 r dx dx T R Ð1 r T Ð1 r { d 2 ( k ) } + r dx dx Ð 2 r dx R dx

d’où :

T H∗ ε min = J min = E { d K2 } Ð r xd

d (n ) ou y (n ) x (n )

Filtre de Wiener

y (n )

+ –

Nota : 1°) L’erreur ε ( n ) donnée par l’équation (39) peut être mise sous la forme suivante :

ε EQM = ε min + ( H Ð H∗ ) T R ( H Ð H∗ )

e (n ) On peut écrire cette relation :

ö ( n ) du filtre de Wiener, comparée à la sortie Figure 7 – Sortie y désirée y (n)

ε = ε min + V T RV V apparaît comme l’écart par rapport au filtre optimal de Wiener.

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite. © Techniques de l’Ingénieur, traité Mesures et Contrôle

R 7 228 − 9

FILTRAGE OPTIMAL _____________________________________________________________________________________________________________________

2°) Reprenons la relation (38) : e ( n ) = d ( n ) Ð X NT ( n ) H et multiplions les deux membres par X N :

Si l’on applique la transformée de Fourier discrète à l’équation (46), le filtre optimal devient :

e ( n ) X N ( n ) = d ( n ) X N ( n ) Ð X N ( n ) X NT ( n ) H E { e ( n ) X N ( n ) } = r dx Ð RH

ϕ dx ( ω ) h opt ( ω ) = -------------------ϕ xx ( ω )

(42)

où :

On constate que si H prend sa valeur optimale R Ð1 r dx la relation (42) devient :

∞

E { e ( n ) XN ( n ) } = O

ϕ dx ( ω ) =

La résolution de l’équation de Wiener-Hopf nécessite l’inversion de la matrice R. Les algorithmes classiques de type Gauss-Jordan ont une complexité arithmétique de l’ordre de N 3. De plus, pour des applications temps réel, le calcul de R Ð1 doit être effectué chaque fois que de nouveaux échantillons du signal deviennent disponibles. Il faut remarquer que la structure particulière de la matrice d’autocorrélation R (à savoir le fait qu’elle soit symétrique et de Tœplitz) permet de réduire le nombre d’opérations à N 2 par l’utilisation de l’algorithme de Levinson. En outre, l’estimation de R et r nécessite la prise en compte d’un grand nombre de réalisations, ce qui n’est pas très réaliste dans la pratique. Nous développerons, dans le paragraphe 4, une méthode de résolution récursive de l’équation de Wiener qui ne nécessite pas la connaissance a priori de la matrice d’autocorrélation R.

2.5.2 Filtre de Wiener à réponse impulsionnelle infinie

ϕx x ( ω ) =

x (n Ð k)

(43)

On cherche le filtre h opt ( n ) à réponse impulsionnelle infinie (il a des valeurs non nulles pour toutes les valeurs entières de n) qui minimise l’erreur quadratique moyenne J (n) en comparant la sortie du filtre y ( n ) avec la sortie désirée d ( n ) :

{ e2 ( n ) }

= E

2.6 Application du filtre de Wiener non causal discret pour le débruitage de la parole

Soit le signal bruité :

x (n) = s (n) + b (n)

(48)

Le filtre de Wiener non causal (47) est donné par la relation :

ϕ sx ( ω ) h opt ( ω ) = ------------------ϕ xx ( ω ) Si l’on suppose que le signal et le bruit ne sont pas corrélés, on obtient successivement :

{ [ d ( n ) Ð y ( n ) ]2 }

∞

 = E  d ( n ) Ð ∑ h opt ( k )  k = Ð∞

r xx ( n ) e Ðj ωn

sont respectivement les densités interspectrale et spectrale de puissance.

k = Ð∞

J (n) = E

∑

k = Ð∞

où s ( n ) représente le signal de parole et b ( n ) le bruit additif.

∞

h opt ( k )

r dx ( n ) e Ðj ωn

2.6.1 Nouvelle expression du filtre

Nous recherchons un filtre numérique de réponse impulsionnelle h opt ( n ) de longueur infinie. Soit x ( n ) le signal d’entrée, alors la sortie y ( n ) du filtre devient :

∑

∑

k = Ð∞ ∞

Elle signifie que l’erreur est non corrélée avec le signal X N ( n ) , on dit aussi qu’elle est orthogonale au signal.

y (n) =

(47)

ϕ xx ( ω ) = ϕ ss ( ω ) + ϕ bb ( ω ) 2

x (n Ð k)

  

(44)

et :

ϕ sx ( ω ) = ϕ ss ( ω )

On notera :

Le filtre optimal de Wiener devient :

r xx ( k ) = E { x ( n ) x ( n Ð k ) }   r xd ( k ) = E { x ( n ) d ( n Ð k ) }   r dd ( k ) = E { d ( n ) d ( n Ð k ) } 

(45)

ϕ ss ( ω ) h opt ( ω ) = ----------------------------------------------ϕ ss ( ω ) + ϕ bb ( ω )

(49)

Cette dernière relation est analogue à la relation (35) que nous avons établie précédemment pour le filtre de Wiener dans le cas continu.

alors la relation (44) devient :

On peut encore écrire : ∞

J ( n ) = r dd ( 0 ) Ð 2

∑

ϕ xx ( ω ) Ð ϕ bb ( ω ) ϕ bb ( ω ) - = 1 Ð -------------------h opt ( ω ) = ----------------------------------------------ϕ xx ( ω ) ϕ xx ( ω )

h opt ( k ) r dx ( k )

k = Ð∞ ∞

+

∑

∞

∑

h opt ( k ) h opt ( , ) r xx ( , Ð k )

k = Ð∞ , = Ð∞

La minimisation de J (n) conduit à l’équation : ∞

∑

h opt ( k ) r xx ( , Ð k ) = r dx ( l )

∀,

(46)

k = Ð∞

qui constitue l’équation de Wiener-Hopf dans le cas de filtres non causaux.

R 7 228 − 10

(50)

2.6.2 Mise en œuvre du filtre La mise en œuvre du filtre nécessite : — que l’on segmente le signal x (n) en blocs à traiter successivement ; — que l’on connaisse ou que l’on estime les densités spectrales du bruit et du signal. Dans la mesure où l’on ne peut pas traiter tous les échantillons simultanément, on segmentera le signal en blocs de N échantillons

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite. © Techniques de l’Ingénieur, traité Mesures et Contrôle

____________________________________________________________________________________________________________________

(N = 256). Ces blocs peuvent être disjoints ou se recouvrir partiellement. On adoptera dans notre exemple un recouvrement de 50 %, c’est-à-dire que le début de chaque bloc va coïncider avec l’échantillon 127 du bloc précédent. Ainsi, pour un nombre total d’échantillons égal à M.N, on aura à traiter 2M − 1 blocs. En outre, en tenant compte de la spécificité du signal de parole, il est possible d’estimer la densité spectrale du bruit. En effet, le signal de parole est une succession de segments qui correspondent à la présence ou à l’absence de signal. Par ailleurs, dans la mesure où l’enregistrement commence par un silence, on peut, en faisant l’hypothèse que seul le bruit est présent dans ce cas, avoir une estimation de sa densité spectrale de puissance. Il existe plusieurs estimateurs de densité spectrale ; nous adopterons ici celui basé sur la méthode du périodogramme qui conduit à l’estimateur suivant :

FILTRAGE OPTIMAL

Amplitude Signal bruité

6 000 4 000 2 000 0 – 2 000 – 4 000 – 6 000

1 ϕ xx , i ( ω ) = ---- X i ( ω ) 2 N où

Xi ( ω )

0

Nous adopterons pour la densité spectrale de puissance du bruit la moyenne des densités spectrales de puissance calculées pour les premiers blocs où il y a seulement du bruit, que nous noterons ϕ bb ( ω ) . Cette procédure est appelée procédure de Bartlett :

ϕ bb

Xi ( ω ) 2

ϕ bb ( ω ) ( ω ) = 1 Ð -----------------------ϕ xx , i ( ω )

2 500

Amplitude

Signal débruité

2 000 0 – 2 000 – 4 000

Ensuite, on filtre le signal bruité :

Sö i ( ω ) = H opt, i ( ω )

2 000

4 000

i=1

Ainsi, dans les blocs où le signal de parole est présent, l’estimation de la densité spectrale nous donne ϕ xx . Dans le cas du bloc i, on aura alors l’expression suivante du filtre de Wiener :

h opt, i

1 500

Figure 9 – Signal de parole bruité

6 000

2M Ð 1

∑

1 000

Échantillons

est la transformée de Fourier du bloc i

1 ( ω ) = -----------------------------N (2M Ð 1)

500

Xi ( ω )

– 6 000

où X i ( ω ) est la transformée de Fourier du signal bruité pour le bloc courant i. La séquence temporelle est la partie réelle de la transforö i ( ω ) . Elle nous donnera le mée de Fourier inverse du spectre S signal de parole non bruité. Dans la figure 8, on présente la densité spectrale de puissance du bruit capté dans une voiture R25 roulant sur autoroute à 130 km/h (signaux aimablement fournis par la société Matra).

0

500

1 000

1 500

2 000

2 500

Échantillons Figure 10 – Signal débruité par le filtre de Wiener

Amplitude 6 000

Signal original

4 000 Amplitude [dB] 0

2 000

– 10

0

– 20

– 2 000

– 30 – 4 000 – 40 – 6 000 – 50 – 60 0

500

1 000

1 500 2 000

2 500 3 000

3 500 4 000

0

500

1 000

1 500

2 000 2 500 Échantillons On remarque bien que le signal de parole est une succession de segments pseudo-périodiques et aléatoires correspondant respectivement aux sons voisés et non voisés

Fréquence [Hz] Figure 8 – Spectre du bruit

Figure 11 – Signal de parole original

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite. © Techniques de l’Ingénieur, traité Mesures et Contrôle

R 7 228 − 11

FILTRAGE OPTIMAL _____________________________________________________________________________________________________________________

m (t ) = b (t ) + s (t )

h

coefficients sont déterminés par les propriétés statistiques des processus.

s (t )

ö (t) Figure 12 – Filtre de Wiener fournissant une estimation s du signal recherché

3.2 Algorithme du filtre de Kalman

2.6.3 Résultats expérimentaux Dans les figures 9, 10 et 11, on donne successivement des enregistrements du signal bruité, du signal obtenu par filtrage de Wiener et, enfin, le signal original dans une application sur un signal de parole capté dans une voiture R25 roulant sur autoroute à 130 km/h.

L’approche de Kalman a consisté à transformer l’équation intégrale (3.1), décrivant le filtre de Wiener dans le cas continu, en équation différentielle. Nous nous limiterons ici à un exposé élémentaire basé sur une méthode algébrique. On adoptera les notations suivantes : les grandeurs scalaires, vectorielles ou matricielles seront respectivement représentées par des lettres minuscules, des lettres soulignées ou des lettres en majuscule.

En comparant les signaux des figures 10 et 11, on constate que le niveau du bruit a été réduit dans la figure 10.

3.2.1 Formulation du problème Soit un signal x ( k ) représenté dans l’espace d’état par les équations suivantes :

3. Filtre de Kalman

x ( k + 1 ) = Φ ( k + 1, k ) x ( k ) + G ( k ) u ( k )

3.1 Propriétés du filtre de Kalman

y (k) = H (k) x (k) + v (k)

La forme la plus générale de l’équation de Wiener-Hopf peut s’écrire :

u ( k ) , processus générateur, et v ( k ) bruit de mesure sont supposés indépendants, blancs et de moyenne nulle :

t

∫h

( t , τ ) k mm ( τ , σ ) d τ = k ms ( t , σ )

(51)

t0

avec t 0 < σ < t ,

E

que nous écrivons plus simplement :

E

t

∫

h ( t , τ ) k m ( τ , σ ) d τ = k ms ( t , σ )

Le filtre de Wiener permet de construire, à partir d’un message observé m (t), une estimation sö ( t ) du signal recherché (figure 12). L’estimation de la sortie est : t

∫h

(t, τ) m (t, τ) dτ

(52)

t0

L’équation (51) est, d’une manière générale, difficile à résoudre dès que l’on aborde des cas non triviaux. Elle est, en outre, peu adaptée au traitement par calculateur à cause de son caractère non récursif. On se propose de développer ici un filtre qui permette d’estimer de manière récursive le signal noyé dans le bruit. Il s’agit du filtre de Kalman, plus approprié au traitement par calculateur. Kalman a introduit ce filtre, dès 1960, à partir de la représentation des systèmes dans l’espace d’état par des équations différentielles matricielles du premier ordre. Ce filtre est basé sur le fait qu’un processus aléatoire peut être modélisé comme étant la sortie d’un système linéaire gouverné par un bruit blanc, alors que, dans le cas du filtre de Wiener, les systèmes sont représentés par des équations de covariance. Au lieu de décrire les systèmes linéaires qui génèrent les messages en termes de réponse impulsionnelle, l’approche de Kalman amène une description par des équations différentielles dont la solution est le signal recherché. Par ailleurs, au lieu de spécifier la solution optimale comme sortie d’un système linéaire gouverné par une équation intégrale, comme dans le cas de l’approche de Wiener, l’estimation optimale donnée par Kalman est ici solution d’une équation différentielle dont les

R 7 228 − 12

E

E

{u (k)}

= 0

E

{v (k)}

= 0

{ u ( k ) vT ( l ) } = 0

∀ k, l

{ u ( k ) uT ( l ) } = Q ( k ) δ ( k , l ) et E

t0

sö ( t ) =

(53)

{x

( 0 ) uT ( k ) } = 0

{ v ( k ) vT ( l ) } = R ( k ) δ ( k , l )

et E

{x

∀k

∀k ( 0 ) vT ( k ) } = 0

(54) ∀k

Q ( k ) et R ( k ) sont respectivement les matrices de covariance du processus générateur et du bruit de mesure, δ ( h , l ) est le symbole de Kronecker. ■ Estimation du vecteur d’état Il s’agit d’estimer le vecteur d’état x ( k ) compte tenu des informations disponibles à l’instant n, postérieur, antérieur ou identique à l’instant k. On considérera les trois cas suivants : 1) k = n : il s’agit, dans ce cas, de déterminer une estimation de l’état, compte tenu de toutes les mesures disponibles à l’instant considéré n. C’est le cas du filtrage ; 2) k < n : on ne prendra en compte ici qu’une partie des mesures disponibles. On fait alors un lissage ou une interpolation ; 3) k > n : il s’agit, dans ce cas, de prédire une estimation du vecteur d’état. C’est ce qu’on appelle prédiction ou extrapolation. Nous noterons pour ces différents cas une telle estimation par : xö ( k ⁄ n ) , c’est-à-dire l’estimation à l’instant k compte tenu des informations disponibles à l’instant n. On définit aussi une telle estimation par la valeur moyenne du vecteur d’état compte tenu des mesures y ( 1 ) , ... y ( n ) . Cette moyenne est exprimée par l’espérance mathématique conditionnelle, notée comme suit :

xö ( k ⁄ n ) = E

{x (k) ⁄ y

( 1 ), ..., y ( n ) }

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite. © Techniques de l’Ingénieur, traité Mesures et Contrôle

(55)

____________________________________________________________________________________________________________________

Le filtre de Kalman est donc un filtre à minimum de variance.

Nota : on peut aussi, de manière plus commode, noter xö ( k ⁄ n ) = E [ x ( k ) ⁄ y ( n ) ] , cette notation condensée est équivalente à celle explicitée par la relation (55).

La matrice de transition Φ ( k , k Ð 1 ) caractérise l’évolution du vecteur d’état. Intuitivement, on peut admettre que l’évolution de l’estimation du vecteur d’état xö ( k ⁄ n ) sera, elle aussi, caractérisée par cette même matrice de transition. En d’autres termes, nous aurons :

xö ( k ⁄ k Ð 1 ) = Φ ( k , k Ð 1 ) xö ( k Ð 1 ⁄ k Ð 1 )

(56)

Nous devrions exprimer le critère J ( k ) en fonction du gain. Il est plus aisé d’opérer sur l’expression de la matrice de covariance de l’erreur P ( k ⁄ k ) plutôt que sur l’expression du critère J ( k ) . ■ Calcul de la matrice de covariance P ( k ⁄ k ) en fonction du gain K (k) Nous allons, à partir des équations d’état et de l’expression de l’estimation, expliciter la matrice de covariance P ( k ⁄ k ) en fonction du gain K ( k ) . Reprenons les relations (53) et (56).

Cette relation peut être obtenue de manière rigoureuse en appliquant la définition donnée par la relation (55) à l’équation (53) soit :

E

{x (k) ⁄ y

x (k) = Φ (k, k Ð 1)

( 1 ), ..., y ( k Ð 1 ) } =

Φ ( k ⁄ k Ð 1 ) E [ x ( k Ð 1 ) ⁄ y ( 1 ) , ..., y ( k Ð 1 ) ]

{u (k Ð 1) ⁄ y

+ G (k Ð 1) u (k Ð 1)

P (k ⁄ k Ð 1) = E

sachant que :

Notre objectif est d’avoir une estimation récursive du vecteur d’état qui, à partir d’une estimation à l’instant k nous fournisse, si une mesure est disponible à l’instant k + 1, une nouvelle estimation, compte tenu de cette dernière mesure. Nous adopterons un estimateur linéaire de la forme :

P (k ⁄ k Ð 1) = Φ (k, k Ð 1) E

(57)

{ ( x ( k Ð 1 ) Ð xö ( k Ð 1 ⁄ k Ð 1 )

( x ( k Ð 1 ) Ð xö ( k Ð 1 ⁄ k Ð 1 ) ) T } Φ T ( k , k Ð 1 ) + G ( k Ð 1 ) E

La quantité H ( k ) xö ( k ⁄ k Ð 1 ) que nous noterons yö ( k ) apparaît ainsi comme une prédiction de la mesure. Elle serait parfaite s’il n’y avait pas de bruit.

{ u ( k Ð 1 ) uT ( k Ð 1 ) } GT ( k Ð 1 )

Cette relation peut se simplifier en tenant compte de la relation (54) et du fait que :

E

{ ( x ( k Ð 1 ) Ð xö ( k Ð 1 ⁄ k Ð 1 ) ) ( x ( k Ð 1 ) Ð xö ( k Ð 1 ⁄ k Ð 1 ) ) T } ∆ = P (k Ð 1 ⁄ k Ð 1)

Nota : on note, dans certains exposés sur le filtre de Kalman, les quantités xö ( k ⁄ k ) et xö ( k ⁄ k Ð 1 ) comme suit :

Il vient alors :

xö ( k ⁄ k Ð 1 ) = xö ( k Ð 1 )

P ( k ⁄ k Ð 1 ) = Φ ( k , k Ð 1 ) P ( k Ð 1 ⁄ k Ð 1 ) ΦT ( k , k Ð 1 )

On devra s’attendre à ce que K ( k ) ne dépende pas de la mesure pour que cette linéarité soit assurée. La relation (57) exprime le fait que la nouvelle estimation du vecteur d’état, à l’instant k, est égale à l’estimation à l’instant k Ð 1 , mise à jour avec un certain poids. Cette mise à jour tient compte de l’écart entre la mesure effective et la mesure prédite. L’optimalité du filtre recherché vient du fait que le poids K (k ) à accorder à cette mise à jour résulte de la minimisation de l’erreur au sens des moindres carrés, entre l’état et son estimation. Le gain K (k ) ainsi obtenu est appelé gain du filtre de Kalman.

+ G ( k Ð 1 ) Q ( k Ð 1 ) GT ( k Ð 1 ) Nous devons, à ce stade du calcul, exprimer P ( k ⁄ k ) en fonction de P ( k ⁄ k Ð 1 ) . À cet effet, reprenons l’expression (58) de l’erreur x÷ ( k ) et remplaçons l’estimée xö ( k ⁄ k ) par son expression donnée par la relation (57) :

x÷ ( k ) = x ( k ) Ð xö ( k ⁄ k ) x÷ ( k ⁄ k ) = x÷ ( k ⁄ k Ð 1 ) + K ( k ) [ y ( k ) Ð H ( k ) xö ( k ⁄ k Ð 1 ) ]

(58) soit :

Il s’agira donc de minimiser la quantité :

J (k) = E

[ ( x ( k ) Ð xö ( k ⁄ k Ð 1 ) ) ( x ( k ) Ð xö ( k ⁄ k Ð 1 ) ) T ]

D’où, en prenant les espérances mathématiques des deux membres de l’équation (60) :

E [ u ( k ) ⁄ y ( 1 ) , ..., y ( k Ð 1 ) ] = 0

Soit x÷ ( k ) = x ( k ) Ð xö ( k ⁄ k ) l’erreur d’estimation.

x÷ ( k ⁄ k ) = x ( k ) Ð xö ( k ⁄ k Ð 1 ) Ð K ( k ) [ H ( k ) x ( k ) + v ( k ) ]

{ x÷ T ( k ⁄ k ) xö ( k ⁄ k ) }

Ð H ( k ) xö ( k ⁄ k Ð 1 ) ] = [ I Ð K ( k ) H ( k ) ] [ x ( k ) Ð xö ( k ⁄ k Ð 1 ) ] Ð K ( k ) v ( k )

Si nous introduisons la matrice

P (k ⁄ k) = E

(61)

{ [ x ( k ) Ð xö ( k ⁄ k ) ] [ x ( k ) Ð xö ( k ⁄ k ) ] T } = E { x÷ ( k ) x÷ T ( k ) }

Or : (59)

appelée matrice de covariance de l’erreur, on peut vérifier que l’on a: trace P ( k ⁄ k ) = J ( k ) .

(60)

D’après la relation (59) définissant P ( k ⁄ k ) , nous aurons :

xö ( k ⁄ k Ð 1 ) = Φ ( k , k Ð 1 ) xö ( k Ð 1 ⁄ k Ð 1 )

et

xö ( k Ð 1 ⁄ k Ð 1 )

x ( k ) Ð xö ( k ⁄ k Ð 1 ) = Φ ( k , k Ð 1 ) [ x ( k Ð 1 ) Ð xö ( k Ð 1 ⁄ k Ð 1 ) ]

( 1 ), ..., y ( k Ð 1 ) }

xö ( k ⁄ k ) = xö ( k ⁄ k Ð 1 ) + K ( k ) [ y ( k ) Ð H ( k ) xö ( k ⁄ k Ð 1 ) ]

u (k Ð 1)

En les soustrayant, il vient :

soit encore :

xö ( k ⁄ k ) = xö ( k )

x (k Ð 1) + G (k Ð 1)

xö ( k ⁄ k Ð 1 ) = Φ ( k , k Ð 1 )

+G (k Ð 1) E

FILTRAGE OPTIMAL

P (k ⁄ k) ∆ = E

{ x÷ ( k ⁄ k ) x÷ T ( k ⁄ k ) }

(62)

Si on substitue l’expression (61) dans la relation (62), il vient, en tenant compte de la relation matricielle ( AB ) T = B T . A T :

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite. © Techniques de l’Ingénieur, traité Mesures et Contrôle

R 7 228 − 13

FILTRAGE OPTIMAL _____________________________________________________________________________________________________________________

P (k ⁄ k) =

Les quantités P . P Ð1 et R . R Ð1 étant égales à la matrice unité, on peut les insérer dans l’expression du gain, soit :

[ I Ð K ( k ) H ( k ) ] E [ ( x ( k ) Ð xö ( k ⁄ k Ð 1 ) ) ( x ( k ) Ð xö ( k ⁄ k Ð 1 ) ) T ] . [ I Ð K ( k ) . H ( k ) ] T + K ( k ) E { v ( k ). v T ( k ) } K T ( k ) soit :

P ( k ⁄ k ) = [ I Ð K ( k ) . H ( k ) ] . P ( k ⁄ k Ð 1 ) . [ I Ð K ( k ) . H ( k ) ]T + K (k) R (k)

KT

(k)

(63)

K ( k ) = P ( k ⁄ k ) P Ð1 ( k ⁄ k ) P ( k ⁄ k Ð 1 ) H T ( k ) R ( k ) R Ð1 ( k )

[ ( H ( k ) P ( k ⁄ k Ð 1 ) H T ( k ) + R ( k ) ) ] Ð1 = P ( k ⁄ k ) P Ð1 ( k ⁄ k ) P ( k ⁄ k Ð 1 ) H T ( k ) R Ð1 [ ( H ( k ) P ( k ⁄ k Ð 1 ) H T ( k ) + R Ð1 ( k ) + I ) ] Ð1 Nous allons remplacer dans cette dernière relation P Ð1 ( k ⁄ k ) par sa valeur donnée par la relation (67), il vient :

K ( k ) = P ( k ⁄ k ) [ P Ð1 ( k ⁄ k Ð 1 ) + H T ( k ) R Ð1 ( k ) H ( k ) ]

■ Expression du gain du filtre de Kalman K (k) Rappelons que nous devons chercher l’expression du gain K (k) qui minimise le critère :

J = E

P ( k ⁄ k Ð 1 ). H T ( k ) R Ð1 ( k ). [ ( H ( k ) P ( k ⁄ k Ð 1 ) H T ( k ) R Ð1 ( k ) + I ) ] Ð1 = P ( k ⁄ k ) [ I + H T ( k ) R Ð1 ( k ) H ( k ) P ( k ⁄ k Ð 1 ). H T ( k ) R Ð1 ( k ) ]

{ x÷ T ( k ⁄ k ) x÷ ( k ⁄ k ) } =

Notons au passage que nous avons : trace P ( k ⁄ k ) = E

{ x÷ T ( k ⁄ k ) x÷ ( k ⁄ k ) } = J

[ [ ( H ( k ) P ( k ⁄ k Ð 1 ) H T ( k ) R Ð1 ( k ) + I ) ] Ð1 ] P ( k ⁄ k ) H T ( k ) R Ð1 ( k ) [ I + H ( k ) P ( k ⁄ k Ð 1 ) H T ( k ) R Ð1 ( k ) ] [ ( H ( k ) P ( k ⁄ k Ð 1 ) H T ( k ) R Ð1 ( k ) + I ) ] Ð1

(64)

La relation (63) est une forme quadratique en K (k). Nous obtiendrons la valeur du gain optimal en résolvant l’équation suivante : ∂ J (k) -------------------- = 0 ∂ K (k) ∂ et en tenant compte de la relation --------- ( trace KAK T ) = 2 KA ∂ K

soit, en définitive :

K ( k ) = P ( k ⁄ k ) H T ( k ) R Ð1 ( k ) La quantité P (k/k) traduit la confiance que nous pouvons avoir dans le modèle adopté.

R traduit le niveau de bruit sur la mesure. Nous pouvons, en outre, écrire, dans le cas scalaire par exemple, l’équation récursive de mise à jour de l’estimation comme suit :

soit : Ð [ I Ð K ( k ) H ( k ) ] P ( k ⁄ k Ð 1 ) HT ( k ) + K ( k ) R ( k ) = 0

xö ( k ) = xö ( k Ð 1 ) + K ( k ) [Terme de correction]

d’où :

Il apparaît ainsi que :

K ( k ) = P ( k ⁄ k Ð 1 ) H T ( k ) [ H ( k ) P ( k ⁄ k Ð 1 ) H T ( k ) + R ( k ) ] Ð1 (65) Dans le cas continu, cette dernière équation prend la forme d’une équation différentielle de Riccati. C’est pour cette raison, et par analogie, que la relation (65) est appelée ici équation de Riccati discrète. Elle peut être réécrite sous la forme suivante :

P (k ⁄ k) = P (k ⁄ k Ð 1) Ð K (k) H (k) P (k ⁄ k Ð 1) =

[I Ð K (k) H (k)] P (k ⁄ k Ð 1)

(66)

La matrice K ( k ) H ( k ) P ( k ⁄ k Ð 1 ) est définie positive (cf. nota ci-dessous) ou, à la rigueur, semi-définie positive. Il s’ensuit que tous ses éléments diagonaux seront soit positifs, soit nuls. Par conséquent, on pourra affirmer que : trace P ( k ⁄ k ) < trace P ( k ⁄ k Ð 1 )

1) à R (k) constant, si P (k) est faible, le gain sera faible. Ainsi, on fera davantage confiance à l’estimation obtenue à partir du modèle. Si P (k) est élevé, ce qui traduit notre faible confiance dans le modèle, K sera élevé et la contribution du terme de correction pondéré par le gain sera plus forte ; 2) à P (k) constant, si R (k) est faible, nous aurons donc des mesures faiblement bruitées. La valeur élevée du gain donnera plus de poids à la mesure. Si R (k) est élevé, le gain sera faible et le poids du deuxième terme sera plus faible.

3.3 Initialisation et mise en œuvre du filtre de Kalman ■ Conditions initiales x (0) et P (0).

et donc le critère J (k) ira en décroissant, ce qui garantit la convergence du filtre. Nota : soit x un vecteur et A une matrice, on dit que la matrice A est définie positive si la forme quadratique x T A x est strictement positive. Dans le cas contraire, on dit qu’elle est définie négative. Une matrice telle que x T A x > 0 est dite semi-définie positive.

Pour pouvoir utiliser l’ensemble des équations récurrentes constituant le filtre de Kalman, on doit choisir les conditions initiales, x ( 0 ) et P ( 0 ⁄ 0 ) :

xö ( 0 ) = E

■ Autre expression du gain K (k)

K ( k ) = P ( k ⁄ k Ð 1 ) H T ( k ) [ H ( k ) P ( k ⁄ k Ð 1 ) H T ( k ) R ( k ) ] Ð1 Si P et R sont inversibles, on aura en appliquant le lemme d’inversion matricielle à la relation (66) :

R 7 228 − 14

(0)}

car, au démarrage du filtre, la seule information que l’on pourrait espérer avoir au mieux est la valeur moyenne du vecteur d’état,

Reprenons l’expression du gain donnée par la relation (65) :

P Ð1 ( k ⁄ k ) = P Ð1 ( k ⁄ k Ð 1 ) + H T ( k ) R Ð1 H ( k )

{x

(67)

et E

{[x

( 0 ) Ð xö ( 0 ) ] [ x ( 0 ) Ð xö ( 0 ) ] T } = P ( 0 ⁄ 0 ) = P ( 0 )

■ Nous résumons ci-après l’ensemble des équations du filtre. Équation du modèle x ( k + 1 ) = Φ ( k + 1 , k ) x ( k ) + G ( k ) u ( k ) Équation d’observation y ( k ) = H ( k ) x ( k ) + v ( k )

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite. © Techniques de l’Ingénieur, traité Mesures et Contrôle

____________________________________________________________________________________________________________________

Information a priori E

{u (k)}

= 0

E

{v (j)}

E

{u (k)

E

{ v ( k ) vT ( j ) } = δ ( k , j ) R ( k ) E

uT

= 0

(j)} = δ (k, j) Q (k)

v (k ) +

u (k ) G

{ u ( k ) v T ( j ) } = 0 ∀k , j

+

Z –1

+

+

x (k ) H

y (k ) +

+

+

Φ (k,k –1)

+

K

–

x (k /k) +

Z –1 Φ (k,k –1)

H (k )

Équation du filtre

x (k /k –1)

xö ( k ⁄ k ) = xö ( k ⁄ k Ð 1 ) + K ( k ) [ y ( k ) Ð H ( k ) xö ( k ⁄ k Ð 1 ) ]

(68)

ou encore xö ( k ⁄ k ) = Φ ( k ⁄ k Ð 1 ) xö ( k Ð 1 ⁄ k Ð 1 ) + K ( k ) [ y ( k ) Ð H ( k ) xö ( k ⁄ k Ð 1 ) ] Expression du gain et

FILTRAGE OPTIMAL

K ( k ) = P ( k ⁄ k Ð 1 ) H T ( k ) [ H ( k ) P ( k ⁄ k Ð 1 ) H T ( k ) + R ( k ) ] Ð1   K ( k ) = P ( k ⁄ k ) H T ( k ) R Ð1 ( k )  (69)

Variance a posteriori P ( k ⁄ k ) = [ I Ð K ( k ) H ( k ) ] P ( k ⁄ k Ð 1 )

Figure 14 – Schéma bloc du filtre de Kalman

Le choix des valeurs initiales est délicat. En effet, un mauvais choix de x (0), c’est-à-dire à la limite une valeur arbitraire, n’est pas catastrophique en ce sens que l’algorithme excité par les mesures apportera les corrections nécessaires. Par contre, le traitement des mesures n’améliore pas la covariance de l’erreur au fur et à mesure de son traitement. Nous pouvons traduire notre ignorance de P (0) en adoptant P (0) = a I, a ayant une valeur élevée (100 à 10 000). ■ Le déroulement des différentes opérations peut être décomposé comme suit : — pour k = 1, on calcule successivement :

Variance a priori

P ( k ⁄ k Ð 1 ) = Φ ( k , k Ð 1 ) P ( k Ð 1 ⁄ k Ð 1 ) ΦT ( k , k Ð 1 ) + G ( k Ð 1 ) Q ( k Ð 1 ) GT ( k Ð 1 )

x (k–1 /k–1)

P (1/0) à partir de la relation (71), en tenant compte de P (0) (70)

K (1) à partir de la relation (69)

Remarques générales :

X (1) par la relation (68)

Si nous prenons les expressions des matrices de covariance d’erreur P ( k ⁄ k Ð 1 ) , appelées aussi erreur a priori et P ( k ⁄ k ) erreur a posteriori d’une part, et du gain K (k) d’autre part :

P (1/1) — et on continue en faisant k = 2, ..., comme suit : k = 1 : P (1/0) ; K (1) ; x (1) ; P (1/1)

P ( k ⁄ k Ð 1 ) = Φ ( k , k Ð 1 ) P ( k Ð 1 ⁄ k Ð 1 ) ΦT ( k , k Ð 1 ) + G ( k Ð 1 ) Q ( k Ð 1 ) GT ( k Ð 1 )

k = 2 : P (2/1) ; K (2) ; x (2) ; P (2/2)

(71)

K ( k ) = P ( k ⁄ k Ð 1 ) H T ( k ) [ H ( k ) P ( k ⁄ k Ð 1 ) H T ( k ) + R ( k ) ] Ð1 P (k ⁄ k) = [I Ð K (k) H (k)] P (k ⁄ k Ð 1)

(72)

il apparaît que l’on peut déterminer ces trois quantités indépendamment de toute mesure et avant même que celle-ci ne soit traitée par le filtre. En effet, les paramètres qui interviennent dans ces trois expressions dépendent soit des paramètres du modèle, soit de la connaissance a priori des statistiques du processus générateur et du bruit de mesure (ces informations font aussi partie du modèle). ■ Ces quantités peuvent être déterminées dans l’ordre où elles sont écrites. On a besoin pour ce faire de la connaissance de P (0/0) que nous noterons P (0).

k = 3 : etc. Ceci permet de représenter ces opérations par l’organigramme de la figure 13. ■ Le schéma bloc du filtre, quant à lui, est représenté par la figure 14.

Exemple : nous allons, sur un exemple scalaire, illustrer le déroulement des différentes opérations de mise en œuvre du filtre de Kalman. Au préalable, reprenons les équations (69), (70), (72) qui se réduisent dans ce cas à :

P (k ⁄ k Ð 1) = P (k Ð 1 ⁄ k Ð 1) + Q K ( k ) = [ Φ 2 P ( k Ð 1 ⁄ k Ð 1 ) + Q ] [ Φ 2 P ( k Ð 1 ⁄ k Ð 1 ) + Q + R ] Ð1 et

Calcul a priori

P ( k ⁄ k ) = R [ Φ 2 P ( k Ð 1 ⁄ k Ð 1 ) + Q ] [ Φ 2 P ( k Ð 1 ⁄ k Ð 1 ) + Q + R ] Ð1

Injection de P (0)

= R.K (k)

P (k /k –1)

Considérons le modèle scalaire avec les valeurs numériques suivantes :

K (k ) k=k+1 Prise en compte des mesures

x (k ) Mise à jour

P (k /k) Calcul a posteriori Figure 13 – Organigramme du filtre de Kalman

Φ = 1 ; G = 1 ; P ( 0 ) = 10 ; Q = 20 ; R = 10 et H = 1. Les trois équations précédentes deviennent :

P ( k ⁄ k Ð 1 ) = P ( k Ð 1 ⁄ k Ð 1 ) + 20 K ( k ) = [ P ( k Ð 1 ⁄ k Ð 1 ) + 20 ] [ P ( k Ð 1 ⁄ k Ð 1 ) + 30 ] Ð1 et P ( k ⁄ k ) = 10 K ( k )

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite. © Techniques de l’Ingénieur, traité Mesures et Contrôle

R 7 228 − 15

FILTRAGE OPTIMAL _____________________________________________________________________________________________________________________

Le tableau suivant résume l’évolution des différents paramètres :

k 0 1 2 3

P (k ⁄ k Ð 1) Ð 30 27,5 27,4

K (k) Ð 0,75 0,74 0,74

P (k ⁄ k) 10 7,5 7,4 7,4

où

u (k)

représente le signal d’excitation qui permet de générer s (k)

Soit le vecteur d’état constitué par les p dernières valeurs du signal s (k) :

x ( k ) = [ s ( k Ð p + 1 ), ..., s ( k ) ] T

On constate que pratiquement, au bout de la troisième itération, k = 3, P (k/k) et le gain n’évoluent plus et l’on peut écrire :

(74)

et l’observation y (k). On peut représenter ce système dans l’espace d’état :

P (k ⁄ k) = P (k Ð 1 ⁄ k Ð 1) = P (k)

x ( k + 1 ) = Φx ( k ) + Gu ( k )

(75)

en combinant alors les deux équations donnant P (k/k) et K (k), on voit que P (k) obéira alors à l’équation du deuxième degré :

y ( k ) = Hx ( k ) + b ( k )

(76)

P k2 + 20 P k Ð 200 = 0

où :

.

Si cette prédiction est parfaite, la correction apportée par la quantité :

.

est la matrice de transition. Les matrices G de commande et H d’observation sont telles que :

xö ( k ⁄ k ) = xö ( k ⁄ k Ð 1 ) + K ( k ) [ y ( k ) Ð H ( k ) xö ( k ⁄ k Ð 1 ) ] y (k) représente la mesure effective et H (k) xö ( k ⁄ k Ð 1 ) la prédiction de la mesure.

.

■ Notion d’innovation

.

Φ=

xö ( k + 1 ⁄ k + 1 ) = xö ( k ⁄ k ) + 0,74 [ y ( k + 1 ) Ð xö ( k ⁄ k ) ] Reprenons l’équation de l’estimateur récursif dans le cas scalaire :

0 1 . . . . 0 . . . . . . . . . . . . . . . . 0 0 1 – ap – ap –1 . . . . – a1 .

qui admet pour solution Pk = 7,4. On peut donc considérer R = 3 comme permettant d’atteindre le régime permanent. L’équation de l’estimateur devient :

G = H T = [ 0, 0, ..., 0, 1 ] T u (k) est le processus générateur et b (k) le bruit additif contaminant le signal. On suppose que u (k) et b (k) sont indépendants, blancs, gaussiens, de moyennes nulles et de variances respectives σ u2 et σ b2 :

ν ( k ) = y ( k ) Ð H ( k ) xö ( k ⁄ k Ð 1 ) E

{u (k)}

= 0

E

{b (k)}

= 0

est nulle. On appelle ν ( k ) l’innovation. On peut montrer que si le filtre est optimal, le processus ν ( k ) est un bruit blanc de valeur moyenne nulle, c’est-à-dire qu’il ne contient plus d’information pouvant enrichir la mise à jour de l’estimation de l’état. Ainsi, en testant la « blancheur » de l’innovation, on peut apprécier le degré d’optimalité et donc les performances du filtre. Nota : les performances du filtre de Kalman peuvent se dégrader et le filtre risque de diverger si la dynamique du bruit − processus générateur − n’est pas gaussienne.

E

{u (k) b (l)}

= 0

E

{u (k) u (l)}

= σ u2 δ ( k , l ) et E

{x

(0) u (l)} = 0

E

{b (k) b (l)}

= σ b2 δ ( k , l ) et E

{x

(0) b (l)} = 0

Il s’agit d’estimer le vecteur d’état x ( k ) compte tenu des informations disponibles à l’instant k. La composante p du vecteur d’état est l’échantillon s (k) du signal de parole. Pour cela, on peut utiliser le filtre de Kalman [équations (68)-(70)] :

3.4 Filtrage de Kalman d’un signal de parole noyé dans un bruit blanc gaussien

xö ( k ⁄ k Ð 1 ) = Φxö ( k Ð 1 ⁄ k Ð 1 ) Nous allons, dans ce qui suit, formaliser l’utilisation du filtre de Kalman en tant que filtre qui permet d’augmenter le rapport signal à bruit. Nous présenterons ensuite son application au débruitage du signal de parole dans le cas du radiotéléphone.

3.4.1 Mise en œuvre du filtre

p

∑ ai s ( k Ð i ) + u ( k ) i=1

R 7 228 − 16

K ( k ) = P ( k ⁄ k Ð 1 ) H T ( k ) [ H ( k ) P ( k ⁄ k Ð 1 ) H T ( k ) + σ b2 ] Ð1 P (k ⁄ k) = [I Ð K (k) H (k)] P (k ⁄ k Ð 1)

On considère le signal de parole s (k) modélisé par un processus autorégressif (AR) d’ordre p :

s (k) = Ð

xö ( k ⁄ k ) = xö ( k ⁄ k Ð 1 ) + K ( k ) [ y ( k ) Ð H ( k ) xö ( k ⁄ k Ð 1 ) ]

(73)

P ( k ⁄ k Ð 1 ) = Φ P ( k Ð 1 ⁄ k Ð 1 ) Φ T + GG T σ u2 Dans le paragraphe suivant, on va estimer les paramètres du modèle autorégressif du signal utilisé pour la mise en œuvre du filtre de Kalman.

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite. © Techniques de l’Ingénieur, traité Mesures et Contrôle

____________________________________________________________________________________________________________________

FILTRAGE OPTIMAL

3.4.2 Estimation de paramètres Amplitude

On va écrire la fonction d’autocorrélation de y (k), que nous noterons ryy, sachant que le signal et le bruit ne sont pas corrélés et que le bruit est un bruit blanc :

r yy ( l ) = E + E

{y (k) y (k Ð l)}

{s (k) s (k Ð l)}

= E

{ b ( k ) b ( k Ð l ) } = r ss ( l ) + r bb ( l ) = r ss ( l ) + δ ( l )

Signal original

6 000 4 000 2 000

σ b2

(77)

La fonction d’autocorrélation de l’observation y est égale à la fonction d’autocorrélation du signal s pour les valeurs de l non nulles.

0 – 2 000 – 4 000

Reprenons l’équation (73) : – 6 000

p

∑ ai s ( k Ð i ) + u ( k )

s (k) = Ð

0

500

1 000

1 500

2 000 2 500 Échantillons

i=1

Multiplions-la par s (k − l ) et prenons l’espérance mathématique ; il vient :

E

{s (k) s (k Ð l)}

Figure 15 – Signal de parole original

= Amplitude

p

  E  ∑ Ð ai s ( ( k Ð i ) s ( k Ð l ) ) + u ( k ) s ( k Ð l )  i = 1 

(78)

Signal bruité

6 000 4 000

où :

E E E

{u (k) s (k Ð l)} {s (k) s (k Ð l)}

{s (k Ð i) s (k Ð l)}

2 000

= δ ( l ) σ u2

0

= r ss ( l )

(79)

= r ss ( l Ð i )

– 4 000

on peut donc réécrire l’expression (78) :

– 6 000

p

r ss ( l ) = Ð

– 2 000

∑ ai rss ( l Ð i ) + δ ( l ) σu2

0

500

1 000

1 500

i=1

On peut voir que le deuxième terme dans la partie droite est non nul uniquement pour l = 0. On peut donc écrire l’équation (79) pour p + 1 < l < 2 p en utilisant l’équation (77), soit :

2 000 2 500 Échantillons

Figure 16 – Signal de parole bruité

p

∑ ai ryy ( l Ð i )

r yy ( l ) = Ð

Amplitude

i=1

Pour calculer les paramètres ai, on écrira :

ryy (1) . . . . ryy (p) . . . . . . . . ryy (p) . . . . ryy (2p –1)

–1

.

.

.

.

.

– ap . . =– . . – a1

Signal débruité

6 000

ryy (p +1) . . . . ryy (2p)

Si l’on dispose de N échantillons du signal y (n), on peut estimer sa fonction d’autocorrélation par :

4 000 2 000 0 – 2 000 – 4 000

r yy

1 ( l ) = ---N

N

∑

y (i) y (i Ð l)

i = l+1

Rappelons, pour les lecteurs qui ne sont pas familiers avec les signaux de parole, que l’on considère que le signal de parole est une succession de segments de signaux quasi périodiques (voyelles) et aléatoires (consonnes).

– 6 000 0

500

1 000

1 500

2 000

2 500

Échantillons

Figure 17 – Signal de parole débruité par filtrage de Kalman

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite. © Techniques de l’Ingénieur, traité Mesures et Contrôle

R 7 228 − 17

FILTRAGE OPTIMAL _____________________________________________________________________________________________________________________

L’estimation de la variance du bruit de mesure sera effectuée pendant les périodes de silence à partir de la relation : 1 σ b2 = ---N

d (n ) X (n )

N

HN

∑ b2 ( i )

–

Filtre adaptatif

i=1

e (n ) e (n ) = d (n) – y (n)

et celle de la variance du bruit générateur par la relation :

σ u2 = r yy ( 0 ) Ð σ b2

+

y (n )

Figure 18 – Filtrage adaptatif

3.4.3 Résultats expérimentaux

4.2 Algorithme du type gradient stochastique ou (LMS)

On considère l’expérimentation suivante où le signal de parole est le signal utile. Ce signal est perturbé par un bruit blanc. Cette expérimentation peut couvrir le cas de la réduction du bruit dans l’habitacle d’une voiture pour le radiotéléphone. Il s’agit, à défaut d’annuler l’effet du bruit, tout au moins de le réduire.

4.2.1 Position du problème

Ainsi, on présente successivement une trame de signal original (figure 15), de signal bruité (figure 16) et, enfin, le signal débruité (figure 17). Nota : on peut considérer le problème de l’estimation d’un signal noyé dans un bruit comme la résolution d’un problème d’estimation conjointe des paramètres du modèle et de l’état. Dans ce cas, on est conduit à la résolution d’un problème de filtrage non linéaire en utilisant le filtre de Kalman étendu (EKF).

Nous allons reprendre les notations que nous avons utilisées au § 2.5. Soit H N le vecteur des paramètres du filtre à estimer : T = [ h ( 1 ), ....., h ( N ) ] HN T ( n ) le vecteur des échantillons du signal d’entrée et soit X N

On peut aussi aborder la résolution du problème en estimant les variances Q et R. L’approche de Mehra sur l’estimation de Q et R a été étendue à l’extraction de signaux noyés dans du bruit blanc ou coloré.

T ( n ) = [ x ( n ), x ( n Ð 1 ), ....., x ( n Ð N + 1 ) ] XN

qui donnent le signal de sortie :

4. Introduction au filtrage adaptatif 4.1 Différents algorithmes On appellera filtre adaptatif un filtre numérique dont les coefficients évoluent en fonction des signaux reçus. Ces coefficients seront estimés par des algorithmes récursifs, au sens d’un certain critère. Les critères qui sont généralement retenus pour l’obtention de ces algorithmes sont du type moindres carrés car ils conduisent aux résultats les plus simples à interpréter. Les filtres peuvent être à réponse impulsionnelle finie (RIF) ou infinie (RII), quant à la structure, elle est soit transverse, soit en treillis. Les familles d’algorithmes qui en découlent ont des complexités arithmétiques différentes et leur comportement dépend du type d’excitation et de l’absence ou de la présence de bruit. On considèrera deux grandes familles d’algorithmes basés sur le gradient stochastique ou LMS (Least Mean Squares) et sur les moindres carrés récursifs type MCR (ou RLS). Dans cette dernière famille, la recherche d’algorithmes à très faible complexité a conduit à développer des algorithmes dits rapides, appelés FTF (Fast Transversal Filters). Le développement d’algorithmes rapides de type FTF qui peuvent opérer en temps réel résulte de la réduction des redondances. Cette réduction de la redondance « fragilise » en quelque sorte ces algorithmes sous l’effet conjugué des arrondis et des troncatures quand ces algorithmes doivent opérer sur des processeurs à virgule fixe. Il en résulte ainsi une instabilité que l’on a essayé de circonscrire durant les dernières années en développant des versions dites stabilisées.

R 7 228 − 18

T (n) H (n) HN N

On représente aussi le filtrage adaptatif par le schéma de la figure 18, que nous avons présenté dans le paragraphe 2.51, où la sortie du filtre adaptatif yö ( n ) est comparée à une sortie désirée d (n) qui joue le rôle de la sortie réelle que l’on cherche à approximer.

4.2.2 Algorithme LMS Nous allons présenter ici une méthode itérative, celle dite du gradient par la méthode de la plus grande pente (steepest descent) qui effectue la recherche du minimum de la surface représentant l’erreur quadratique moyenne (EQM) en suivant la direction et le sens opposé du gradient. Cette approche permet l’obtention d’une solution récursive de l’équation de Wiener dans le cas discret [relation (41)] :

H∗ = R Ð1 r dx La méthode du gradient permet de construire cette solution récursive par la relation suivante : α H ( n ) = H ( n Ð 1 ) Ð --- grad ε EQM 2 où :

ε EQM = E

{ d ( n Ð 1 ) Ð H NT ( n Ð 1 ) X N ( n Ð 1 ) } 2

et α est un paramètre dit gain ou pas, qui permet de contrôler la convergence de l’algorithme du gradient. Soit :

H (n) = H (n Ð 1) + α E

{X (n Ð 1) e (n Ð 1)}

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite. © Techniques de l’Ingénieur, traité Mesures et Contrôle

(80)

____________________________________________________________________________________________________________________

E

Cet algorithme est d’un usage limité car il est difficile de connaître { X ( n Ð 1 ) e ( n Ð 1 ) } que l’on va donc remplacer :

E

{X (n Ð 1) e (n Ð 1)}

= X (n Ð 1) e (n Ð 1)

L’algorithme (80) devient :

H (n) = H (n Ð 1) + α X (n Ð 1) e (n Ð 1)

Car, même dans ce dernier cas, c’est-à-dire dans le cas où :

E

(81)

La valeur moyenne est ainsi remplacée par sa valeur instantanée.

{ X ( n ) XT ( m ) }

E

{ X ( n ) XT ( n Ð 1 ) }

=

où λ max est la plus grande valeur propre de la matrice d’autocorrélation du signal d’entrée. Dans le cas d’une séquence blanche de puissance σ x2 , cette condition devient : 2 α < ----------Nσ x2

n≠m

0

0

σ b2

0

0

σ b2

...

0

0

0

σ b2

0

On voit ainsi que, même avec une séquence blanche, on ne peut pas satisfaire « l’hypothèse d’indépendance ». Ceci montre l’aspect un peu « irréaliste » de cette hypothèse. Reprenons l’équation de mise à jour des coefficients du filtre sous la forme : T X (n)] X (n) H ( n + 1 ) = H ( n ) + α [ d ( n ) Ð HN

qui peut être réécrite, en notant que T X ( n ) = ( H T X ( n ) )T = XT ( n ) H HN : N N

On montrera aussi que la vitesse de convergence, proportionnelle au pas, est inversement proportionnelle à la dispersion des valeurs propres et qu’elle est indépendante des conditions initiales. On peut démontrer qu’il existe une valeur de α qui maximalise la vitesse moyenne de convergence du LMS, cette valeur est : 2 α = -------------------------------λ min + λ max

H ( n + 1 ) = ( I Ð α X ( n ) X T ( n ) ) H ( n ) + αd ( n ) X ( n ) et, en prenant les valeurs moyennes des deux membres :

E [ H ( n + 1 ) ] = ( I Ð α R ) E ( H ( n ) ) + α . R . H∗ en notant R = E [ X ( n ) X (n) .

Dans le cas d’une séquence blanche gaussienne, le pas maximalisant la vitesse de convergence est :

α opt

n=m

= 0

0

On va montrer que la condition qui assure la convergence est : 0 < α < 2 ⁄ λ max

= σ b2

La matrice d’autocorrélation définie par la relation (81) donne par exemple pour m = n Ð 1 :

(82)

Ce sont les différentes itérations successives pour effectuer la mise à jour de H ( n ) qui assurent la réalisation de la moyenne et la convergence en moyenne de l’algorithme.

FILTRAGE OPTIMAL

XT

(83)

( n ) ] la matrice d’autocorrélation de

Définissons le vecteur erreur U ( n ) , différence entre la valeur moyenne de H ( n ) et la valeur optimale H∗ soit :

1 = ----------Nσ x2

U ( n ) = E ( H ( n ) ) Ð H∗ Si nous retranchons H∗ de chacun des deux membres de l’équa-

■ Erreur quadratique résiduelle À l’inverse du gradient déterministe, le LMS à pas constant ne permet pas d’atteindre le minimum de la surface d’EQM, à cause du bruit introduit par les fluctuations du gradient estimé. Le vecteur H ( n ) des paramètres suit des fluctuations aléatoires vis-à-vis de la trajectoire obtenue pour l’utilisation du gradient déterministe. Ces fluctuations sont à moyenne nulle et leur variance, proportionnelle au pas, reste bornée. ■ Convergence en moyenne du LMS

tion (83), il vient :

E ( H ( n + 1 ) ) Ð H∗ = ( I Ð αR ) E ( H ( n ) ) Ð ( I Ð αR ) H∗ ainsi U ( n + 1 ) = E ( H ( n + 1 ) ) Ð H∗ d′où U ( n + 1 ) = ( I Ð αR ) E ( H ( n ) Ð H∗ ) soit U ( n + 1 ) = ( I Ð αR ) U ( n )

(84)

Reprenons la relation (81) :

H (n) = H (n Ð 1) + α e (n Ð 1) X (n Ð 1) La démonstration de la convergence est basée sur l’hypothèse d’indépendance de H ( n ) et de X ( n ) . En fait, compte tenu de la nature récursive de l’équation de mise à jour des coefficients du filtre, H ( n ) dépend en fait, non seulement de X ( n Ð 1 ) , mais aussi de X ( n Ð 2 ) , X ( n Ð 3 ) , ... D’où le caractère contraignant de cette hypothèse qui nécessite que les vecteurs d’entrée { X ( n ) } soient non corrélés, c’est-à-dire :

E [ X ( n ) XT ( m ) ] = 0

n≠m

C’est une condition encore plus forte que celle qui nécessite que l’entrée soit une séquence blanche.

Cette dernière équation est intéressante dans le sens où elle nous fournit une expression récursive d’évolution de l’erreur d’estimation des paramètres du filtre. Ainsi, démontrer que le LMS converge en moyenne vers la solution optimale H∗ se ramène à démontrer que l’erreur U ( n ) converge vers zéro. Pour ce faire, nous allons effectuer un découplage des équations de mise à jour en utilisant la décomposition suivante que nous établirons en annexe (en fin d’article) :

R = QΛQ T

(85)

où les colonnes de la matrice Q sont les vecteurs propres de R :

Q = [ Q 1 -------- Q N ]

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite. © Techniques de l’Ingénieur, traité Mesures et Contrôle

R 7 228 − 19

FILTRAGE OPTIMAL _____________________________________________________________________________________________________________________

et Λ la matrice des valeurs propres notée comme suit :

Λ =

λ1

... 0

0

λ2

or λ i > 0 d’où trace (R) > λ max

2 2 ------------ > -------------N λ max ∑ λi

avec Q T Q = I où I est la matrice identité.

i=1

λN

0

soit :

On peut donc remplacer la condition (88) par la condition plus restrictive :

Reprenons l’équation (84) et substituons-y l’équation (85) :

2 0 < α < -------------N

U ( n + 1 ) = ( I Ð α Q Λ QT ) U ( n )

∑ λi

(86)

L’analyse de la convergence va être effectuée en tirant avantage de la relation (85). Pour ce faire, on va considérer un vecteur U ′ ( n ) obtenu par la transformation linéaire suivante :

i=1

Par ailleurs, R étant une matrice de Tœplitz, on pourra écrire : N

∑ λi = trace ( R ) = N r ( 0 ) = N E [ x 2 ( n ) ] i=1

U ′ ( n ) = QT U ( n ) d’où

Ainsi, si nous prémultiplions les deux membres de l’équation (86) par la matrice Q T et si nous utilisons la relation Q T Q = I , il vient :

2 2 0 < α < ----------------------------------- = ------------N . E ( x2 ( n ) ) N σ x2

Q T U ( n + 1 ) = [ Q T Ð αQ T QΛ ] U ( n )

2 0 < α < ---------------------------------------------------------------------------------------------------N ( puissance du signal d ′entrée )

et, en utilisant la propriété Q T Q = I , il vient :

Il s’agit d’une condition certes plus restrictive sur la stabilité du filtre, mais plus facile à mettre en œuvre.

Q T U ( n + 1 ) = [ I Ð αΛ ] Q T U ( n ) , soit : U ′ ( n + 1 ) = [ I Ð αΛ ] U ′ ( n )

(87)

et, compte tenu du fait que Λ est une matrice diagonale, on pourra écrire pour chaque composante i du vecteur U ′ ( n + 1 ) , les relations suivantes :

U i′ ( n + 1 ) = ( 1 Ð αλ i ) U i′ ( n )

■ Dispersion des valeurs propres Pour les processus stationnaires, pourvu qu’il y ait indépendance, l’algorithme LMS converge en moyenne vers la valeur optimale. Mais cette convergence ne se fait pas de manière uniforme. Ceci est un problème majeur pour le LMS. Reprenons l’équation (87), soit :

U i′ ( n + 1 ) = ( 1 Ð αλ i ) U i′ ( n )

Le caractère récursif de cette relation permet d’écrire :

U i′ ( n ) = ( 1 Ð αλ i ) n U i′ ( 0 ) Ainsi, l’algorithme va converger si la limite de l’erreur U i′ ( n ) tend vers zéro quand n devient infini, soit :

(89)

Si la condition 0 < α < 2 ⁄ λ max est satisfaite, tous les coefficients du vecteur erreur vont décroître en amplitude. Néanmoins, le taux de décroissance de chaque élément − ou mode − dépend de la quantité 1 Ð αλ i . Ceci explique que certains modes vont converger plus rapidement que d’autres. Ainsi, le LMS converge de manière non uniforme vers la valeur optimale.

lim n → ∞ ( U i′ ( n ) ) → 0 ceci a lieu si 1 Ð αλ i < 1 Ainsi, le vecteur U ′ converge vers zéro et le filtre converge vers H∗ .

Ce phénomène de convergence non uniforme est dû à la disparité des valeurs propres.

Les valeurs propres λ i sont, comme nous le verrons dans l’annexe, réelles et positives, car la matrice R est symétrique définie positive donc :

■ Constante de temps de convergence

Ð 1 < 1 Ð αλ i < 1

U i′ ( n ) = ( 1 Ð αλ i ) n U i′ ( 0 )

soit :

Notons l’instant considéré n = t i . La constante de temps se défi-

2 2 α > 0 ; α < ---- ; 0 < α < ---λj λj

j = 0.... L Ð 1

Cette condition doit être satisfaite pour toutes les valeurs de λ , ainsi la borne supérieure sera : 2 0 < α < ------------λ max où

λ max

Reprenons la relation :

nit par le temps nécessaire pour que l’amplitude soit réduite à U i′ ( 0 ) ----------------- , c’est-à-dire : e

U i′ ( 0 ) U i′ ( t i ) = ( 1 Ð αλ i ) ti U i′ ( 0 ) = ---------------e

(88)

est la plus grande valeur propre de la matrice R

Soit, en prenant les logarithmes des deux derniers termes de la relation (90) :

Nota : considérons la trace de la matrice R :

t i ln ( 1 Ð αλ i ) = Ð 1

N

trace ( R ) =

∑ λi i=1

R 7 228 − 20

(90)

où

ln

désigne le logarithme népérien

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite. © Techniques de l’Ingénieur, traité Mesures et Contrôle

____________________________________________________________________________________________________________________

m(n) = s (n) + b0 (n)

+ + –

e (n )

Amplitude 6 000

y (n )

b1 (n )

FILTRAGE OPTIMAL

4 000

H 2 000

Figure 19 – Schéma bloc du système annulateur de bruit

0 – 2 000

Considérons que λ i ,, 1 , alors :

– 4 000

ln ( 1 Ð αλ i ) ≈ αλ i – 6 000 0

1 d’où t i = --------αλ i C’est une quantité qui dépendra des valeurs propres. En fait, t i sera « imposée » par la plus faible parmi les valeurs propres.

4.3 Application à l’annulation du bruit

500

1 000

1 500

2 000

2 500

Échantillons Figure 20 – Signal de parole non bruité

Amplitude 6 000

Considérons le message reçu :

m ( n ) = s ( n ) + b 0 ( n ) et supposons que l’on dispose d’une source de bruit b 1 ( n ) corrélé à b 0 ( n ) (figure 19). Nous allons montrer, dans ce cas, que l’utilisation d’un filtre adaptatif permettra, à partir du bruit mesuré b 1 ( n ) , de disposer d’une estimation du bruit b 0 ( n ) . Cette estimation est notée y (n). Ainsi, pour obtenir une estimation du signal s (n), il suffira de soustraire l’estimation du bruit b 1 ( n ) du message m (n). Soit e (n) = s (n) + b 0 ( n ) − y (n) un signal d’erreur.

E

{ e 2 ( n ) } = E { s 2 ( n ) } + E { ( b0 ( n ) Ð y ( n ) ) 2 } +2 E

{ s ( n ) ( b0 ( n ) Ð y ( n ) ) }

Nous supposerons que s (n) n’est corrélé ni avec b 0 ( n ) ni avec b 1 ( n ) . Le signal s (n) n’est pas corrélé non plus avec y (n). Ainsi, l’expression de l’erreur E ( e 2 ( n ) ) devient :

E

4 000 2 000 0 – 2 000 – 4 000 – 6 000

0

500

1 000

1 500

2 000

2 500

Échantillons Figure 21 – Signal de parole bruité (rapport signal à bruit égal à 7,35 dB)

{ e 2 ( n ) } = E { s 2 ( n ) } + E { ( b0 ( n ) Ð y ( n ) ) 2 }

On voit que la valeur minimale de E ( e 2 ( n ) ) est obtenue quand :

Amplitude 6 000

y ( n ) = b0 ( n ) alors e ( n ) = s ( n ) est le résultat désiré. ■ Conditions expérimentales de mise en œuvre des simulations. Pour satisfaire la condition sur la corrélation des bruits b 1 ( n ) et b 0 ( n ) , le bruit b 0 ( n ) sera généré par le processus suivant :

b 0 ( n ) = 0,5 b ( n Ð 1 ) + 0,35 b ( n Ð 2 ) Ð 0,3 b ( n Ð 3 ) Ð 0,2 b ( n Ð 4 ) + 0,1 b ( n Ð 5 ) Ð 0,2 b ( n Ð 6 ) + 0,1 b ( n Ð 7 ) Ð 0,1 b ( n Ð 8 ) + 0,1 b ( n Ð 9 ) On prendra N = 10 le nombre d’échantillons du signal d’entrée et α = 10 Ð8 le gain de l’algorithme défini dans la relation (82). Les résultats obtenus avec cet annuleur de bruit sont consignés sur les trois figures ci-après : — un signal de parole non bruité, figure 20; — un signal de parole bruité, figure 21, avec un rapport signal à bruit de 7,35 dB ; — le signal de parole débruité, avec un rapport signal à bruit de 17,85 dB (figure 22).

4 000 2 000 0 – 2 000 – 4 000 – 6 000

0

500

1 000

1 500

2 000

2 500

Échantillons Figure 22 – Signal de parole débruité par filtrage adaptatif (rapport signal à bruit égal à 17,85 dB)

Nota : 1°) Les propriétés de convergence du LMS se dégradent très rapidement dans le cas d’entrées corrélées ; en plus de l’alternative d’utilisation des algorithmes de type moindres carrés, il existe plusieurs approches basées sur les transformations dans le domaine

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite. © Techniques de l’Ingénieur, traité Mesures et Contrôle

R 7 228 − 21

FILTRAGE OPTIMAL _____________________________________________________________________________________________________________________

fréquentiel. Ceci est motivé par l’existence de l’algorithme de la transformée de Fourier FFT qui conduit, en général, à une double amélioration : rapidité de convergence et diminution de la complexité. On peut aussi citer l’utilisation d’autres transformations, de type Hartley par exemple, ou opérant sur des décompositions sur des bases d’ondelettes.

1 Compte tenu de la symétrie de R, avec R = R T, les vecteurs propres qui correspondent à des valeurs propres distinctes sont orthogonaux :

2°) Filtres d’ordre adaptatifs Dans le cas de bruits impulsionnels, la combinaison de filtres non linéaires de type filtres d’ordre avec des filtres adaptatifs permet de réduire l’effet des bruits impulsionnels, aussi bien pour l’estimation de paramètres de modèles de signaux bruités que pour le rehaussement de signaux perturbés par des bruits impulsionnels. 3°) L’algorithme LMS n’était jusqu’à présent pas considéré comme un algorithme optimal. En effet, l’approximation faite dans la relation 4.4 en prenant :

T Q = 0 pour toute paire de vecteurs m, n. Qm n

En effet, soit λ 1 et λ 2 deux valeurs propres distinctes :

RQ 1 = λ 1 Q 1

E {X (n) e (n)} = X (n) e (n)

RQ 2 = λ 2 Q 2

fait que l’algorithme ne découle pas de l’utilisation d’un critère d’optimalité exact mais d’un critère approché. Ainsi, on le considère comme un algorithme sous-optimal. La raison pour laquelle nous l’avons introduit dans la famille des algorithmes optimaux résulte de notre interprétation des travaux récents (1996) menés par Th. Kailath et ses collègues qui ont permis de montrer que le LMS est optimal au sens du critère dit « H infini ».

(93)

En prenant les transposées et en postmultipliant (93) par Q 2 , il vient :

Q 1T R = λ 1 Q 1T

5. Annexe. Quelques propriétés de la matrice R

Q 1T RQ 2 = λ 1 Q 1T Q 2 En postmultipliant (A.5) par Q 1T , on aura :

Reprenons l’équation (cf. § 2.51) :

Q 1T RQ 2 = λ 2 Q 1T Q 2

ε = ε min + ( H Ð H∗ ) T R ( H Ð H∗ )

Or R T = R d’où :

= ( ε min + V TRV )

λ 1 Q 1T Q 2 = λ 2 Q 1T Q 2

Les performances de la surface représentant l’erreur dépendent de la matrice R, en particulier de ses valeurs et vecteurs propres.

Or λ 1 ≠ λ 2 d’où Q 1T Q 2 = 0

Soit λ i et Q i les valeurs propres et vecteurs propres de R :

RQ i = λ i Q i avec λ i scalaire Q i vecteur colonne

(91)

Ainsi, les vecteurs propres correspondant aux valeurs propres λ 1 et λ 2 sont orthogonaux.

RQ = QΛ

soit : ( R Ð λ i I ) Q i = O Cet ensemble d’équations homogènes n’a de solution non nulle que si nous avons : det ( R Ð λ i I ) = O Cette dernière relation dont les solutions seront notées λ 1 , .... λ N est appelée équation caractéristique. On peut alors écrire, compte tenu de (91), et en définissant la matrice Λ des valeurs propres :

λ1 0 0 Λ =

2 La matrice R est réelle, toutes ses valeurs propres doivent être réelles. On peut esquisser une démonstration par l’absurde. Supposons que λ 1 soit complexe. L’équation caractéristique de R est un polynôme de degré N. Si λ 1 est complexe, alors λ 1* , qui est la valeur conjuguée, doit être aussi une valeur propre.

λ 1 complexe entraîne que Q 1 complexe

0 λ2 .. .

λ 1* → Q 1

0 ... λ N

On ne peut plus alors vérifier la relation d’orthogonalité λ 1 Q 1T Q 2 = λ 2 Q 1T Q 2 . On en déduit l’impossibilité d’avoir une valeur

λ1 0 0 R [ Q 1 , Q 2 , ... Q N ] = [ Q 1 , Q 2 , ... Q N ]

0 λ2 .. .

0 .. .

propre complexe. On peut normaliser les vecteurs propres pour qu’ils soient normés à 1. Dans ce cas, la matrice Q est dite orthonormée :

0 ... λ N

qui peut être réécrite :

QQ T = I avec Q Ð1 = Q T

RQ = QΛ ou encore R = QΛQ Ð1 R est symétrique définie positive.

(92) La relation Q Ð1 = Q T signifie qu’il y aura toujours un inverse pour Q. Aussi la relation (92) devient :

On appelle cette dernière relation la forme normale de R.

Q est la matrice des vecteurs propres.

R 7 228 − 22

R = QΛQ T

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite. © Techniques de l’Ingénieur, traité Mesures et Contrôle

____________________________________________________________________________________________________________________

FILTRAGE OPTIMAL

Références bibliographiques [1]

WIENER (N.). – Extrapolation, interpolation and smoothing of stationary time series. MIT Press, 1949. Cambridge, Ma.

[2]

LEE (Y.W.). – Statistical theory of communication. J. Wiley, 1960.

[3]

VAN TREES (H.L.). – Detection, estimation and modulation theory. Tome 1, J. Wiley, 1968.

[4]

THOMAS (J.B.). – An introduction to statistical communication theory. J. Wiley, 1969.

[5]

PAPOULIS (A.). – Probability, random variables and stochastic processes. Mc Graw Hill, 1965.

[6]

KALMAN (R.E.), BUCY (R.S.). – New results in linear filtering and prediction theory. Trans. ASME, Series D, Journal of Basic Eng. vol-38, p. 95-101, 1960.

[7]

KALMAN (R.E.). – A new approach to linear filtering and prediction problems. ASME, Journal of Basic Eng. vol-82-D, p. 34-45, 1960.

[8]

NAJIM (M.). – Modélisation et identification en traitement du signal. Masson, Paris, 1988.

[9]

SORENSON (H.W.). – Kalman Filtering, Theory and Applications. IEEE Press, New York, 1985.

[10]

BELLANGER (M.). – Analyse du signal et filtrage numérique adaptatif. Masson, 1989.

[11]

GABREA (M.), MANDRIDAKE (E.), NAJIM (M.). – A single microphone noise canceller

based on adaptive Kalman Filter. Proceedings EUSIPCO’96, sept. 96, Trieste, Italie. [12]

[13]

[14]

DENTINO (M.), Mc COOL (J.) et WIDROW (B.) . – Adaptive filtering in the frequency domain. Proc. IEEE, vol. 66, pp. 1658-1659, déc. 1978. FERRARA (E.R.). – Fast implementation of LMS adaptive filters. IEEE Trans. on Acous. Speech and Signal Proc. vol. ASSP-28, pp. 474-475, april 1980. MANSOUR (P.), GRAY Jr (A.H.). – Unconstrained frequency domain adaptive filter. IEEE Trans. on Acous. Speech and Signal Proc, vol. ASSP 30, pp. 726-734, Oct. 1982.

[15]

FLORIAN (S.), BERSHAD (N.J.). – A weighted normalized frequency domain LMS adaptive algorithm. IEEE Trans. on ASSP, vol. ASSP36, pp. 1 002-1 007, july 1988.

[16]

ATTALLAH (S.), NAJIM (M.). – The influence of the regularity and the subband coding decomposition structure on the convergence behaviour of the normalized wavelet transform. IC ASSP 1996.

[17]

HAWEEL (T.I.), CLARKSON (P.M.). – A class of order statistic LMS algorithms. IEEE Trans. on S.P. vol SP-40 n° 1, p. 44-51, jan 1992. Cette approche généralise celle du filtre median LMS développée dans : P.M. CLARKSON and T.I. HAWEEL : A median LMS algorithms. Electron Letters, vol. 25, pp. 520-522, 1989.

[18]

OTTAVIANI (R.), SETTINERI (R.) et NAJIM (M.). – Order statistics fast Kalman filter. IEEE ISCAS 96 - mai 1996 - Atlanta.

[19]

NASSER ELDIN (S.) et NAJIM (M.). – Sample restoration based on nonlinear fast Kalman filtering. Actes de la Conférence Intern. Workshop on Theory of Sampling, Aveiro, June 16/19, 1997. Portugal.

[20]

MEHRA (R.K.). – On the identification of variances and adaptive Kalman filter. IEEE Trans. on Automatic Control vol. AC-15, n° 2, pp 175-84, April 1970.

[21]

MEHRA (R.K.). – On line identification of linear dynamic systems with applications to Kalman filtering. IEEE Trans. on Automatic Control vol. AC-16, n° 1, p. 12-21, February 1971.

[22]

HASSIBI (B.), SAYED (A.H.), KAILATH (Th). – H infinity optimality of the LMS algorithm. Trans. on Signal Proces. vol. 44, n° 2, pp. 267280, 1996.

[23]

WIDROW (B.), STEARNS (M.). – Adaptive digital signal processing. Prentice Hall, NJ, 1985.

[24]

THERRIEN (Ch. W.). – Discrete random signal and statistical signal processing. Prentice Hall, 1992.

[25]

MACCHI (O.). – Adaptive processing. The least means square approach with applications in transmission. J. Wiley 1994.

Toute reproduction sans autorisation du Centre français d’exploitation du droit de copie est strictement interdite. © Techniques de l’Ingénieur, traité Mesures et Contrôle

R 7 228 − 23