Try   HackMD

Auteur: Rémi Maubanc
Professeur: M. Zitouni

Optimisation & Filtrage

Introduction

L'object du cours est d'apprendre à synthétiser des filtres numériques, analogiques et adaptatifs (des filtres numériques qui ont des coeficients qui bougent automatiquement en fonction de la situation).
Role du filtrage:

  • Mise en forme du spectre: à partir d'un gabarit déterministe (passe-bas, passe-haut, rejecteur de bandes). Comme pour les bruits ou atténuations créés par le passage dans un canal de transmision
  • Indentification d'un système: Pour représenter un système inconnu à l'aide d'un filtre numérique linéaire à partir uniquement des signaux d'entrée et de sortie.
    Image Not Showing Possible Reasons
    • The image file may be corrupted
    • The server hosting the image is unavailable
    • The image path is incorrect
    • The image format is not supported
    Learn More →
  • Réduction du bruit d'un message: Pour réduire le bruit pour peut lisser/moyenner un signal pour éliminer le bruit. Lorsque le spectre du signal utile et celui du bruit sont superposables. La méthode classique du filtrage ne marche pas. La méthode utilisée est l'optimisation d'une fonction coût par la méthode des moindres carrés : Filtre de Wiener
    Image Not Showing Possible Reasons
    • The image file may be corrupted
    • The server hosting the image is unavailable
    • The image path is incorrect
    • The image format is not supported
    Learn More →

Filtres numériques

Rappel

  • Un filtre analogique:

    y(t)=x(t)h(t)=+x(t).h(tτ)dτ
    Y(jω)=H(jω).X(jω)
    (Produit simple)
    Image Not Showing Possible Reasons
    • The image file may be corrupted
    • The server hosting the image is unavailable
    • The image path is incorrect
    • The image format is not supported
    Learn More →

  • Un filtre numérique: c'est un algorithme qui effectue une opération de récurrence. Trouver un filtre numérique c'est trouver les coeficients de

    h(n).
    y(n)=x(n)h(t)
    (Convolution numérique)
    =k=0N1x(k)h(nk)

    Image Not Showing Possible Reasons
    • The image file may be corrupted
    • The server hosting the image is unavailable
    • The image path is incorrect
    • The image format is not supported
    Learn More →

    Te
    est le temps échantillonnage avec
    1Te=Fefm

Pour caractériser et calculer les filtres numériques, on doit utiliser la transformée en

Z.

Transformée en
Z

Enchantillonnage idéal

Le rôle de l'échantionnage est de remplacer un signal analogique donc continu par une succession de points.

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

Peigne de Dirac :

n=+δ(tnTe)
Pour échantillonner un signal, il faut le multiplier par un peigne de Dirac.
xe(t)=x(t).n=δ(tnTe)

xe(t)=n=x(nTe).δ(tnTe)

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

Effet de la Transformée de Fourier sur un signal échantillonné:

xe(t).
TE(xe(t))=TF(x(t).n=δ(tnTe))

TE(xe(t))=TF(x(t))TF(n=δ(tnTe))

TF(x(t))=X(f)1Ten=δ(fnTe)

Xe(f)=TF(xe(t))=X(f)1Ten=infinfδ(fnTe)
Xe(f)=1Ten=infinfX(fnTe)

La TF montre l'effet de l'échantillonnage sur le spectre d'un signal. C'est la périodisation du spectre. Regardons l'effet de la transformée de Laplace sur un signal échantillonné.

TL(xe(t))=+xe(t).eftdt
TL(xe(t))=+(n=x(nTe).δ(tnTe)).eftdt

TL(xe(t))=n=(x(nTe)+δ(tnTe).eftdt)

TL(xe(t))=n=x(nTe).+δ(tnTe).efnTedt

TL(xe(t))=n=x(nTe).efnTe.+δ(tnTe)dt

TL(xe(t))=n=x(nTe).efnTe

Sachant que
Te=1s
par conversion et
x(t)=0
pour
t<0
(signal causal)
TL(xe(t))=n=0xnepn

On pose
ep=z

TL(xe(t))=n=0xnzn=X(z)

Donc la transformée en

Z n'est rien d'autre que la Transformée de Laplace d'un signal échantillonné.

Correspondance entre la TZ et le plan de Laplace

z=epTe;
Te=1sec
par convention
p est un complexe donc
p=σ+jω

z=e(σ+jω)Te=eσTe.ejωTe

|z|=eσTe
;
arg(z)=ωTe

Rappel:

ejωTe=cos(ωTe)+jsin(ωTe)
|ejωTe|=cos2(ωTe)+sin2(ωTe)=1=1

