TAS2022
TAS2022
1
Traitement Avancé du Signal Master RT M1 2022/2023
Programme
1. Transformée en Z
2. Structures, fonctions de transfert, stabilité et implémentation des filtres numériques (RIF et RII)
2
Traitement Avancé du Signal Master RT M1 2022/2023
1. Méthodes paramétriques
1. Dualité temps-fréquence
5. Transformée de Wigner-Ville
6. Analyse Temps-Echelle,
3
Traitement Avancé du Signal Master RT M1 2022/2023
Amplitude
Onde acoustique : délivré par un micro (parole, musique, …) 0
Images, Vidéos 0.954 0.956 0.958 0.96 0.962 0.964 0.966 0.968 0.97 0.972 0.974
Temps (s)
Remarque: Tout signal physique comporte une composante aléatoire (perturbation externe, bruit, erreur de
mesure, etc …). Le bruit est défini comme tout phénomène perturbateur gênant la perception ou
l’interprétation d’un signal, par analogie avec les nuisances acoustiques (interférence, bruit de fond, etc.). La
différentiation entre le signal et le bruit est artificielle et dépend de l’intérêt de l’utilisateur : les ondes
électromagnétiques d’origine galactique sont du bruit pour un ingénieur des télécommunications par
satellites et un signal utile pour les radioastronomes.
Les fonctions du traitement du signal peuvent se diviser en deux catégories : l’élaboration des signaux
(incorporation des informations) et l’interprétation des signaux (extraction des informations). Les principales
fonctions intégrées dans ces deux parties sont les suivantes [1]:
4
Traitement Avancé du Signal Master RT M1 2022/2023
1. Théorie des systèmes linéaires et invariants dans le temps (SLIT) discrets et continus
Un système linéaire est un modèle de système qui applique un opérateur linéaire à un signal d'entrée.
C'est une abstraction mathématique très utile en automatique, traitement du signal, mécanique et
télécommunications. Les systèmes linéaires sont ainsi fréquemment utilisés pour décrire un système non
linéaire en ignorant les petites non-linéarités.
Un système est continu si à une entrée continue x(t), il fournit une sortie continue y(t). Un système
discret fera correspondre à une suite d'entrées discrètes x(n) une suite de sorties discrètes y(n).
- Le système sera dit linéaire si quand on applique une entrée k x(t) (ou k.x(n) ) , la sortie sera k.y(t) (ou k.y(n)
). Si deux entréesx1(t) et x2(t) engendrent deux sorties y1(t) et y2(t) alors x1(t) + x2(t) engendrera y1(t) + y2(t).
(De même, dans le cas discret : si deux entrées x1(n) et x2(n) engendrent deux sorties y1(n) et y2(n) ) alors x1(n) +
x2(n) engendrera y1(n) + y2(n))
Exemple
3
x1+x2
y1+y2
2
-1
-2
-3
0 0.002 0.004 0.006 0.008 0.01 0.012 0.014
- S’il y a invariance dans le temps, une translation de l'entrée x(t)x(t-τ) (oux(n)x(n-m)) se traduira par une
même translation dans le temps de la sortie y(t)y(t-τ).(ouy(n)y(n-m)).
x(t-τ) y(t-τ)
Système
x(n-m) y(n-m)
Si le système est invariant, cela implique que le système réagit de la même façon quel que soit l’instant auquel
nous appliquons ses excitations. Cette propriété exprime que la caractéristique du système ne dépend pas de
l’origine du temps, on parle encore de stationnarité.
5
Traitement Avancé du Signal Master RT M1 2022/2023
Convolution
Si les hypothèses de linéarité et d'invariance temporelle sont vérifiées, on peut caractériser le système par sa
réponse impulsionnelle h(t) ( ou h(n) ).
x(t) y(t)
h(t) (ou h(n) )
x(n) y(n)
On peut en déduire l'effet d'une entrée quelconque sous la forme d'une convolution. Cette dernière est
l’opération de traitement de signal la plus fondamentale. Elle indique que la valeur du signal de sortie à
l’instant t (ou n )est obtenue par la sommation (intégrale) pondérée des valeurs passées du signal d'excitation
x(t) ( ou x(n) ). La fonction de pondération est précisément la réponse impulsionnelleh(t) ( ou h(n) ):
y (t ) = x (t ) ∗ h (t ) = x (τ ) h (t − τ ) dτ = x (t − τ ) h (τ ) dτ
∞ ∞
y ( n) = h ( n ) ∗ x( n) = h ( m) x ( n − m ) =
m = −∞
x ( m) h ( n − m )
m = −∞
La réponse impulsionnelle h(t) (ou h(n))est le signal qu'on obtient en sortie y(t)=h(t)( ou y(n)=h(n)) si on
applique en entrée une impulsion "de Dirac'' x(t)=δ(t) ( ou x(n)=δ(n) ). Ainsi, le Dirac est l'élément neutre de
l'opération de convolution:
x (t ) ∗ δ (t ) = x (t )
δ ( n) ∗ x ( n) = x ( n)
x(t )δ (t − t o ) = x(t o )δ (t − t o )
+∞
x ( t ) ∗ δ (t − t 0 ) = x (t − t 0 )
−∞
x(t )δ (t − t ) dt = x(t
o 0 )
Le calcul de la convolution consiste donc à calculer la somme du produitx(τ)h(t -τ)(ou x(m)h(n -m)) .
Le signal h(t -τ) (ou h(n-m) ) est simplement le signal initialh(τ) ( ou h(m) ), retourné dans le temps pour
donner h(-τ) ( ou h(-m)) puis translaté de t ( ou n ). En calculant alors l’ensemble des produits obtenus en
faisant « glisser » h, c’est-à-dire pour tous les décalages de t ( oun ) , on obtient le produit de convolution
pour tout t ( ou n ).
h(-m)
0 N-1 n 0 n 0 m
6
Traitement Avancé du Signal Master RT M1 2022/2023
On distingue 3 cas :
x(m)
h(n-m)
t +T / 2
1 1
h (t ) = Π (t )
T T
y (t ) =
T x (τ ) dτ
t −T / 2
N/2
1
h(n) =
1
Π ( n ) y(n) =
N + 1 N +1
x(n + m)
N + 1 m =− N / 2
Calculer y(n)=x(n)*h(n)
Y(n)={2, 3, 5, 10, 3, 9}
Remarques
Si on applique à un SLIT une entrée sinusoïdale réelle ou complexe de fréquence f0, alors, la sortie sera une
sinusoïde dont l'amplitude et la phase pourront être modifiées mais qui conservera la même forme (une
sinusoïde) et la même fréquence f0. On dit que les sinusoïdes sont les fonctions propres des SLIT.
Un système linéaire invariant est un système dont le comportement dans le temps, peut-être décrit par :
M i N i
- soit par une équation différentielle linéaire à coefficients constants: a i d y ( t ) = bi d x ( t )
i=0 dt i i=0 dt i
M N
- soit par une équation aux différences, pour le cas discret: a i y ( n − i ) = bi x ( n − i ) ,
i=0 i=0
7
Traitement Avancé du Signal Master RT M1 2022/2023
2. Stabilité et causalité
Une contrainte importante pour la formalisation de nombreux problèmes est de respecter la notion de
causalité (les effets ne peuvent pas précéder la cause). Dans le cas des SLIT, cette causalité se traduit par le
fait que pour: h(t) = 0 pour t<0 ( ou h(n) = 0 pour n<0 ).
t +∞
x(t) = 0, t < t0 alors y(t) = 0, t < t0 h(t ) = 0, t < 0 , y (t ) = −∞ x(τ )h(t − τ ) dτ = 0 x(t − τ )h(τ ) dτ
n
- si h et x sont causaux y ( n) = h ( n − m) x ( m )
m =0
Remarque:
Nous pouvons envisager mémoriser les signaux d’entrée et faire un traitement de ceux-ci en temps différé,
les systèmes utilisés ne sont plus alors nécessairement causaux car pour élaborer la sortie à l’instant ti( ou ni
), nous disposons en mémoire des entrées aux instants suivants. C’est souvent le cas en traitement d’image,
Une autre notion fondamentale est la stabilité des systèmes. La propriété de stabilité des systèmes
Système stable
bouclés est non seulement une performance mais une exigence pour le bon fonctionnement d’une boucle
7
6.5
d’asservissement ou de régulation. Une boucle instable est une boucle inutilisable. La définition la plus
6
5.5
4.5
3.5
3
0 5 10 15 20 25 30 35 40 45 50
Système instable
10
h2(t) 3
) 1
0
0 5 10 15 20 25 30 35 40 45 50
On dit qu'un système est stable si, en lui appliquant une entrée bornée quelconque, la sortie reste
ou
Remarque : Pour le reste, les définitions seront données pour le cas discret.
3. Energie et puissance
Toute transmission d’information s’accompagne de transferts d’énergie. En effet, les signaux continus
ou discrets sont essentiellement caractérisés par l’énergie ou la puissance qu’ils véhiculent. Ce sont les seules
grandeurs physiques auxquelles sont sensibles les détecteurs. Beaucoup de capteurs physiques mesurent
une énergie ou une quantité quadratique. Par exemple, les capteurs optiques mesurent une intensité, les
compteurs d’électricité mesurent une énergie, etc. Compte tenu de la définition fondamentale, l’énergie du
signal entre les instants t et t+dt est : |x(t)|2 dt (puissance instantanée multipliée par le temps).
+∞
2
Soit un signal x(n) à temps discret, tel que x(n) existe et converge. Alors le signal est dit à énergie
−∞
finie et la valeur de cette somme est appelée énergie du signal :
+∞
E x = x(n)
2
−∞
Exemples:
x(n) = Rect(n/N) énergie finie. x(n) = a (constante) et x(n) = A sin(2πf0n) ne sont pas à énergie finie
Pour un signal périodique, cette somme ne converge pas. On peut néanmoins définir la puissance d'un signal
x(n) périodique de période N par :
1 N / 2−1 1 N −1
Px = Px =
2 2
x ( n) ou x(n)
N −N / 2 2.N − N
Dans le cas général, on parle de signaux à puissance moyenne finie définie par:
1 N / 2−1 1 N −1
Px = limN →∞ Px = limN →∞
2 2
x(n) ou x(n)
N −N / 2 2.N − N
Exemples:
Il existe des signaux ni périodiques, ni d'énergie finie, pour lesquels la puissance ne peut être définie, comme
par exemple la rampe x(n)=n. Il s'agit là de définitions mathématiques, en pratique, un signal mesuré ne l'est
jamais sur un intervalle de temps infini. On peut commencer à visualiser un signal à un instant qu'on prendra
comme origine des temps, et dans ce cas on arrêtera son examen au bout d'un temps Tobs :
Nobs
Ex = x ( n)
2
n =0
9
Traitement Avancé du Signal Master RT M1 2022/2023
Remarque
Le calcul de l'énergie ou la puissance permet d'obtenir une première caractérisation du signal. Par
ailleurs, la théorie du signal a largement développé des méthodes d’étude basées sur la corrélation pour
caractériser le comportement temporel du signal.
Exercice d'application :
4. Corrélation et auto-corrélation
La fonction de corrélation permet de mesurer le degré de ressemblance entre deux signaux en fonction
d’un décalage. Considérons x(n) et y (n) deux signaux d'énergie finie, la fonction d'intercorrélation Rx,y(k) est
∞
définie par: Rxy (k ) = x(n) y
n = −∞
*
(n − k )
1 N / 2 −1
Pour des signaux à puissance moyenne finie, elle vaut : Rxy (k ) = x(n) y* (n − k )
N n=−N / 2
Exemples
Soient un signal aléatoire et sa version décalée de 50s. On remarque que les signaux se ressemblent le plus
quand y(n) est décalé de 50 secondes.
Intercorrélation maximum à 50 s
3 3 200
x(n) y(n)
2
2
150
1
1
0 100
0
-1
50
-1
-2
0
-3 -2
-4
0 50 100 150 200 250 -3 -50
0 50 100 150 200 250 -250 -200 -150 -100 -50 0 50 100 150 200 250
Pour l'auto-corrélation, on remplace y(n) par x(n) on obtient l'expression de l'auto-corrélation pour les
∞
R
signaux à énergie finie: xx (k ) =
x(n) x* (n − k)
n=−∞
L'auto-corrélation permet de détecter des régularités, des profils répétés dans un signal comme un
10
Traitement Avancé du Signal Master RT M1 2022/2023
Propriétés :
Auto-corrélation des signaux périodiques : Le calcul sur une seule période suffit. L’auto-corrélation
d’un signal périodique est-elle même périodique. Par définition, le signal périodique ressemble parfaitement
à lui-même, décalé d’une ou plusieurs périodes.
x(n) périodique Autocorrélation de x(n)
- signaux périodiques 1 1
1 N / 2 −1 0.8 0.8
Rx (k ) = x(n) x* (n − k )
N n =− N / 2 0.6 0.6
0.4 0.4
0.2 0.2
Extraction d'un signal noyé dans du bruit, mesure d'un temps ou retard, détection d'un signal périodique
(Voir TP n° 1).L'exemple ci-dessous illustre l'auto-corrélation d'un signal sinusoïdal d'amplitude 1 noyé dans
du bruit Gaussien de variance 1.
4 200
sinusoide bruitée Autocorrélation de x(n)
3 150
2
100
1
50
0
0
-1
-50
-2
-100
-3 -0.015 -0.01 -0.005 0 0.005 0.01 0.015
0 0.002 0.004 0.006 0.008 0.01 0.012 0.014
La corrélation est largement utilisée dans les systèmes radar. Ainsi, pour détecter un avion, on envoie une
impulsion, puis on reçoit une version retardée, atténuée et bruitée de cette impulsion . L'intercorrélation du
signal reçu et émis présentera un pic à l'instant correspondant au retard.
11
Traitement Avancé du Signal Master RT M1 2022/2023
signal émis
0.8
0.6
0.4
0.2
0
0 50 100 150 200 250 300 350 400 450 500
1.5
signal reçu: bruité et retardé de 50S 10
Intercorrélation du signal reçu et émis avec pic en 50s
1
8
0.5 6
4
0
2
-0.5
0
-1
-2
-1.5
0 50 100 150 200 250 300 350 400 450 500 -4
0 50 100 150 200 250 300 350 400 450 500
Remarque: La notion de bruit est relative, elle dépend du contexte. Le rapport signal/bruit désigne la qualité
de la transmission d'une information par rapport aux parasites. Il est défini par: SNRdb=20 Log(PS/PB)
3 1.5
SNR=0 db SNR=20 db
2 1
1 0.5
0 0
-1
-0.5
-2
-1
-3
0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 -1.5
0 0.002 0.004 0.006 0.008 0.01 0.012 0.014
- Si les ai sont ≠ de 0, le système est dit récursif (RII), il est non récursif s'il ne dépend que des x(n-i) (RIF)
K
- Si le système est à réponse impulsionnelle de durée finie (RIF), alors : y ( n) =
h(m) x (n − m)
m =0
Dans ce cas, le système numérique est une fenêtre centrée sur les K plus récents échantillons.
+∞
- Si le système est à réponse impulsionnelle de durée infinie (RII) : y ( n) = h (m ) x (n − m )
m =0
Dans ce cas, il est nécessaire de connaître tous les échantillons présents et passés, le système à une mémoire
de longueur infinie.
12
Traitement Avancé du Signal Master RT M1 2022/2023
avec comme réponse impulsionnelle h(n)=δ(n)+a1δ (n-1)+a2δ (n-2)+.......+akδ (n-k) qui, on peut le constater,est
bien finie.
Exemple 2y(n)=x(n)+a1y(n-1) est l'équation aux différence finies d'un filtre RII
y(n)=x(n)+a1x(n-1)+a12x(n-2))+a13x(n-3)+a14x(n-3)+.....+ a1my(n-m)
En poursuivant le procédé à l'infini y(n) dépend d'une infinité de x(n-k) ce qui en fait un filtre RII.
Exemple d'application
Les séquences x(n) (réel) et y(n) représentent respectivement l’entrée et la sortie d’un système discret. Pour
chaque cas, identifiez celles représentant
a) des systèmes linéaires, b) des systèmes causals, c) des systèmes invariants aux translations de n,
d) des systèmes assurément ou possiblement stables (en fonctions des constantes)
1. y(n) = x(n) + bx(n-1) 2. y(n) = x(n) + bx(n+1) 6. y(n) = bx(n) b : constante réelle
3. y(n) = nx(n) 5. y(n) = x(n)en 7. y(n) = |x(n)|
4. y(n) = x(n) sin(2πf0n)
a) Tous sauf 6 et 7 b) Tous sauf 2 c) Les systèmes 1, 2, 6 et 7 d) 1 (b finie), 2 (b finie), 4, 6 (b finie) et 7.
13
Traitement Avancé du Signal Master RT M1 2022/2023
La transformée de Fourier est un outil précieux d'analyse et de traitement des signaux. Cependant,
dans certains problèmes (comme le filtrage numérique), les limites de la TF sont vite atteintes. La transformée
en Z, qui s'applique aux signaux discrets, généralise la TF et permet de dépasser ces limites [10]. Elle est tout-
à-fait analogue à la transformée de Laplace, mais plus facile à utiliser. Ce type de transformée permet de
décrire aisément les signaux à temps discret et la réponse des systèmes linéaires invariants soumis à des
entrées diverses. C'est un outil qui permet de calculer la réponse impulsionnelle d’un système linéaire
invariant décrit par une équation aux différences finies. Elle permet l'interprétation directe des
caractéristiques des signaux et des filtres dans le domaine des fréquences [9].
1. Transformée en Z et propriétés
+∞
• Définition : La TZ est la généralisation de la TFTD ( X ( f ) =
n=−∞
x(nT )e
e
−2π j f nTe
).
En effet, comme cette transformation est une série infinie, elle n’existera que pour lesvaleurs de z pour
lesquelles cette série converge.
Exemples :
- x(n) = δ(n)X(z) = 1,
L'ensemble des valeurs de la variable complexe z pour lesquelles la série converge est appelée Région De
Convergence (RDC): +∞
RDC = z ∈ C / x ( n).z − n < +∞
n = −∞
De façon générale, on montre que la RDC est un anneau de convergence centré sur l'origine défini par :
14
Traitement Avancé du Signal Master RT M1 2022/2023
r2 Re(z)
RDC
- x(n)=0 pour n<n0r2= +∞,
Exemples
+∞
a n si n ≥ 0 z
- Soit a > 0, x ( n ) = X(z)= a n
z −n =
z−a
, convergente pour z > a .
0 si non n =0
0 si n ≥ 0 z
- Soitb > 0, y ( n ) = Y(z)= − b z n −n
=
z −b
, convergente pour z < b .
− b si n < 0
n
n<0
a n si n ≥ 0
Soient a >0 ,b> 0, w ( z ) =
z z
- W(z)= − , convergente pour b > z > a .
b si n < 0
n z−a z−b
Les propriétés qui sont les plus utilisées sont résumées comme suit :
Exemples
15
Traitement Avancé du Signal Master RT M1 2022/2023
1 jw 0 n
x ( n ) = cos( wo n )U ( n ) = (e + e − jw 0 n )U ( n )
2
1 1 1 1 − cos( w0 ) z −1
X (z) = + = avec z >1
2 1 − e jw 0 z −1 1 − e − jw0 z −1 1 − 2 cos( w0 ) z −1 + z −2
Les systèmes linéaires invariants décrits par une équation aux différences finies possèdent une
transformée en Z rationnelle c'est ainsi que celles-ci vont s’écrire comme le rapport de deux polynômes en
z-1.
M N M N
i =0
i
−i
Y ( z) = bi .z −i X ( z)
i =0
b .z i
−i
b N ( z ) b0 M − N ∏ (z − z ) = K z
N
∏ (z − z )
N
H ( z) = i =0
= 0 z M −N = i =1 M −N i =1
i i
z
∏ (z − p ) ∏ (z − p )
M M M
a0 D ( z ) a0
a .z
i =0
i
−i
i =1 i i =1 i
On appelle zéros, les valeurs de z pour lesquelles H(z)=0 et on appelle pôles, les valeurs de z pour
lesquelles H(z) est infini (annule le dénominateur). C'est ainsi que H(z) possède N zéros (zi), M pôles (pi). Si
M>N, elle possède (M-N) zéros en 0, sinon (N-M) pôles en 0.
Ainsi, la position de ses pôles et de ses zéros ( +le facteur d’amplitude K=b0/a0) va nous fournir une
description complète de H(z) (par conséquent de h(n) et H(f)) donc du comportement du système. H(z) peut
donc être représentée sous la forme d’un cercle modélisant la position des pôle set des zéros dans le plan
complexe.
Exemple
1 Im(z)
3z − 2
H ( z) =
( z − 1)(z + 0.5)
-1
X X
Un zéro en 2/3 et deux pôles p1= -0.5 et p2= 1 1 Re(z)
Remarques
-1
- Dans la plupart des systèmes, les ai et le bi sont réels les pôles et les zéros sont soient réels soient des
paires de complexes conjuguées.
- Rappelons que le rayon d'un système causal se trouve à l'extérieur d'un cercle. Par ailleurs, s'il est stable :
+∞
h ( n ) < ∞ , puisque
n
H ( z) =
h( n).z − n , il suffit donc que z=1 fasse partie de la RDC.
n = −∞
16
Traitement Avancé du Signal Master RT M1 2022/2023
- Pour un système causal et stable, tous les pôles sont à l’intérieur du cercle unité (|pi|<1, ∀i). Le domaine
de convergence ne peut contenir de pôles puisque la TZ ne converge pas aux pôles. S’il est anti-causal, il sera
stable si les pôles sont à l’extérieur du cercle unité.
N
- Si le filtre est non-récursif H ( z ) = bi .z . Un filtre RIF a tous ses pôles à l’origine et sera donc toujours
−i
i =0
stable.
20
15 0.8
0.6
10
0.4
200 0.2 1
5
0 0.9
150
0 -0.2 0.8
-5
-0.6 0.6
50
-0.8 0.5
-10
0 2 4 6 8 10 12 14 16 18
0
-1 0.4
0 2 4 6 8 10 12 14 16 18
0.3
-50
0.2
-100
0.1
-150
0 5 10 15 20 25 1 Im(z) 0
0 5 10 15 20 25
X X 300
X 250
200
1 X XX X X X 150
-1 1 100
Re(z)
0.8
0.6 50
0.4 X 0
0.2
X
0
X -100
-50
-0.2 0 2 4 6 8 10 12 14 16 18
-0.4
-0.6
-1 1
1
-0.8
1 0.9
-1 0.8
0 5 10 15 20 25 0.8
0.8
0.6
0.7
0.6
0.4 0.6
0.4
0.2 0.5
0 0.2 0.4
0.3
-0.2 0
0.2
-0.4
-0.2
0.1
-0.6
-0.4
0
0 5 10 15 20 25
-0.8
0 5 10 15 20 25
-0.6
0 2 4 6 8 10 12 14 16 18
- A un pôle pi simple ou multiple va correspondre une réponse impulsionnelle qui converge si |pi|<1. Elle
divergera dans le cas contraire, soit si |pi|>1 .
- Sachant qu'à chaque pôle complexe est associéun pôle conjugué cela donnera une réponse impulsionnelle
h(n) oscillante(cosinus ou sinus) amortiesi|pi=1,2|< 1 ou divergentesi |pi=1,2|> 1.
- Dans un système à phase minimale, tous les zéros sont à l’intérieur du cercle unité (|zi|<1, ∀i).
On suppose que le cercle unité (|z|=1) ∈ RDC de X (z). On restreint le calcul de X(z) au cercle unité
en posant z = e 2πj f Te.
Lorsqu'un zéro est placé sur un point donné du plan en z, la réponse fréquentielle sera de 0 au point
considéré. Un pôle quant à lui produira un pic au point correspondant. Plus les pôles ou les zéros sont
proches du cercle unité, plus ils influencent la réponse en fréquence [11].
17
Traitement Avancé du Signal Master RT M1 2022/2023
– un zéro sur le cercle unité introduit une annulation du module pour la fréquence correspondant
– Un zéro au voisinage du cercle unité introduit une atténuation dans le module de la réponse en fréquence.
Atténuation d’autant plus importante que le zéro est proche du cercle unité.
– Un pôle sur le cercle unité introduit une résonance infinie dans le module de lar éponse en fréquence pour
la fréquence correspondante.
– Un pôle au voisinage du cercle unité introduit une résonance d’autant plus importante dans le module de
la réponse en fréquence que le pôle est proche du cercle unité.
1) Sur la figure ci-dessous le cercle complet correspond à une fréquence d'échantillonnage fe. Des pôles
proches du cercle unité sont à l'origine de larges pics tandis que des zéros proches ou sur le cercle unité
produisent des minima. Ce tracé nous permettra d'identifier la nature du filtre. On peut aussi avoir une idée
sur son comportement général : passe-bas, passe-haut ou passe-bande, connaître sa ou ses fréquences de
coupure.
f=fe/4
2) 1 10
0.8 f=fe/6 9
0.6
X 8
0.4 7
Partie Imaginaire
0.2
2 60° 6
f=fe/2 0
2 f=0
5
-0.2
4
-0.4
3
-0.6
X 2
-0.8
1
-1
0
-1 -0.5 0 0.5 1 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
z 2 − 2z + 1 j 72°
3) H ( z) = Zéros double en z=-1, pôles p12 = ±0.6e
z − 0.371z + 0.36
2
- Un zéro double en z = 1 ⇒ |H(f)| = 0 pour f = 0 -Des pôles proches du cercle unité ⇒ maxima.
3
1
0.8
2.5
0.6
X
0.4
2
Partie Imaginaire
0.2
2
0
1.5
-0.2
-0.4 1
-0.6 X
-0.8 0.5
-1
-1 -0.5 0 0.5 1 0
Partie Réelle 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
Remarque : Puisque les coefficients du filtre sont réels, les pôles et zéros sont réels (sur l'axe des réels) ou
paires de complexes conjugués.
18
Traitement Avancé du Signal Master RT M1 2022/2023
1. La relation générale de la transformée en z inverse est donnée par l’équation donnée par l'intégrale de
1 n −1
2.π . j
Cauchy : x(n) = X ( z).z .dz , où C est un contour fermé parcouru dans le sens inverse des aiguilles
C
Re s[ z n −1
X ( z)] z = pi =
1 d m−1
(m − 1)! dz m−1
[ ]
(z − pi )m z n−1 X ( z) z = pi
2. Transformée inverse par division polynômiale :Il est possible de calculer la transformée en Z inverse selon
les puissances croissantes de z-1(système causal) ou selon les puissances décroissantes de z (système anti-
−1
causal). X ( z ) = C z − n TZ
→x ( n) = C n n
n
3. L’idée générale de cette approche consiste à trouver pour une fonction X(z) complexe un développement
en fonctions en Z plus simples et pour lesquelles une transformée inverse est connue:
X ( z ) = X i ( z ) TZ
→x(n) = xi (n)
−1
i i
où les Xi(z) sont des fonctions dont les TZ-1 sont connues (Voir page 38).
Un filtre numérique est constitué d'un groupement de circuits logiques astreints à un processus de
calcul (ou algorithme) qui confère à ce filtre une fonction déterminée (passe-bas, passe-haut, passe-bande,
réjecteur de bande, intégrateur[y(n)=(x(n)+x(n-1))/2], différentiateur[y(n)=(x(n)-x(n-1))/2], ...). Il faut
souligner que certains filtres ne sont pas conçus pour arrêter une fréquence, mais pour modifier légèrement
le gain à différentes fréquences, comme les égaliseurs. Ce sont tous des systèmes linéaires, discrets, invariants
dans le temps et unidimensionnels. De plus, pour qu’ils soient physiquement réalisables, il faut qu’ils soient
nécessairement causaux.
19
Traitement Avancé du Signal Master RT M1 2022/2023
Les filtres représentés ci-dessus sont idéaux. Dans un cas réel il 1.4
Amplitude
passantes et atténuées ne sont également pas idéales, elles 0.6
0
0 0.2 0.4 0.6 0.8 1 1.2 1.4
Bande passante Bande de transition Bande atténuée
Idéalement, il est souhaitable qu'un filtre possède une phase linéaire dans la bande passante. Une phase
linéaire assurera un même déphasage pour toutes les fréquences (pas de distorsion). Les filtres FIR peuvent
générer des filtres à phase linéaire à la condition que la réponse impulsionnelle soit symétrique.
Et la dérivée de cette dernière par rapport à f fournit le ‘retard de groupe ‘, défini donc par :
dφ ( f )
β =−
df
et qui correspond au retard subi par le signal après être passé par un filtre. Si la phase est linéaire (filtres RIF
symétrique), le retard est constant et le signal à la sortie aura donc une distorsion minimale puisque l’effet
de la phase sur le signal sera un simple décalage temporel (primordial dans un système audio)
Exemple: y ( n ) =
1
(x( n) + 2 x ( n − 1) + x (n − 2) )
4
−2π j f
h(n) =
1
(δ ( n ) + 2δ ( n − 1) + δ ( n − 2) ) et H( f ) = e cos2 (πf )
4
20
Traitement Avancé du Signal Master RT M1 2022/2023
1 0 2
0.9 1.8
-0.5
0.8 1.6
0.7 -1 1.4
0.6
Phase de H(f)
-1.5
0.5 1
-2 0.8
0.4
0.3 0.6
-2.5
0.2 0.4
-3
0.1 0.2
0 -3.5 0
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
fréquence (Hz) fréquence (Hz) fréquence (Hz)
Pour un filtre RII, le retard de groupe n'est pas constant. On souhaite avoir un minimum de retard en
s'orientant faire des filtres à minimum de phase. Un filtre est dit à minimum de phase si lui et son inverse
sont tous deux stables et causaux. Ce qui implique que tous les pôles et zéros de la fonction de transfert H(z)
sont à l'intérieur du disque unité. Un filtre à phase minimale présente une réponse impulsionnelle
(coefficients de H(Z)) à tendance décroissante car ses zéros sont en module inférieurs ou égaux à l'unité.
y(n) est obtenu plutôt à partir des échantillons “récents” du signal d'entrée qu'à partir des échantillons
“anciens”. Il faut noter que tout filtre causal et stable peut être décomposé en un produit de filtres à phase
minimale par une cellule passe-tout.
• RII ou RIF ?
Les filtres RII, on l'avantage qu'ils sont efficaces. Avec très peu de pôles et zéros on peut assurer la
plupart des réponses fréquentielles dont on peut avoir besoin dans les applications audio. Cependant, le
filtre étant récursif, les erreurs de précision numérique deviennent une question d'importance, car ils peuvent
s'amplifier et devenir hors de contrôle, d'abord dans la forme de bruit, mais éventuellement dans la forme
d'instabilité. La forme de la réponse impulsionnelle n'est pas facile à déterminer, non plus, car elle est définie
indirectement par les pôles et zéros de H(f).
Par contre, les filtres RIF n'ont jamais des problèmes d'instabilité, car la sortie n'est qu'une somme
finie d’échantillons de l'entrée. Cependant, quand la réponse impulsionnelle est longue, le numéro
d'opérations peut devenir un facteur décisif quand il faut choisir entre RIF ou RII. Un autre avantage des RIF
est le retard de groupe constant, qui permet d'avoir une distorsion de phase minimale sur le signal traité [13].
L’application d’un filtre numérique implique le calcul de la sortie y(n) à l’instant t=nTe à partir des sorties et
entrées précédentes plus la valeur courante de l’entrée.
M N
a y(n − i) = b x(n − i)
i =0
i
i =0
i En prenant a0=1, on obtient
21
Traitement Avancé du Signal Master RT M1 2022/2023
…
b0 b1 b2 bN
..
….. y(n)
Un filtre numérique est généralement constitué des éléments suivants : un ou plusieurs organes de
retard (ce sont des registres à décalage jouant le rôle de mémoires retardées), pilotés par une horloge de
période; des opérateurs arithmétiques (additionneurs et multiplieurs); des registres fournissant les
coefficients de pondération du filtre [18]. Ainsi, à chaque top d’horloge Te, les valeurs des registres subissent
un décalage permettant de calculer la nouvelle sortie. Cette structure peut être aussi réalisée par logiciel.
La synthèse d’un filtre est un ensemble de processus qui débute par la définition des caractéristiques
du filtre, jusqu’à sa réalisation informatique et/ou électronique, en passant par la détermination de ses
coefficients. Pour synthétiser un filtre numérique, on considère connu le gabarit du filtre analogique et on
cherche un système numérique caractérisée par une fonction de transfert H(z) à insérer dans le circuit ci-
dessus permettant de satisfaire le gabarit analogique.
x(t) y(t)
Système analogique
Ha(p)
La détermination de la fonction de transfert d’un filtre numérique, par une méthode directe, n’est pas
toujours très simple. Par contre, le problème qui consiste à transformer un filtre analogique en un filtre
numérique est relativement simple. De ce fait, de nombreuses méthodes sont proposées pour concevoir un
filtre numérique à partir du filtre analogique équivalent. Dans tous les cas, la synthèse d’un filtre numérique
est une approximation d’un filtre analogique idéal équivalent. Il est nécessaire de contraindre un certain
nombre de paramètres.
22
Traitement Avancé du Signal Master RT M1 2022/2023
Les fonctions modèles utilisées pour la synthèse des filtres sont soit la réponse impulsionnelle soit la
réponse en fréquence (celle-ci est préférée) de filtres analogiques connus. Si l'on emploie la réponse
impulsionnelle, les éléments h(n) de la réponse impulsionnelle numérique sont obtenus en calculant h(t), la
réponse impulsionnelle du filtre analogique, aux instants t=nTe.
0.015
0
(droit et avec une droite de transition perpendiculaire). Il s'ensuit
-0.005
que les filtres qui vont pouvoir être réellement synthétisés n'ont
-0.01
pas de réponse fréquentielle correspondant à la fonction porte, 0 5
t(s)
10 15
2 10
1.8
8
1.6
6
Réponse impulsionnelle h(n)
1.4
1.2
4
H(f) idéal
0.8
TF 2
0.6 0
0.4
Troncature
-2
0.2
0 -4
-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 -10 -8 -6 -4 -2 0 2 4 6 8 10
frequence (Hz) t(s)
2 10
1.8
8
1.6
1.4 6
Réponse impulsionnelle h(n)
1.2
4
H(f) réel
0.8 TF -1 2
0.6
0
0.4
-2
0.2
0 -4
-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 -10 -8 -6 -4 -2 0 2 4 6 8 10
frequence (Hz) t(s)
Comme l'illustre la figure suivante, la troncature (multiplication par une porte de largeur NTe) du sinc
dans le domaine temporel se traduira par une convolution dans le domaine fréquentiel du filtre idéal avec
un sinc s'annulant tous les NTe. Pour de grandes valeurs de N, les sinc dans la bande passante se
compenseront les uns les autres mais autour des points de discontinuité (fréquence de coupure), les
ondulations restent apparentes.
FGE, USTHB [ assiakourgli@[Link] / [Link]
23
Traitement Avancé du Signal Master RT M1 2022/2023
0.06
0.04
0.02
-0.02
-5 -4 -3 -2 -1 0 1 2 3 4 5
0.06
0.04
0.02
-0.02
-5 -4 -3 -2 -1 0 1 2 3 4 5
On peut observer que les différences vis-à-vis du filtre idéal (soit la fonction porte) sont principalement
les ondulations dans la bande passante et dans la bande atténuée ainsi que la largeur de la transition.
C'est ainsi que les spécifications du filtre vont être définies par un gabarit fréquentiel linéaire ou en dB
(décibels). Ce gabarit indique la ou les fréquences de coupure, la largeur de la bande de transition minimale
souhaitée, le maximum d’ondulation de la bande passante et de la bande atténuée, la fréquence
d'échantillonnage et éventuellement l'ordre maximal permis.
En pratique, plus les fréquences fa et fp sont proches, plus l'ordre du filtre devra être élevé. Pour un filtre idéal,
ces valeurs seraient confondues.
L'emploi des filtres RIF peut se révéler attrayant eu égard à ses nombreux avantages : stabilité
inconditionnelle (Tous les pôles sont en 0), phase linéaire possible. Néanmoins, ils présentent l'inconvénient
de nécessiter un plus grand nombre de coefficients que les filtres RII pour obtenir les mêmes caractéristiques
24
Traitement Avancé du Signal Master RT M1 2022/2023
fréquentielles à cause de l'absence de pôles hors 0. Ainsi, toute fonction de filtrage numérique stable et
causale peut être approchée par la fonction de transfert d'un filtre RIF.
Rappelons que la sortie d'un filtre RIF va s'exprimer comme une combinaison linéaire d'un ensemble fini
d'éléments d'entrée :
N −1 N −1 N −1
y(n) = bi x(n − i) d’où H ( z ) = bn .z −n
H( f ) = b .e n
− 2 π j f nTe
i=0 n =0 n =0
Ainsi, les coefficients de pondération ne sont rien d'autre que les valeurs de la réponse impulsionnelle du
filtre. Ces coefficients constituent les coefficients du développement en série de Fourier de la fonction de
transfert H(f) (voir TFTD chapitre 1)
Du fait qu'un filtre RIF possède une fonction de transfert polynomiale (non rationnelle), il ne peut être obtenu
par transposition d'un filtre continu. Les deux méthodes les plus utilisées pour l’approximation des filtres
numériques RIF sont alors:
o Développement par série de Fourier : cette série est ensuite tronquée par des fonctions fenêtres pour
limiter la réponse impulsionnelle. Les coefficients de Fourier coïncident avec les échantillons de la
réponse impulsionnelle du Filtre.
o Echantillonnage de la réponse fréquentielle : Cette méthode fait appel à la TFD. Celle-ci est appliquée
aux coefficients recherchés bi pour obtenir une suite fréquentielle qui corresponde à la réponse
fréquentielle du filtre.
Il existe d'autres méthodes telles les méthodes d’optimisation qui sont basées sur la minimisation d’un critère
d’erreur entre la courbe réelle et le filtre idéal [9].
1. A partir du gabarit idéal du filtre, on détermine les coefficients du filtre en limitant le calcul à N valeurs
réparties symétriquement autour de n=0. Puis, on calcule de la TFTD inverse du filtre idéal qui nous
permettra de retrouver les échantillons de la réponse impulsionnelle soient les coefficients du filtre :
1 fe / 2
H ( f )e 2 π j f n Te df N impair ( Filtre de type I )
f e − fe / 2
h( n) = fe / 2
1 H ( f )e − j π f Te e 2 π j f n Te df N pair ( Filtre de type II )
f e − f/ 2
e
Remarques :
• Le filtre de type I est le plus employé (C’est celui que nous développons dans ce qui suit) hormis pour
les différenciateurs.
25
Traitement Avancé du Signal Master RT M1 2022/2023
• Le filtre de type II possède un zéro en -1 (Nombre de zéros est impair donc répartition de zéros en paires
conjuguées + un zéro en fe/2 puisque une somme de cosinus) , il ne peut donc être employé pour un
passe-haut ou un coupe-bande.
• Il existe des filtres de type III (N impair) et IV (N pair) permettant d’obtenir une réponse impulsionnelle
anti-symétrique avec déphasage linéaire également (somme de sinus au lieu de cosinus). Ils possèdent
tous deux des zéros en 1 (sin en f=0 est nulle), leur utilisation n’est donc pas adaptée aux filtres passe-
bas. En outre le filtre de type III possède également un 0 en -1, on ne peut l’employer que pour un passe-
bande.
2. Cette méthode produit une série infinie de coefficients, on limite, alors la réponse impulsionnelle à N
échantillons (troncature). Sachant que la troncature induit des ondulations, on peut faire appel aux fenêtres
de pondération pour les atténuer. Ainsi, la réponse impulsionnelle idéale h(n) sera multipliée par la fenêtre
discrète wN(n) de longueur N:
3. Il ne reste plus qu'à décaler la réponse impulsionnelle h'(n) pour avoir une solution causale.
Remarque : Pour le choix de fc il faudra faire attention aux fréquences de coupure à prendre en compte. Afin
d'avoir de bons résultats lors de la synthèse, ce ne sont pas les fréquences de coupure du filtre idéal qu'il faut
utiliser mais il faut déplacer celles-ci afin de les centrer dans la zone de transition. Pour un passe-bas,
l'augmenter de la demi zone de transition (∆f), soit fp+∆f et pour un passe-haut, la diminuer de la demi zone
de transition. Pour un passe-bande, diminuer la première fréquence de coupure de la demi zone de transition
et augmenter la seconde de la demi zone de transition, pour un rejecteur de bande, on fera l'inverse.
|H(f)|
Exemple :
fe / 2
1 2 π j f nT e
h (n) =
fe H ( f )e
− fe / 2
df
]
fc
1 1
e2π j f nTe df = e2π j f nTe
fc
On pose fc=(fp+fa)/2 h(n) =
f e − fc 2πjnfeTe
− fc
f
On normalise les fréquences par rapport à fe/2 (ou fe) c.a.d que l’on remplace partout fc par fc/(fe/2), on obtient
alors :
Tout aussi facilement, on peut déterminer les réponses impulsionnelles d’un passe-haut (1-H(f)), d’un
passe-bande (différence de 2 passe-bas) et d’un coupe-bande donnés.
FGE, USTHB [ assiakourgli@[Link] / [Link]
26
Traitement Avancé du Signal Master RT M1 2022/2023
Les valeurs de la réponse impulsionnelle idéale h(n) sont données au tableau suivant. Les fréquences fc
indiquées dans ce tableau (fréquences de coupure désirées) s'expriment également en fréquences normalisées
(divisées parfe/2).
h'N(n)=h(n).w(n)
N −1
1 si n ≤
si w( n) = 2 W ( f ) = sin( N πf ) H'N ( f ) = H( f ) *W( f )
0 ailleurs sin( πf )
La troncature temporelle introduit des ondulations et induit une zone de transition moins rapide
déterminée par la largeur du lobe principal. Un compromis est à faire entre la raideur et l'amplitude des
ondulations. Notons que cette méthode donne des ondulations de même amplitude dans la bande passante
et dans la bande atténuée. 1.4
Passe-bande avec Fen Rectangulaire (N=51)
Passe-bande avec Fen Hanning (N=51)
Remarque: Si N augmente, l'étendue des oscillations 1.2
Passe-bande avec Fen Rectangulaire (N=91)
Pour diminuer les oscillations : on utilise les fenêtres de pondération qui permettent d’obtenir de forte
atténuation des oscillations mais cela se fait au détriment de la largeur de la bande de transition qui devient
plus grande.
27
Traitement Avancé du Signal Master RT M1 2022/2023
largeur de la zone de transition ∆f (spécifiée aussi au départ) et du type de la fenêtre w(n),on déterminera la
longueur de la réponse impulsionnelle N.
Exemple : On veut synthétiser un filtre passe-bas de fréquence de coupure fc = fe/10 avec ∆f =fe/5 et une
ondulation en bande atténuée > 50 db (voir TP n°3) 1.5
sin( π n / 5 )
a- h ( n ) = 0 .2 0.5
πn / 5
16
Imaginary Part
16
0 X
b- On choisit w(n) comme étant la fenêtre de Hamming
sin( π n / 5 )
d- on calcule les valeurs h(n)= 0 .2 pour -8≤n≤8 -1.5
πn / 5
-1.5 -1 -0.5 0 0.5 1 1.5
Real Part
N -8 -7 -6 -5 -4 -3 -2 -1 0
h(n) -0,0399 -0,0456 -0,0329 0 0,0493 0,1064 0,1597 0,1974 0,2110
N 1 2 3 4 5 6 7 8
h(n) 0,1974 0,1597 0,1064 0,0493 0 -0,0329 -0,0456 -0,0399
0.3
Réponse impusionnelle h(n)
0.25 Réponse impusionnelle h'(n)
0.15
sin(π n / 5 ) 2πn
h' N (n) = 0.2 0.54 + 0.46 cos( ) 0.1
πn / 5 16 0.05
0
f- On translate le résultat de 8 échantillons.
-0.05
0 5 10 15 20
28
Traitement Avancé du Signal Master RT M1 2022/2023
N 0 1 2 3 4 5 6 7 8
h(n) -0,0399 -0,0456 -0,0329 0 0,0493 0,1064 0,1597 0,1974 0,2110
h'N(n) -0,0031 -0,0050 -0,0067 0 0,0255 0,0730 0,1325 0,1826 0,2023
N 9 10 11 12 13 14 15 16
h(n) 0,1974 0,1597 0,1064 0,0493 0 -0,0329 -0,0456 -0,0399
h'N(n) 0,1826 0,1325 0,0730 0,0255 0 -0,006 -0,005 -0,0030
1.4
1.5
Passe-bas avec Fen Rectangulaire (N=17)
Passe-bas avec Fen Hamming (N=17)
1.2
Passe-bas idéal
1
1
0.5
0.8
Imaginary Part
16
16
0
X
0.6
-0.5
0.4
-1
0.2
-1.5
0
0 200 400 600 800 1000 1200 1400 1600 1800 2000 -1 -0.5 0 0.5 1 1.5 2 2.5
Real Part
fc1=400 fe=4000Hz
On peut remarquer l'emplacement de zéros autour du premier 0 en 1 permet d'atténuer les lobes secondaire
et de maintenir une réponse cste autour du zéro.
La réalisation concrète d'un filtre numérique consistera en fait à matérialiser l'algorithme de calcul pour
la structure retenue. On aura la possibilité de travailler : Soit en logique câblée (assemblage d'organes
logiques, tels que portes, mémoires, etc ...),soit en logique programmée (organisation autour d'un processeur
de traitement du signal (DSP) ou, même, utilisation d'un microprocesseur(micro-ordinateur) standard).
29
Traitement Avancé du Signal Master RT M1 2022/2023
Le principe de définition du cahier des charges d'un filtre récursif (RII) se déroule de la même manière
que pour un filtre non récursif (RIF) mais la synthèse se fera de différente manière. L'intérêt d'employer un
filtre RII réside principalement dans la possibilité d'obtenir une bande de transition étroite pour un ordre
raisonnable bien qu'il présente d'une part, un risque d'instabilité du à une grande sensibilité numérique des
coefficients mais que l'on peut toutefois contrôler en déterminant une structure mieux adaptée, et d'autre
part, une variation de phase fortement non linéaire.
Rappelons que pour un filtre RII, la sortie s'exprime comme une combinaison linéaire d'un ensemble
fini d'éléments d'entrées et de sortie: N
N M
−i
b .z i
y ( n ) = bi x ( n − i ) − H ( z) = i=0
ai y ( n − i ) d’où M
i =0 i =1 1 + ai . z − i
i =1
La méthode directe consiste à placer des pôles et des zéros aux fréquences utiles mais la plus courante
(indirecte dite aussi de transposition) est l’utilisation des méthodes de synthèse des filtres analogiques
aboutissant à une fonction H(p) correspondant aux spécifications. Une fonction permettant le passage du
plan p dans le plan z (p = fct(z)) est ensuite utilisée pour obtenir H(z). Cette fonction doit maintenir la stabilité
du filtre analogique et maintenir, au mieux, les caractéristiques de la réponse fréquentielle H(f) du filtre
numérique.
Pour concevoir un filtre analogique, on peut employer des filtres passifs obtenus par combinaison de
résistances, de condensateurs et/ou de bobines ou encore utiliser des filtres actifs comportant un élément
amplificateur (transistor, AO, etc.) qui permet donc de modifier les amplitudes des signaux. Ainsi, un filtre
passif de base (de premier ordre) peut être composé d'une cellule RC ou d'une cellule RL.
- les filtres passifs (L,C, quartz, etc ;) sont utilisés pour les hautes fréquences
- les filtres actifs (R,C, ALI) utilisent des amplificateurs linéaires integrés (ALI) limités à 1Mhz
- les filtres à capacité commutés (R et C intégrés, ALI, Interrupteur commandé MOS) ont des fréquences
(<10 Mhz) programmables
FGE, USTHB [ assiakourgli@[Link] / [Link]
30
Traitement Avancé du Signal Master RT M1 2022/2023
Ci-dessous, sont donnés quelques exemples de structure pour des filtres passifs passe-bas de premier et
second ordre suivis de filtres actifs passe-bas de premier et second ordre aussi.
Pour obtenir des filtres d'ordre N plus élevé, on emploiera des filtres du premier et du deuxième ordremis
en cascade. Les pentes asymptotiques seront proportionnelles au nombre de cellules (+N ↑ et +∆f↓). A noter
également que les filtres passe-bas ou passe-haut peuvent avoir un nombre entier d'ordre (1, 2, 3...) tandis
que les filtres passe-bande ou coupe-bande ne peuvent qu'avoir un ordre pair (2, 4, 6, ...) car ils sont formés
de paires de cellules : 2 cellules RC ou une cellule RC et une cellule RL.
Pour synthétiser des filtres analogiques répondant à un gabarit donné, on choisira parmi un ensemble de
filtres analogiques testés et éprouvés donc connus pour leurs propriétés en termes de pente d’atténuation et
d’ondulation dans la bande passante et atténuée telles que :
Remarque : Le filtre de Bessel est intéressant pour filtrer les signaux large bande (modulations HF haut
débit : PSK, 8-PSK, OFDM…) en préservant les phases. Il n’a pas intérêt dans le domaine du filtrage
numérique.
31
Traitement Avancé du Signal Master RT M1 2022/2023
1. Un filtre de Butterworth est caractérisé par le fait que la réponse d’amplitude est maximalement plate
dans la bande passante et monotonement décroissante à partir d’une certaine fréquence (fréquence de
coupure).
L’amplitude d’un filtre analogique Passe-bas de Butterworth d’ordre N est définie par l’expression :
H (ω ) = H (ω ) =
1 ou 1
2N 2N
ω ω
1 + 1 + ε .
2
ωc ωc
ωc est la pulsations de coupure limites et ε est un paramètre de conception qui fixe la région de tolérance
dans la bande passanteδ1=1/(1+ε2)1/[Link] filtre analogique aura N pôles placés sur lecercle unité (à partie réelle
négative).
L’ordre N est déterminé par la zone de transition (ωa-ωp)et les ondulations permises en bande passante AP=20
log(1+δ1) et bande atténuée (Aa=-logδ2)
10 Aa Ap
1 1
10 log (10 − 1) /(10 − 1) log 2 − 1 / 2 − 1
Aa 10
log10 − 1
= δ2 δ1
N ≥=
N≥ 2 log(ω a / ω p ) 2 log(ω a / ω p )
2 log(ωa / ω p )
Remarque : Plus N est élevé, plus grande est la sélectivité (zone de transition ∆f étroite) mais un filtre RII
de fait de sa récursivité a une période d’initialisation avant le régime permanent qui augmente avec l’ordre
du filtre. En outre, son déphasage augmente avec l’ordre du filtre.
Déterminer le filtre et N.
N ≥=
( )
log (10 4 ) / 0.2589
= 4.80 ≈ 5
2 log(3 / 1)
2. Un filtre de Tchebychev est caractérisé généralement par une ondulation équilibrée dans la bande
coupée. La forme analytique du module de réponse fréquentielle d’un filtre analogique de Tchebychev (ou
Chebyshev) d’ordre N est donnée par :
32
Traitement Avancé du Signal Master RT M1 2022/2023
H (ω ) =
1
ω
1 + ε 2 .T N2
ωc
L’ordre N est déterminé par la zone de transition (ωa-ωp) et les ondulations permises en bande passante
AP=20 log(1+δ1) et bande atténuée (Aa=-logδ2) *
Aa
10
Ap
10 10 Aa 10
log 2.10 20 / 10 − 1
Aa Ap Ap
log 10 − 1
10 − 1 + 10 10 − 1
10 − 1 − 1
N≥ N≥
log ω a / ω p + (ω / ω p ) − 1
2
log ω a / ω p +
(ω a / ω p ) − 1
2
a
Exemple
Aa
10
Ap
log 2.10 20 . / 10 − 1
. 2.100 / 0.5089
N≥ = = 3.39 ≈ 4
log ω a / ω p +
(ω a / ω p ) − 1
2
3+ 8
Dans la pratique, on utilise trois valeurs d'ondulation, ε = 0.1 dB, 0.5 dB et 1 dB. Pour ces 3 valeurs HN(p) est
donnée :
33
Traitement Avancé du Signal Master RT M1 2022/2023
On rappelle que les fonctions de transfert HN(p) des filtres analogiques polynômiaux (Butterworth,
Tchebychev,Bessel, Cauer, etc.) sont données pour des fréquences de coupure normalisées et uniquement
pour desfiltres passe-bas [2].
1 p ω
Pulsation centrale
Passe-bas p = p / ωA Passe-bande p = + A ω A = ω A1ω A2
B ωA p
−1 Largeur de Bande
1 p ω A
Passe-haut p =ωA / p Coupe-bande p =
B ω
+
p
B = (ωA2 −ωA1)
A
A partir de l'expression de HN(p) normalisée, on dénormalise en remplaçant p par les valeurs données
au tableau précédent permettant d'aboutir à une fonction de transfert dénormalisée H(p).
La Méthode de la transformation bilinéaire a pour objectif de faire coïncider au mieux les domaines
analogique et numérique. Les filtres qui en sont dérivés sont plus stables que sont obtenus à travers l'emploi
de la méthode de la variance impulsionnelle. Mais, en contrepartie, elle introduit une distorsion sur l'axe des
fréquences.
On sait qu’un filtre analogique est caractérisé par sa fonction de transfert H(p) et un filtre numérique est
défini par sa fonction de transfert H(z) avec z= e2̟jfTe= epTe. Ce qui implique que p=Ln(z)/Te. Dans les deux cas,
cette transformation entraîne une relation non linéaire entre les fréquences fA du domaine analogique et les
fréquences fN du domaine numérique.
Pour conserver le caractère polynomial des fonctions de transfert, une approximation de ln(z) par les
séries de Laurent est employée :
1 − z −1 1 1 − z −1 3 1 1 − z −1 5
ln( z ) ≈ 2 + + + .........
1 + z −1 3 1 + z −1 5 1 + z −1
En ne conservant que le terme du premier ordre, on définit la transformation bilinéaire :
ln( z ) 2 1 − z −1 2 z − 1
p= ≈ =
Te Te 1 + z −1 Te z + 1
jω N Te
Sachant que p = jω = ln( z ) = 2 e − 1 2 e jω N Te / 2 (e jω N Te / 2 − e − jω N Te / 2 ) 2 j sin(ω N Te / 2)
jω T
= jω T / 2 jω T / 2 − jω T / 2
=
+1 +e Te cos(ω N Te / 2)
A
Te Te e N e
Te e N e
(e N e N e
)
Ainsi, on peut montrer qu’on obtient la correspondance entre fréquence analogique fA(ou pulsation
analogique ωA) et la fréquence numérique fN (ou ωN) par :
2 ω N Te 2 π f N
ωA = tg = tg
Te 2 Te f e
Cette équation montre qu'il n’y a pas égalité entre pulsation
analogique et pulsation discrète et que la relation les liant n’est pas fe π fN
fA = tg
non plus linéaire puisqu'il y a distorsion des fréquences, y compris π fe
du retard de groupe.
34
Traitement Avancé du Signal Master RT M1 2022/2023
Pour faire la synthèse d’un filtre numérique par la transformation bilinéaire, on procède comme suit :
- On définit les caractéristiques souhaitées du filtre numérique (fréq d’échantillonnage, de coupure, etc.)
2 π fN
- On calcule les pulsations analogiques ωA correspondant aux pulsations numériques ω A = tg
Te f e
- On détermine le gabarit du filtre analogique HN(p) normalisée d'ordre n (Chebyshev, Butterworth, etc.) qui
servira de modèle au filtre numérique et on écrit la fonction de transfert dénormalisée H(p) de ce filtre
analogique (qu’il faut recalculer en fonction de ωA).
2 1 − z −1
- On applique la transformation bilinéaireà H(p) en remplaçant p = , ce qui donne la fonction H(z).
Te 1 + z −1
Exemple 1: On désire concevoir un filtre passe-bas numérique de premier ordre à partir d’une fonction de
transfert d'un filtre RC dans le domaine continu (HN(p)=1/(1+p)). La fréquence de coupure désirée est fN=30
Hz et la fréquence d’échantillonnage est fe=150 Hz.
2
0.7265
2 π 2 1 Te
On calcule d'abord ω A = tg = 0.7265puis H ( p ) = =
Te 5 Te (1 + p ) p = p / ω A p + 2 0.7265
Te
2
0.7265
Te 0.7265( z + 1) 0.7265( z + 1)
H ( z ) = H ( p ) p = 2 z −1 = = =
Te z +1 2 z −1 2 ( z − 1) + 0.7265( z + 1) 1.7265 z − 0.2735
+ 0.7265
Te z + 1 Te
2 π fe
On peut remarquer que f A = tg = 0.7265 = 34.69 ≠ f N = 30
Te (2π ) 5 π
Exemple 2 : On désire réaliser un filtre numérique passe-bas du second ordre avec les caractéristiques
suivantes [18]: fréquence de coupure fN= 500 Hz (gain de 1), fréquence d'échantillonnage fe=5 kHz,
2 ωN Te 2 π 2
ωA = tg = tg = 0.325
Te 2 Te 10 Te
4
0.0114
0.1075ω A
2
0.1075 Te2
H N ( p) = H ( p) = H N ( p = p / ω A ) = 2 =
p + 0.1075 p + 1
2
p + 0.1075 pω A + ω A
2
2 4
p 2 + 0.035 p + 2 0.1056
Te Te
4
0.0114
Te2 0.0114
H ( z ) = H ( p ) p = 2 z −1 = =
z −1 z −1
2 2
Te z +1 2 z −1 2 2 z −1 4
+ 0.035 + 2 0.1056 + 0.035 + 0.1056
Te z + 1 Te Te z + 1 Te z +1 z +1
0.0114 ( z + 1) 0.0114 ( z + 1)
2 2
H (z) = =
( z − 1) + 0.035( z − 1)( z + 1) + 0.1056 (z + 1) 1.14 z − 1.7888 z + 1.0706
2 2 2
35
Traitement Avancé du Signal Master RT M1 2022/2023
z −1
z +1
Pour un filtre numérique passe-bas ou passe-haut, on peut déterminer H(z) directement à partir de HN(p) en
π fN
utilisant la pulsation normalisée Ω c = tg et remplacer p dans HN(p) comme suit:
f
e
1 z −1
- pour un passe-bas par p = soit H ( z ) = H N ( p ) 1 z −1
Ωc z +1 c
p=
Ω z +1
z +1
- pour un passe-haut par p = Ω c soit H ( z ) = H N ( p ) z +1
z −1 p=Ω
c
z −1
- pour un passe bande ou un rejecteur de bande, la détermination de ωA s'obtient par la pulsation centrale
ωA = ωA1ωA2
Remarque: Puisqu'il n'est pas possible de faire coïncider précisément la réponse fréquentielle du filtre
analogique (sur[0, +∞[ et numérique (sur [0, fe/2] la transformation bilinéaire aura pour objet de faire
coïncider au mieux les deux réponses. Cette méthode donne d'assez bons résultats à condition que les
fréquences caractéristiques ne soient pas trop proches de la demi-fréquence d'échantillonnage.
Cette approche peut aussi permettre d’améliorer les caractéristiques d’implantation d’un algorithme
en adoptant pour chacune de ses étapes un échantillonnage que l’on appelle critique dans le sens où la
fréquence de traitement est choisie égale à la fréquence de Nyquist. Ainsi la chaîne de traitement
correspondant à l’émetteur bande de base d’un modulateur DPQSK où l’on a trois fréquences de traitement :
la fréquence bit, la fréquence symbole et la fréquence d’échantillonnage [20]
Dans les systèmes considérés aux chapitres précédents, n’était considéré qu’une seule fréquence ou
cadence d’échantillonnage fe =1/Te. Mais il arrive que la fréquence maximale fmax soit diminuée ( filtrage passe-
bas ) ou augmentée (le cas d’une modulation). Dans les deux cas, il est plus judicieux d’adapter la fréquence
d’échantillonnage afin de minimiser le temps de calcul. Il s’agira d’une décimation ou sous-échantillonnage
FGE, USTHB [ assiakourgli@[Link] / [Link]
36
Traitement Avancé du Signal Master RT M1 2022/2023
Exemple : Supposons que l’on veuille synthétiser un filtre passe-bas avec les spécifications
suivantes :fp=100Hz fa=300Hz Aa>50 db fe=20kHz.
L’emploi de la fenêtre de Hamming nous donne N=331, soit un nombre de coefficients important. Or
le filtre occupe une petite portion de la bande passante (passe-bas avec fc=1% fe). On peut donc sous-
échantillonner et puis sur-échantillonner sans perte d’informations.
Rappelons que lorsqu’on souhaite échantillonner un signal, nous devons respecter le théorème de
Shannon à savoir prendre une fréquence d’échantillonnage fe/ 2 >fmax. Si ce signal a subi des traitements qui
ont eu pour résultat, entres autres, la modification de la fréquence maximale. On parlera alors de sous-
échantillonnage ou de sur-échantillonnage.
x(n) D xD (m )
1 0 2
Exemple : Le signal suivant est le résultat de la somme de 2 sinusoïdes de fréquences 100 et 300Hz. On
le décime de 2 et de 6. On remarque que lorsque fe=fe/6 il y a repliement puisque la condition de Shanon n’est
plus respectée.
37
Traitement Avancé du Signal Master RT M1 2022/2023
2 0.5
0.4
1
0.3
0
0.2
-1
0.1
-2 0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 -1000 -800 -600 -400 -200 0 200 400 600 800 1000
2 0.5
0.4
1
0.3
0
0.2
-1
0.1
-2 0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 -500 -400 -300 -200 -100 0 100 200 300 400 500
2 0.8
1 0.6
0 0.4
-1 0.2
-2 0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 -200 -150 -100 -50 0 50 100 150 200
0.4
1
0.3
0
0.2
-1
0.1
-2 0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 -1000 -800 -600 -400 -200 0 200 400 600 800 1000
2 0.5
0.4
1
0.3
0
0.2
-1
0.1
-2 0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 -500 -400 -300 -200 -100 0 100 200 300 400 500
1 0.5
0.4
0.5
0.3
0
0.2
-0.5
0.1
-1 0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 -200 -150 -100 -50 0 50 100 150 200
38
Traitement Avancé du Signal Master RT M1 2022/2023
Rappels :
∑ = ∑ En discret ∑ # = ∑$
! &
" '
telle que =
$ %
∑ ∑$
/012&
, - . %
' 3 /$
$
Fonction de transfert
∑$ -∑ - 3 3
/012 4
' '
$ . %
, ∑$ )5$. * avec 5$ =
4 /01
' '
$ . %
Réponse en fréquence
( .∑ # , 6 6 ∗ $ ∑$ )6 *
8 8
% $
, 6 ∑$ % )6 * Périodisation 9:
8 8 ;:
$ $ <
Exemples:
1
, B )5=% > * C )5= *C )5= *C )5=? > *D
4 4 4 4
4
> >
o Décimation de 2
, ∑. )5 . / * avec 5 , B ) /* C ) *D Y f H 6 C 6 )6 *I
4 /01 4 4
8
/ /
%
x(n) D xD (m)
n 2 . D DD …2D 2
39
Traitement Avancé du Signal Master RT M1 2022/2023
Théoriquement, si le signal d’origine a été bien échantillonné, la reconstruction parfaite est possible,
donc le sur-échantillonnage aussi. Cependant, la formule de reconstruction (avec les sinus cardinaux) est
compliquée à mettre en pratique. On fait appel à des idées plus simples [25]:
Le signal d’origine n’a pas de contenu fréquentiel au-delà de fe/2. Le processus d’interpolation est
susceptible d’en ajouter. On le fait donc suivre d’un filtrage passe-bas à la coupure fe/2 pour supprimer le
contenu ajouté. C’est la raison pour laquelle la méthode d’interpolation n’a pas beaucoup d’importance.
J J J
0.4
1
0.3
0
0.2
-1
0.1
-2 0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 -1000 -800 -600 -400 -200 0 200 400 600 800 1000
2 0.25
0.2
1
0.15
0
0.1
-1
0.05
-2 0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 -2000 -1500 -1000 -500 0 500 1000 1500 2000
2 0.08
1 0.06
0 0.04
-1 0.02
-2 0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 -4000 -3000 -2000 -1000 0 1000 2000 3000 4000
40
Traitement Avancé du Signal Master RT M1 2022/2023
Après sur-échantillonnage, la forme du spectre reste la même aucune énergie n’a été ajoutée au signal),
mais la fréquence d’échantillonnage est maintenant de Dfe. Le sur-échantillonnage fait apparaître ce que l’on
appelle des spectres miroirs et pour obtenir une signal dont la forme temporelle correspondrait au signal
x(n) que l’on aurait échantillonné à une fréquence Dfe.
Il faut supprimer ces spectres miroirs et donc filtrer par un filtres passe-bas le signal sur-échantillonné
y(m) avec un filtre ayant une bande atténuée à partir de fe/2
2 0.5
0.4
1
0.3
0
0.2
-1
0.1
-2 0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 -1000 -800 -600 -400 -200 0 200 400 600 800 1000
2 0.5
0.4
1
0.3
0
0.2
-1
0.1
-2 0
0 0.02 0.04 0.06 0.08 0.1 0.12 -2000 -1500 -1000 -500 0 500 1000 1500 2000
2 0.4
1 0.3
0 0.2
-1 0.1
-2 0
0 0.02 0.04 0.06 0.08 0.1 0.12 -4000 -3000 -2000 -1000 0 1000 2000 3000 4000
Nous pouvons observer l'effet l'interpolation suivie du filtrage passe-bas anti-miroirs (équivalent à
un 0 padding dans le domaine fréquentiel)
Exemples:
o Interpolation de 4 Y z X z=
41
Traitement Avancé du Signal Master RT M1 2022/2023
o Interpolation de 2 Y z X z
Remarque :
La décomposition polyphases intervient dans les bancs de filtres, elle permet une représentation simple
pour le sur et sous échantillonnage [21]. Quand un signal subit une décimation d’un facteur 2, cela revient à
séparer les échantillons avec un indice pair de ceux avec un indice impair, ou encore, à séparer deux phases
distinctes ayant des délais (aussi appelé écarts de phase) différents dans le vecteur [22]. Si on considère une
décimation d’un facteur M, le principe est exactement le même mais avec M phases distinctes.
Supposons que l’on veuille changer numériquement la fréquence d’échantillonnage d’un signal x(n) et
passer d’une fréquence d’échantillonnagef1 à une fréquence d’échantillonnage f2 où f2/f1=L/M. Cela conduit
donc à sur-échantillonner le signal x(n) d’un facteur L puis à le décimer d’un facteur M.
o Le sur-échantillonnage par L doit être suivie d’un filtre passe-bas coupant à fe/2L pour
supprimer les spectres dus au sur-échantillonnage où fe est la fréquence d’échantillonnage du
signal de sortie.
o La décimation par M doit être précédée d’un filtrage anti-repliement ayant une fréquence de
coupure de fe/2M, si fe est la fréquence d’échantillonnage du signal d’entrée.
On ne conservera alors qu’un seul des deux filtres passe-bas : celui dont la fréquence de coupure est la
plus basse
↑R HL ( ) HM( ) ↓#
y(n)=xL/M(n)
x(n)
Remarque : Si Let M sont premiers entre eux, décimation et interpolation sont commutatifs, ci-dessous
un exemple :de conversion où L = 5 et M = 2 [31]
Q% % que Q% R S# 1
En outre, d'après la théorie des nombres, si L et M sont premiers entre eux alors ils existent deux entiers
TU V W$ TU V U$
Exemple : On souhaite sur une réduction de fréquence d’échantillonnage d’un signal audio d’un
facteur 2/3 (passage de 48kHz à 32kHz). A noter que décimation et interpolation ne sont pas commutatives
sauf si L et M sont premiers entre eux [20].
42
Traitement Avancé du Signal Master RT M1 2022/2023
Remarques
Un filtre de décimation/interpolation à bande passante très étroite est complexe à mettre en œuvre (Ordre
élevé)
. Il peut être est remplacé par la mise en cascade de filtres plus simples. Les spécifications des filtres
individuels sont considérablement assouplis puisque la spécification globale du filtre est partagée entre
plusieurs filtres d'ordre inférieur.
Ainsi, pour un facteur de décimation ou interpolation, il est plus judicieux d'opérer par paliers (multi-étages).
Si on doit décimer de 15, cela nécessite un filtre d'ordre élevé. On peut alors procéder en décimant par 3 puis
par 5. Ce qui nous permettra d'employer 2 filtres de plus faible ordre.
43
Traitement Avancé du Signal Master RT M1 2022/2023
az − 1
1. Soit H(z) la fonction de transfert d’un SLIT causal avec : H ( z ) = avec a réel.
z−a
Déterminer les valeurs de a pour lesquelles H(z) correspond à un système stable. Prendre une valeur de a=0.5.
Représenter alors les pôles et zéros de la fonction, la région de convergence. Donner et tracer|H(f)|. Ce filtre
est-il à minimum de phase?
2. On suppose que le tracé des pôles et des zéros de ce système est le suivant :
- Est-ce un filtre RIF ou RII ? (Justifier votre réponse)
- Donner l’allure approximative de H(f)
- Déterminer H(z) puis déterminer et tracer h(n) 3
X
- A partir de h(n), étudier la stabilité, la causalité et l’invariance de ce filtre.
- Calculer et tracer H(f) pour au moins 3 valeurs
- Ce filtre possède-t-il un retard de groupe constant (justifier)
- Déterminer sa réponse pour une entrée échelon x(n)=U(n).
3. On considère que l'équation aux récurrences du système suivant est donnée comme suit :
y(n)= x(n) - x(n-2) - 0.878 y(n-2)
- Etudier la causalité et l'invariance de ce système
- Est-ce un filtre RIF ou RII ?
- Déterminer H(z) et donner le tracé des pôles et zéros, en déduire le rôle de ce filtre :
- Donner les allures approximatives de h(n) et |H(f)|
- Déterminer la réponse impulsionnelle h(n) et tracer la pour les 3 premières valeurs
- On suppose que fe= 6 kHz, quelle sera la sortie du filtre si l'on donne en entrée : un signal bruité par une
sinusoïde de 500 Hz puis un signal composé de 2 sinusoïdes l'une de 1500 Hz et l'autre de 2500Hz
4. Obtenir les coefficients d'un filtre à RIF passe-bas par la méthode de fenêtrage pour obtenir les
spécifications suivantes :
- Fréquence de coupure idéale : fc=1.75 kHz
- Largeur de transition : ∆f=0.5 kHz
- Atténuation en bande atténuée : A=-20 log10(δ)> 51 dB avec δ =min(δ1,δ2)
- Fréquence d'échantillonnage : fe=8 kHz
5. On souhaite approcher un filtre idéal passe-haut par un filtre à réponse impulsionnelle finie, synthétisé
par la méthode du fenêtrage. Ce filtre doit répondre aux spécifications suivantes :
Fréquence de coupure fc=2 kHz
Largeur de transition : ∆f=0.5 kHz
Atténuation en bande atténuée : A=-20 log10(δ)> 40 dB avec δ =min(δ1,δ2)
Fréquence d'échantillonnage : fe=8 kHz
- Déterminer l’expression mathématique exacte de h’(n).
- Calculer h’(0) et tracer approximativement h’(n).
- Quel est l’intérêt du fenêtrage et quel est son inconvénient ?
- Quel est l’inconvénient de cette technique de synthèse des filtres?
- Citer un avantage et un inconvénient de la synthèse par des filtres RIF de même que par les RII.
44
Traitement Avancé du Signal Master RT M1 2022/2023
7. On considère le filtre passe-bas analogique normalisé à un pôle unique HN(p)=1/(p+1). On considère que
la fréquence de coupure normalisée à -3db se produit pour fc=fe /10 . Par transformation bilinéaire, trouver le
filtre numérique passe-bas équivalent H(z).
8. Soit la fonction de transfert passe bas dénormalisée H(p) suivante : (p+0.1)/((p+0.1)2+ωa2). Sachant que
ωa=4, utiliser la méthode de la transformée bilinéaire pour transformer ce filtre en numérique. Le filtre
numérique aura sa résonnance pour ωN = π/([Link]).
10. Soit le signal suivant, on veut le décimer par 2 puis l'interpoler par 2 également
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
- Donner le tracé du signal décimé puis interpolé
- On considère que la fréquence du signal original est 10kHz, donner la fréquence à la sortie du décimateur
puis celle de l'interpolateur.
11. On considère le signal x(n) dont la TFTD est donnée comme ci-contre ,
↓2 ↓2
- Préciser le rôle de H(z) dans chaque cas
H(z)
x(n) y(n) x(n) y(n)
↓4 H(z) ↓4
x(n) y(n) x(n) y(n)
↑2 ↑2 H( )
x(n) y(n) x(n) y(n)
↑4 ↑4
-
H( )
y(n)
assiakourgli@[Link] / [Link]
FGE, USTHB [ x(n)
x(n) y(n)
45
Traitement Avancé du Signal Master RT M1 2022/2023
12. Avant de convertir un signal numérique x(n) (ci-dessous signal et sa TFTD),en signal analogique, on
décide de l’interpoler par D.
x(n) D y(n)
|X(f)|
2 0.5
0.45
1.5
0.4
1
0.35
0.5
0.3
0 0.25
0.2
-0.5
0.15
-1
0.1
-1.5
0.05
-2 0
0 0.002 0.004 0.006 0.008 0.01 -2000 -1500 -1000 -500 0 500 1000 1500 200
- Tracer le signal et sa TFD dans le ca ou l'on fait suivre l'interpolation d'un filtre passe-bas à fe/2
13. Un signal audio est enregistré a été transmis à une fréquence de 30 kHz, à la réception avant de le convertir
en analogique, on opère un changement de fréquence telle que la nouvelle fréquence soit 45Hz.
Solutions
Les coefficients du filtre sont symétriques, il suffira par conséquent de calculer les valeur de h(0)à h(26).
46
Traitement Avancé du Signal Master RT M1 2022/2023
5. ∆f=0.5 fc= 0.125 h( n) = − f c sin(π n f c ) pour k≠0 et h(0)=1-fc A>40 Hanning N=51
πnf c
6. y(n)=b3(x(n-3)+x(n-4))+b2(x(n-2)+x(n-5))+b1(x(n-1)+x(n-6))+b0(x(n)+x(n-7))
Type II , pas adéquat pour un passe-haut (un zéro en -1)
7. ω A = 2 tg π = 0.65 H ( p ) = H N ( p ) =
ωA
=
0 .65 / Te
Te 10 Te p /ωA
p +ωA p + 0 .65 / Te
−1
2 1− z 0.65/ Te 0.245(1+ z −1 )
p= H ( z) = H ( p) 2 1−z −1 = =
Te 1+ z −1 p=
Te 1+ z−1 2 1− z −1 1− 0.509z −1
+ 0.65/ T
Te 1+ z −1
e
1 − z −1
4 + 0.1
2 π 1 − z −1 1 + z −1 0.125 + 0.006z −1 − 0.118z −2
8. ωa = 4 = tg Te = 0.5 p = 4 H ( z ) = =
Te 4 1 + z −1 1 − z −1
2
1 + 0.0006z −1 − 0.950z −2
4 −1
+ 0.1 + 16
1+ z
9.
-600 -450 -300 -150 -75 -50 50 75 150 300 450 600 Hz
47
Traitement Avancé du Signal Master RT M1 2022/2023
H(z) ↓2
x(n) y(n)
-75 -50 50 75 Hz
↓4
Période de -37,5 à 37,5 Hz
H(z)
x(n) y(n)
-37,5 37,5
↑2 H( )
x(n) y(n)
-300 -50 50 30
FGE, USTHB [ assiakourgli@[Link] / [Link]
48
Traitement Avancé du Signal Master RT M1 2022/2023
↑4 H( )
x(n) y(n)
12. ’
’
f1 = 30 fe = 90 f2= fe /2=3 f1/2=45
↑3 H3( ) H2( ) ↓2
x(n) y(n)=x3/2(n)
’ ’
fc= f1/2= fe /6=15 fc= fe /4=22.5
49
Traitement Avancé du Signal Master RT M1 2022/2023
Exercice supplémentaires
0.35
h(n) 0.25
0.2
- Déterminer les pôles et zéros de ce filtre puis donner leur tracé. En déduire un tracé approximative de |H(f)|
- Calculer et tracer |H(f)| puis en déduire le tracé du module de la TFD pour N=6.
3. Au cours de la transmission d'un signal numérique (échantillonné à une fréquence de 2,5 kHz) , il a été
affecté d'un bruit localisé entre les bandes de fréquence 350 Hz et 550 Hz. On veut éliminer le bruit par
l'emploi d'un filtre RIF possédant une possédant une bande de transition ∆f=100 Hz. Concevoir ce filtre par
la méthode du fenêtrage. On souhaite une atténuation en bande atténuée A=-20 log10(δ)> 20 dB.
- Tracer H(f) idéal
1.4
- Déterminer h(n) et l'ordre N du filtre 6
Imaginary Part
0.8
14
fréquence de 0 à fe/2 et des pôles et zéros d’un 0
0.6
filtre numérique synthétisé par la méthode des -2
fenêtres sont donnés comme suit : 0.4
50
Traitement Avancé du Signal Master RT M1 2022/2023
7. Soit le signal suivant :h(n) = sin(2*pi*150*n) + sin(2*pi*350*n) et le module de la TFD dont les tracés
respectifs sont donnés ci-dessus : (voir exo 9)
- Déterminer Te.
- Quel est l’intérêt de l’interpolation ?
- Tracer le signal (les 20 premières valeurs) et la TF obtenue pour une interpolation de 2 puis de
4(sans filtrage passe-bas)
- Pourquoi emploie-t-on un filtre passe-bas lors l’interpolation, où le place-t-on ?
- Quel est le type du filtre de Tchebeychev employé dans ce cas et pourquoi ?
- Donner la décomposition polyphases finale pour une interpolation de 2 en donnant l’expression
des filtres polyphases
x(n) 3 y(n)
-fmax fmax=fe/5
- Tracer la TFTD obtenue après décimation (fe=30kHz)
- On veut placer un filtre anti-repliement, ou doit-on le placer ? (justifier), tracer la TFTD à nouveau
- Donner la décomposition polyphases initiale et finale en donnant l’expression de chaque filtre.
9. Soit le signal suivant :h(n) = sin(2*pi*150*n) + sin(2*pi*350*n) et le module de la TFD dont les tracés
respectifs sont donnés ci-après :
2 0.7
1.5 0.6
1
0.5
0.5
0.4
0
0.3
-0.5
0.2
-1
-1.5 0.1
-2 0
0 0.005 0.01 0.015 0.02 -1000 -500 0 500 1000
- Tracer le signal et la TF obtenue pour une décimation de 2 puis de 4(avec filtrage passe-bas)
- Donner la décomposition polyphases finale pour une décimation de 4 en donnant l’expression
des filtres polyphases
- Quel est l’intérêt de la décomposition polyphases
- Pourquoi emploie-t-on un filtre passe-bas lors de la décimation
51
Traitement Avancé du Signal Master RT M1 2022/2023
But du TP : Dans ce TP, on teste deux méthodes différentes pour synthétiser un filtre numérique : un filtre
RIF par la méthode des fenêtres et un filtre RII par la méthode de la transformation bilinéaire.
1. Rappels
Un filtre de réponse impulsionnelle finie (RIF) possède une fonction de transfert polynomiale. Il ne
peut pas être obtenu par transposition d’un filtre continu, comme cela est fait pour les filtres RII. Les filtres
RIF présentent l’inconvénient de nécessiter un grand nombre de coefficients pour obtenir les mêmes
caractéristiques fréquentielles. Mais par contre, ils sont inconditionnellement stables. On peut synthétiser des
filtres RIF à phase linéaire, c’est-à-dire à temps de propagation de groupe constant.
L’intérêt des filtres récursifs (RII) est leur faible coût en calcul. Avec très peu de pôles et zéros on peut
assurer la plupart des réponses fréquentielles dont on peut avoir besoin dans les différentes applications. Les
inconvénients des filtres récursifs sont : leur non-linéarité en phase linéaire (phase linéaire : temps de
propagation constant pour tout fréquence), et leur instabilité numérique. En effet, le filtre étant rétroactivé,
les erreurs de précision numérique deviennent une question d'importance, car elles peuvent s'amplifier et
devenir dehors contrôle, d'abord dans la forme de bruit, mais éventuellement dans la forme d'instabilité. A
noter que les filtres RII peuvent être conçus par des méthodes semblables à ceux utilisé pour les filtres
analogiques.
1.4
2. Démo
# -*- coding: utf-8 -*-
52
Traitement Avancé du Signal Master RT M1 2022/2023
[Link]()
#Effet du fenêtrage
from [Link] import windows
[Link](4)
N=41;
A=[Link]([[Link](N), [Link](N),[Link](N) ]);
i=1;
for Fen in A:
n1 = [Link](-(N//2),N//2+1)
b1 = fcn*[Link](n1*fcn)*Fen;
f,Hf1= [Link](b1,a,L,fs=fe);
[Link](2,3,i); [Link](b1); [Link]('Filtre * Fen %d'%i);
[Link](212); [Link](f, [Link](Hf1),label='Fen %d'%i)
[Link]('Module du Filtre'); [Link](True); [Link]('Fréquence (Hz)');[Link]('Amplitude')
[Link]()
i+=1
#Effet du fenêtrage (affichage en db)
[Link](5)
fc=0.2;fe=2; fcn=2*fc/fe; N=41;
a = [Link]([1]);
L = 256;
A=[Link]([[Link](N), [Link](N),[Link](N) ]);i=1;
for Fen in A:
n2 = [Link](-(N//2),N//2+1)
b2= fcn*[Link](n2*fcn)*Fen;
f,Hf1= [Link](b2,a,L,fs=fe);
[Link](2,3,i); [Link](b2); [Link]('Filtre * Fen %d'%i);
[Link](212); [Link](f, 20*np.log10([Link](Hf1)),label='Fen %d'%i)
[Link]('Module du Filtre'); [Link](True); [Link]('Fréquence (Hz)');[Link]('Amplitude')
[Link]()
i+=1
# Transformation bilinéaire
fp=195;fa= 205; fe=1000;
FGE, USTHB [ assiakourgli@[Link] / [Link]
53
Traitement Avancé du Signal Master RT M1 2022/2023
att_p=1; att_a=40;
wpn = fp*2*[Link];
wan = fa*2*[Link];
wpa = 2*fe*[Link]([Link]*fp/fe) #Compensation de pulsation
[Link](7);
z,p=zplane(Bz1,Az1);
3. A faire en TP
# # t = [Link](0,N)*Te;
# # f0 = 1500; b = 20*[Link](2*[Link]*f0*t);
# # Xb =x +b ;
# # [Link](np.int8(Xb), fe)
"""
----------------------A faire--------------------------------------
54
Traitement Avancé du Signal Master RT M1 2022/2023
55
Traitement Avancé du Signal Master RT M1 2022/2023
[Link](2)
freq = [Link](-NF/2,NF/2)*fe/NF;
[Link](321); [Link](t, x); [Link](322); [Link](freq, [Link](TFx));
x1 = [Link](x,2*N); t1 = [Link](0.0, N*Te, Te/2); TFx1 = [Link](x1,NF);
TFx1 = [Link](TFx1); freq = [Link](-NF/2,NF/2)*2*fe/NF;
[Link](323); [Link](t1, x1); [Link](324); [Link](freq, [Link](TFx1));
x2 = [Link](x,4*N); t2 = [Link](0.0, N*Te, Te/4); TFx2 = [Link](x2,NF);
TFx2 = [Link](TFx2); freq = [Link](-NF/2,NF/2)*4*fe/NF;
[Link](325); [Link](t2, x2); [Link](326); [Link](freq, [Link](TFx2));
[Link](4)
freq = [Link](-NF/2,NF/2)*fe/NF;
[Link](321); [Link](t, x); [Link](322); [Link](freq, [Link](TFx));
x1 = [Link](x,2); t1 = t[::2]; TFx1 = [Link](x1,NF); TFx1 = [Link](TFx1);
freq = [Link](-NF/2,NF/2)*fe/(2*NF);
[Link](323); [Link](t1, x1); [Link](324); [Link](freq, [Link](TFx1));
x2 = [Link](x,6); t2 = t[::6]; TFx2 = [Link](x2,NF); TFx2 = [Link](TFx2);
freq = [Link](-NF/2,NF/2)*fe/(6*NF);
[Link](325); [Link](t2, x2); [Link](326); [Link](freq, [Link](TFx2));
56
Traitement Avancé du Signal Master RT M1 2022/2023
#----------------------A faire--------------------------------------
#1
# Tester le sous échantillonnage avec et sans filtrage pase-bas
# Afficher le signal et sa TFD dans chaque cas et zoomer sur une meme zone d'une durée de 50 ms et
commenter
# Ecouter le signa obtenu à chaque fois
# Perçoit-on une déformation du son ? A quoi est-elle due ?
# Jusqu'à quel facteur de décimation peut-on aller sans déformation
# Que se passet-il si on garde la même fe (originale) à l'écoute?
#2
# Tester le sur échantillonnage avec et sans filtrage passe-bas
# Afficher le signal et sa TFD dans chaque cas et zoomer sur une meme zone d'une durée de 50 ms ou
20ms et commenter
# Ecouter le signa obtenu à chaque fois ([Link] avec la nouvelle fe en multipliant par facteur d'Interp)
# Perçoit-on une déformation du son ? pourquoi?
# Que se passet-il si vous garde la même fe à l'écoute?
#3
# On veut paser à un signal dont la fréquence d'échantillonnage est 1.5*fe
# Proposer 2 solutions
# Laquelle vous permet de conserver le plus les valeurs originelles ?
57
Traitement Avancé du Signal Master RT M1 2022/2023
Les processus aléatoires (stochastiques) décrivent l'évolution d’une grandeur aléatoire en fonction du
temps (ou de l’espace). On peut définir un processus stochastique comme étant une famille de variables
aléatoires indexées par le temps définies sur un même espace de probabilités Ω. Un processus aléatoire ne
possède pas de représentations temporelles analytiques. Chaque signal aléatoire observé représente une
réalisation particulière de ce processus. Le signal de parole, le signal radar, l'électrocardiogramme, l'électro-
encéphalogramme, les signaux sismiques en sont des exemples. Rappelons qu'un signal de réception -
constitué d'un signal informatif (aléatoire ou déterministe), d'un signal aléatoire d'interférence et de bruit lié
au canal de transmission- est aléatoire par nature.
Exemple 1 : Si l'on prend plusieurs résistances identiques (de même valeur) et que l'on mesure la tension. On
trouvera une valeur non nulle due à l'agitation thermique des électrons libres dans la résistance. Les tensions
mesurées au cours du temps sur l'ensemble des résistances fourniront plusieurs signaux aléatoires tous
différents qui constituent le processus aléatoires.
18
0
16
R U .................... 0
14
0
12
0
10
p( x; t i ) = p ( xi )
0 10 20 30 40 50 60 70 80 90 100
la densité est ti tj
- Deux instants ti et tj permettent de définir deux variables aléatoires xiet xj on peut définir des densités de
probabilités conjointes par: ( )
p x; t i , t j = p ( xi , x j )
Exemple 2 : Signal sinusoïdal à phase aléatoire X (t , ϕ ) = a cos(ω t + ϕ ) avec une phase φ uniformément
1
0.6
0.4
- φ est à valeur réelle continue 0.2
-0.4
-0.6
- Un signal particulier X(t,φi) est déterministe -0.8
58
Traitement Avancé du Signal Master RT M1 2022/2023
1. Espérances mathématiques
Statistiques de 1er ordre: Elles permettent de décrire le signal à un instant donné. Le processus aléatoire devient
une simple variable aléatoire que l'on peut décrire à l'aide de moments à condition que sa probabilité soit
connue ∀ ti.
+∞
- moyenne statistique: µ x (t i ) = E{x(t i )} = E{xi } = x(t ) ⋅ p(x, t )dx
i i
{ }
−∞
- variance: σ x2 (t i ) = E x(t i ) − µ x (t i )2 = E {xi2 }− µ x (t i )2
+∞
{ }
Où E xi = x(t i ) ⋅ p ( x, t i ) dx est la valeur quadratique moyenne (on parle aussi de puissance) et
2 2
σestx (t i )
dit écart-type−∞
Exemple 1: Soit le processus stochastique x(t) défini par :x(t) = a + b.t où aet b sont des v.a. dont les probabilités
sont connues, alors :
µ X (t ) = E[a + b t i )] = µ a + µ b t
Reprenons l’exemple 2 : Pour un instant donné tk, on peut calculer des moments statistiques de la variable
aléatoire X(tk,u)
- Moyenne :
2π 2π
1
µ X (t ) = E[ X (t )] = f ϕ (ϕ )acos(ω t + ϕ )dϕ = a cos(ω t + ϕ )dϕ = 0
0
2π 0
-Variance :
2π
a2
σ (t ) = E[ X (t )] − µ (t ) = E[ X (t )] =
2
X
2 2
X
2
f ϕ (ϕ )a cos (ωt + ϕ )dϕ =
2 2
0
2
Supposons, maintenant, que ϕ est déterministe et que a est aléatoire de densité p(a), de moyenne µa et
variance σa2, alors :
On dit que le processus est connu à 2 instants t1 et t2, si ∀t1, t2, la probabilité conjointe est connue. Dans la
mesure où l'espérance fait intervenir le produit de deux variables aléatoires (2 instants), on parlera de
statistiques d'ordre 2. on peut estimer les statistiques suivantes:
59
Traitement Avancé du Signal Master RT M1 2022/2023
t − τ. 2
On peut étendre ces notions pour mesurer ce lien entre deux processus aléatoires par :
3. Processus stationnaires
Beaucoup de processus aléatoires observés en pratique ont des propriétés statistiques qui ne dépendent pas
du temps ou l'observation est faite. On dit que ce sont des processus stationnaires.
On désigne par stationnaire les processus dont les caractéristiques statistiques sont indépendantes de l'origine du
temps. Un processus aléatoire est dit stationnaire au sens strict lorsque toutes ces caractéristiques statistiques
c'est à dire tous ses moments à tout ordre sont indépendants de l'origine du temps.
⇔ p(x, ti ) = p( x );
µ x (ti ) = µ x
Rxy (t1, t2 ) = Rxy (τ ) avec τ = t2 − t1
L
60
Traitement Avancé du Signal Master RT M1 2022/2023
La stationnarité ne signifie pas pour autant que le processus est indépendant du temps mais plutôt que
ses propriétés statistiques ne dépendent pas du moment auquel on commence à les estimer. Ainsi, la
température est un processus non stationnaire alors que le lancer de dé est un processus stationnaire.
On peut considérer que de nombreux phénomènes sont approximativement stationnaires sur des durées
d'observation finies. Comme nous n'avons pas toujours accès à p(x,ti), on se contentera de la stationnarité au
sens large(SSL) lorsque seul µx et Rx sont indépendants de t.
Ainsi :
- Un processus stationnaire est dit stationnaire d'ordre 1 si sa moyenne et sa variance sont constantes et
donc indépendantes de tout décalage temporelle :
et µ x (t ) = µ x = cste σ x2 (t i ) = σ x2
- Un processus stationnaire est dit stationnaire d'ordre 2 si ses statistiques d'ordre 2 ne dépendent que de
l'écart temporel entre les deux instants t1 et t2:
R x (t1 , t 2 ) = R x (τ ) avec τ = t 2 − t1
Exemple d'application:
Montrer que le processus x(t) = [Link](2̟t) + [Link](2̟t) où a, b sont des v.a décorrélées de moyenne nulle et
de variance unité est stationnaire au sens large.
3.5
2.5
1.5
Cas discret :
61
Traitement Avancé du Signal Master RT M1 2022/2023
R x (0 ) Rx (1) R x ( 2) .............. Rx ( N − 1)
R (1) R x ( 0) Rx (1) ............. Rx ( N − 2)
x
Rx (n1 , n2 ) = Rx ( 2)
Rx (1) R x ( 0) ............. Rx ( N − 3)
............ ............. ............. ............
Rx ( N − 1) R x ( N − 2) Rx ( N − 3) ............. Rx (0)
Rx (n1 − n2 ) = Rx (k)
Ainsi, la matrice de corrélation d’un processus stationnaire discret est une matrice Toeplitz carrée.
Une matrice carrée est dite Toeplitz si tous les éléments d’une même diagonale ou sous-diagonale sont égaux.
On voit directement que c’est le cas ici. Cette propriété est directement liée à la propriété de stationnarité (au
sens large) du processus.
Puissance et DSP : Les personnes qui conversent (dans un café ou en cours) génèrent un signal aléatoire qui
selon son volume global possèdera une puissance. Si X(t) est SSL, on peut calculer l’espérance de la puissance
instantanée par:
PX = E ( x(t ) 2 ) = R x (0 ) = µ x + σ x
2 2
Quant au spectre pour un signal aléatoire, il faut noter que chaque réalisation fournira un spectre
différent (Voir exemple ci-dessous)
62
Traitement Avancé du Signal Master RT M1 2022/2023
Or, nous avons déjà un outil statistique qui contient une information unique lorsque le processus est
considéré SSL: c’est la fonction de corrélation statistique dont la TF fournira une information sur la
distribution fréquentielle de la puissance moyenne du signal.
On définit, alors, la densité spectrale de puissance d’un signal aléatoire X(t) stationnaire au sens large
comme la fonction de la fréquence f donnée par la TF de la fonction de corrélation statistique du signal:
+∞
S x ( f ) = R x (τ ) e -2πjfτ dτ
−∞ +∞
La puissance moyenne totale du processus est : Px = S x ( f ) df = R x (0)
−∞
Comme le spectre représente une moyenne sur l'ensemble des réalisations possibles du processus,
une réalisation particulière peut toujours avoir un spectre de puissance différent de Sx(f).
Bruit blanc : Un bruit blanc est un signal aléatoire stationnaire au second ordre dont la densité spectrale de
puissance est constante sur tout l'axe des fréquences.
Il est défini par : R x (τ ) = σ x2δ (τ ) ce qui implique que sa moyenne sera nulle (µx=0) et qu'il est décorrélé.
Par analogie avec la lumière blanche, on dénomme un tel signal bruit blanc. Il existe donc différents
type de bruit blanc en fonction de la variable aléatoire décrivant le bruit. Les plus répandus sont le bruit
blanc gaussien où la variable aléatoire suit une loi gaussienne et le bruit blanc uniforme où la variable
aléatoire suit une loi uniforme.
63
Traitement Avancé du Signal Master RT M1 2022/2023
4. Processus ergodiques
Il n'est pas toujours possible de réaliser un nombre suffisant de mesures pour établir les propriétés
statistiques d'un processus alé[Link] est plus facile de lancer un dé 1000 fois que de réquisitionner 1000
personnes pour lancer l000 dés. Tout comme pour caractériser le bruit thermique, on réaliserait plusieurs
mesures en usant de la même résistance au lieu d'en employer 1000. Cela signifie que l'on assimilera les
résultats obtenus sur une réalisation à ceux obtenus pour un instant donné tipour différentes réalisations.
Un processus aléatoire stationnaire est dit ergodique lorsque les valeurs moyennes statistiques et temporelles
sont identiques. µ x = µ x ; Rx (τ ) = Rx (τ )L
T
1 2
où x = lim x(t ) ⋅ dt et Rxy (τ ) = x(t ) ⋅ y * (t − τ )
T → +∞ T T
− 2
D'un point de vue pratique lorsqu'un phénomène aléatoire est ergodique et stationnaire, on peut
mesurer avec un seul appareil fiable à partir de n'importe quel instant.
T T
1 2 1 2
µ x = lim x(t ) ⋅ dt µ x = µ x Rx (τ ) = lim x(t ) x(t − τ ) ⋅ dt Rx (τ ) = Rx (τ )
T → +∞ T T → +∞ T
−T 2
−T 2
T −T / 2 T −T / 2
2
Il est donc ergodique d'ordre 1
Si le processus est ergodique Sx(f) représente le spectre de puissance de n'importe qu'elle réalisation xi(t):
+∞
Px = S x ( f ) df = lim
1 T /2
−∞ T → +∞ T
−T / 2
xi (t ) 2 dt Où xi(t) est une réalisation de x(t).
Dans le cas où le signal x(t) est ergodique, stationnaire, on détermine Sx(f) par le calcul de la DSP temporel
sur sur T:
S x ( f ) = TF {R x (τ )} = lim
1 2
XT ( f ) Où XT(f) est la TF de x(t) limité à l'intervalle [0 T]
T →∞ T
5. Exemples de Processus
Il existe de nombreuses classes de processus particuliers : les processus de Markov (et notamment les chaines
de Markov quand t est discret), les martingales, les processus Gaussiens, les processus de Poisson.
• Poisson
64
Traitement Avancé du Signal Master RT M1 2022/2023
La distribution de Poisson est destinée à modéliser la fréquence des événements au cours d'un intervalle de
temps fixe. Étant donné un intervalle de temps de longueur, t, et le taux d'arrivée des événements, lambda,
c d.
le nombre attendu d'événements au cours de cet intervalle de temps est :
Pour que le processus de Poisson soit réalisable, les hypothèses suivantes doivent être remplies :
1. La probabilité de succès sur une très courte période de temps est le paramètre lambda multiplié par
cette période de temps.
2. La probabilité que plus d'un événement de succès se produise dans l'intervalle de temps défini n'est
pas significative. En d'autres termes, la probabilité que plus d'une expérience réussisse dans un intervalle
de temps fixe est très faible, et donc non importante ou non significative.
3. La probabilité qu'un événement réussi se produise pendant un intervalle de temps défini ne dépend
pas de ce qui s'est passé précédemment. Autrement dit, chaque expérience réussie est indépendante de
l'expérience précédente.
Le processus de Poisson est connu en statistique comme un processus stochastique qui essaie d'enregistrer
des événements très improbables en temps continu.
Exemples: panne, accident nucléaire, défaut de fabrication dans une chaine. Il est également utilisé en
Télécommunications pour estimer le nombre d'appels dans un laps de temps donné.
• Markovien
L’axiome de Markov traduit que la probabilité de n’importe quel comportement futur, le présent étant connu,
n’est pas modifié par toute connaissance supplémentaire du passé.
Un processus est markovien, si pour tout instant u, pour toute valeur donnée Xu=x, la probabilité pour que
le processus prenne la valeur y à un instant quelconque t ultérieur (t > u), ne dépend pas des valeurs prises
par le processus avant l'instant u. Le processus est dit sans mémoire.
Les chaînes de Markov sont des processus markoviens qui évoluent dans le temps discret selon des
transitions probabilistes. Elles vont donc décrire un ensemble de phénomènes aléatoires selon le principe des
processus stochastiques.
Exemples : la hauteur du niveau d'un barrage à un instant t, le cours du pétrole à la fin de la journée
• Gaussien
Les processus gaussiens jouent un rôle crucial dans la théorie des processus stochastiques du fait que :
- Certains processus stochastiques peuvent être approchés par des processus gaussiens,
- Théorème central limite et approximations : Chaque fois qu’un phénomène peut être considéré comme la
résultante d’un grand nombre de causes aléatoires indépendantes (ex: notes obtenues à l'examen, poids d'un
enfant), on peut résumer que ce phénomène suit une loi de distribution normale: c’est le théorème central
limite. Ainsi, la distribution statistique de la somme de n V.A. indépendantes et de même loi tend vers la loi
Normale quand n tend vers l’infini.
FGE, USTHB [ assiakourgli@[Link] / [Link]
65
Traitement Avancé du Signal Master RT M1 2022/2023
Un processus (Xt) réel est gaussien si : ∀n ∀, t1;…….; tn, le vecteur (Xt1=x1 ;……;Xtn=xn) est gaussien.
n
Y= α X k k = αT X
k=1
Rappelons que toutes les probabilités marginales d'un processus gaussien sont gaussiennes et que toute
combinaison linéaire de marginales d'un processus gaussien est gaussienne.
1
1
( )
− X −µx T CX −1 X −µx( )
p( x1, x2 ,....,xn ) = e 2
(2π )n det(CX )
X1 µ X 1 X 1 µ X 1
. . . .
Où X = . , µX = . , XC = . − .
. . . .
X n µ Xn X n µ Xn
Exemple : Un mouvement brownien (standard) réel est un processus gaussien centré à trajectoires
continues. C'est un processus à accroissements indépendants stationnaires et gaussiens.
Les quantités E{x} ,E{x2}, ,Sxx(f),... sont impossibles à calculer sur un ordinateur car elles nécessiteraient un
nombre de points infini. Sur calculateur, on ne dispose que d'une séquence discrète et finie de N points. En
réalité, on calcule des estimées de ces grandeurs en supposant généralement le signal stationnaire et
ergodique. On remplace alors le calcul des moyennes statistiques par des moyennes temporelles.
Nous calculerons la DSP sur un nombre de point finis. Ainsi, la transformée de Fourier calculée sur la
séquence finie sera donc convoluée avec le spectre de la fenêtre rectangulaire c'est-à-dire un sinus cardinal.
Rappelons que les propriétés spectrales du sinus cardinal ne sont pas bien adaptées à l'analyse spectrale du
signal (lobes secondaires importants). Pour y remédier, on emploie des fenêtres de pondérations (Hamming,
Hanning, Kaiser,etc.)
N −1
1
Y ( f ) avec Y ( f ) = yk e − 2πjfk
2
La méthode du périodogramme à prendre la TF d'une réalisation : S xx ( f ) =
ˆ
N k =0
N −1 N −1 N −1 N −1
1 1
Sˆ xx ( f ) =
N
yl e −2πjfl ym e 2πjfm
l =0 m =0 N
y y
l =0 m=0
l m e −2πjf ( l − m )
( )
N −1 N −1 N −1
1
D’où =
N
y y
k =0 l =k
l l −k e −2πjfk = Rˆ yy (k )e − 2πjfk = TF Rˆ yy (k )
k =0
{ } N −1 N −1 N −1
1
E Sˆ xx ( f ) = E Rˆ yy (k )e − 2πjfk = E y y l l −k e − 2πjfk
k =0 N k =0 l = k
66
Traitement Avancé du Signal Master RT M1 2022/2023
{ } 1 N −1 N −1 N −k
N −1
E Sˆxx ( f ) == E{yl yl −k }e−2πjfk = Ryy (k ) e−2πjfk = S yy ( f ) * N sin c( fN )2
N k =0 l = k k =0 N
Afin de diminuer la variance de cet estimateur, on peut utiliser un périodogramme moyenné. Cela consiste
à séparer le signal en K tranches (de longueur N/K), à calculer le périodogramme sur chaque tranche et à faire
la moyenne. Du fait des K moyennages, la variance est presque divisée par K : néanmoins, les tranches étant
plus courtes, la résolution diminue(∆f=fe/N→fe/K). En pratique, un taux de recouvrement de 40% donne de
bons résultats. Au delà, l’indépendance entre les tranches n’est plus respectée.
[
Sˆ xx ( f ) = TF Rˆ xx (k ) ]
A partir d'un échantillon (x0, x1, …..,xN-1) indépendants et identiquement distribuées, plusieurs estimateurs
sont alors possible:
N −1 N −1
1 1
Rˆ xx (k ) =
N
x ( n) x * ( n − k ) Rˆ xx (k ) =
N −k
x( n) x * ( n − k )
n=k
n =k
Le premier estimateur est biaisé car on ne tient pas compte du 0.5
unbiased
biased
0.35
0.2
0
-5 -4 -3 -2 -1 0 1 2 3 4 5
Frequence (Hz)
c) Dans cette méthode, on recherche un modèle AR, MA ou ARMA pour la séquence xk. Dans le cas d'un
modèle auto-régressifxk = a1xk-1+....+arxk-r+uk 2
1
Sˆ xx ( f ) = r
u0
1 + ak e − i 2πfk
k =1
Un signal aléatoire est amené à être transmis, analysé, transformé, etc. Conserve-t-il son caractère
aléatoire? sa stationnarité ? Que deviennent sa moyenne et son auto-corrélation statistiques lors d’un filtrage
linéaire. On examinera la transformation des caractéristiques du signal dans le domaine fréquentiel qui
67
Traitement Avancé du Signal Master RT M1 2022/2023
permettra d'aborder la notion de filtre formeur. Ces connaissances nous permettront d'envisager une
application directe qui est le filtrage optimal et adapté.
Soit un système linéaire et invariant dans le temps (SLIT) défini par sa x(t) h(t) y(t)
réponse impulsionnelle h(t) ou sa fonction de transfert H(f) :
Rappelons que la réponse du système linéaire et invariant à un signal quelconque déterministe est donnée
par :
y (t ) = x (t ) * h (t ) = x (τ ) h(t − τ ) dτ
Cette expression implique que si x(t) est un signal aléatoire, le signal en sortie y(t) est forcément un signal
aléatoire puisque la sortie est une somme pondérée de l'entrée. Il faudra donc caractériser y(t) de manière
statistique, de même que pour x(t), en étudiant les grandeurs statistiques déjà vues au chapitre précédent.
Supposons que x(t) est un processus aléatoire, la moyenne µy(t) du signal de sortie (aléatoire aussi) peut se
calculer par :
µ y ( t ) = E {y ( t )} = E {x (t ) * h (t )} = E { h (τ ) x (t − τ ) d τ }
= h (τ ) E {x ( t − τ )}d τ = h (τ ) µ x ( t − τ ) = h ( t ) *µ x ( t )
{ } { } {
Ry (t1 , t2 ) = E y(t1 ) y* (t2 ) = E x(t1 ) * h(t 1).x* (t2 ) * h* (t 2 ) = E h(τ1 ) x(t1 −τ1 )dτ1 h* (τ 2 ) x* (t2 −τ 2 )dτ 2 }
= h(τ )h (τ )E{x(t −τ ) x (t
1
*
2 1 1
*
2 }
−τ 2 ) dτ1dτ 2
-
µ y (t ) = h(τ ) µ x (t − τ )dτ = h(τ ) µ x dτ =µ x h(τ )dτ = µ x H (0) ( f = 0)
Ainsi, si le signal d’entrée est stationnaire au sens large (SSL) alors le signal de sortie sera aussi SSL
68
Traitement Avancé du Signal Master RT M1 2022/2023
Densité spectrale
En appliquant la transformée de Fourier aux deux membres de l’équation précédente, la densité spectrale de
puissance est obtenue :
Ainsi, la densité spectrale de puissance du signal de sortie est égale à la densité spectrale de puissance du
signal d’entrée multipliée par le carré du module de la réponse en fréquences du système (la phase de H(f)
n'intervient pas). C'est une propriété très importante à la base de nombreuses applications dont notamment
la notion de filtre formeur et le filtrage optimal.
Par ailleurs, si l’on souhaite calculer la corrélation du signal de sortie en éludant la lourdeur de calcul
inhérente au produit de convolution, on calculera la densité spectrale de puissance du signal de sortie et on
prendre la transformée de Fourier inverse.
( )
R y (τ ) = TF −1 H ( f ) S x ( f ) = H ( f ) S x ( f )e 2πjfτ df
2 2
{ }
La puissance moyenne du signal de sortie est alors obtenue par : E y (t ) 2 = R y (0) = H ( f ) S x ( f ) df
2
Soient y1(t) la sortie d'un système h1(t) dont l'entrée est un signal aléatoire x1(t) et y2(t) la sortie d'un système
h2(t) dont l'entrée est un signal aléatoire x2(t). La formule de l'interférence permet de relier l’intercorrélation
entre les sorties de deux filtres à celles des entrées de ces filtres:
S y1 y 2 ( f ) = H 1 ( f ) S x1 x 2 ( f ) H 2* ( f )
et R yx (τ ) = R x (τ ) * h (τ ) S yx ( f ) = H ( f ) S x ( f )
τ
−
Exemple 1 : Soit le processus stochastique x(t) SSL de corrélation statistique: R x (τ ) = σ e 2 θ
et le système
linéaire et invariant de réponse impulsionnelle h(t)=5e-2tU(t). On obtient alors:
5
- µ y (t ) = µ x H (0) = 0 × = 0×5/ 2 = 0
2 + 2πjf f =0
5
2
2 −τ 25. 2σ 2θ
S
- y ( f ) = H ( f ) S
2
( f ) = . TF σ e θ = .
x
2 + 2πjf 4 + (4π 2 f 2 ) 1 + 4π 2 f 2θ 2
Exemple 2 : Considérons par exemple le cas où le signal en entrée est un bruit blanc b(t). Sa DSP est donc
une fonction constante. Alors, en sortie, on aura un signal tel que :Sy(f) = constante|H(f)|2
69
Traitement Avancé du Signal Master RT M1 2022/2023
H( f ) = Sy( f ) /σb
2 2
On pourra donc décrire le signal aléatoire par les paramètres du filtre et la variance du bruit blanc. En
effet, si on fait passer le bruit blanc dans un filtre linéaire et stationnaire à paramètres ajustables et si on
obtient le signal désiré (à modéliser) à la sortie du filtre, alors on peut dire que toute l’information spectrale
est contenue dans le filtre représenté par ses coefficients (Voir exo 1 du TP n°3)
40
Filte Formeur
DSP du signal de sortie
30
20
Amplitude (dB)
10
-10
-20
-30
0 2000 4000 6000 8000 10000 12000
Frequence (Hz)
Remarque: Ainsi le bruit blanc joue pour l'aléatoire l'équivalent de la distribution de Dirac pour le
déterministe.
Rappelons que le bruit blanc n'a pas d'existence physique car il serait de puissance infinie. Une
approximation du bruit blanc à bande limitée appelé aussi bruit blanc coloré est défini par :
σ2 f <B
Sb ( f ) = b
0 ailleurs
La transmission d'un signal s'accompagne de distorsions (dues aux milieux de transmission) qu'il serait
souhaitable d'éliminer ou au moins d'atténuer avant tout traitement ultérieur. Si l'on connait le signal
(déterministe) d'origine, on parlera de détection. Dans le cas contraire ou si le signal est aléatoire on
emploiera le terme estimation. Il existe de nombreuses approches de détection, on peut, entre autres procéder
par filtrage qui permettra de rehausser le signal noyé dans le bruit.
70
Traitement Avancé du Signal Master RT M1 2022/2023
Considérons un signal déterministe x(t) supposé connu, dont on souhaite tester la présence possible
dans une observation s(t). Le bruit d'observation est supposé quant à lui stationnaire SSL de densité
spectrale Sb(f). On cherche un filtre H(f) qui maximise le SNR à un instant précis T0.
On suppose donc que le signal utile x(t) est noyé dans un bruit b(t) stationnaire SSL additif, d'où :
s(t) = x(t)+b(t)
On filtre le signal par un filtre linéaire dont la réponse impulsionnelle est h(t). A la sortie du filtre, on
obtient un signal :
Puis( x2 (T0 ))
A la sortie du filtre et à l’instant T0, le SNR s’écrit : SNR (T0 ) =
Puis(b2 (T0 ))
Le but étant de trouver h(t), on va alors exprimer le SNR en fonction de h(t) (ou H(f)), ainsi :
2 2
2πjfT0
Puis ( x2 (T0 )) = x2 (T0 ) = TF −1 ( X 2 ( f )) = X ( f ) H ( f )e
2
df
- Le signal b2(t) est issu du filtrage d'un signal aléatoire b(t) SSL donc il est aussi SSL, c'est ainsi que sa
puissance pourra se formuler comme suit:
{ }
Puis (b2 (T0 )) = E b2 (T0 ) 2 = Rb 2 (τ = 0) = TF −1 ( S b 2 ( f ))
τ =0
= S b ( f ) H ( f ) 2 df
Ainsi, on arrive à exprimer le SNR comme suit : 2
a( f )b ( f )df
*
2
2πjfT0
Puis ( x2 (T0 )) X ( f ) H ( f )e df SNR (T0 ) =
SNR(T0 ) = = a( f )a ( f ) df
*
S ( f ) H( f )
2
Puis (b2 (T0 )) b df
Avec a( f ) = Sb ( f ) H ( f ) et b ( f ) = X * ( f ) e −2 πjfT 0 / S b ( f )
Pour trouver H(f) qui maximise ce rapport, on fait appel à l’inégalité de cauchy-schwartz:
a( f )b ( f )df ≤ a( f )a * ( f ) df b( f )b ( f ) df
* *
a( f )b ( f )df
*
71
Traitement Avancé du Signal Master RT M1 2022/2023
avec égalité (soit le maximum) lorsque a(f)=kb(f) Cette dernière nous fournit l'expression suivante du
filtre optimal H(f) :
H ( f ) optimal = k . X * ( f )e −2πjfT 0 / S b ( f )
2
X( f )
L'expression du SNR devient : SNR(T0 )max = S ( f ) df
b
(Ce maximum est indépendant de T0)
Cas particulier : Par ailleurs, si le bruit b(n) est blanc, on parle de filtre adapté:
= k /σ x * (T 0 − t )
− 2 π jfT 0 2
= k /σ
2
H ( f ) adaptél b X *
( f )e donnant h ( t ) adaptél b
Ex
SNR(T0 ) max =
σ b2
x(t) h(t)=x(T0-t)
La réponse impulsionnelle du filtre représente le
signal utile x(t) renversé et décalé de [Link] filtrage adapté
revient à effectuer l'inter-corrélation entre l'observation
et le signal à détecter. T0
Cette réponse n’est pas causale ce qui ne permet pas de l'appliquer en temps réel. Cependant, ce filtre
peut-être appliqué à un signal sauvegardé.
Exemple 2 : En sonar ou en radar, on cherche à localiser une « cible ». Cette cible peut être par exemple un
bateau ou un avion. Pour cela, on procède de la façon suivante : on émet un signal x(t), qui parcourt la
distance d jusqu’à la cible, sur laquelle il sera réfléchi en direction d’un récepteur. Le récepteur reçoit alors le
signal y(t)bruité, atténué (de a) et retardé de TAR (voir TP n°3)
On utilisera alors comme filtre adapté le filtre adapté à x(t), soit h(t) = x*(T0−t)
La sortie de ce filtre sera l'intercorrélation des signaux y(t) et x(t). On obtient alors :
Sachant que l’autocorrélation est maximale en 0, l'intercorrélation sera maximale pour t = T0+TAR ce qui nous
permettra de déterminerTAR.
72
Traitement Avancé du Signal Master RT M1 2022/2023
1 1
0.5 0.5
0 0
-0.5 -0.5
-1 -1
0 50 100 150 200 250 300 0 50 100 150 200 250 300
Signal émis Filtre adapté
10 0.2
5 0.1
0 0
-5 -0.1
-10 -0.2
0 100 200 300 400 500 600 700 800 900 1000 0 200 400 600 800 1000 1200 1400
Signal reçu Signal reçu filtré
Remarques
- Le choix du signal 'x'est important : Il est souhaitable que le max de sa fonction d'autocorrélationRxx soit
facile à identifier (pic bien visible comme pour la fonction triangle)
-Dans le cas où le signal est sinusoïdal, le filtre adapté est un passe bande idéal centré sur la fréquence du
signal (détection synchrone).
On utilise ce signal pour déterminer la distance d'un objet. Le signal reçu par le recepteur aprèsréflexion sur
l'objet est :y(t)=αx(t-TD) + b(t) où α est une constante réelle positive etb(t) est un bruit blanc de densité
spectrale de puissance σ2. On veut maximiser le rapport signal sur bruit
1. Déterminer la réponse impulsionnelle du filtre h(t) telle que h(t)2dt=1 (prendre T0=T).
2. Donner le rapport signal sur bruit après filtrage.
3. Exprimer l'inter-corrélation temporelle Ryx(t-T0) en fonction de l'auto-corrélation de x(t)et de l'inter-
corrélation des signaux b(t) et x(t)
4. Comment déterminer le temps TD?
5. Pourquoi ne peut-on pas utiliser le filtrage optimal pour un signal déterministe dont on ignore
l'expression?
Il fait partie de la famille des estimateurs linéaires à variance minimale qui consistent à chercher une estimée
n
ŷ de y qui soit une fonction linéaire des observations, c'est-à-dire: yˆ = h x
i =1
i i = hT x
Dans le cas de l'estimation en moyenne quadratique, le filtre h doit alors être déterminé de telle sorte que la
variance de l'erreur d'estimation E{(y-ŷ)2} soit minimale :
{ } {(
E ( y − yˆ ) = E y − hT x
2
) (y − h x )}= E{y y }− E{y
T T T T
} { } {
h T x − E h xT y + E h x T hT x }
FGE, USTHB [ assiakourgli@[Link] / [Link]
73
Traitement Avancé du Signal Master RT M1 2022/2023
On commence par dériver par rapport à h puis on annule la dérivée, ce qui nous donne
Dans de nombreuses applications, les signaux temporels sont entachés d’un bruit que l'on souhaite
supprimer ou du moins réduire. Comme le signal utile aléatoire occupe les mêmes bandes de fréquences que
le signal parasite, on ne peut recourir au filtrage classique. Le filtre de Wiener apporte une solution à ce
problème lorsque le processus est stationnaire. Les filtres de Wiener sont dits optimum au sens du critère de
l'erreur quadratique moyenne entre leur sortie et une sortie désirée. .
On cherche donc un filtre qui minimiser a l'erreur. Il est pratique de chercher à minimiser e2 (n) car c'est une
fonction quadratique facilement dérivable. Par ailleurs, étant donné que les signaux sont aléatoires, la
fonction coût à minimiser est l'erreur quadratique moyenne (MSE) définie par : ξ (n) = E(e
2
(n))
Si on suppose que le filtre recherché H est un filtre RIF de longueur N (d'ordre N-1), on peut en calculer les
coefficients par résolution d'un système linéaire d'équations. h = [b0 b1 ... bN −1 ]
T
N −1
ˆ(n)peut alors s'écrire : yˆ (n) =
Le signal estimé y b x(n − i)
i =0
i
2
Rappelons que l’on cherche à minimiser ξ (n) = E{e (n)} = E y(n) − bi x(n − i )
N −1
2
i =0
Pour en obtenir le minimum, il suffit de chercher de dériver et d'annuler la fonction coût par rapport aux
variables bi de la réponse impulsionnelle du filtre. La dérivée de la fonction coût par rapport au jème point de
la réponse impulsionnelle est donnée par :
74
Traitement Avancé du Signal Master RT M1 2022/2023
∂ξ ∂ 2 ∂e(n) ∂
∂b j
= E { }
e (n) = E 2e(n) = E 2e(n) {− b j x(n − j)}
∂b j ∂b j ∂b j
∂ξ N −1
= − E{2e(n) x(n − j )} = − E 2 y (n) − bi x(n − i) x(n − j )
∂b j i =0
En faisant l'hypothèse que les signaux x(n) et y(n) sont stationnaires, et en annulant la dérivée, on trouve :
N −1
R yx ( j ) = bi Rxx ( j − i)
i =0
Ce qui pour les différentes valeurs de j, nous donne le système d’équations suivantes à résoudre :
Il faut noter que l'obtention des coefficients du filtre repose sur la connaissance de la fonction
d'autocorrélation du signal d'entrée et de l'intercorrélation entre les signaux d'entrée et de sortie désirée.
Exemple 1 : On suppose que l’observation x(n)=y(n)+bb(n)et que le bruit additif bb(n) est centré et non corrélé
au signal. Simplifions les équations de Wiener-Hopf en conséquences:
Rxx(k)=E{(y(n)+bb(n))(y(n-k)+bb(n-k))}= Ryy(k)+Rbb(k)
Le signal à estimer y(n)a pour la fonction d’autocorrélation Ry(k)=αk0<α<1. Il est décorrélé du bruit blanc
bb(n) de variance σb2. Cherchons h(n)tel que H(z) = b0 + b1z-1
1 + σ b2 α b0 1
. = H (z) =
1
[(1 + σ 2
− α 2 ) + ασ b2 z −1 ]
α 1 + σ b2 b1 α (1 + σ )
2 2
b −α 2 b
75
Traitement Avancé du Signal Master RT M1 2022/2023
Exemple 3: Comme illustration du filtrage de Wiener, on prend la mesure de l’activité cardiaque d’un fœtus
à l’aide d’un électrocardiogramme (ECG) pris au niveau de l’abdomen de la mère (le signal x(n)) et qui va
naturellement être perturbé par l’ECG de celle-ci auquel se rajoute le bruit thermique des électrodes et des
équipements électroniques. Pour retrouver l'ecg du foetus, on réalise une deuxième mesure fournissant
l’ECG de la mère (le signal bb). On peut employer alors le filtre de Wiener pour estimer le signal y(n)
représentant l'ECG du foetus.
10 500
ecg mère fft ecg mère
400
0
300
-10
0 500 1000 1500 2000 2500 3000 200
5
100
ecg foetus
0 0
0 500 1000 1500 2000 2500 3000
-5
0 500 1000 1500 2000 2500 3000
500
5
fft ecg foetus
bruit 400
0
300
-5
0 500 1000 1500 2000 2500 3000 200
10 100
f-ecg+m-ecg+bruit
0 0
0 500 1000 1500 2000 2500 3000
-10
0 500 1000 1500 2000 2500 3000
On peut observer que les TF des deux signaux ecg (mère et foetus) occupent la même plage de fréquences. A
partir des auto-corrélations de l'ECG mesuré sur l'abdomen x(n) et celui de la mère bb(n), on retrouve les
paramètres bi du filtre qu'on appliquera à x(n) pour obtenir une estimée de y(n) soit l'ECG foetal.
10
observation
Remarque
0
-5
0 500 1000 1500 2000 2500 3000
76
Traitement Avancé du Signal Master RT M1 2022/2023
1. Soient X(t) et Y(t) deux signaux aléatoires indépendants et stationnaires au sens large (SSL) tels que
−a τ
RX(τ) = 9 e et RY(τ) = b δ(τ) +c . On considère Z(t) = X(t) + Y(t),
- Quelles conditions doit-on imposer sur les constantes a, b et sur c (Justifier)
- Calculer les statistiques de premier ordre de Z(t)
- Calculer Rz(t,τ). Z(t) est-il SSL?
- Peut-on calculer sa DSP Sz(f)? si oui calculer la.
2. Soit le processus aléatoire Z(t) = [Link](2πf0t) - Ysin(2πf0t), où X, Y sont des variables aléatoires gaussiennes
indépendantes à moyennes nulles et de variances σx2 = σy2 = σ2 ; f0 est une constante.
a) Calculer la valeur moyenne E{Z(t)} et la variance σz2.
b) Calculer la fonction d’autocorrélation Rz(t1,t2) et examiner si le processus est SSL.
3. L’entrée x(t) du circuit est un bruit blanc de fonction d’autocorrélation Rx(τ) = σx2δ(τ)
- Déterminer la densité spectrale de la sortie y(t), notée Sy(f)
- Déterminer la fonction d’autocorrélation de la sortie y(t), notée Ry(τ) ainsi C
que sa puissance R
- Remplacer R par une self L et C par une résistance et reprendre les
questions
4. Les signaux x(n) et y(n) ont été obtenus en filtrant, au moyen d'un filtre à réponse impulsionnelle finie, un
bruit blanc b(n) gaussien centré de variance σ2.
A]L'équation de filtrage est la suivante : x(n) = 2.b(n) + 0.5.b(n -1) -0.2.b(n - 2) + 0.1.b(n - 3)
• Calculez, en fonction de σ2, les coefficients d'autocorrélation d'ordre 0,1,2,3 du signal x(n) .
• On notera Rxx (0), Rxx (1), Rxx (2), Rxx (3) ces coefficients.
• La répartition des niveaux d'amplitude du signal x(n) est elle gaussienne (sans justifier)?
B]y(n) = b(n) - b(n - 2)
• Donnez la TZ de la réponse impulsionnelle du filtre qui a permis d'obtenir y(n) à partir de b(n).
• Placez les zéros de ce filtre sur un cercle unité, quelles sont la ou les fréquence(s) coupées par ce filtre ?
• Tracez approximativement sa réponse en fréquence.
C]On considère maintenant le coefficient d'intercorrélation Rxy(k), entre les signaux x(n) et y(n) donné par
Rxy(k) = E[x(n) y(n - k)*].
Calculez Rxy (0), Rxy (1), Rxy (2), Rxy(-1), Rxy(-2)
5. Soit X(t) = A + b(t) un signal réel aléatoire, où A est une constante réelle et b(t) est un bruit blanc de densité
spectrale de puissance σb2, et soit un filtre moyenneur de réponse impulsionnelle:
1 Y(t)
h (t ) = Π (t − T / 2 ) X(t)
T T h(t)
- Exprimer l’autocorrélation statistique du signal d’entrée Rxx(t,τ). X(t) est-il SSL ?
- Déterminer la moyenne statistique de y(t).
77
Traitement Avancé du Signal Master RT M1 2022/2023
8. On considère un processus de la forme X(n) = θ + W(n) où W(n) est un processus gaussien de moyenne 0
et de variance 1 tel que W(n) et W(j) sont indépendants si n≠j. On suppose que θ suit une loi N (0,σ2)
indépendante de W(n) n∈Z.
- Caractériser le filtre de Wiener permettant d’estimer θ à partir de Xn, Xn−1, ..., Xn−N+1.
- Calculer les coefficients du filtre.
9. On désire dé-bruiter un signal de parole z(n) corrompu par du bruit additif b(n) indépendant du signal
sonore et ce, par filtrage de Wiener. On suppose connues quelques valeurs d’autocorrélation pour les deux
signaux telles que :
RZZ[0] = 1.5; RZZ[1] = 0.5; RZZ[2] = 0.25; RZZ[3] = 0.125 ; RZZ[4] = 0.0625
Rbb[0] = 1; Rbb[1] = 0.25; Rbb[2] = 0.0625; Rbb[3] = 0.015625
-Pourquoi ne peut-on pas employer un filtre adapté ou un filtre moyenneur ?
-Donner les équations de Wiener-Hopf permettant d’estimer z(n).
-Déterminer le filtre de Wiener d’ordre 2 permettant de retrouver le signal utile zˆ(n) puis exprimer H(z).
78
Traitement Avancé du Signal Master RT M1 2022/2023
Solutions
1. Voir interro 1 2016/2017
2. µ z (t) = 0 σ z2 (t) = σ 2 Rz (τ ) = σ 2 cos(2πf 0τ ) SSL
3. Rxx(t,τ)=A2+σb2δ(τ)=fct(τ) µx(t)=ASSL H(f)=sinc(fT)e-πjfT µy(t)=A filtre moyenneur
4. H(f)=1/(1+2πjfRC), Sy(f)= σx2/(1+4(πfRC)2 Ry(τ)=σx2/2RC e-|τ|/RC, P=Ry(0)=σx2/2RC, interro1 16/17
[Link](0)=(b02+b12+b22+b32).σ2 Rx(1)=(b0b1+b1b2+b2b3).σ2 Rx(2)=(b0 b2+ b1 b3).σ2 Rx(3)=(b0 b3).σ2Rx(k≥4)=0
H(z)=(z2-1)/z2 Rxy(0)=2.2 Rxy(1)=0.4 Rxy(2)=-0.2 Rxy(3)=0.1 Rxy(-1)=-0.5 Rxy(-2)=-2
6. x(t) Gaussien de moyenne s(t) et de variance 1. x(t) non stationnaire Sb’(f)=|H(f)| Ss’(f)=|H(f) |sinc2(f) h(t)=
2 2
7. Interro1 14/15 8. Rxx(0)=1+σ2 Rxx(k>0)= σ2 Rθx(k)=σ2 bi=σ2/(Nσ2+1) 9. µz=0, µb=0, signal (parole) aléatoire,
10. µs=0, µb=0 , Sb(f)=0.8. Quand signal utile et bruit occupent même plage de fréquences. b0=0.8/2.6 et b1=0.
Exercices supplémentaires
1. Soit le signal aléatoire défini par :x(t) =A1ej2πf1t + A2e j2πf2t où A1 et A2 sont deux v.a gaussiennes,
décorrelées, centrées et de variance σ2
- Donner la d.d.p conjointe fA1,A2(A1,A2)
- Déterminer la moyenne statistique du signal x(t) ainsi que sa fonction d’autocorrélation et la D.S.P
- Déterminer la d.d.p de x(t)
- Le processus x(t) est-il stationnaire au 2éme ordre ?
2. Le signal aléatoire x(t) réel SSL possède une fonction d’autocorrélation de la forme Rx(τ) = σx2e-β|τ|.
Un autre signal est lié à x(t) par l’équation déterministe suivante : y(t) = ax(t) + b, où a et b sont des
constantes données.
- Quelle est la fonction d’autocorrélation de y(t) ? En déduire µy et σy2
- Quelle est la fonction d’intercorrélation et de covariance de x(t) et y(t) ?
- Calculer le coefficient de corrélation. Ce résultat était-il prévisible, pourquoi ?
3. Parmi les fonctions d'auto-corrélation satitistiques suivantes lesquelles peuvent être celles d'un processus
aléatoire réel?
Rx1(τ)= Λ 2 (τ ) − 2 Rx2(τ)= − Λ 2 (τ ) + 2 −2τ
Rx4(τ)= e U(τ )
−2τ Rx5(τ)=δ(τ-1)+δ(τ+1)
Rx3(τ)= e
4. On considère le signal x(n) = g(n)cos(2πk0n/N+ϕ) défini pour n∈ [0,N-1] , et où g(n) est une fonction
aléatoire SSL, indépendante de ϕ v.a uniforme ente 0 et 2π .
- Calculez l’autocorrélation de x(n)
- Déduisez en sa densité spectrale de puissance, en fonction de la DSP de g, Sg(f).
5. Soit le signal aléatoire x(t) SSL dont la DSP est donnée par Sx(f)= σ 2 B Π ( f )
B
79
Traitement Avancé du Signal Master RT M1 2022/2023
6. On dispose d’un signal reçu qui est la version bruitée, retardée et atténuée d'un signal d’intérêt s(n). Le
bruit b(n) est supposé blanc gaussien de variance σ2. Le problème est de déterminer l’amplitude A et le retard
n0dans le signal reçu x(n) = A.s(n -n0) + b(n). Sachant que le rapport signal-à-bruit est maximum en sortie du
filtre adapté de réponse h(n) = s(-n).
- Vérifier que la sortie y(n) de ce filtre s’exprime comme la somme de deux fonctions de corrélation.
- Calculer la variance σb'2du bruit de sortie.
- On prend pour s(n) une impulsion rectangulaire de largeur L, tracer alors un exemple de signal reçu et
de sortie du filtre adapté.
h(t)
[Link] considère un processus de la forme X(n) = b(n) + αb(n-1) + W(n) où W(n) est un processus gaussien
de moyenne 0 et de variance σ2 tel que W(n) et W(j) sont indépendants si n≠j. On suppose que b(n) est une
variable aléatoire uniforme à valeurs dans {−1,1}, indépendante de(W(n) )n∈Z et que de même, b(n) et
b(j) sont indépendants si n ≠ j (P(b(n) = 1) = P(b(n) = −1) =1/2.
- Construire un filtre d'ordre 3 qui estime b(n) .
10. On considère un problème d'estimation d'un signal s(n) bruité et ayant subi un écho.
Le signal observé est x(n) = s(n) + 0.5s(n −1) + b(n) .
On suppose connue l'autocorrélation du signal utile et qu'elle a pour expression 0.5|k| et l'on suppose que le
signal utile est décorrélé du bruit dont l'autocorrélation est Rbb(k)=0.25|k|.
1. Déterminer les moyennes statistiques de x(n) et b(n)
2. Donner les équations de Wiener-Hopf permettant d’estimer s(n)
)
3. Déterminer le filtre de Wiener d’ordre 2 permettant de retrouver le signal utile s(n).
)
4. Exprimer s(n)et commenter
80
Traitement Avancé du Signal Master RT M1 2022/2023
Solutions
1. Voir interro 1 2014/2015 [Link](τ)=σx2e-β|τ| +b2 Sy(f)= 2 βσx2/(β2+4(πf)2)+ b2δ(f) Rxy(τ)=a Rx(τ)
3. Voir interro 1 2015/2016 [Link](m)=Rg(m)cos(2πk0m/N)/2 Sx(f)=[Sg(f-k0)+Sg(f+k0)]/4
5. Rx(τ)=σ2B2sinc(Bτ),µx=0 ,σx2=σ2B2, x(t) SSL+h(t) SLITy(t) SSL, Sy(f)= σ 2 B Π ( f ) , Rx(τ)=σ2ABsinc(Aτ)
A
1/ 2
6.y(n)=[Link](n-n0)+ Rbs(n) σb'2=Rb'(0) Sb’(f)=σb2|H(f)|2Rb'(0)=σb2 H () f
2
df
−1 / 2
81
Traitement Avancé du Signal Master RT M1 2022/2023
Buts : Manipuler des signaux aléatoires, les caractériser grâce à leurs moments d'ordre 1 (moyenne,
variance) et d'ordre 2 (autocorrélation, covariance). Aborder et acquérir les notions de stationnarité et
d'ergodicité. Estimation spectrale. Aborder la notion de Filtre Formeur.
1. Démo
# -*- coding: utf-8 -*-
import numpy as np; import [Link] as plt;
import [Link] as sp
N_R=500; N_va=1000; A=1; f0=.002; X=[];
t = [Link](0,N_va-1,N_va)
phi=[Link](0,2*[Link],N_R);
[Link](1);[Link](t,X[0,:]);[Link](t,X[1,:]);[Link](t,X[2,:]);[Link](t,X[3,:]);
"""Sationnarité d'ordre 1"""
moy=[Link](X,0); var=[Link](X,0);
[Link](2); [Link](211); [Link](moy); [Link](212); [Link](var);
"""Sationnarité d'ordre 2"""
Cx=[Link](X, rowvar=False);
[Link](3); [Link](Cx);
[Link](4); [Link](Cx[0,:]); [Link](Cx[20,:]);[Link](Cx[40,:]);
"""Ergodisme d'ordre 1"""
X1=X[1,:];
moy_t=[Link](X1); var_t=[Link](X1);
"""Ergodisme d'ordre 2"""
Cx_t=[Link](X1-moy_t,X1-moy_t,'full')/N_va
Rx_t=[Link](X1,X1,'full')/N_va
[Link](5)
[Link](Cx_t,label='Cov temp'); [Link](Rx_t,label='Cor temp');
[Link]()
# A = [Link](-1,1,N_R);
# for i in range (N_R):
# X=[Link](X,A[i]*[Link](2*[Link]*f0*t))
#X=[Link](N_R,N_va)
# mux=0; varx=1;
# X=mux+[Link](varx)*[Link](N_R,N_va)
82
Traitement Avancé du Signal Master RT M1 2022/2023
NF=2048;
TFb = [Link](bb,NF); TFb=[Link](TFb);
Sbf = [Link](TFb)**2/N; ff = [Link](-NF/2, NF/2-1, NF)*fe/NF;
s=[Link](Sbf)
[Link](313); [Link](ff,Sbf); [Link](True);[Link]("DSP du bruit blanc");
2. A faire
Télécharger le fichier ‘vous avez du courrier en [Link]’ et le placer dans le même répertoire que
votre programme commençant comme suit :
import numpy as np; import [Link] as plt;
from [Link] import wavfile as wf; import winsound ;
fname = '[Link]';
[Link](fname, winsound.SND_FILENAME)
fe, X = [Link](fname);
Te=1/fe; N=len(X); t = [Link](0, N-1, N)*Te;
[Link](1);[Link](211); [Link](t,X); [Link](True);
[Link]('temps'); [Link]('Amp'); [Link]('Phrase');
N1=21000;N2=21500; Y=X[N1:N2] ; ty=[Link](N1, N2-1, N2-N1)*Te;
[Link](212); [Link](ty,Y); [Link](True);
[Link]('temps'); [Link]('Amp'); [Link]('morceau de la phrase');
#zz=np.int8(Y) ; [Link]("[Link]", fe, zz);
83
Traitement Avancé du Signal Master RT M1 2022/2023
#[Link]("[Link]",winsound.SND_FILENAME) ;
1. Le signal étudié est un processus aléatoire ou une réalisation d’un signal aléatoire?
2. Calculer la DSP. Est-elle d’origine statistique ou temporelle?
3. Calculer et visualiser la moyenne et la variance temporelles pour différentes parties du signal qui se
chevauchent pour des durées de
N=1000,500, 250, 125. Pour quelle valeur de N peut-on considérer le processus stationnaire?
1. Comparer les 4 techniques en faisant varier A1, A2 et A3 et commenter.
2. Rapprocher les fréquences et observer.
3. Commenter les quatre techniques d’estimation de la DSP en faisant varier N.
4. Quelle est l'avantage et l'inconvénient du moyennage?
84
Traitement Avancé du Signal Master RT M1 2022/2023
1. Démo
# -*- coding: utf-8 -*-
import numpy as np; import [Link] as plt
" Filtrage adapté "
N = 500; var = 0.16 ; mu = 0;
x=[Link](10); Retard = 100
y = [Link](var) * [Link](N) + mu
y[Retard:Retard+len(x)]= y[Retard:Retard+len(x)] + x
#Retard=300; y[Retard:Retard+len(x)]= y[Retard:Retard+len(x)] + x
Ryx = [Link](y,x,"full"); Ryx=Ryx[len(x)-1:]
[Link](1)
[Link](311); [Link](x); [Link]('signal émis')
[Link](312); [Link](y); [Link]('signal reçu')
[Link](313); [Link](Ryx); [Link]('Intercorrélation entre signal émis et signal reçu')
import [Link] as la; import [Link] as sp;
fecg=[]; mecg=[];
2. A faire
FGE, USTHB [ assiakourgli@[Link] / [Link]
85
Traitement Avancé du Signal Master RT M1 2022/2023
86
Traitement Avancé du Signal Master RT M1 2022/2023
Introduction
Nous avons vu précédemment qu'un signal aléatoire peut être modélisé (synthétisé) comme la réponse
d'un filtre linéaire à une excitation sous forme de bruit blanc tel que H ( f ) = S x ( f ) / σ b . Ce filtre formeur
2 2
H(f) est aussi dit processus générateur. Ces paramètres associés à la variance du bruit σb2 constituent le
modèle mathématique correspondant au signal aléatoire. Le concept de processus générateur de signal a été
particulièrement développé et appliqué avec des filtres numériques. Ce sont les modèles de signaux les plus
utilisés en traitement statistique du signal (estimation, prédiction...). Selon la nature du filtre, on peut obtenir
différents modèles de signaux (AR,MA, ARMA, etc.)
Exemple : Lors d’un appel par GSM (téléphone portable), le portable qui fait office d'un micro- ordinateur
regroupant différentes fonctionnalités dont l'analyse, la synthèse, le codage, etc. va nous permettre de
modéliser la parole (aléatoire) en opérant un codage LPC par tranches de dizaines de ms. Elle consiste à
retrouver les paramètres du filtre formeur h(n) pour chaque tranche y(n) enregistrée et analysée. Ce sont ces
paramètres (ai et bj) qui seront transmis pour produire un signal de synthèse approchant le signal original
y(n).
x(n) Bruit Blanc y(n) signal connu
h(n) à modéliser
à déterminer
b j z−j
Sy( f ) = σb H( f )
j=0 2 2
H ( z) = M
a z
i=0
i
−i
Ci-dessous la phrase "vous avez du courrier en attente" échantillonnée à une fréquence fe =22050 (45531
échantillons) et un zoom sur une voyelle de durée 25 ms (500 échantillons).
87
Traitement Avancé du Signal Master RT M1 2022/2023
40
phrase: Vous avez du courrier en attente Filter formeur AR d ordre 20
0.5
Amplitude (dB)
20 Spectre de la voyelle
Amplitude
0 0
-20
-0.5
-40
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 2000 4000 6000 8000 10000 12000
Temps (s) Frequence (Hz)
40
Voyelle sur 25 ms Spectre du son synthétisé
0.5
Amplitude (dB)
20
Amplitude
0 0
-20
-0.5
-40
0.954 0.956 0.958 0.96 0.962 0.964 0.966 0.968 0.97 0.972 0.974 0 2000 4000 6000 8000 10000 12000
Temps (s) Frequence (Hz)
En comparant l'allure du filtre formeur H(f) à celle de Sy(f), on note que l'on retrouve l'allure générale du
spectre de la voyelle, notamment les fréquences dont la puissance est maximale. Sachant que le spectre de la
voyelle comporte 500 valeurs et que le filtre H(f) équivalent est obtenu à partir de 20 coefficients, il vaut
mieux transmettre 21 coefficients (20 ai + variance du bruit) que 500 valeurs.
C'est ainsi que les modèles autorégressifs sont d'un emploi de plus en plus répandu en traitement du signal
: codage et transmission par prédiction linéaire, synthèse de parole, reconnaissance, etc.
Remarque : Le signal de parole est un processus aléatoire non-stationnaire à long terme, mais il est considéré
comme stationnaire dans des fenêtres temporelles d’analyse de l’ordre de 20 à 30ms. Cette propriété de
stationnarité à court terme permet donc une analyse et modélisation progressive du signal de parole Pour
éviter toutes pertes d'information, on veillera à prendre des fenêtres chevauchantes.
Les signaux autorégressifs sont obtenus par passage d'un bruit blanc dans un filtre purement récursif.
Ce filtre est donc de réponse impulsionnelle infinie.
N
H ( z ) = 1 / 1 + a i z −i
i =1
N
A partir de H(z), on peut déterminer l'équation aux différences : y (n ) = x ( n) − a y (n − i)
i =1
i
Cela signifie que le signal y(n) est supposé être prédictible en fonction d'un certain nombre de ses
valeurs antérieures.
- µx = E{x(n)} = 0 - {
R xx ( 0 ) = E x ( n ) 2 = σ } 2
- Rxx(k) = E{x(n)x(n − k)} = 0 pour k ≠ 0
88
Traitement Avancé du Signal Master RT M1 2022/2023
Calculons alors R yy (k )
N
R yy (k ) = E{y(n) y(n − k )} = E{x(n) y(n − k )} − ai E{y(n − i) y(n − k )}
i =1
N N
= E x(n) x(n − k − k ' ) h(k ' ) − ai R yy (k − i) = E{x(n) x(n − k − k ' )}h(k ' ) − ai R yy (k − i)
k' i =1 k' i =1
N N
= Rxx (k + k ' ) h(k ' ) − ai R yy (k − i) = Rxx (k )h(0) + Rxx (k + 1)h(1) + ....... − ai R yy (k − i )
k' i =1 i =1
N N
-Si k=0, Ryy (0) = Rxx (0).h(0) − ai Ryy (−i) = σ 2 .1 − ai Ryy (i)
i =1 i =1
N N
- Pour k=1 à N, Ryy (k ) = 0 − ai Ryy (k − i) = −ai Ryy (k − i)
i =1 i=1
On peut utiliser une forme matricielle : R. a = S (Rappelons que pour un signal réel Ryy(k)=Ryy(-k)
R est la matrice d’autocorrélation dont le terme général rij ne dépend que de la différence i-j (Matrice
de Toeplitz). La résolution de ces équations dites de Yule-Walker permet de connaître les paramètres du filtre
et la variance du bruit blanc.
On peut reformuler cette matrice sous la forme suivante (On l'emploiera en TP) :
C B 4
6444444447 44444444 8 647 8
R yy (0) R yy (1) ...... R yy ( P − 1) a1 − R yy (1)
R (1) R yy (0) ...... R yy ( P − 2) a 2 − R yy (2)
yy
. =
...... . .
R yy ( P − 1) R yy ( P − 2) ...... R yy (0) a P − R yy ( P )
P
σ 2 = Ryy (0) + ai Ryy (i)
i =1
Remarques : Nous ne disposons pas d'un processus aléatoire mais d'une seule réalisation soit y(n), il n'est pas
donc pas possible de calculer l'auto-corrélation statistique Ryy(k). Cette dernière sera remplacée par
l'autocorrélation temporelle en faisant l'hypothèse que le processus est ergodique (voir exo 1 du TP n°4).
Il existe divers algorithmes (Burg, Levinson) qui permettent d'estimer assez rapidement les ai et σ sans
passer par l'inversion matricielle. Tout comme il est possible de déterminer l'ordre N adéquat (critère AIC).
89
Traitement Avancé du Signal Master RT M1 2022/2023
Exemple d'application
• (a1=-Rx(1)/Rx(0) σ2=Rx(0)(1-a12) )
• (Rx(0)= σ2/(1-a12) Rx(1)=-a1σ2/(1-a12)Rx(k)=(-a1)kσ2/(1-a12) )
2. Refaire le même travail pour un modèled’ordre 2 tel que : x(n)=-a1.x(n-1)-a2.x(n-2)+b(n)
Réponses :
• a1=Rx(1)[Rx(2)-Rx(0)]/[Rx(0)2-Rx(1)2] a2=[Rx(1)2-Rx(0)Rx(2)]/[Rx(0)2-Rx(1)2]
• σ2=Rx(0)+Rx(1)2[Rx(2)-Rx(0)]/[Rx(0)2-Rx(1)2]+Rx(2)[Rx(1)2-Rx(0)Rx(2)]/[Rx(0)2-Rx(1)2]
• Rx(1)=-a1Rx(0)/(a2+1) Rx(2)=(-a2+a12/(1+a2))Rx(0)Rx(0)=(1+a2)σ2/(1+a2- a12–a12-a23+a2a12)
-Algorithme de Levinson
C’est un algorithme qui permet de résoudre tout système du type Ax = b avec A Toeplitz, donc en particulier
les équations normales de Yule-Walker :
Le principe d l'algorithme de Levinson est de trouver de façon récursive les paramètres d’ordre k en fonction
des paramètres d’ordre k − 1.
- Initialisation : i 1 [ 1
jkk
jkk %
l 1 |i | mnn 0
∑2p4
i. [ [ [)
jkk . rs4 o2p4 q jkk . q
t2p4
/
i. u i. u C i. [ i. [ u vwxy u 1, … … … . , [ 1
l. 1 |i. [ | l.
90
Traitement Avancé du Signal Master RT M1 2022/2023
Les signaux à moyenne mobile sont obtenus par passage d'un bruit blanc dans un filtre purement
transverse.
M
Ce filtre est aussi appelé filtre à réponse impulsionnelle finie : H ( z) = b z
i =0
i
−i
Le signal y(n) est supposé pouvoir s'écrire comme une combinaison linéaire d'échantillons décorrélés
entre eux, ce qui peut se formaliser comme une combinaison linéaire d'échantillons d'un bruit blanc x(n).
M
On a donc : y(n) = b x(n − i)
i =0
i
M M
et µ y = E{y(n)} = bi µ x = µ x bi = µ x .H f (0)
i =0 i =0
On cherche les paramètres du filtre qui génèrent y(t) à partir de x(t), bruit blanc centré :
M M
R yy ( k ) = E {y ( n ) y ( n − k )} = E bi x ( n − i ). b j x ( n − j − k ) .
i =0 j =0
M M M M
R yy (k ) = bi . b j E{x(n − i).x(n − j − k )} = bi . b j Rxx ( j + k − i)
i =0 j =0 i =0 j =0
• Si j+k≠i R yy ( k ) = 0
M −k
• Sinon R yy ( k ) = σ b j + k .b j
2
j =0
Le problème est non linéaire en fonction des coefficients, il faut un algorithme de programmation non
linéaire pour obtenir bi à partir des Ryy (k). Cependant, l'algorithme de Durbin permet d'approcher la solution
optimale avec de bons résultats. Le principe de cet algorithme consiste à identifier le modèle MA d'ordre M
avec un modèle AR d'ordre N>>M (TP n°4). En effet, tout modèle MA peut être identifié à un modèle AR
M ∞
d'ordre infini: bi z −i = 1/ ai z −i
i =0 i =0
Exemple
• Calculer µx.
• En supposant que x(n) est connu, déterminer les paramètres du modèle.
• Connaissant les paramètres du modèle, déterminer les Rx(k)
FGE, USTHB [ assiakourgli@[Link] / [Link]
91
Traitement Avancé du Signal Master RT M1 2022/2023
• Calculer µx.
• Connaissant les paramètres du modèle, déterminer les Rx(k) Réponses
µx =0 Rx(0)=(1+ b12).σ2 Rx(1)=b1.σ2 Rx(k)=0 pour k≥2
b1=(Rx(0)±sqrt(Rx(0)2-4Rx(1)2)/2Rx(1) σ2=2/(Rx(0)±sqrt(Rx(0)2-4Rx(1)2)
Remarque : Il est très important de remarquer que nous ne disposons que d'une seule réalisation du signal
aléatoire à modéliser y(n), de ce fait l'auto-corrélation statistique Ryy(k) est obtenu de l'auto-corrélation
temporelle en considérant le processus ergodique.
4. Modèle ARMA
Les signaux ARMA sont obtenus par passage d'un bruit blanc dans un filtre récursif appelé aussi filtre
à réponse impulsionnelle infinie (R.I.I). Ces signaux sont une combinaison des signaux AR etMA. La fonction
de transfert du filtre présente un numérateur et un dénominateur:
M N M
N
Soit y(n) = b x(n − i) − a y(n − i)
i i
H ( z) = b z
i =0
i
−i
1 + a i z −i
i =1
i =0 i =1
N M −k
La corrélation statistique de y(n) s’écrit alors : R yy ( k ) = − ai R yy ( k − i ) + σ 2 b j + k .b j
i =1 j =0
La modélisation ARMA peut se décomposer en une modélisation AR suivie d'une modélisation MA.
Le modèle AR présente une simplicité de calcul par rapport aux modèles MA et ARMA du fait où les
coefficients AR sont solutions d’un système linéaire d’équations. Alors que la détermination des coefficients
MA et ARMA requiert la résolution d’équations non linéaires. Cependant, le modèle ARMA permet de
modéliser aussi bien les minima que les maxima de la DSP et est donc moins restrictif que le modèle AR.
Les applications des modèles AR, MA, ARMA sont nombreuses, entre autres :
- la modélisation et la prédiction de série temporelle dite séries chronologiques Une série chronologique
est une suite formée d’observations au cours du temps que l'on cherche à modéliser pour la prédiction de
données futures. Ainsi, en finance, cela permet de modéliser le cours des devises ou du pétrole. Alors qu'en
météorologie, les cela permet de faire des prévisions sur la température ou les précipitations. Dans chacun
des cas, on essaiera à partir d'un échantillon de données de construire le meilleur modèle qui s'ajuste ces
données.
- l'estimation du spectre d’un signal aléatoire, etc. Cette dernière application est basée sur
l'identification des paramètres du modèle considéré : Le modèle AR est bien adapté aux signaux composés
92
Traitement Avancé du Signal Master RT M1 2022/2023
de raies pures dans du bruit blanc. Alors que le modèle MA est bien adapté aux signaux dont la puissance
est nulle dans certaines bandes de fréquences.
Exemple : Reprenons à nouveau l'exemple de parole, pour lequel, nous avons testé les trois variantes
40 35
Spectre de la voyelle Filter formeur AR d ordre 20
30
30
25
20
20
Amplitude (dB)
Amplitude (dB)
10 15
0 10
5
-10
0
-20
-5
-30 -10
0 2000 4000 6000 8000 10000 12000 0 2000 4000 6000 8000 10000 12000
Frequence (Hz) Frequence (Hz)
30 25
Filter formeur MA d ordre 20 Modèle ARMA (6,8)
20
20
15
10
Amplitude (dB)
Amplitude (dB)
10
0
5
-10
0
-20
-5
-30 -10
0 2000 4000 6000 8000 10000 12000 0 2000 4000 6000 8000 10000 12000
Frequence (Hz) Frequence (Hz)
Rappelons que pour un ordre du filtre faible, le modèle ne sera pas en mesure de capturer toutes les
informations pertinentes du signal, à contrario s’il est trop élevé, il modélisera le bruit. Rappelons qu'avec
le périodogramme (Algorithmes rapides : FFT) et ses variantes :
Nous n'irons pas plus loin, dans le cadre de ce cours, le but n'étant que d'introduire les notions de
modélisation et de prédiction linéaire.
Dans le filtrage de Wiener, le critère d’optimalité est stochastique : on désire minimiser la moyenne de
l’erreur au carré. Cela requiert des statistiques de second ordre des processus (moyennes et corrélations). En
93
Traitement Avancé du Signal Master RT M1 2022/2023
pratique, la matrice de corrélation m{{ [ et l'inter-corrélation mn{ [ ne sont pas connues et doivent donc
être estimés à partir des données observées (temporelles) en supposant les processus ergodiques.
Le second problème est lié à l'hypothèse même de stationnarité qui n'est pas vérifiée à long terme. Le filtrage
de Wiener devient inadéquat pour les situations dans lesquelles le signal ou le bruit sont non stationnaires.
Dans de telles situations le filtre optimal doit être variable dans le temps. La solution à ce problème est fournie
par le filtrage adaptatif.
y(n) -
Le filtrage adaptatif comporte une mise à jour récursive des paramètres (coefficients) du filtre. L’algorithme
part de conditions initiales prédéterminées et modifie de façon récursive les coefficients du filtre pour
s’adapter au processus. Pour cela, les coefficients de la réponse impulsionnelle du filtre sont adaptés en
fonction de l'erreur par une boucle de retour.
- L’algorithme du gradient stochastique (LMS) est une version déterministe dans laquelle les grandeurs
statistiques impliquées sont remplacées par des valeurs instantanées.
- L'algorithme des moindres carrés (RLS) qui est un algorithme récursif du LMS basé sur la
minimisation de la somme des erreurs quadratiques sur un laps de temps donné.
- [
A l’instant n, ces paramètres sont notés h(n) = b0 (n) b1 (n) ... bN −1 (n) et]T
- les entrées du filtre forment le vecteur X | , 1 ,…… }C1 ~ .
- A chaque instant n, la sortie du filtre est un signal (• que l’on souhaite aussi proche que possible
du signal original y(n).
On rappelle que la solution optimale donnée par la solution des équations de Wiener-Hopf est : m{{ ℎ mn{
2
i =0
L’algorithme procède en trois étapes, une étape d’initialisation et deux étapes qui sont itérées :
94
Traitement Avancé du Signal Master RT M1 2022/2023
L’idée de cet algorithme est que les ajustements en sens opposé au gradient (descente du gradient) vont
conduire h(n) vers la valeur hopt qui est associée à l’erreur minimale ξmin. Si la dérivée est positive alors on
va diminuer ℎ de manière à se déplacer vers le minimum de la fonction de coût. Le terme correctif sera
donc choisi négatif et vice-versa
Cet algorithme est aussi appelé algorithme du gradient déterministe puisque les quantités mn{ et m{{ sont
parfaitement connues et déterminées. On parle aussi de filtrage de Wiener adaptatif.
tend vers la solution optimale de l’équation de Wiener-Hopf, sans inverser la matrice m{{ . Néanmoins, à
La méthode de descente (gradient) est un algorithme qui estime, par itérations successives, une solution qui
chaque itération, il faut calculer les corrélations statistiques m{{ et mn{ que nous avons remplacé par des
corrélations temporelles sous hypothèse d'ergodicité.
gradient de ξ (n) soit ∇‚ ƒ 2mn{ C 2m{{ ℎ en remplaçant les moyennes statistiques (ergodiques :
L'algorithme adaptatif connu sous le nom de Least-Mean-Square (LMS) repose sur la simplification du
mŠ{{ X X
mŠn{ ( X
A noter que minimiser l’erreur quadratique instantanée est moins efficace et moins précis que corriger
l’erreur quadratique moyenne.
Soit ℎŠ C1 ℎŠ CcX
95
Traitement Avancé du Signal Master RT M1 2022/2023
L'algorithme LMS converge en moyenne, c’est-à-dire que, pour n∞ , ℎŠ hopt des équations de Wiener-
Hopf. A noter que l’algorithme de descente du gradient tend vers l’erreur minimale ξmin lorsque n tend vers
l'infini et ceci de façon exponentielle régulière, en raison de l’estimation exacte des m{{ et mn{ . Au contraire,
Si 1 , ξ est somme des erreurs quadratiques; Si d … 1, les erreurs passées sont pondérées avec un facteur
(d'oubli) qui décroît exponentiellement.
- ξ ∑. d .
où e k y k (• k y k ℎ X [
m{{† ℎ mn{†
m{{† ∑[ 1d X [ X [ d ∑[ 1d X [ X [ dm{{† 1 CX [ X [
[ 1 [ 1
Avec
mn{† ∑[ 1d ( [ X [ d ∑[ 1d ( [ X [ C( [ X [ dmn{† C( [ X [
[ 1 [ 1
et
- ℎŠ m{{† mn{†
Remarque : Le filtrage des moindres carrés se place en déterministe. Alors que le filtre de Wiener utilise la
méthode des moindres carrés en stochastique.
96
Traitement Avancé du Signal Master RT M1 2022/2023
Exemples d'application
Le filtrage adaptatif (au sens de Wiener, au sens des moindres carrés, le filtrage de Kalman, le filtrage
particulaire) est utilisé pour divers objectifs entre autres : la modélisation, la prédiction, l'identification des
systèmes, l'annulation d'interférence (écho, bruit), etc. permettent la reconstitution d’un signal en milieu bruit
- Annulation d’un brouilleur tel un signal perturbé par un brouilleur sinusoïdal de fréquence connue
mais dont on ignore l’amplitude A et la phase φ.
- Kalman : Le filtre de Kalman est employé pour estimer des paramètres d'un système évoluant dans le temps
à partir de mesures bruités. Il est particulièrement adapté pour le tracking du fait qu'il permet non seulement
de prédire les paramètres mais aussi de rectification les erreurs de mesure et du modèle. Il opère de façon
similaire aux filtres adaptatifs. La prédiction se fait en employant l'estimation précédente. Le calcul d'erreur
permet, ensuite, de faire la mise à jour de cette prédiction grâce aux nouvelles mesures.
A noter que la version étendue (Filtre de Kalman étendu) permet de prendre en compte une modélisation
non linéaire.
- Particulaire (méthode de Monté-Carlo séquentielle) : Une alternative aux filtres de Kalman étendu. E filtre
a pour but d'estimer la densité postérieure des variables qu'on cherche à estimer (dites variables d'état) en
fonction des observations (variables d'observation) en supposant que les variables d'état constituent une
chaine de Markov de premier ordre.
Il est utilisé particulièrement pour des applications de navigation : suivi de cible, guidage de missile,
robotique, etc.
97
Traitement Avancé du Signal Master RT M1 2022/2023
3. Soit un processus AR défini par : y(n)=-a1y(n-1) – a2 y(n-2) + x(n) où x(n) bb décorrélé de variance 1
• Calculer µy(n) puis sans calcul, expliquer pourquoi y(n) est SSL.
• Montrer que pour k>0, Ry(k)=-a1Ry(k-1)-a2 Ry(k-2)
• Déterminer a1et a2
4. On veut modéliser le signal y(n) dont l'autocorrélation statistique est donnée à la figure ci-contre.
Ry(k). On suppose que l'autocorrélation statistique du signal d'entrée est 4δ(k)
o Déterminer les moments statistiques d'ordre 1 du signal à modéliser.
2
o S'agit-il d'un modèle AR ou MA? Justifier
o Déterminer l'ordre du modèle 1 k
• En supposant que 2 coefficients sont égaux, à partir
des valeurs de Ry(k), déduire que l'un des coefficients est nul -5 -4 -3 -2 -1 1 2 3 4 5
puis déterminer les coefficients de ce modèle puis donner l'équation aux différences liant y(n) et x(n)
Solutions
1. Modèle AR (Ry(k)≠0 ) x(n) entrée bb, y(n) signal aléatoire à modéliser h(n) filtre formeur (modèle
math)entrée bb + ergodisme a1= - α σx2=1-α2
[Link](k)=σx2δ(k) Sx(f)=σx2 MA d'ordre 2 Ry(0)=1+b12+b22Ry(1)=b1(1+b2) Ry(2)= b2 Ry(k)= 0 pour k≥0
3. µy(n)=0 entrée bb SSL sortie SSL Interro1 14/15 4. Voir examen 17/18 5. Examen 16/17
98
Traitement Avancé du Signal Master RT M1 2022/2023
Exercices supplémentaires
[Link] considère un signal aléatoire stationnaire x(n) et l'on suppose connu ses coefficients d'autocorrélation :
R(0) = 3σ2 , R(1) = 2σ2 , R(2) = σ2 , R(3) = 0
- On cherche le filtre MA d'ordre 3 de ce signal. Identifiez les paramètres du filtre.
- Si le signal x(n) a été obtenu par filtrage d'un bruit blanc gaussien de variance σ2 par un filtre à réponse
impulsionnelle finie de fonction de transfert H(z) = 1+ a.z-1 + b.z-2, en déduire les valeurs de a et b.
2. Soit un filtre formeur dont l'équation aux différences est y(n)=0.5 (x(n)+x(n-1)+x(n-2)+x(n-3))
• Expliquer la notion de filtre formeur puis identifier ce modèle linéaire AR ou MA
• Donner la moyenne, l'autocorrélation et la DSP de son entrée x(n).
• Calculer et tracer Ryy(k) puis déterminer Syy(f)
3. Les signaux x(n) et y(n) ont été obtenus en filtrant, au moyen d'un filtre à réponse impulsionnelle finie, un
bruit blanc b(n) gaussien centré de variance σ2.
x(n) = 2.b(n) + 0.5.b(n -1) -0.2.b(n - 2) + 0.1.b(n - 3)y(n) = b(n) - b(n - 2)
• Identifier les 2 modèles puis pour chacun, calculer et tracer les coefficients d'autocorrélation
• x(n) et y(n) sont-ils Gaussiens (justifier)
• Calculer les intercorrélations Rxy(k), et Ryx(k) et commenter
• Donner quelques applications des modèles AR,MA, ARMA
4. Soit un filtre formeur dont l'équation aux différences est y(n)= - αy(n-1) - β y(n-2) + x(n)
• Identifier l'ordre du modèle linéaire AR.
• Déterminer les paramètres du modèles, on suppose que Ryy(k)=2*0.5|k|
5. Soit un filtre formeur dont l'équation aux différences est y(n)= αy(n-1) + x(n)
• Identifier l'ordre du modèle linéaire AR.
• Déterminer la moyenne de y(n) et montrer que Ryy(k)=αk Ryy(0), déduire une condition sur α.
• Montrer que Ryy(0)=Rxx(0)/(1-α2)
• Tracer Ryy(k) (Prendre Rxx(0)=σ2=1).
• Citer 2 applications concrètes des modèles AR.
6. Soit un filtre formeur dont l'équation aux différences est y(n)= 0.25 y(n-1) - 0.25 y(n-2) + x(n)
• Identifier ce modèle linéaire AR ou MA (Justifier)
• Déterminer la moyenne de y(n) et donner l'expression de Ryy(k).
• Calculer et tracer Ryy(k) (Prendre Rxx(0)=σ2=1).
Solutions
1. b0=1 b1=1 et b2=1 etpar identification, on trouve a=1 et b=1
2. MA, µx=0 et Sx(f)=σ2Ry(0)=σ2, Ry(1)=Ry(-1)=0.75 σ2, Ry(-2)=Ry(2)=0.5 σ2, Ry(-3)=Ry(3)=0.25 σ2,
Ry(k>3)=[Link](f)=σ2 (1+1.5 cos(4πf)+ cos(6πf) + 0.5 cos(8πf))
3. Interro2 15/16 4. Ratt 15/16 5. Examen 15/16 6. Ratt 14/15
99
Traitement Avancé du Signal Master RT M1 2022/2023
Ce TP a pour objectif :
- Modéliser un signal aléatoire par les modèles AR et MA (par approximation AR) et d'en extraire les
informations utiles.
- Tester une méthode de filtrage adaptatif (LMS)
Exercice 1: Télécharger le fichier ‘vous avez du courrier en [Link]’ et le placer dans le même répertoire
que votre programme puis sélectionner une partie entre 21000 et 21500
I. Démo
# -*- coding: utf-8 -*-
# fname = '[Link]';
# #[Link](fname, winsound.SND_FILENAME)
# fe, S = [Link](fname);
# Te=1/fe; N=len(S); t = [Link](0, N-1, N)*Te;
# [Link](1);[Link](211); [Link](t,S); [Link](True);
# [Link]('temps'); [Link]('Amp'); [Link]('Phrase');
# L=500;N1=21000;N2=N1+L; y=S[N1:N2] ; ty=[Link](N1, N2-1, N2-N1)*Te;
# [Link](212); [Link](ty,y); [Link](True);
# [Link]('temps'); [Link]('Amp'); [Link]('morceau de la phrase');
100
Traitement Avancé du Signal Master RT M1 2022/2023
fecg=[]; mecg=[];
""" Lecture des fihiers """
with open('[Link]') as f:
for i in f: [Link](float(i))
fecg = [Link](fecg)
with open('[Link]') as f:
for i in f : [Link](float(i))
mecg = [Link](mecg)
""" Rajout du bruit """
N=len(mecg); var = 0.01;
bb=(var**0.5)*[Link](N);
x=fecg+bb+mecg; #x=[Link](x)
y=fecg; #y=[Link](y)
# fe=10000 ;Te=1/fe;
# f0= 100;
101
Traitement Avancé du Signal Master RT M1 2022/2023
# t = [Link](0,N)*Te
# y = 0.5*[Link](2*[Link]*f0*t)
# x = y + bb
# Ncoef=5; mu=0.2;
# hp = [Link](Ncoef) # Initialisation des coeffcients du filtre
# L = len(x) # Longueur du signal
# yest = [Link](N) # y(n) estimé
# En = [Link](N) # Erreur
# for i in range(Ncoef,L-1):
# X= [Link](x[i-Ncoef+1:i+1]); #Prendre X=[x(n),x(n-1),...,x(n-Ncoef+1)))]
# yest[i] = [Link](hp, X)
# En[i] = x[i] - yest[i]
# hp = hp + mu * En[i] * X
# Ncoef=3; mu=0.6;
# hs = [Link](Ncoef) # Initialisation des coeffcients du filtre
# yest = [Link](N) # y(n) estimé
# En = [Link](N) # Erreur
# for i in range(Ncoef,N-1):
# X= [Link](x[i-Ncoef+1:i+1]); #Prendre X=[x(n),x(n-1),...,x(n-Ncoef+1)))]
# yest[i] = [Link](hs, X)
# En[i] = y[i] - yest[i]
# hs = hs + mu * En[i] * X
# print(hs)
# [Link](3); [Link](y,'b', label= 'Sortie'); [Link](yest,'r', label= 'Sortie estimée')
# [Link](En,'g', label= 'Erreur instantanée'); [Link](); [Link]()
# def ModeleAR(y,P):
102
Traitement Avancé du Signal Master RT M1 2022/2023
# y=[Link](y)
# R = [Link](y,y,mode='full')[len(y)-1:]
# C = [Link](R[0:P])
# B = -R[1:P+1]
# a = [Link](C).dot(B)
# a = [Link](1,a)
# sigma_carre = sum(a*R[0:P+1]);
# return a,sigma_carre
# fe,x = [Link]('[Link]')
# [Link](x, fe)
103
Traitement Avancé du Signal Master RT M1 2022/2023
La plupart des signaux du monde réel ne sont pas stationnaires, et c’est justement dans l’évolution de
leurs caractéristiques (statistiques, fréquentielles, temporelles, spatiales) que réside l’essentiel de
l’information qu’ils contiennent. Les signaux vocaux et les images en sont des exemples courants [27].
Rappelons que l’analyse de Fourier permet une caractérisation globale du signal (on intègre de -∞ à +∞), on
perd toute localisation temporelle ou spatiale, l’idéal est de faire appel à une transformation qui nous apporte
l’information sur le contenu fréquentiel tout en préservant la localisation (temporelle ou spatiale) afin
d’obtenir une représentation temps/fréquence ou espace/échelle du signal. C’est ainsi que deux théories ont
été élaborées, la transformée de Fourier à fenêtre puis la transformée continue par ondelettes.
Les ondelettes, famille de fonctions déduites d’une même fonction, appelée ondelette mère, par
opérations de translations et dilatations, ont trouvé, de par la puissance de leur théorie, des applications dans
de nombreux domaines aussi variés que les mathématiques (analyse, probabilités, fractales), le traitement du
signal (compression, astronomie, sismique), la physique (mécanique quantique, turbulence).
−∞
Où h(n) est une fenêtre de pondération (Rectangulaire, Bartlett, Hanning, Hamming, Gaussienne etc.)
+∞ 2
x(τ )h (τ − t )e dτ
− j 2πfτ
S x (t , f ) = *
−∞
Il faut relever n’existe pas qu’une seule TFCT puisqu’elle dépend de : La durée de la fenêtre (choisie
pour que le signal soit supposé stationnaire sur cette durée), la forme de la fenêtre (compromis largeur-
hauteur des lobes), le taux de recouvrement entre les fenêtres.
Le rôle de la fenêtre h(t)(dont l’énergie doit valoir 1) est de découper un voisinage de longueur L du point
t, dans lequel le contenu fréquentiel est analysé. On conçoit qu’il y a un compromis entre la longueur L de
104
Traitement Avancé du Signal Master RT M1 2022/2023
h(t) qui représente la résolution temporelle, qui induit une résolution fréquentielle en fe/L, et la capacité de la
TFCT à suivre des modulations plus ou moins rapides. Ces 2 résolutions évoluent en inverse l’une de l’autre.
Il a été montré (principe d’incertitude d’Heisenberg) que la fenêtre qui réalise le meilleur compromis temps-
fréquence est la fenêtre gaussienne.
Exemple : Analyse d’un signal chirp allant de f1 puis f2 et inversement (voir TP n°6)
0.5 0.5
0 0
-0.5 -0.5
-1 -1
0 0.5 1 1.5 2 0 0.5 1 1.5 2
FFT de y1 FFT de y2
30 30
20 20
10 10
0 0
-100 -50 0 50 100 -100 -50 0 50 100
Comme illustré ci-haut, on perd toute localisation temporelle puisqu’on obtient la même TF pour les
deux signaux. Pour y remédier on calcule la TFCT en effectuant la TFD pour différents intervalles et
recouvrements
TFCT avec Tfen=100 Recouv=20%, NF=1024 TFCT avec Tfen=100 Recouv=20%, NF=1024
1.4 1.4
1.2 1.2
1 1
Time
Time
0.8 0.8
0.6 0.6
0.4 0.4
0 20 40 60 80 100 0 20 40 60 80 100
Frequency (Hz) Frequency (Hz)
TFCT avec Tfen=10 Recouv=20%, NF=1024 TFCT avec Tfen=10 Recouv=20%, NF=1024
1.5 1.5
Time
Time
1 1
0.5 0.5
0 20 40 60 80 100 0 20 40 60 80 100
Frequency (Hz) Frequency (Hz)
Les TFDs sur chaque fenêtre glissante fournissent le spectrogramme qui permet d’adapter la TF à la
caractérisation des signaux non stationnaires. On obtient, alors, une représentation temps-fréquence
permettant de localiser la distribution de l’énergie simultanément en temps et fréquence. Rappelons que la
FGE, USTHB [ assiakourgli@[Link] / [Link]
105
Traitement Avancé du Signal Master RT M1 2022/2023
Où m{ ™ š ∗
™ › š C ™/2 ∗
™/2 ›
L'expression C ™/2 ∗ ™/2 peut-être alors considérée comme une corrélation instantanée, pour
laquelle on peut définir une densité spectrale instantanée plus communément appelée Transformée temps-
fréquence de Wigner-Ville :
5œ{ , 6 • C ™/2 ∗
™/2 8ž
›™
Cette transformée conserve bon nombre des propriétes de la TF sauf qu'elle ne conserve pas la propriété de
linéarité :
∑Ÿ
. . 5œ{ , 6 ∑Ÿ
. 5œ{2 , 6 C ∑Ÿ
., . 5œ{2 {‡ , 6
Le second terme montre que cette transformée crée des fréquences (interférences) qui n'ont pas d'existence
réelle. En outre, il fait apparaitre des énergies négatives engeandrant des difficultés d'interprétation. Une
version pseudo lissée permet de résoudre partiellement quelques uns de ces problèmes.
Exemples :
- 8U ¡
5œ{ , 6 6 6%
-( cos 2¥6%
/01¦U ! p/01¦U !
1
5œn , 6 § •H 8U ¡
C 8U ¡
IH
8U ¡
C 8U ¡
I 8ž
›™©
4
/̈ /̈ /̈ /̈
1
5œn , 6 §• 8U ž 8ž
›™ C • = 8U ¡ 8ž
›™ C • = 8U ¡ 8ž
›™ C • 8U ž 8ž
›™©
4
1
5œn , 6 ) 6 6% C = 8U ¡
6 C = 8U ¡
6 C 6 6% *
4
106
Traitement Avancé du Signal Master RT M1 2022/2023
5œn , 6 =
K 6 6% C 6 C 6% L C 2cos 4¥6% 6
1 -
¯ªo,« ¯ • °ª - 3° › ‖ª‖
i i
Exemple : 1.5
1 ³u 0 ´ ´
Haar a=0.5, b=0.5
Haar a=1, b=2
Haar a=3, b=4
La fonction de Haar ª ² 1
1
³u ´ ´1
Amplitude de l ondelette de Haar
0.5
0 iuQQ xy³ 0
³u - ´ ´-C
o
√o
Ainsi ªo,« ª) * ²
-0.5
¡ «
³u -C ´ ´-Ci
o
√o o
√o
-1
107
Traitement Avancé du Signal Master RT M1 2022/2023
Exemples d’ondelettes
!/
ª )1 *
¡/
/µ/
t√ t/
- L’ondelette chapeau mexicain (Paul à l’ordre 2)
!/
ª /µ/ 8¡
t√
- L’ondelette de Morlet
L’ondelette de Morlet est une sinusoïde complexe modulée par une gaussienne. L’ondelette de Paul
décroît plus vite que celle de Morlet et autorise des localisations en temps plus précises. L’ondelette est
obtenue par dérivée d’une gaussienne et permet des localisations temporelles d’une qualité légèrement
inférieure à celle de Paul [25]
En décomposant le signal x(t) sur cette famille, on obtient ainsi les coefficients d’ondelettes5 {,¶ i, -
qui caractérisent le coefficient de la décomposition du signal x(t) dans cette base, soit :
1 - ∗
5 i, - • ªo,«
∗
› • ª- 3 ›
{,¶
√i i
Pour un facteur d’échelle assez grand, la représentation des coefficients d’ondelettes en fonction de b,
la position, donne une représentation de “la forme générale de la fonction”. Par contre un facteur d’échelle
faible correspond à une représentation des singularités.
La transformée en ondelettes est un opérateur linéaire, invariant par translation, et par dilatation. Quelle
que soit l’échelle et quel que soit l’endroit, l’analyse du signal se fait avec la même fonction. La transformée
en ondelettes d’un signal n’est pas unique, elle dépend de l’ondelette mère utilisée. En effet, l’ondelette mère
ψ (t) devra avoir une bonne localisation (nulle en dehors d’un certain intervalle), et devra être oscillante le
nombre de moments nuls correspond au nombre d’oscillations).
On peut montrer que si la fonction analysante (l’ondelette) est correctement choisie, la transformation en
ondelettes est inversible. Le signal x(t) peut être reconstruit après double intégration suivant le facteur
d’échelle a et le paramètre de translation b :
1
• • 5 i, - ªo,«
∗
›i ›-
i {,¶
108
Traitement Avancé du Signal Master RT M1 2022/2023
Exemple d’application : Analyse du signal chirp avec l’ondelette de Haar (Voir TP n°6)
Absolute Values of Ca,b Coefficients for a = 1 2 3 4 5 ...
Les coefficients de l’ondelette sont 127
représentés en fonction du temps où les 120
113
valeurs les plus élevées sont de couleur 106
99
claire. On note qu’un coefficient à une
92
amplitude d’autant plus grande que 85
78
l’ondelette ressemble au signal sur la
scales a
71
portion analysée. Lorsque la fenêtre est 64
57
étroite (onde étroite), on observe les 50
hautes fréquences et lorsqu’elle est large, 43
36
ce sont les coefficients des basses 29
fréquences qui sont élevées. 22
15
8
1
By scale Values of Ca,b Coefficients for a = 1 2 3 4 5 ...
50 100 150 200 250 300 350 4
Toutefois bien que la transformation time (or space) b
0
l’information engendrée par la -2
13 300
d’ondelettes discrètes dyadiques. 11
9 200
250
7 150
5 100
3 50
scales a
1
time (or space) b
i 2 - [2
{ K2 , [2 L 2/ š ªK2 [L ›
1
Avec ces valeurs de a et b, l’équation devient :
Si la fonction x(t) est discrétisée, en supposant une période d’échantillonnage égale à 1, pour des raisons
{ K2 , [2 L 2/ ∑ ªK2 [L
1
de simplicité, l’équation s’écrit alors :
109
Traitement Avancé du Signal Master RT M1 2022/2023
L’ondelette doit être très compacte ainsi, ses coefficients doivent être, pour la plupart, proche de zéro.
Cette condition dépend de quelques paramètres comme la régularité de la fonction, le nombre de moments
nuls et la taille de son support. La compacité du support impose que les filtres (h0) qui vont engendrer la
fonction d’échelle (donnant une approximation grossière du signal) et l’ondelette (h1 fournissant les détails)
soient à réponse impulsionnelle finie.
Cette condition impose d’utiliser des filtres miroirs conjugués, Les 1.5
Filtre passe-bas et passe-haut de décomposition
La transformée en ondelettes discrète sur des bases orthonormées se ramène à des opérations de filtrage
numérique suivies de sous échantillonnage. La reconstruction est parfaite et s’effectue également par des
filtrages numériques précédés de sur-échantillonnage [27]
110
Traitement Avancé du Signal Master RT M1 2022/2023
i | ~ J ℎ% |2 [~i |[~
La décomposition suivant l’algorithme de Mallat est obtenue [27]:
.
› | ~ J ℎ |2 [~i |[~
.
Si on appelle i^ le signal i
Remarque
et › le signal ›
après interpolation (insertion d'un 0 entre les 2 échantillons successifs)
^
après interpolation alors [27]:
i| ~ J“6% | [~i^ |[~ C 6 | [~›^ |[~”
.
puisque les trois autres filtres se calculent à partir de ce dernier (filtres de reconstruction identiques aux filtres
d’analyse mais retournés dans le temps).
\ -% - C- -? ? ?
-? C - - C -% ?
˜% -% C - C- C -? ? ?
-? C - C- C -% ?
˜ -% - C- -? ?
-% - C- -? ?
Ainsi, considérons par exemple l’ondelette de Daubechies à 6 coefficients telle que la fonction d’échelle h
111
Traitement Avancé du Signal Master RT M1 2022/2023
soit :
ℎ% = [0.2352 0.5706 0.3252 -0.0955 -0.0604 0.0249],
ℎ est obtenu en retournant ℎ dans le temps et en inversant le signe des coefficients impairs, ce qui nous
donne :
ℎ = [-0.0249 -0.0604 0.0955 0.3252 -0.5706 0.2352]
Filtres de reconstruction identiques aux filtres d’analyse mais retournés dans le temps):
f0= [0.0249 -0.0604 -0.0955 0.3252 0.5706 0.2352],
f1= [0.2352 -0.5706 0.3252 0.0955 -0.0604 -0.0249 ]
Ci – dessous sont donnés les coeffcients d’ondelettes de Daubechies 3 avec le tracé des pôles et zéros ainsi
que des réponses en fréquence (Voir TP n°6)
Filtre passe-bas de décomposition Filtre passe-haut de décomposition Filtre passe-bas de décomposition Filtre passe-haut de décomposition
1 1
1
Imaginary Part 1
Imaginary Part
0.5 0.5 0.5
3 5 5 3
0 0
0 0
-0.5
-1
-0.5 -0.5 -1
0 2 4 6 0 2 4 6 -1 10 2 -1 -0.5 0 0.5 1
Real Part Real Part
Filtre passe-bas de Reconstruction Filtre passe-haut de Reconstruction
1 1 Filtre passe-bas de Reconstruction Filtre passe-haut de Reconstruction
1
1
Imaginary Part
Imaginary Part
0.5 0.5 0.5
3 5 5 3
0 0
0 0
-0.5
-1
-0.5 -0.5
0 2 4 6 0 2 4 6 -1
-1 -0.5 0 0.5 1 -2 -1 0 1
Real Part Real Part
Filtre passe-bas et passe-haut de décomposition Filtre passe-bas et passe-haut de Reconstruction
1.5 1.5
1 1
0.5 0.5
0 0
0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 0.5
Comme les coefficients ne sont pas symétriques, le déphasage ne sera pas linéaire. Rappelons que la
seule ondelette orthogonale qui soit symétrique est l’ondelette de Haar.
112
Traitement Avancé du Signal Master RT M1 2022/2023
Exemple d’application 1 : On suppose que le signal à analser est une rampe et que l’ondelette est celle de
Haar.
Décomposition
Niveau 0 0,1,2,3,4,5,6,7
,, , ,
¼ ½ ?
√ √ √ √ √ √ √ √
Niveau 1
, , ,
. .
,, , ,
¼ ½ ?
√ √ √ √ √ √ √ √
Niveau 1
, , ,
, 0, 0, , 0, , 0,
¼ ½ ?
√ √ √ √ √ √ √ √
0, , 0, , 0, ,0 ,0, ,0
1 1 5 5 9 9 13 13
, , , , , , ,
2 2 2 2 2 2 2 2
, , , , , , ,
Niveau 0 0,1,2,3,4,5,6,7
Exemple d’application 2 :
On suppose que le filtrage passe-bas ℎ% est un opérateur de moyennage entre 2 coefficients et que le passe-
haut ℎ est un opérateur qui encode la demi-différence entre 2 coefficients [cours 2013] tels que ℎ% ={0.5, 0.5}
, ℎ ={-0.5, 0.5} , 6% =2*{0.5, 0.5} , 6 =2*{0.5, -0.5}
113
Traitement Avancé du Signal Master RT M1 2022/2023
L’opérateur ℎ% donne une représentation approximative alors que ℎ donne une représentation plus
Exemple 3 Voici le résultat de la décomposition du signal ecg, ne sont gardés à la fin que la dernière
approximation et tous les niveaux de détails grâce auxquels on pourra reconstruire le signal original (TP n°6)
ecg original ecg reconstruit
1000 1000
0 0
-1000 -1000
0 500 1000 1500 2000 0 500 1000 1500 2000
Coefficients Détail N 1 Coefficients Détail N 2
20 100
0 0
-20 -100
0 500 1000 1500 0 200 400 600
Coefficients Détail N 3 Coefficients d Approximation N 3
500 5000
0 0
-500 -5000
0 100 200 300 0 100 200 300
114
Traitement Avancé du Signal Master RT M1 2022/2023
Remarques :
• Les filtres de Daubechies conduisent à des filtres à réponse impulsionnelle finie (RIF) mais ne sont pas
linéaires en phase.
• On peut construire des bases d’ondelettes conduisant à des filtres RIF linéaires en phase telles les B
splines. Il s’agit des bases biorthogonales où familles de décomposition et reconstruction sont
différentes et orthogonales uniquement entre elles [29]:
\% 6 ˜%∗ 6 C \% 6 C 6 /2 ˜%∗ 6 C 6 /2 2 et \ 6 ˜ ∗ 6 C \ 6 C 6 /2 ˜ ∗ 6 C 6 /2 2
\% 6 ˜ ∗ 6 C \% 6 C 6 /2 ˜ ∗ 6 C 6 /2 0 et \ 6 ˜%∗ 6 C \ 6 C 6 /2 ˜%∗ 6 C 6 /2 0
h1(n) = -(-1)n f0 (L - 1 - n)
f1(n) = (-1)n h0 (L - 1 - n) Filtre passe-bas de décomposition Filtre passe-haut de décomposition
1
Filtre passe-bas de décomposition Filtre passe-haut de décomposition
Imaginary Part
Imaginary Part
1 1 1
0.5
0.5 3 7 5 3
0.5 0 0
0 -0.5
0 -1
-0.5 -1
-1 0 1 2 -1 -0.5 0 0.5 1
-0.5 -1
0 2 4 6 8 0 2 4 6 8 Real Part Real Part
Filtre passe-bas de Reconstruction Filtre passe-haut de Reconstruction
Filtre passe-bas de Reconstruction Filtre passe-haut de Reconstruction
0.6 1 1
Imaginary Part
Imaginary Part
1
0.5 0.5
0.4
3 5 7 3
0 0 0
0.2
-0.5 -0.5
-1
0 -1 -1
0 2 4 6 8 0 2 4 6 8 -1 -0.5 0 0.5 1 -2 -1 0 1
Real Part Real Part
115
Traitement Avancé du Signal Master RT M1 2022/2023
On notera la symétrie des réponses impusionnelles des filtres passe-bas et l'anti-symétrie des filtres
passe-haut d'analyse et de reconstruction conduisant dons à un déphasage linéaire. On notera également la
dualité entre les filtres d'analyse et reconstruction.
Remarque : Le lifting scheme permet une implantation très simple des décompositions en ondelettes
et de leurs opérations inverses en employant la factorisation polyphases [21]
Les ondelettes sont bien adaptées au débruitage des signaux. Le principe est simple, on décompose le
signal bruité et l’on force à zéro les coefficients qui ne passent pas un seuil. Puis l’on reconstruit le signal
(Voir TP n°6). De la même façon on peut compresser le signal.
Pour une image, on appliquera deux filtres le long des lignes l’un passe-haut (H), l’autre passe-bas (L),
ensuite, on ne considère qu’une colonne sur deux et on applique à nouveau les 2 filtres. En ne considérant
qu’une ligne sur deux, on obtient 4 images de taille réduite de moitié. En réitérant ce processus sur l'image
d’approximation (LL), on obtient une pyramide d’images.
Cet exemple donne le résultat d’une compression jpeg et jpeg2000 (Ondelettes) avec un bpp de 0.2
116
Traitement Avancé du Signal Master RT M1 2022/2023
84 ¡
vwxy … %
1. Déterminer la TFCT du signal suivant (h est une fenêtre porte de largeur T):
À 8/ ¡
vwxy ® %
1 ³u 0 ´ ´
3. Soient les ondelettes suivantes :
- La fonction de Haar ª ² 1 ³u ´ ´1
0 iuQQ xy³
!/
- L’ondelette chapeau mexicain ª Â 1 /
4. Soit le filtre passe bas de décomposition suivant : ℎ% {-0.1294, 0.2241, 0.8365, 0.4830}, déterminer le
filtre passe-haut de décomposition et les filtres de reconstruction
8. On décompose un signal x(n) par l'ondelette de Haar ℎ% “ , ” suivant le schéma ci-dessous. Les
√ √
8 premières valeurs du signal x(n) sont ={4, 4, 6, 8, 8, 10, 9, 11}
↓2 ↓2
i i
H0(z) H0(z)
X(z)
↓2 ↓2
›
H1(z) H1(z)
›
FGE, USTHB [ assiakourgli@[Link] / [Link]
117
Traitement Avancé du Signal Master RT M1 2022/2023
↓2 ↓2
›
›
H (z) H1(z)
1
Solutions
³u ÈK 6 6 L vwxy C …
1.
8 84 ¡
⎧ %
⎪ ³u ÈK 6 6 L 8 8/ ¡
vwxy ®
⎪
%
C 2* ³u È -) C 2* 6
C
) 6 3 C
8 84 ) 0 *
⎨ 0
2 4
0
⎪
C * ³u È -)2 C * 6
C
⎪) 2 6 3 vwxy … …
8 8/ ) 0 C *
2 4
0 0 %
⎩
2. 5œ{ , 6 |Á | 6 6 C |Á | 6 6 C 2Á Á cos 2¥ 6 6 )6 *
84 8/
TFdeψ(t)
0.4 0.06
0.4
ψ(t)
ψ(t)
0 0.05
0.3 0.2
0.04
-0.5 0.03
0
0.2
0.02
-0.2
-1 0.01
0.1
-0.4 0
-10 -5 0 5 10 -2 -1 0 1 2
-1.5 0 t f
0 0.5 1 -1000 -500 0 500 1000
118
Traitement Avancé du Signal Master RT M1 2022/2023
Reconstruction
Energie
Exercices Supplémentaires
1. Wigner-Ville
√ √
- Donner le schéma de reconstruction permettant de passer d’un niveau à un autre.
- Reconstruire le signal à partir du niveau 1.
- Donner le signal décomposé aux niveaux 2 puis 3.
- Vérifier que l’énergie se conserve à tous les niveaux (1, 2 et 3)
- Proposer un moyen de compresser ce signal.
119
Traitement Avancé du Signal Master RT M1 2022/2023
- L’énergie se conserve-t-elle ?
- Reconstruire le signal à partir du niveau 2.
4. A la transmission d'un signal ecg, il a été affecté d'un bruit issu de la combinaison de 3 sinusoïdes (voir
figure ci-dessous : signal bruité + Module de sa TFD).
En analysant le signal par les ondelettes discrètes dyadiques, nous avons obtenu les graphes suivant:
- Identifier pour chacun des 3 signaux : le signal (approximé ou détail), le niveau de décomposition
- Pour chacun des 3 signaux, esquisser le module de la TFD
- Proposer une solution pour supprimer le bruit (explication détaillée)
Solutions :
2. examen 17/18
3. Ratt 17/18
4. Examen 20/21
120
Traitement Avancé du Signal Master RT M1 2022/2023
1. Démo
# -*- coding: utf-8 -*-
import numpy as np; import [Link] as plt;import [Link] as sp;
NF = 1024; Te=0.005; fe=1/Te;N=round(2/Te);
t = [Link](0, 2, N);
x1= [Link](t, f0=50, t1=2, f1=1, method='linear');
x2= [Link](t, f0=1, t1=2, f1=50, method='linear');
""" TFD"""
TFx1 = [Link](x1,NF); TFx1 = [Link](TFx1);
TFx2 = [Link](x2,NF); TFx2 = [Link](TFx2); freq = [Link](-NF/2,NF/2)*fe/NF;
[Link](1)
[Link](221); [Link](t, x1); [Link]('Chirp de 50 à 1 Hz'); [Link](222); [Link](freq, [Link](TFx1));
[Link](223); [Link](t, x2); [Link]('Chirp de 1 à 50 Hz'); [Link](224); [Link](freq, [Link](TFx2));
"""Spectrogramme : TFCT"""
f, tt, Sxx1 = [Link](x1, fe, nperseg=100, noverlap=20);
f, tt, Sxx2 = [Link](x2, fe, nperseg=100, noverlap=20);
[Link](2)
[Link](221); [Link](t, x1); [Link](222); [Link](tt, f, Sxx1,shading='auto')
[Link]('Frequence [Hz]'); [Link]('Temps [sec]'); [Link]('Spectrogramme Tfen=100, TRec=20');
[Link](223); [Link](t, x2); [Link](224); [Link](tt, f, Sxx2,shading='auto')
[Link]('Frequence [Hz]'); [Link]('Temps [sec]'); [Link]('Spectrogramme Tfen=100, TRec=20');
f, tt, Sxx1 = [Link](x1, fe, nperseg=20, noverlap=4);
f, tt, Sxx2 = [Link](x2, fe, nperseg=20, noverlap=4);
[Link](3)
[Link](221); [Link](t, x1); [Link](222); [Link](tt, f, Sxx1,shading='auto')
[Link]('Frequence [Hz]'); [Link]('Temps [sec]'); [Link]('Spectrogramme Tfen=20, TRec=4');
[Link](223); [Link](t, x2); [Link](224); [Link](tt, f, Sxx2,shading='auto')
[Link]('Frequence [Hz]'); [Link]('Temps [sec]'); [Link]('Spectrogramme Tfen=20, TRec=4');
"""Wigner-Ville """
from [Link] import WignerVilleDistribution
WVx1=WignerVilleDistribution(x1,timestamps=t)
WVx2=WignerVilleDistribution(x2,timestamps=t)
tfr_x1, t_x1, f_x1 = [Link]()
tfr_x2, t_x2, f_x2 = [Link]()
[Link](4)
[Link](221); [Link](t, x1); [Link](222); [Link](t_x1, f_x1*fe, [Link](tfr_x1),shading='auto')
[Link]('Frequence [Hz]'); [Link]('Temps [sec]'); [Link]('Wigner-Ville');
[Link](223); [Link](t, x2); [Link](224); [Link](t_x2, f_x2*fe, [Link](tfr_x2),shading='auto')
[Link]('Frequence [Hz]'); [Link]('Temps [sec]'); [Link]('Wigner-Ville');
121
Traitement Avancé du Signal Master RT M1 2022/2023
2. A faire
122
Traitement Avancé du Signal Master RT M1 2022/2023
A l'aide de pywt.dwt2, décomposer et reconstruire une image sur plusieurs niveaux, mettre des détails
à zéros et commenter.
Supposons que l'on veuille connaître le poids 'p' ou la taille 'T' que devrait avoir un enfant de 4 ans.
moyenne sur ces mesures nous fournira une estimation v̂ du poids une estimation Š de la taille. Cet
Pour ce faire, on prendra un échantillon d'enfants de 4 ans que l'on pèsera et dont on mesurera la taille. Une
échantillon devra bien sûr être le plus grand possible et le plus hétérogène possible (plusieurs ethnies).
L'estimation de la durée moyenne d'une communication téléphonique est très utile aux opérateurs de
téléphonie pour ajuster les prix et offres en conséquence ou pour la gestion du trafic. Autre exemple : les
intentions de vote. Dans chacun des cas, un échantillon est considéré pour représenter la population.
1. Estimateur
Dans de nombreuses situations, on ne dispose pas directement de mesures sur la variable d'intérêt
mais que d'une observation liée à notre variable inconnue. Le but des techniques d’estimation est d’utiliser
123
Traitement Avancé du Signal Master RT M1 2022/2023
les observations pour extraire de l’information sur la grandeur d’intérêt. Ceci doit fait de la meilleure façon
qui soit et implique plusieurs choix dont la relation mathématique entre l'observation et la variable à estimer.
Un estimateur est donc une valeur Š calculée sur un échantillon Y tiré au hasard (dit aussi
observations), ce qui fait de Š une variable aléatoire possédant une moyenne et une variance.
Exemple
Considérons que l’on veuille estimer la valeur d’une résistance R (TP n°5). On procède comme suit : On
réalise N expériences indépendantes où l’on mesure le courant I(k) passant par la résistance et la tension U(k)
à ses bornes. Ces mesures de courant et de tension sont toutes bruitées par un bruit dont on ignore la densité
de probabilité.
Variable à estimer : X= R
Bien que X et Y ne soient pas supposés avoir une densité de probabilité, on peut construire différents
estimateur de R en fonctions des observations U(k) et I(k)
1.7
N
U ( k ) I (k )
R1
R2
1.6 R3
Estimateur 1 : Rˆ1 = k =1
N 1.5
I
k =1
2
(k ) 1.4
1.3
N
1 U (k )
Estimateur 2 : Rˆ 2 =
N
k =1 I ( k )
1.2
1.1
N
U (k ) 1
Estimateur 3 : Rˆ 3 = k =1
N
0.9
I (k )
k =1
0.8
0.7
0 10 20 30 40 50 60 70 80 90 100
Ces figures montrent l’évolution de l’estimée d’une résistance R=1 en fonction du nombre de mesures
pour U=I=1. Ci-contre une estimation pour N=100 et ci-bas pour N=1000
On remarque que le premier estimateur donne des valeurs comprises entre 0.9 et 1. Le deuxième
fournit une estimation proche de 1 et 1.2 Ils sont donc tous deux fortement biaisés. Le troisième estimateur
est quant à lui non biaisé.
Ainsi, pour pouvoir déterminer d’une façon constructive une règle d’estimation, il faut définir un
critère qui évalue la qualité des résultats, et définir l’estimée comme l’application de Y en X qui optimise ce
critère (minimiser une erreur ou une distance, dite aussi fonction de coût). Deux mesures de qualité de
l'estimation sont largement utilisées, il s'agit du biais et de la variance.
Le biais est la moyenne de l’écart et la variance est la puissance de l’écart (mesure les fluctuations de
l’estimateur autour de la valeur moyenne).
Le biais indique la valeur moyenne de l’erreur d’estimation, trois cas sont possibles:
• E [ x̂ } =x ou µx pour toutes les valeurs possibles du paramètre. On dit alors que l’estimée est non-biaisée ;
• E [ x̂ } =x+b ou µx+b. Si b est indépendant de x, dans ce cas l’estimateur a un biais constant et connu, qui
peut toujours être éliminé ;
125
Traitement Avancé du Signal Master RT M1 2022/2023
La variance doit être aussi petite que possible, de façon à que l’estimée soit concentrée autour de la
vraie valeur du paramètre.
Plus b et σ2 sont faibles, meilleur est l’estimateur (estimateur non biaisée à variance minimale).
Malheureusement, la diminution de l’un provoque l’augmentation de l’autre. Dans ces conditions, un
estimateur biaisé pourrait être préférable, si le biais est faible, à un estimateur non biaisé mais de très grande
variance. C'est le compromis biais-variance , on peut démontrer que b 2=E { x̂ −x}2 - σ2 où E { x̂ −x}2 est l'erreur
quadratique moyenne.
Spectre réel
Spectreestimé
Xˆ ( f )
X(f )
Variance
Biais
On dit qu’un estimateur est consistant s’il tend vers la vraie valeur du paramètre quand le nombre
d’observations tend vers infini (comportement asymptotique) :
lim bN = lim σ N = 0
N → +∞ N → +∞
Si l’estimateur est non biaisé, nous avons une borne inférieure de la variance (borne de Cramer Rao). Si
la variance atteint la borne de CR, on dit que l’estimateur est efficace. Autrement dit, un estimateur non-biaisé
est appelé efficace si sa variance dans est plus petite que celle de n’importe quelautre estimateur non-biaisé.
Remarque: Dans ce qui suit, on parlera de variables xi iid (indépendantes et identiquement distribuées).
L'exemple classique est celui du lancé de dé ou de pièce de monnaie : Les variables aléatoires représentent
chaque résultat des lancers (0 pour face et 1 pour pile) suivent toutes la même loi de Bernoulli. Bien que les
lancers soient successifs, le résultats d'un lancé donné ne dépend pas des résultats précédents et n'aura
aucune influence sur les prochains lancés (pas de lien de dépendance)
Ce sont des variables aléatoires qui suivent toutes la même loi de probabilité (donc même moyenne et même
variance) et sont indépendantes.
- E{xi}=m ∀i
- E{(xi-m)2}=σ2
Exemple d'application
Supposons que l'on observe N échantillons indépendants Y1,Y2, …….., Yn d'une variables aléatoire Y de
moyenne m= E{Yi} et de variance σ2=E{Yi2}-m2.
FGE, USTHB [ assiakourgli@[Link] / [Link]
126
Traitement Avancé du Signal Master RT M1 2022/2023
1. On désire estimer X=m (déterministe) en supposant σ2 connue. Pour cela, on étudie les deux estimateurs
suivants :
1 n 2 n
X1 = Yi
n i =1
X2 = [Link]
n ( n + 1) i =1
- Biais de X1= E {Xˆ 1 }− X=0 - Variance de X1=σ2/n
2. On désire estimer X= σ2 en supposant m connue. Pour cela, on étudie les deux estimateurs suivants :
1 n 1 n
X1 =
n i=1
(Yi − m) 2 et X1 =
n i =1
(Yi − µY ) 2
Les quantités E{x} ,E{x}2, ,Sxx(f),... sont impossibles à calculer sur un ordinateur car elles nécessiteraient un
nombre de points infini. Sur calculateur, on ne dispose que d'une séquence discrète et finie de N points. En
réalité, on calcule des estimées de ces grandeurs en supposant généralement le signal stationnaire et
ergodique. On remplace alors le calcul des moyennes statistiques par des moyennes temporelles.
[Link] de Moyenne
Soit N échantillons (x0, x1, …..,xN-1) indépendants et identiquement distribuées (même loi avec même
paramètre) d’un signal aléatoire stationnaire. On définit l’estimée m̂ de E{x}
N −1
1
mˆ =
N
x
k =0
k
On remarque donc qu'on remplace la moyenne statistique par la moyenne temporelle sur la séquence finie.
Cet estimateur est non biaisé. En effet, du fait que le signal soit iid ,∀k;E [xk} = m, alors :
N −1
E {mˆ } =
1
E{x } = N .N .m = m
1
k
b{mˆ } = 0
N k =0
1 N −1 1 N −1
2
2
mˆ {
σ = var{mˆ } = E (mˆ − E{mˆ }) 2
} = E xk − E xk
N k =0 N k =0
1 N −1 1 N −1 1 N −1
2
2
= E xk − E{xk } = 2 E ( xk − E{xk })
N k =0 N k =0 N k =0
127
Traitement Avancé du Signal Master RT M1 2022/2023
E{( xk − E{xk }) 2 } = 2
1 N −1 1 N −1
2 σ
σ m2ˆ = 2
=σ2 / N
N k =0 N i =0
Cet estimateur est consistant, puisque pour N→∞, biais et variance tendent vers 0.
B. Estimateurs de variance
Selon que l’on connaisse ou non la moyenne, on utilisera l’un des deux estimateurs suivants :
N −1
1 1 N −1
σˆ 1 = (x − m) 2 σˆ 2 = ( x k − mˆ ) 2
2 2
N − 1 k =0
k
N k =0
{ } E{(x }
N −1
− m) = σ 2
1 bσˆ 2 = 0
E σˆ 1 =
2 2
- Pour le premier, on a : k 1
N k =0
E {σˆ } = E {( x } { }
N −1
1 N −1
− mˆ ) = E ( xk − m + m − mˆ )
1
(N − 1)
2 2 2
- Pour le deuxième, on a : 2
k =0
k
(N − 1) k =0
1 N −1 1 N −1 1 N −1
{ }
2 2
1 N −1
E σˆ 2 =
(N − 1)
E xk − m + m − x j = E (xk − m) − N ( x j − m)
2
k =0 N k =0 ( N − 1) k =0 k =0
Sachant que E {( x k − m )} = 0 , on a :
1 N −1
{ } { } { }
N −1 2
1
E ( xk − m ) + ( x j − m) − E ( xk − m )
2
E σˆ 2 =
2 2 2
(N − 1) k =0
N k =0 N
E{(x }
1 N −1 2 N −1
2σ 2
1 N −1 2 σ 2 2σ 2 Nσ 2 N − 1
− m)
1
(N − 1)
= σ + 2 2
− = .σ + − = =σ
2
( N − 1) i=0 N ( N − 1) N
i
i =0 N i =0 N N
Les deux estimateurs ne sont pas biaisés puisque on obtient la vraie valeur à chaque fois.
C. Estimateurs de la corrélation
Dans le cas discret, la fonction d'autocorrélation d'un signal aléatoire x; supposé ergodique, est définie par :
1 n=+ N
Rxx (k ) = E{x(n) x * (n − k )} = E{x(n) x * (n − k )} = lim x(n) x * (n − k )
N →∞ 2 N + 1 n=− N
128
Traitement Avancé du Signal Master RT M1 2022/2023
A partir d'un échantillon (x0, x1, …..,xN-1) indépendants et identiquement distribuées, plusieurs estimateurs
sont alors possible:
N −1 N −1
1 1
Rˆ xx (k ) =
N
x ( n) x * ( n − k ) Rˆ xx (k ) =
N −k
x( n) x * ( n − k )
n=k
n =k
Le premier estimateur est biaisé car on ne tient pas compte du nombre d'échantillons disponibles qui
varie avec le pas k. Le deuxième en tient compte, il est non biaisé. En effet :
{
E Rˆ xx (k ) = } 1 N −1
N − k n=k
E{x( n) x * ( n − k )} =
1 N −1
Rxx (k ) = Rxx (k )
N − k n=k
On peut démontrer que la variance de cet estimateur tend vers 0 quand N tend vers l'infini. Cet estimateur
est donc consistant.
On observe par contre que lorsque k approche du nombre d'échantillons N, la variance de l'estimateur
non biaisé devient excessive (voir exemple pour un bruit blanc et une sinusoïde).Alors que celle de
l'estimateur biaisé varie beaucoup moins. C'est une des raisons pour lesquelles cet estimateur est souvent
utilisé par la suite, malgré son biais.
a.Méthode du périodogramme
Supposons que l'on dispose d'une séquence de N points du signal aléatoire : x [ x0,x1,...,xN-1}.
Soit donc la séquence yk obtenue à partir de la séquence xk pondérée par la fenêtre fk : yk=[Link]
La transformée de Fourier calculée sur la séquence finie sera donc convoluée avec le spectre de la fenêtre
rectangulaire c'est-à-dire un sinus cardinal. Rappelons que les propriétés spectrales du sinus cardinal ne sont
pas bien adaptées à l'analyse spectrale du signal (lobes secondaires importants). Pour y remédier, on emploie
des fenêtres de pondérations (Hamming, Hanning, Kaiser,etc.)
( )
N −1 N −1 N −1
1
=
N
y y
k =0 l =k
l l −k e −2πjfk = Rˆ yy (k )e − 2πjfk = TF Rˆ yy (k )
k =0
FGE, USTHB [ assiakourgli@[Link] / [Link]
129
Traitement Avancé du Signal Master RT M1 2022/2023
{ } N −1 N −1 N −1
1
E Sˆ xx ( f ) = E Rˆ yy (k )e − 2πjfk = E y y l l −k e − 2πjfk
k =0 N k =0 l = k
{ } 1 N −1 N −1 N −k
N −1
E Sˆxx ( f ) == E{yl yl −k }e−2πjfk = Ryy (k ) e−2πjfk = S yy ( f ) * N sin c( fN ) 2
N k =0 l = k k =0 N
Afinde diminuer la variance de cet estimateur, on peut utiliser un périodogramme moyenné. Cela consiste à
séparer le signal en K tranches (de longueur N/K), à calculer le périodogramme sur chaque tranche et à faire
la [Link] fait des K moyennages, la variance est presque divisée par K : néanmoins, les tranches étant
plus courtes, la résolution diminue(∆f=fe/N→fe/K). En pratique, un taux de recouvrement de 40% donne de
bons résultats. Au delà, l’indépendance entre les tranches n’est plus respectée.
[
Sˆ xx ( f ) = TF Rˆ xx (k ) ] 0.5
unbiased
biased
0.4
0.25
0.15
0.1
2
0.05
xx r
u0
séquence xk. Dans le cas d'un modèle auto- 1+ a e
k =1
k
− i 2πfk
régressifxk = a1xk-1+....+arxk-r+uk
1. Estimateurs
Le but des techniques d’estimation est d’utiliser les observations X (aléatoires) pour extraire de l’information
sur la grandeur d’intérêtY. Suivant que Y soit déterministe ou aléatoire, on fera appel à différentes
techniques. Ainsi :
Pour estimer une variable certaine, dans le cas le plus général, on souhaite que l’estimateur soit non
biaisé, on cherche donc l’estimateur à minimum de variance dont le plus connue est le maximum de
130
Traitement Avancé du Signal Master RT M1 2022/2023
vraisemblance (MV) basée sur l'emploi de la probabilité des valeurs observées p(X/Y). Si on autorise un
biais, on cherche alors à minimiser l’erreur quadratique moyenne.
Lorsque Y la variable aléatoire à estimer est aléatoire, on aura recours aux estimateurs Bayesiens.
Plusieurs cas de figures sont alors possibles :
- Lorsqu'on suppose p(Y) connue, on utilisera en fonction de la fonction de coût à minimiser, différents
estimateurs tels celui du maximum à postériori (MAP).
- Si seules sont connues les moments d’ordre 1 et 2 de p(Y)etp(X), on utilise alors souvent l’estimateur linéaire
non-biaiséà variance [Link] filtre de Wiener en est un exemple.
Dans le cas où on ne possède aucune information statistique sur X et Y et que la seule dont on dispose
est que X est une mesure bruitée de fct(Y), on adoptera l'estimateur des moindres carrés pour estimer Y qui
s'affranchit de tout cadre probabiliste.
1. Les éléments d’une population possèdent un caractère X qui suit une loi de probabilité dont la densité est
donnée par :
³u 0 … … Ð
.{ Ï
vÎ À Î>
0 ³u w
Où k =cste et θ>0 est le paramètre inconnu.
- Déterminer la constante k.
N
1
- Montrer que l’estimateur
N
X
i =1
i de θest biaisé et en déduire un non biaisé.
2. Soit (X1,X2,…….XN) avec N>10 un échantillon aléatoire issu d’une population suivant une loi de Bernoulli
de paramètre p. p(x)=p.δ(x-1)+(1-p).δ(x) . Considérons les trois estimateurs du paramètre p :
) N ) 1 N /2 ) 1 10
p1 = X i − 1 / N
i =1
p2 = X 2i
N / 2 i =1
p3 = X i
10 i =1
- Comparer les trois estimateurs (consistance)
3. Soit (X1,X2,…….XN) une population de moyenneµet de variance σ2. On considère les deux estimateurs de la
moyenne :
1 N N N
µ1 = Xi
N i =1
µ 2 = ai X i / ai
i =1 i =1
ai ≥ 1
4. Soient (X1,X2,…….XN), N variables aléatoires indépendantes suivant des lois de poisson de paramètre λ.
Soient les estimateurs suivants du paramètre λsuivant :
N
x i
λˆ2 =
x1 + x N
λ̂1 = i =1
2
N
FGE, USTHB [ assiakourgli@[Link] / [Link]
131
Traitement Avancé du Signal Master RT M1 2022/2023
5. On considère un processus de la forme X(n) = θ + W(n) où W(n) est un processus gaussien de moyenne 0
et de variance 1 tel que W(n) et W(j) sont indépendants si n≠j. On suppose que θ suit une loi N (0,σ2)
indépendante de W(n) n∈Z.
- Caractériser le filtre de Wiener permettant d’estimer θ à partir de Xn, Xn−1, ..., Xn−p+1.
- Calculer les coefficients du filtre.
6. On considère un processus de la forme X(n) = b(n) + αb(n-1) + W(n) où W(n) est un processus gaussien
de moyenne 0 et de variance σ2 tel que W(n) et W(j) sont indépendants si n≠j. On suppose que b(n) est une
variable aléatoire uniforme à valeurs dans {−1,1}, indépendante de(W(n) )n∈Z et que de même, b(n) et
b(j) sont indépendants si n ≠ j (P(b(n) = 1) = P(b(n) = −1) =1/2.
- Construire un filtre qui estime b(n) .
7. On désire dé-bruiter un signal de parole z(n) corrompu par du bruit additif b(n) indépendant du signal
sonore et ce, par filtrage de Wiener. On suppose connues quelques valeurs d’autocorrélation pour les deux
signaux telles que :
RZZ[0] = 1.5; RZZ[1] = 0.5; RZZ[2] = 0.25; RZZ[3] = 0.125 ; RZZ[4] = 0.0625
Rbb[0] = 1; Rbb[1] = 0.25; Rbb[2] = 0.0625; Rbb[3] = 0.015625
-Pourquoi ne peut-on pas employer un filtre adapté ou un filtre moyenneur ?
-Donner les équations de Wiener-Hopf permettant d’estimer z(n).
ˆ(n).
-Déterminer le filtre de Wiener d’ordre 2 permettant de retrouver le signal utile z
-Exprimer H(z).
Solutions
5 N
1. k=4 b=(4θ/5)-θ θ= xi
4N i=1
2. p(x)=pδ(x-1)+(1-p)δ(x) µx=p σx2=p(1-p)
- p1 b=1/N σ =p(1-p)/N
2 (=0 N→∞) consistant
- p2 b=0 σ2=2.p(1-p)/N (=0 N→∞) consistant
- p3 b=0 σ2=p(1-p)/10 (≠0 N→∞) non consistant
p1meilleur estimateur pour N→∞
2
n
n
3.µ1 voir cours µ2 : b=0 ,σ2=
i =1
ai σ / ai
2 2
i =1
FGE, USTHB [ assiakourgli@[Link] / [Link]
132
Traitement Avancé du Signal Master RT M1 2022/2023
4. bλ1=0 bλ2=0 σλ12=σ2/N σλ22=σ2/2 Le premier (variance plus petite pour N>2)
5. Rxx(0)=1+σ2 Rxx(k>0)= σ2 Rθx(k)=σ2 bi=σ2/(Nσ2+1)
6. Rxx(1:3)=(1+α2+σ2, α, 0) Rbx(1:3) =(1, 0, 0)
7. µz=0, µb=0, signal (parole) aléatoire,
8. µs=0, µb=0 , Sb(f)=0.8. Quand signal utile et bruit occupent même plage de fréquences.
b0=0.8/2.6 et b1=0.
Exercices supplémentaires
1. Soit X une variable aléatoire qui suit une loi Binomiale B(16; p), où p est un paramètreinconnu. On se
propose d'estimer p a l'aide d'un échantillon de taille n de X : (x1; .....;xN). On se propose d'étudier les
estimateurs suivants :
1 N 1 N
1 1 N
pˆ 1 = Xi
N i=1
pˆ 2 =
16 N
Xi
i =1
pˆ 3 = X1 +
16
Xi
N − 1 i =2
- Lesquels sont sans biais, consistants.
- Lequel doit-on employer de préférence?
On rappelle que pour une loi Binomiale B(n; p) de paramètre n et p, moyenne et variance sont les suivants
m = np σ 2 = np (1 − p )
2. Soit D la variable aléatoire représentant la durée d'une communication téléphonique. On suppose que la
[Link] , l ,.......lŸ les N variances obtenues à partir de N échantillons aléatoires i.i.d. de tailles respectives
n1, n2, ......nN.
On considère l'estimateur de la variance suivant :l• 4 t4 / t/Ñ ……….………… Ò tÒ
/ / /
4 4 …….. Ò
- Est-ce un estimateur biaisé?
- S'il est biaisé, en proposer un non biaisé.
4. Soient les variables aléatoires X1 et X2 i.i.d. de moyenne µ et de variance σ2. Comparer les deux estimateurs
de µ suivants:
- ĉ ĉ
Ó4 Ó/ oÓ4 «Ó/
o «
avec a et b réels.
[Link] (x1, x2, .........xN) des échantillons indépendants correspondant à une population dans la moyenne est
µ et de variance σ2 . On considère les 2 estimateurs suivants de la moyenne:
Ô Ÿ ∑Ÿ
c q q c
Ô i C 1 i où 0≤a≤1
• Calculer le biais et la variance de chaque estimateur.
• Etudier leur consistance.
• Lequel vous semble préférable (Justifier) ?
6. On considère un processus de la forme X(n) = S(n) − 2S(n-1) + S(n-2) + W(n) où W(n) est un processus
gaussien de moyenne 0 et de variance σ2 tel que W(n) et W(j) sont indépendants si n≠j. On suppose que
133
Traitement Avancé du Signal Master RT M1 2022/2023
S(n) suit une loi N (0,1) indépendante de (W(n) ) n∈Z et que de même, S(n) et S(j) sont indépendants si n ≠
j.
- Donner l’équation de Wiener–Hopf permettant de calculer les coefficients du filtre(anti causal) de Wiener
d’ordre 3 permettant d’estimer S(n-2) à partir de X(n) , X(n-1) et X(n-2) .
-Vérifier que si σ = 0, la solution est : S(n-2) = -( X(n) + 3X(n-1) + X(n-2) )/5
[Link] considère une observation x(n)=s(n)+w(n) où w(n) est un bruit blanc de variance σ2 indépendant de
x(n). On suppose aussi que le signal utile x(n) est centrée et réduit et que les x(n) sont indépendants.
- Dans quel cas, utilise-t-on le filtre de Wiener ?
- Déterminer le filtre de Wiener d’ordre 3 permettant de retrouver le signal utile sˆ(n).
- Donner l’expression de sˆ(n)en fonction de x(n).
8. On considère un problème d'estimation d'un signal θ bruité.
Le signal observé est x(n) =θ +b(n).
On suppose que θ suit une loi uniforme sur [−θ0, θ0] et qu'elle est décorrélée du bruit b(n) qui possède une
DSP qui vaut σ2.
• Que signifie θ et b(n) sont décorrélés ? quel est le lien avec l'indépendance?
• Calculer la moyenne et la variance de b(n). Est-il SSL?
• Déterminer le filtre de Wiener d’ordre N permettant de retrouver le signal utile θ.
• Exprimer H(z) puis θ.
2. -΋ 0 l΋
Î/
?Ÿ
134
Traitement Avancé du Signal Master RT M1 2022/2023
l l•′
Ÿt / Ø/
t
3.-t×/
4 …….. Ò
Ò
4 .
&4 Ñ&4 Ñ ……..&Ò
Buts : -Aborder les notions d'estimation ainsi que les propriétés des estimateurs (biais, variance) à travers un
exemple
- étudier les estimateurs de la moyenne, de la variance, de la corrélation et de la DSP.
- Application du filtrage de Wiener
Exercice 1: On reprend les 3 estimateurs de la valeur d'une résistance R donnés comme suit :
N N
U ( k ) I (k ) 1 N
U (k ) U (k )
Rˆ1 = k =1
N
Rˆ 2 =
N
k =1 I ( k )
Rˆ 3 = k =1
N
I
k =1
2
(k ) I (k )
k =1
Le programme suivant permet de simuler des courants et tensions bruitées et d'étudier les 3 estimateurs.
clc;clear all; close all;
randn('state', sum(100*clock));N=1000; it=1:N;R=1;
for k=2:N
U=1+0.22*randn(1,k); I=1+0.22*randn(1,k);
R1(k)=sum(U.*I)/sum(I.^2); R2(k)=sum(U./I)/k; R3(k)=sum(U)/sum(I);
end
figure; plot(it(1:100),R1(1:100),'-.r',it(1:100),R2(1:100),':g',it(1:100),R3(1:100),'-b');
legend('R1','R2','R3');
figure; plot(it,R1,'-.r',it,R2,':g',it,R3,'-b');legend('R1','R2','R3');
for k=2:N
biais_R1(k)=mean(R1(1:k))-R;biais_R2(k)=mean(R2(1:k))-R;biais_R3(k)=mean(R3(1:k))-R;
end
figure ;subplot(311); plot(biais_R1); grid; legend('Biais de R1');
subplot(312); plot(biais_R2); grid; legend('Biais de R2');
subplot(313); plot(biais_R3);grid; legend('Biais de R3');
for k=2:N
var_R1(k)=var(R1(1:k));var_R2(k)=var(R2(1:k));var_R3(k)=var(R3(1:k));
135
Traitement Avancé du Signal Master RT M1 2022/2023
end
figure ;subplot(311); plot(var_R1); legend('Var R1');
subplot(312); plot(var_R2); legend('Var R2');subplot(313); plot(var_R3);legend('Var R3');
Exercice 2: Utiliser le programme ci-dessus pour construire les estimateurs de la moyenne et de la variance
N −1 N −1
1 1 1 N −1
mˆ = xk σˆ 1 = (x − m) 2
2
σˆ 2 = ( x k − mˆ ) 2
2
k
N k =0 N k =0 N − 1 k =0
1. Prendre m=1 et σ2=2. Etudier le biais, la variance et la consistance de l'estimateur de la moyenne
2. Quelle est la différence entre les deux estimateurs de la variance ?
[Link] les 2 estimateurs de la variance sur les 100ères valeurs. Quel est le biais de chacun d’eux ?
4. Que remarque-t-on lorsque N augmente ? Les 3 estimateurs sont-ils efficaces?
Exercice 3 :Les deux estimateurs classiques de la fonction d'autocorrélation sont donnés par :
1 N −k 1 N −k
Rˆ kk (k ) = Rˆ kk (−k ) = x(n) x(n + k ) et Rˆ kk (k ) = Rˆ kk (−k ) = x(n) x(n + k)
N i =1 N − k i =1
Exercice 4 : Il existe différentes méthodes d'estimation de la densité spectrale de puissance. Les méthodes les
-
plus couramment utilisées sont :
Périodogramme simple: Sˆ xx ( f ) =
1
XT ( f )
2
Corrélogramme: TF de la corrélation xx
Sˆ ( f ) = TF Rˆ xx (k ) [{ }]
N
- Périodogramme moyenné et Périodogramme moyenné avec recouvrement (overlapping)
136
Traitement Avancé du Signal Master RT M1 2022/2023
Exercice 5: L’extraction d’ECG du fœtus se fait à partir de plusieurs électrodes placées en différents endroits
du ventre de la mère ; mais les enregistrements obtenus sont des mélanges de l’ECG du fœtus noté (FECG)
et celui de la mère noté (MECG) auquel se rajoute le bruit thermique des électrodes et des équipements
électroniques d'où l'emploi du filtre de Wiener.
137
Traitement Avancé du Signal Master RT M1 2022/2023
Buts : -Aborder les notions d'estimation ainsi que les propriétés des estimateurs (biais, variance) à travers un
exemple
- étudier les estimateurs de la moyenne, de la variance, de la corrélation et de la DSP.
- Application du filtrage de Wiener
Exercice 1: On reprend les 3 estimateurs de la valeur d'une résistance R donnés comme suit :
N N
U ( k ) I (k ) 1 N
U (k ) U (k )
Rˆ1 = k =1
N
Rˆ 2 =
N
k =1 I ( k )
Rˆ 3 = k =1
N
I
k =1
2
(k ) I (k )
k =1
Le programme suivant permet de simuler des courants et tensions bruitées et d'étudier les 3 estimateurs.
import numpy as np; import [Link] as plt;
N=1000;R1=[];R2=[];R3=[];R=1;
for k in range(2,N):
U=1+0.22*[Link](k);I=1+0.22*[Link](k);
R1=[Link](R1,sum(U*I)/sum(I**2));
R2=[Link](R2,sum(U/I)/k);
R3=[Link](R3,sum(U)/sum(I));
[Link](1);
[Link](311); [Link](R1); [Link](True); [Link]('R1');
[Link](312); [Link](R2); [Link](True);[Link]('R2');
[Link](313); [Link](R3); [Link](True);[Link]('R3');
bR1=[];bR2=[];bR3=[];VarR1=[];VarR2=[];VarR3=[];
for k in range(2,len(R1)):
bR1=[Link](bR1, [Link](R1[:k])-R);
bR2=[Link](bR2, [Link](R2[:k])-R);
bR3=[Link](bR3, [Link](R3[:k])-R);
VarR1=[Link](VarR1, [Link](R1[:k]));
VarR2=[Link](VarR2, [Link](R2[:k]));
FGE, USTHB [ assiakourgli@[Link] / [Link]
138
Traitement Avancé du Signal Master RT M1 2022/2023
VarR3=[Link](VarR3, [Link](R3[:k]));
[Link](2);
[Link](321); [Link](bR1); [Link](True); [Link]('Biais R1');
[Link](323); [Link](bR2); [Link](True);[Link]('Biais R2');
[Link](325); [Link](bR3); [Link](True);[Link]('Biais R3');
[Link](322); [Link](VarR1); [Link](True); [Link]('Variance de R1');
[Link](324); [Link](VarR2); [Link](True);[Link]('Variance de R2');
[Link](326); [Link](VarR3); [Link](True);[Link]('Variance de R3');
Exercice 3 : Il existe différentes méthodes d'estimation de la densité spectrale de puissance. Les méthodes les
plus couramment utilisées sont :
- Périodogramme simple: Sˆ xx ( f ) =
1
N
XT ( f )
2
[{
Corrélogramme: TF de la corrélation Sˆ xx ( f ) = TF Rˆ xx (k ) }]
- Périodogramme moyenné et Périodogramme moyenné avec recouvrement (overlapping)
139
Traitement Avancé du Signal Master RT M1 2022/2023
Exercice 4 : L’extraction d’ECG du fœtus se fait à partir de plusieurs électrodes placées en différents endroits
du ventre de la mère ; mais les enregistrements obtenus sont des mélanges de l’ECG du fœtus noté (FECG)
et celui de la mère noté (MECG) auquel se rajoute le bruit thermique des électrodes et des équipements
électroniques d'où l'emploi du filtre de Wiener.
140
Traitement Avancé du Signal Master RT M1 2022/2023
141
Traitement Avancé du Signal Master RT M1 2022/2023
[Link]
[Link]
[Link]
142
Traitement Avancé du Signal Master RT M1 2022/2023
Cours_TS_3A
[Link]. signal pour bien expliquer Somme XTX le produit matriciel Permettant de Rxx(k) normale
143