Estimation spectrale

Yves Goussard

GBM6103A

2 octobre 2013

Yves Goussard (GBM6103A) Estimation spectrale 2 octobre 2013 1 / 28

Traitement de signaux biomédicaux

Mesure

Extraction
Tracé Stockage Analyse Transformation
d’informations

Échantillonnage Représentation Modélisation

Filtrage
Quantification Fréquence Estimation
·

Yves Goussard (GBM6103A) Estimation spectrale 2 octobre 2013 2 / 28

Estimation spectrale

Objectif : représentation d’un signal dans le domaine spectral

Cadre de référence
Cadre déterministe : les outils ont été présentés
Cadre aléatoire
une seule réalisation de disponible
problème d’estimation

Fil conducteur
Présentation des outils d’estimation du spectre d’un signal aléatoire
Comparaison avec le point de vue déterministe

Yves Goussard (GBM6103A) Estimation spectrale 2 octobre 2013 3 / 28

Plan

1 Cadre de travail

2 Rappels : densité spectrale de puissance

3 Analyse spectrale classique (non paramétrique)

4 Analyse spectrale moderne (paramétrique)

5 Exemple : analyse spectrale du signal Doppler

Yves Goussard (GBM6103A) Estimation spectrale 2 octobre 2013 4 / 28

Cadre de travail

Hypothèses générales
Signaux à valeurs réelles
Signaux à temps discret – quantités réduites (n, νr )
N échantillons disponibles : xn ; 0 ≤ n ≤ N − 1

Hypothèses possibles sur le signal x

1 xn déterministe à temps discret
2 xn partie de la réalisation d’un signal aléatoire Xn ; n ∈ Z stationnaire
du 2e ordre

Effet des deux types d’hypothèses sur l’évaluation des caractéristiques

fréquentielles de x ?

Yves Goussard (GBM6103A) Estimation spectrale 2 octobre 2013 5 / 28

Rappels : densité spectrale (1)

Cadre déterministe
F
xn ←→ X (νr )
Densité spectrale d’énergie : |X (νr )|2 = F(γx (n))
γx (n) : fonction d’autocorrélation du signal déterministe
X
γx (n) = xp xp−n
p∈Z

Cadre aléatoire
Z
On voudrait écrire Xn = u(νr ) e 2iπnνr d νr
Difficulté : intégrale stochastique

Yves Goussard (GBM6103A) Estimation spectrale 2 octobre 2013 6 / 28

Rappels : densité spectrale (2)

Cadre aléatoire : résultats admis

L’intégrale stochastique a un sens
Le calcul direct de u(νr ) dépasse le cadre du cours
On a E |u(νr )|2 = Γx (νr ) densité spectrale de puissance de Xn
 

Γx (νr ) = F(γx (n)) ; γx (n) : fonction d’autocorrélation de Xn

γx (n) = E [Xp Xp−n ] indépendant de p (stationnarité de Xn )

En résumé...
Dans les deux cas (déterministe et aléatoire)

Spectre : transformée de Fourier la fonction d’autocorrélation

La nature de la fonction d’autocorrélation diffère selon le cadre adopté

Yves Goussard (GBM6103A) Estimation spectrale 2 octobre 2013 7 / 28

Analyse spectrale non paramétrique

Approche
Xn : signal aléatoire stationnaire du 2e ordre
Objectif : estimation de la fonction d’autocorrélation à partir des N
échantillons observés

estimation F
xn ←→ bx (n) ←→ b
γ Γx (νr )

Difficultés
Impossible d’appliquer la démarche générale pour l’estimation
Nécessité de recourir à des estimateurs empiriques

Yves Goussard (GBM6103A) Estimation spectrale 2 octobre 2013 8 / 28

Estimateurs empiriques de la fonction d’autocorrélation (1)

Forme des estimateurs