Un système analogique est stable si et seulement si les parties réelles de ses pôles sont négatives.

Un système numérique est stable ssi tous ses poles se trouvent à l'intérieur du cercle unité.

z=eσTe.ejωTe;
Te=1sec

|z|=eσ=1
si
σ=0

|z|=eσ<1
si
σ<0

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

Correspondance entre la TZ et la TF

On sait que

X(z)=n=0+xnzn lorsque
σ=0
.
X(z)=n=0+xnejωn
car
z=eσ=ejω

X(z)=n=0+xne2πjnf

Par ailleurs, on sait que
TF(xe(t))=TF(n=0x(nTe)δ(tnTe))

TF(xe(t))=n=0x(nTe)TF(δ(tnTe))

TF(xe(t))=n=0x(nTe)e2πjnfTe
;
Te=1sec

X(z)=n=0xne2πjnf=Xe(f)=xne2πjnf

TZ=TF
lorsque
p=jω
σ=0

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

Propriétés de T
Z

δ(k)TZ1

u(k)={0ailleurs1sik0

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

TZ(u(k))=n=0zk=limn>01(z1)n1z1

C'est une suite géométrique de raison

q=z1
TZ(u(k))=1z1
si
|z1|<1
donc
|z|<1

δ(ki)>(TZ)zi

K.u(k)TZz(z1)2
K2.u(k)TZz(z+1)(z1)z

αK.u(k)TZzza=11z1a

  • Retard :
    x(ni)TZziX(z)
  • Convolution :
    x(k)y(k)TZX(z).Y(z)
  • Valeur initiale :
    x(0+)=limk0x(t)=limzX(z)
  • Valeur finale :
    x()=limk0x(k)=limz1z1zX(z)=x()

Transformée en
Z
comme fraction rationnelle

Le calcul d'un filtre permet d'obtenir comme fonction de transfert :

H(z)=b0+b1.z1+...+bM.zMa+a1.z1+...+aN.zN
Si le filtre est stable, on peut décomposer la fraction
H(z)
en élements simples :
H(z)=(1z1.z1)(1z2.z1)...(1zM.z1)(1p1z1).(1p2.z1)...(1pN.z1)

Comme dans les systèmes physiques réels
M<N
alors :
H(z)=A11p1.z1+A21p2.z1+...+AN1pN.z1

Le filtre

H(z) converge si chaque filtre élémentaire
Ai1pi.z1
converge.
Or
A1pi.z1=limn1(pi.z1)n1pi.z1|pi.z1|<1

Le filtre élementaire converge lorsque
|pi.z1|<1|z|>pi

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

Transformée en
Z
inverse par la méthode des Résidus

Ai1pi.z1=TZ(pin)x(n)=pin

TZ1
par la méthode des Résidus

x(n)=pôles intérieurs résidus zn1.X(z)

Exemple :

X(z)=z2(za)(z1)
2 pôles
2 résidus
x(n)=limza(za).zn1.X(z)+limz1(z1)zn1X(z)

x(n)=limza(za).zn1.z2(za)(z1)+limz1(z1).zn1.z2(za)(z1)

x(n)=an+1a1+11a=an+11a1

Mathlab

B=[0 0 1] %coeff du numérateur
A=[a a+1 1]
%coeff du dénominateur
h=impz(B,A)
%c'est la réponse impulsionnelle du filtre
zplane(B,A)
% tracé de la position des pôles et zéros
P=roots(A)

Z=roots(B)

Exercice

Calculer la

TZ1 de
X(z)=z2z23z+2
. En déduire
x(0)
,
x(1)
,
x(2)
et
x(3)
.
On utilise la méthode des résidus. On remarque que 1 et 2 sont des pôles :
X(z)=z2(z1)(z2)x(n)=limz1(z1)zn1(z2(z1)(z2))=1x(n)=limz2(z2)zn1(z2(z1)(z2))=2n+11x(n)=2n+11=x(0)=1;x(1)=3;x(2)=7;x(3)=15

Equations aux différences et transformation en
Z

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

  1. Le système discret est dit Linéaire lorque
    Rep(ax1(n)+bx2(n))=a.Rep(x1(n))+b.Rep(x2(n))
  2. Invariance : Un SLD est dit invariant ssi
    Rep(x(n))=y(n)Rep(x(na))=y(na)

    The same excitation produces the same results
  3. La stabilité : Un SLD est dit stable si sa réponse impulsionnelle est convergente
    n=0h(n)<
  4. Equation aux différences

Un filtre est souvent définit par une fraction rationnelle :

H(z)=Y(z)X(z)=b0+b1z1+b2z2+...+bMzM1+a1z1+a2z2+...+aNzNY(z)+a1z1Y(z)+...+aNzNY(z)=b0X(z)+b1z1X(z)+...+bMzMX(z)TZ1y(n)=i=0Mbix(ni)j=1Najy(nj)

Discussion

1er cas :
Si

aj=0
j
alors $y_{n} \sum_{i=0}^{M}b_{i}x(n-i). Nous avons alors un filtre FIR (Finite Impulse Response) et ce filtre effectue une moyenne glissante du signal d'entrée. C'est un filtre MA (Moving Average)

Les coeficients de la réponse impulsionnelle sont les

bi.
H(z)=b0+b1z1+b2z2+...+bMzM

2e cas :
Si
aj0
, deux situations sont possibles.
H(z)=b01+a1z1+a2z2+...+aNzNy(n)=b0x(n)y=1Majy(nj)

C'est un filtre AR (Auto-Regressif).
C'est aussi un filtre IIR (Infinite Impulsse Response) si :
H(z)=b0+b1z1+...+bMzM1+a1z1+...+aNzN

C'est un filtre ARMA (Auto-Regressive with Moving Average).
5. Stabilité (détail)
Quelque soit le filtre, on peut le décomposer en fractions simples.
H(z)=A11p1z1+A21p2z2+...+AN1pNzN

Le filtre est stable si chaque filtre d'ordre 1 qui le compose est stable or

Ai1piz1=limn1(piz1)n1piz1 avec
|piz1|<1|z|>|pi|

Et la réponse impulsionnelle est convergente ssi
n=0h(n)<n=0h(n)=n=0i=1N(pi)n<|pi|<1

Un SLD est stable si tous ses pôles sont inclus dans un cercle unité

Filtrage de Wiener, Filtrage adaptatif

Introduction

Généralement le signal enregistré est

x(t)=xw(t)+n(t) avec :

  • xw(t)
    : Wanted signal
  • n(t)
    : Noise signal

Lorsque le spectre du signal utile

xw(t) est supperposable au spectre du bruit
n(t)
, les filtres classiques (passe-bas, passe-haut, stop et rejecteur de bande) sont inopérant.
Il faut utiliser une méthode récursive d'approximation fondée sur la minimisation (de l'énergie) de l'erreur.

Position du Problème

Les signaux à analyser issus d'un capteur ou d'un tout autre dispositif électronique peuvent s'écrire comme suit :

x(t)=xu(t)+u(t)
xu(t)
est le signal utile et
u(t)
le bruit de mesure ou d'enregistrement.

  • Lorsque les spectres de
    xu(t)
    et de
    n(t)
    sont séparés (non-superposables), un filtre classique peut les séparer. Ce qu'on désigne par filtre classique, un système analogique ou numérique qui garde la partie spéctrale de
    x(t)
    qui correspond à celle de
    xu(t)
    et élimine celle de
    n(t)
    . On obtient alors des filtres passe-bas, passe-haut, passe-bande, rejecteur de bandes, etc
  • Lorsque les spectres sont superposables, on doit estimer l'un des deux signaux. Dans la méthode de Wiener, l'estimation se fait à l'aide d'un filtre linéaire en minimisant l'énergie de l'erreur d'estimation. Cette méthode est efficace pour les signaux stationnaires.
    Elle nécessite beaucoup de calcul et ne peut se faire en temps réel.
    C'est pourquoi, on lui prefère la technique du filtrage adaptatif : les coeficients du filtre linéaire d'estimation sont variables au cours du temps.

Exemple

Définition de quelques estimateurs statistiques

On considère l'enregistrement d'un signal sur

N valeurs (
N
)
moyenne=composante continue=limN1Nn=0N1xn=μx

Puissance moyenne

Px=μx2=limN1Nn=0N1xn2

La variance mesure la puissance moyenne des fluctuations autour de la composante continue :

σx=limN1Nn=0N1(xnμx)2σx=limN1Nn=0N1(xn2+μx22.xn.μx)σx=μx2+μx22.μx2σx=μx2μx2
La puissance des fluctuations = puissance totale - puissance de la composante continue.

Dans le cas où le nombre de variable

N n'est pas infiniment grand, on parle d'espérence mathématique :
x¯=1T(T)x(t)dtE(x)=Dx.p(x).dxE{x}=1Nn=0N1xnμx

Moyenne statistique : moyenne théoriquement espérée
E{x2}=1Nn=0N1xn2μx2