γx (n) = E [Xp Xp−n ]
Propriété d’ergodisme :
N−1
u 1 X
γ
bx (n) = xp xp−|n| estimateur non biaisé
N − |n|
p=|n|
N−1
1 X
bxb (n) =
γ xp xp−|n| estimateur biaisé
N
p=|n|

Yves Goussard (GBM6103A) Estimation spectrale 2 octobre 2013 9 / 28

Estimateurs empiriques de la fonction d’autocorrélation (2)

Biais des estimateurs

Si |n| ≤ N − 1
γxu (n)] = γx (n) =⇒ estimateur non biaisé
E [b
 b   
E γbx (n) = 1 − |n| N γx (n) =⇒ estimateur biaisé

Variance des estimateurs

Complexe à calculer, même sous hypothèse gaussienne

Remarques
Estimateur biaisé souvent préférable (meilleur compromis biais -
variance)
On peut écrire
γxu (n)] = RectN (n) γx (n)
 b 
E [b E γ
bx (n) = TriN (n) γx (n)

Yves Goussard (GBM6103A) Estimation spectrale 2 octobre 2013 10 / 28

Corrélogramme

Définition
N−1
X
Γx (νr ) = F(b
b γx (n)) = bx (n) e −2iπnνr
γ
n=−N+1

Biais du corrélogramme
h i
Linéarité de F et E [·] : E b
Γx (νr ) = F(E [b
γx (n)])
Cas non biaisé : h i
γxu (n)] = RectN (n) γx (n) =⇒ E b
E [b Γux (νr ) = RECTN (νr ) ∗ Γx (νr )
Cas biaisé
 b : h i
E γbx (n) = TriN (n) γx (n) =⇒ E bΓbx (νr ) = TRIN (νr ) ∗ Γx (νr )
sin(2N − 1)πνr 1 sin Nπνr 2
 
RECTN (νr ) = TRIN (νr ) =
sin πνr N sin πνr

Yves Goussard (GBM6103A) Estimation spectrale 2 octobre 2013 11 / 28

Biais du corrélogramme

Allure de RECTN (νr ) et TRIN (νr )

TRANSFORMÉE DE FOURIER DES FENÊTRES, N = 7
1
Fenêtre rectangulaire
Fenêtre triangulaire

0.8

0.6

0.4
Amplitude

0.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
Fréquence réduite

Yves Goussard (GBM6103A) Estimation spectrale 2 octobre 2013 12 / 28

Caractéristiques du corrélogramme

Principales propriétés
Cas non biaisé : bΓux (νr ) peut prendre des valeurs négatives !
Γb (νr ) est toujours positif
Cas biaisé : b x
Γux (νr ) > résolution de b
Résolution de b Γbx (νr )
Dans tous les cas, le biais du corrélogramme tend vers 0 lorsque N
tend vers l’infini

Limitations du corrélogramme
Possibilité d’estimées négatives
Deux opérations en cascade −→ étude difficile

Yves Goussard (GBM6103A) Estimation spectrale 2 octobre 2013 13 / 28

Périodogramme (1)

Objectifs généraux
Correction des défauts du corrélogramme
Lien plus direct avec le cadre déterministe

Définition
N−1 2
1 X
Px (νr ) = xn e −2iπnνr
N
n=0

Lien avec le corrélogramme

Γbx (νr )
On a Px (νr ) = b
Le périodogramme a les mêmes caractéristiques que le corrélogramme
sous hypothèse biaisée

Yves Goussard (GBM6103A) Estimation spectrale 2 octobre 2013 14 / 28

Périodogramme (2)
Caractéristiques statistiques
Biais : E [Px (νr )] = TRIN (νr ) ∗ Γx (νr )
Variance : calcul difficile. Sous hypothèses gaussiennes :
Var [Px (νr )] ≥ |E [Px (νr )]|2 ⇒ quand N → ∞, Var [Px (νr )] ≥ Γx (νr )

Estimateur non convergent

Illustration du comportement du périodogramme

Voir la démonstration demo_periodo.m

Avantages et limitations
Simple à employer
Caractéristiques statistiques déficientes
Possibilités d’extension (autres fenêtres sur la fonction de corrélation)
qui ne rendent pas l’estimateur convergent
Yves Goussard (GBM6103A) Estimation spectrale 2 octobre 2013 15 / 28
Périodogramme moyenné (1)

Démarche
Séparation des N échantillons en L blocs de K échantillons
Hypothèse de décorrélation entre les blocs
Calcul d’un périodogramme sur chaque bloc
Moyenne des L périodogrammes obtenus

Biais du périodogramme moyenné

Sur chaque bloc : E [Px (νr )] = TRIK (νr ) ∗ Γx (νr )
En prenant la moyenne : E [Pxm (νr )] = TRIK (νr ) ∗ Γx (νr )

Perte de résolution par rapport au périodogramme simple

Yves Goussard (GBM6103A) Estimation spectrale 2 octobre 2013 16 / 28

Périodogramme moyenné (2)

Variance du périodogramme moyenné

Hypothèse de décorrélation entre blocs : Var [Pxm (νr )] = L1 Var [Px (νr )]
Si K suffisamment grand, Var [Px (νr )] bornée

Estimateur convergent

Illustration du comportement du périodogramme moyenné

Voir la démonstration demo_periodo.m

Extensions possibles
Fenêtre différente sur chaque bloc
Chevauchement des blocs (fonction spectrum de Matlab)
À N donné, réglage du compromis biais - variance

Yves Goussard (GBM6103A) Estimation spectrale 2 octobre 2013 17 / 28

Périodogramme lissé (1)

Démarche
Périodogramme simple considéré comme un signal bruité
Réduction du bruit par filtrage : Px` (νr ) = Px (νr ) ∗ W (νr )
Choix de W (νr ) et lien avec les résultats précédents ?

Hypothèses
w (n) paire, à support compact [−K + 1 , K − 1] avec K  N
Z
W (νr ) normalisée : |W (νr )|2 d νr = 1

Biais du périodogramme lissé

On montre que E Px` (νr ) ≈ W (νr ) ∗ Γx (νr )
 

faible largeur fréquentielle
Biais faible ↔ W :
grande largeur temporelle
Yves Goussard (GBM6103A) Estimation spectrale 2 octobre 2013 18 / 28
Périodogramme lissé (2)

Variance du périodogramme lissé

Étude très délicate
Sous hypothèses restrictives : Var Px` (νr ) ≈ |W (νr )|2 ∗ Γ2x (νr )
 

 W (νr ) petit et uniformément réparti
Variance faible ↔ grande largeur fréquentielle
faible largeur temporelle


Convergence assurée

À N fixé, nécessité d’un compromis biais - variance

Illustration du comportement du périodogramme lissé

Voir la démonstration demo_periodo.m

Yves Goussard (GBM6103A) Estimation spectrale 2 octobre 2013 19 / 28

Exemple de fenêtre de lissage du périodogramme

Fenêtre de Hamming
FENÊTRE DE HAMMING
1

0.8
Amplitude

0.6

0.4

0.2

0
0 10 20 30 40 50 60 70
Nombre d’échantillons

2
10

0
10
Amplitude

−2
10

−4
10
−0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5
Fréquence réduite

Yves Goussard (GBM6103A) Estimation spectrale 2 octobre 2013 20 / 28

Analyse spectrale non paramétrique
Conclusions

Variantes du périodogramme
Algorithmes efficaces (FFT)
Grande variance à petit nombre d’échantillons
En général, utiles si N > 128 ou 256

Méthodes recommandées
Méthode de Welch (périodogramme moyenné)
Périodogramme lissé
Attention au compromis biais - variance

Yves Goussard (GBM6103A) Estimation spectrale 2 octobre 2013 21 / 28

Analyse spectrale paramétrique
·

Rappel : modélisation des signaux stationnaires du 2e ordre

Si Xn suffisamment régulier :

n Xn
filtre h
Bruit générateur Signal aléatoire

n bruit blanc de variance σ2

Γx (νr ) = σ2 |H(νr )|2

·
Principe de l’analyse spectrale paramétrique
Paramétrer le filtre h
Estimer les paramètres de h à partir des échantillons observés de Xn
En déduire H(νr ), puis le spectre de Xn

Yves Goussard (GBM6103A) Estimation spectrale 2 octobre 2013 22 / 28

Analyse spectrale AR

Paramétrisation du filtre h
h : filtre tout pôles
XP
Xn = ap Xn−p + n
p=1
Représentation généralement compacte, simple à manipuler et à
estimer

Spectre AR
Une fois les paramètres estimés :
b2
σ
Γx (νr ) =
b
PP 2
1− ap e −2iπpνr
p=1 b

Yves Goussard (GBM6103A) Estimation spectrale 2 octobre 2013 23 / 28

Estimation des paramètres AR
Utilisation des coefficients de corrélation (1)

Équations de Yule-Walker
P
X P
X
Xn = ap Xn−p + n =⇒ γ(k) = ap γ(p − k) ; k ≥ 1
p=1 p=1
En faisant varier k de 1 à P
    
1 γ(1) γ(2) · · · ··· γ(P) σ2 γ(0)

 0 γ(0) γ(1) · · · ··· γ(P − 1) 
 a1  
  γ(1) 

 .. .. .. ..  ..   .. 
 . γ(1) . . .  .   . 
.. .. .. .. = ..
    
 .. .. ..  
 . . . . . .  .   . 
.. .. .. ..
    
.. ..
. .
    
 . . γ(1)  .   . 
0 γ(P − 1) ··· ··· γ(1) γ(0) aP γ(P)

Yves Goussard (GBM6103A) Estimation spectrale 2 octobre 2013 24 / 28

Analyse spectrale AR
Utilisation des coefficients de corrélation (2)

Procédure
Choix de l’ordre P
Estimation des coefficients de corrélation {γ(0), . . . , γ(P)}
σ2 , b
Résolution de l’équation de Yule-Walker → {b aP }
a1 , . . . , b
Évaluation du spectre

Avantage
Algorithmes simples

Inconvénients
Estimateur des coefficients de corrélation ?
Étapes multiples, caractéristiques globales de l’estimateur spectral ?
Limitation pratique : P ≤ N/2
Yves Goussard (GBM6103A) Estimation spectrale 2 octobre 2013 25 / 28
Analyse spectrale AR
Utilisation directe des échantillons observés (1)

Démarche
Écriture de la relation entrée-sortie sous la forme
x = Xa + 
x, X : vecteur et matrice formés à partie des échantillons observés
a : vecteur des coefficients {a1 , . . . , ap } du modèle
 : vecteur des échantillons du bruit générateur (considéré comme une
incertitude)

Estimation des coefficients

Approche générale sous hypothèses linéaire et gaussienne
−1
σ2

t
MAP : ba = X X + 2I X tx
σa

−1 t
a = X tX
MV : b X x
Yves Goussard (GBM6103A) Estimation spectrale 2 octobre 2013 26 / 28
Analyse spectrale AR
Utilisation directe des échantillons observés (2)

Caractéristiques de la méthode
Bon comportement général (biais, variance)
Estimateur MV : choisir P ≤ N/2
Estimateur MAP : on peut aller jusqu’à P = N

En pratique
Choisir l’ordre P du modèle
Choisir le type d’estimateur (rapport P/N, largeur de bande de Xn ...)
Estimer a
En déduire le spectre de Xn . σ2 s’estime par la puissance des résidus

Illustration du comportement de la méthode

Voir la démonstration demo_AnSpPar.m
Yves Goussard (GBM6103A) Estimation spectrale 2 octobre 2013 27 / 28
Exemple d’application

analyse spectrale du signal Doppler

Voir la démonstration demo_ASDoppler.m

Yves Goussard (GBM6103A) Estimation spectrale 2 octobre 2013 28 / 28