Autocorrélation statistique = puissance moyenne estimée

Ecriture vectorielle

x=[x0x1...xN1];XNT=[x0,x1,...,xN1]Px=μx2=1N.XNT.XN

Filtrage de Wiener classique

De manière générale le problème peut être schématique comme suit:

  • x(t)
    et
    y(t)
    sont des signaux connus et accessibles
  • $\widehat{y}_{v}(t) est une estimation de
    yv(t)
  • e(t)
    est inacessible : c'est le signal recherché
  • ϵ(t)
    est l'écart entre le signal mesuré
    y(t)
    et l'estimation
    y^(t)

Hypothèse de Wiever

Wiener admet l'hypothèse que le precessus inconnu peut être modélisé par un filtre linéraire à moyenne glissante.
Filtre FiR

H(z)=i=0P1ωizi
wi
sont les coefficients du filtre FIZ
P
: ordre du filtre
H(z)
: fonction de transfert
H(z)=Yv^(z)X(z)=Yv^(z)=H(z).X(z)Yv^(z)=i=0P1ωiziX(z)TZ1Yv^(n)=i=0P1ωixni

Yv^(n)
est une moyenne pondérée du signal d'entrée.
ωi
sont les coefficients de la pondération. Nous devons trouver les coefficients
wk
qui rendent l'estimée
Yv^(t)
la plus porche possible de
Yv
.
Pour cela, on doit minimiser l'erreur quadratique moyenne (MSE : Mean Square Error).
J(ϵ)=1Nn=0N1ϵn2=1Nn=0N1(YnYn^)2

Equation de Wiener-Hopf

J(ϵ)=1Nn=0N1ϵn2=1Nn=0N1(ynYn^)2J(ϵ)=1Nn=0N1(ynk=0P1wkxnk)2Fonction coût

Pour trouver les

{wk} qui minimisent la fonction coût de manière optimale, il faut que
J(ϵ)ωi=0
pour
i[0,P1]

J(ϵ)=1Nn=0N1(ynk=0P1wkxnk)2J(ϵ)ωi=1Nn=0N1(ynk=0P1wkxnk)(xni)=01Nn=0N1ynxni=1Nn=0N1k=0P1wk.xnk.xni1Nn=0N1yn.xni=k=0P1wk.1Nn=0N1xnk.xniRYX(i)=k=0P1wk.RXX(ik)

Avec:

  • RYX(i)
    : Intercorrelation numérique des signaux
    y
    et
    x
    .
  • RWW(ik)
    : Autocorrelation de
    x
    calculé au point
    ik

Ce calcul est à fair pour

i=0..P1
i=0RXY(0)=w0RXX(0)+w1RXX(1)+w2RXX(2)+...+wP1RXX(P+1)i=1RXY(1)=w0RXX(1)+w1RXX(0)+w2RXX(1)+...+wP1RXX(P+2)......i=P1RXY(P1)=w0RXX(P1)+...+wP1RXX(0)

On sait que

RXX(i)=RXX(i) car l'autocorrélation est une fonction paire.

[RXX(0)RXX(1)RXX(2)...RXX(P1)RXX(1)RXX(0)RXX(1)...RXX(P2)...............RXX(P1)RXX(P2)RXX(P3)...RXX(0)][ω0ω1...ωP1]=[RXY(0)RXY(1)...RXY(P1)]

RXX.W=RXY où:

  • RXX
    est une matrice positive carré de Toeplitz à diagonale constante. La diagonale
    =RXX(0)=
    énergie du signal.
  • W
    est le vecteur du filtre
  • RXY
    est le vecteur intercorrélation

Wopt=RXX1.RXY C'est le filtre optimale de Wiener Hopf.

Mathlab

a=[1,2,3]A=Toeplitz(a)A=[123213321]

Remarques:

  • Cette méthode nécessite beaucoup de calcul:
    • Autocorrélation
      p2
      fois
    • Intercorrélation
      p
      fois
    • Inversion de matrices
      Le calcul ne peut pas se faire en temps réel.
  • Lorsque les signaux
    x,y
    ne sont pas stationnaires le calcul de
    RXX
    et
    RXY
    doit être fait à chaque pas d'échantillonnage.

On lui préfère le Filtre de Wiener évolutif.

Filtre de Wiener évolutif

Dans ce filtre, l'auto et l'inter corrélation se calcule au fil de l'eau.
Explication
On sait que :

RXY(k)=1Nn=0N1ynx(nk)
où la somme représente le calcul d'une moyenne.

Pour obtenir la moyenne d'un signal, on utilise un filtre passe-bas (filtre RC).
Son expression est la suivante :

H(z)=1λ1λz1 0<λ<1

H(z)=S(z)E(z)S(z)=H(z).E(z)S(z)=(1λ1λz1).E(z)S(z)λz1S(z)=(1λ).E(z)s(n)=(1λ).e(n)+λs(n1)
RYX(k)=(1λ)ynxnk+λRYX(kY)

Avec cette méthode, les résultats sont légèrement moins bons, mais ils demeurent acceptables.

Le filtrage adaptatif

Un filtre adaptatif est un filtre numérique dont les coeficients se modifient en fonction de l'environnement (Entrées/Sorties).

Algorithme récursif des moindres carrés (RLMS)

On veut éviter de résoudre l'équation de Wiener Hopf

W=RXX1.RXY
La fonction coût qu'on a minimisé est
J(ϵ)=1Nn=0N1ϵn2

ϵn2=(ynyn^)2=(ynk=0P1wkxnk)

De manière heuristique, on part d'une valeur initiale de

wk(n).
wk(n)
est le coeficient calculé au pas d'échantillonnage
nTe
.
Wn=Wn1+terme dadaptation

Selon la valeur initiale de
W
(à droite ou à gauche) de
Woptimale
, le terme d'adaptation est positif ou négatif.
W(n)=W(n1)μ2.ϵn2ω

En utilisant la notation matricielle :

yn^=k=0P1Wkxnk=XnT.WXnT=[xnxn1...xnp+1] et W=[w0w1...wp1]ϵ2=(ynXnT.W)2ϵ2w=2(ynXnTW).XnWn=Wn1+μ(ynXnTWn1).Xn

  • Wn
    est le vecteur filtre calculé au rang
    n
  • Wn1
    est le vecteur filtre calculé au rang
    n1
  • (ynXnTWn1)
    est l'erreur à priori

Cet algorithme s'appelle l'Algorithme du Gradient ou RLMS (Recursive Least Mean Square)

wn peut être initialisé à zéro.

Le vecteur

Xn est un vecteur entrgistré
yn
est un échantillon.
Il nous reste à déterminer
μ
: pas d'adaptation.

  • Si
    μ
    est très grand, l'algorithme oscille autour de la valeur optimale.
  • Si
    μ
    est trop petit, la convergence est trop lente.

Algorithme récursif normalisé

L'étude de la convergence de l'algorithme a été faite par Widrow en 1975 et elle aboutie à :

0<μ<2P.σx2
Avec :

  • P
    l'ordre du filtre
  • σx2
    la puissance du signal d'entrée
    Pour les situations pratiques on choisit :
    μ=μ0P.σx2  avec  μ0=0,01...0,1

    L'algorithme devient :
    wn=wn1+μ0P.σx2.(ynk=0P1wkxnk).Xn

    P.σx2
    est le gain d'adaptation normalisé.
    C'est pour cela que l'algorithme s'appelle NLHS (Normised Least Mean Square)

Calcul de la puissance

σx2 de manière itérative:
σx2=Pwiss(x)=1Nk=0N1xk2

La moyenne d'une grandeur peut se calculer par un filtre passe-bas dont la fonction de transfet
H(z)
est :
H(z)=1λ1λz1=S(z)E(z)S(z)=(1λ)E(z)+λz1S(z)s(n)=(1λ)e(n)+λs(n1)

appliqué au calcul de la puissance, on obtient:
σx2=Px=(1λ)xn2+λPx

Application en travaux pratiques

Chapitre 2 : Spécification des filtres

Un filtre a deux rôles, modifier l'amplitude spectrale d'un signal et modifier la phase en fonction de la fréquence :

H(jω)=Y(jω)X(jω) : fonction de transfert
pour concevoir les filtres analogiques, on travaille avec la fonction d'affaiblissement :
A(jω)=20log10|1H(jω)|=20log10|X(jω)Y(jω)|

La phase
Φ=arg(H(jω))

Exemple de gabarits (passe-bandes)

[fp1,fp2]: Bande passante
[0,fp1][fa2,[
: Bandes atténuées
[fa1,fp1][fp2,fa2]
: Bandes de transition
Amax
: Atténuation maximale autorisée en bande passante
Amin
: Atténuation minimale demandée en bande atténuée

Par contre le gabarit d'un filtre numérique se trace avec la fonction de transfert

H:

Gabarit d'un passe-bande numérique
Parfois, il est nécessaire de préciser le gabarit de la phase; on considère:

Φ=ωτ+cste
τ
est le temps de propagation de groupe

Différents types de filtres

Passe-bas