Signaux aléatoires

Processus ARMA

Simon Leglaive

CentraleSupélec

Dans les épisodes précédents…

• Un processus aléatoire est une collection de variables aléatoires indexées par le temps. C’est une fonction de deux variables : l’aléa et le temps.
• Un processus aléatoire est complètement décrit par sa loi temporelle, et partiellement décrit par ses propriétés au second ordre.
• Un processus peut être stationnaire (au sens strict ou au sens large), ce qui simplifie grandement sa description statistique et son traitement.
• La propriété d’ergodicité ajoutée à celle de stationnarité permet de remplacer les moyennes statistiques par des moyennes temporelles.
• La notion de mesure / densité spectrale de puissance permet de caractériser un processus SSL dans le domaine spectral.
• On peut étendre la modélisation des systèmes linéaires invariants dans le temps au cas d’une entrée sous forme de signal aléatoire SSL.

Au programme

• Les processus autoregressifs à moyenne ajustée (ou moyenne mobile) (ARMA, de l’anglais autoregressive moving average) sont définis par le filtrage d’un bruit blanc gaussien, avec un filtre de fonction de transfert rationnelle.

• Ils constituent une classe importante de processus SSL, très répandue pour modéliser les séries temporelles dans des domaines variés (finance, énergie, audio, etc.)

• Nous allons étudier les processus ARMA sous plusieurs points de vue, en tant que modèles paramétriques de :

  • signaux temporels ;
  • fonctions d’autocovariance ;
  • densités spectrales de puissance.

Processus \(MA(q)\)

Définition

On dit que le processus \(\{X(t)\}_t\) est à moyenne ajustée d’ordre \(q\), ou \(MA(q)\), s’il est défini par : \[ \begin{aligned} X(t) &= W(t) + b_1 W(t-1) + b_2 W(t-2) + \cdots + b_q W(t-q) \end{aligned} \] avec \(\{W(t)\}_t \sim BB(0, \sigma^2)\).

Le terme « Moyenne ajustée » est une traduction malheureuse de « moving average », permettant de garder l’acronyme MA, mais qui signifie en réalité « moyenne mobile »

• Un processus MA résulte donc du filtrage d’un bruit blanc par un filtre de réponse impulsionnelle finie (RIF) d’ordre \(q\) :

  \[ X(t) = [h \star W](t) = \sum_{k = 0}^q b_k W(t-k) , \]

  avec \(h(k) = b_k\) pour \(k \in \{0, \dots, q\}\) et \(b_0 = 1\).

• La fonction de transfert du filtre ne contient que des zéros : \[\begin{aligned} B(z) & = 1 + b_1 z^{-1} + b_2 z^{-2} + \cdots + b_q z^{-q} \\ & = (1 - z_1 z^{-1})(1 - z_2 z^{-1}) \cdots (1 - z_q z^{-1}). \end{aligned}\]

Caractéristiques au second ordre

En utilisant les résultats sur le filtrage des processus SSL, on montre que (voir exercice TD) :

• La moyenne est nulle : \(m_X = 0\).

• La fonction d’autocovariance s’écrit

  \[ R_{XX}(k) = \begin{cases} \sigma^2 \displaystyle\sum\limits_{j=0}^{q - |k|} b_j b_{j + |k|} \,\,& \text{si}\,\, |k| \le q \\[.5cm] 0 &\text{ sinon} \end{cases}. \]

• Un processus MA permet donc de modéliser un signal temporel comme le résultat du filtrage RIF d’un bruit blanc.

• Si on réfléchit en terme de statistiques temporelles au second ordre, un processus MA permet de modéliser un signal centré dont la mémoire est finie, c’est-à-dire que la fonction de covariance \(R_{XX}(k)\) est nulle pour tout \(k\) supérieur à un certain ordre, l’ordre du modèle MA.

• On a même la propriété suivante :

  Soit \(R_{XX}(k)\) une fonction d’autocovariance telle que \(R_{XX}(k) = 0\) pour \(k > q\). Alors, il existe un bruit blanc \(\{W(t)\}_t\) et un polynôme \(B(z) = \sum_{k = 0}^q b_k z^{-k}\) de degré inférieur ou égal à \(q\), tels que \(R_{XX}(k)\) soit la fonction d’autocovariance d’un processus \(MA(q)\).

  C’est assez fort, et très utile pour modéliser tout signal aléatoire dont la fonction d’autocovariance serait à support compact.

Caractéristiques spectrales

• La fonction d’autocovariance étant à support compact, celle-ci est sommable et le processus possède une DSP.

• Toujours en utilisant les résultats sur le filtrage des processus SSL, on montre que la DSP s’écrit (voir exercice TD) :

  \[ S_{XX}(\nu) = \sigma^2 \underbrace{ \left\vert\sum_{k = 0}^q b_k e^{-\imath 2 \pi \nu k} \right\vert^2 }_{|\hat{h}(\nu)|^2}, \]

  où

  • \(\sigma^2\) est la DSP du bruit blanc ;
  • \(\hat{h}(\nu)\) est la TFTD de la réponse impulsionnelle du filtre, aussi appelée réponse en fréquence.

• La DSP présente \(q\) anti-résonances au voisinage des fréquences \(\nu_k\), \(k \in \{1,...,q\}\), correspondant aux arguments des racines \(z_k\) du polynôme \[B(z) = \sum_{k = 0}^q b_k z^{-k} = b_0\prod_{k = 1}^q (1 - z_k z^{-1}),\] c’est à dire : \(\nu_k = \displaystyle \frac{1}{2\pi} \arg(z_k)\).

• Autrement dit, à tout zéro de la fonction de transfert est associée une anti-résonance.

• On peut ainsi se représenter le filtrage RIF d’un bruit blanc dans le domaine spectral de la façon suivante :

  On part de la DSP plate du bruit blanc \(S_{WW}(\nu) = \sigma^2\), qu’on vient « modeler » en ajoutant des anti-résonances, autant que l’ordre du modèle.

Illustration

Trajectoire, fonction d’autocovariance et DSP d’un processus \(MA(1)\) avec \(\sigma^2 = 1\) et \(b_1 = -0.9\).

Ecoutons l’entrée et la sortie du filtre. A quoi vous attendez-vous ?

\(W(t)\) :

\(X(t)\) :

Processus \(AR(p)\)

Définition

• Soit l’équation récurrente suivante : \[ X(t) + a_1 X(t-1) + \cdots + a_p X(t-p) = W(t), \tag{1}\] où \(\{W(t)\}_t \sim BB(0, \sigma^2)\).

• On peut montrer que cette équation admet une unique solution SSL si et seulement si le polynôme \[ A(z) = 1 + a_1 z^{-1} + \cdots + a_p z^{-p} \]

  n’a pas de racines sur le cercle unité, c’est-à-dire \(A(z) \neq 0\) pour \(|z| = 1\).

• Cette solution SSL est appelée processus autorégressif d’ordre \(p\), ou \(AR(p)\).

Expression de la solution

• Un processus \(AR(p)\), unique solution SSL de l’Equation 1, admet pour expression :

  \[ X(t) = [h \star W](t) = \sum_{k = -\infty}^{+\infty} h(k) W(t-k), \]

  où les coefficients \(\{h(k)\}_{k\in \mathbb{Z}}\) sont les coefficients du développement en \(z\) de la fonction de transfert \(H(z) = 1/A(z)\) : \[ H(z) = \frac{1}{A(z)} = \sum_{k = -\infty}^{+\infty} h(k) z^{-k},\]

  avec \(h(0) = 1\) et \(\sum_{k = -\infty}^{+\infty} |h(k)| < \infty\) (équivalent à \(A(z) \neq 0\) pour \(|z| = 1\)).

• Un processus AR résulte donc du filtrage d’un bruit blanc par un filtre stable de réponse impulsionnelle infinie \(\{h(k)\}_{k\in \mathbb{Z}}\).

Pour être plus précis, la fonction de transfert \(H(z) = 1/A(z)\) est analytique dans un anneau contenant le cercle unité, car \(A(z) \neq 0\) pour \(|z| = 1\). Elle admet donc un développement en série de Laurent dans cet anneau, et les coefficients \(h(k)\) sont précisément les coefficients de cette série.

Expression de la solution causale et stable

• En particulier, si les racines du polynôme \(A(z)\) sont toutes à l’intérieur du cercle unité, càd \(A(z) \neq 0\) pour \(|z| \ge 1\), le filtre est causal et stable. On a \(h(k) = 0\) pour \(k < 0\) et : \[ H(z) = \frac{1}{A(z)} = \sum_{k =0}^{+\infty} h(k) z^{-k},\]

• Dans ce cas, le processus \(AR(p)\) s’exprime causalement en fonction de \(W(t)\) :

  \[ X(t) = \sum_{k = 0}^{+\infty} h(k) W(t-k).\]

• On remarque alors qu’un processus AR, dans le cas où \(A(z) \neq 0\) pour \(|z| \ge 1\), admet une représentation sous forme \(MA(\infty)\).

• Nous nous intéresserons dans la suite uniquement à ce type de processus AR.

Identification des coefficients de la réponse impulsionnelle

• Comment identifier les coefficients \(\{h(k)\}_{k\in \mathbb{Z}_+}\) de la réponse impulsionnelle du filtre?

• On rappelle que

  • \(A(z) = 1 + a_1 z^{-1} + \cdots + a_p z^{-p}\) ;
  • \(H(z) = 1/A(z) = \sum_{k =0}^{+\infty} h(k) z^{-k}\).

• On remarque : \[ \begin{aligned} 1 &= (1 + a_1 z^{-1} + a_2 z^{-2} + \cdots + a_p z^{-p}) (h(0) + h(1) z^{-1} + h(2) z^{-2} + \cdots) \\ &= h(0) + z^{-1}\big(h(1) + a_1 h(0)\big) + z^{-2}\big(h(2) + a_1 h(1) + a_2 h(0)\big) + \cdots \end{aligned} \]

• Cela nous permet par identification de calculer récursivement la suite \(\{h(k)\}_{k\in \mathbb{Z}_+}\) à partir des coefficients de l’équation récurrente du modèle AR :

  \[ h(0) = 1, \qquad h(1) = - a_1 h(0) = -a_1, \qquad h(2) = -a_1 h(1) - a_2 h(0) = a_1^2 - a_2, \qquad \cdots\]

Caractéristiques au second ordre

A partir des coefficients de la réponse impulsionnelle infinie du modèle AR, il suffit d’utiliser les résultats sur le filtrage des processus SSL pour montrer que (voir exercice TD) :

• La moyenne est nulle et la fonction d’autocovariance s’écrit

  \[ R_{XX}(k) = \sigma^2 \displaystyle\sum\limits_{j=0}^{+\infty} h(j) h(j + |k|). \]

• Nous montrerons aussi dans le cadre d’un travail préparatoire de TP que la fonction d’autocovariance vérifie la récurrence suivante pour \(k \ge 1\) : \[ \begin{aligned} R_{XX}(k) + \sum_{j=1}^p a_j R_{XX}(k-j) = 0 \end{aligned} \]

  et \(R_{XX}(0) + \sum\limits_{j=1}^p a_j R_{XX}(j) = \sigma^2\).

Caractéristiques spectrales

• Toujours en utilisant les résultats sur le filtrage des processus SSL, on montre que la DSP s’écrit (voir exercice TD) :

  \[ S_{XX}(\nu) = \frac{\sigma^2}{\left\vert 1 + \sum\limits_{k=1}^p a_k e^{-\imath 2 \pi \nu k} \right\vert^2} \]

• La DSP présente \(p\) résonances au voisinage des fréquences \(\nu_k\), \(k \in \{1,...,q\}\), correspondant aux arguments des racines \(z_k\) du polynôme \[A(z) = 1 + \sum_{k = 1}^p a_k z^{-k} = \prod_{k = 1}^p (1 - z_k z^{-1}),\] c’est à dire : \(\nu_k = \displaystyle \frac{1}{2\pi} \arg(z_k)\).

• Autrement dit, à tout pôle de la fonction de transfert est associée une résonance.

• On peut ainsi se représenter le filtrage RII d’un bruit blanc dans le domaine spectral de la façon suivante :

  On part de la DSP plate du bruit blanc \(S_{WW}(\nu) = \sigma^2\), qu’on vient « modeler » en ajoutant des résonances, autant que l’ordre du modèle.

Processus \(AR(1)\) (voir exercice TD)

• On considère le processus \(AR(1)\) défini par : \[ X(t) + a X(t-1) = W(t), \] où \(|a|<1\) et où \(\{W(t)\}_t \sim BB(0, \sigma^2)\).

• On en déduit : \(X(t) = \sum\limits_{k=0}^{+\infty} (-1)^k a^k W(t-k).\)

• On obtient les expressions de la fonction d’autocovariance et de la DSP suivantes : \[ R_{XX}(k) = \sigma^2 \frac{(-1)^k a^{|k|}}{1-a^2}, \qquad S_{XX}(\nu) = \frac{\sigma^2}{\left\vert 1 + a e^{-\imath 2 \pi \nu} \right\vert^2}. \]

Le modèle ne fait intervenir que deux instants consécutifs, donc intuitivement on aurait pu s’attendre à ce que la fonction d’autocovariance s’annule pour \(|k| > 1\). En réalité, celle-ci tend vers \(0\) quand \(|k|\) augmente car \(|a| <1\). Cela vient du fait que la corrélation entre \(X(t)\) et \(X(t-2)\) est non nulle car \(X(t)\) dépend de \(X(t-1)\) et \(X(t-1)\) dépend de \(X(t-2)\), et ainsi de suite.

Illustration

Trajectoire, fonction d’autocovariance et DSP d’un processus \(AR(1)\) avec \(\sigma^2 = 1\) et \(a = [0.7, -0.5, -0.9]\).

• Lorsque \(a < 0\), le processus est positivement corrélé, dans le sens où tous les coefficients de la fonction d’autocovariance sont positifs.
• Les exemples de trajectoire précédents montrent que :

  • Des valeurs de \(a\) positives conduisent à des trajectoires très oscillantes, où une valeur positive à tendance à être suivie par une valeur négative.

    Cela est cohérent avec l’allure de la DSP, concentrée dans les hautes fréquences, et avec la fonction d’autocovariance qui présente des anti-corrélations fortes.

  • Des valeurs de \(a\) proches de \(1\) correspondent à des trajectoires très lisses, où les valeurs se suivent.

    Cela est également cohérent avec l’allure de la DSP, concentrée dans les basses fréquences, et avec la fonction d’autocovariance, qui décroît lentement avec le décalage \(k\).

Sinusoïde bruitée vs processus \(AR(2)\)

• Sinusoïde bruitée de fréquence \(\nu_0\) et processus \(AR(2)\) dont la résonance se situe à la fréquence \(\nu_0\).
• La sinusoı̈de bruitée présente de fortes irrégularités d’amplitude, mais on voit que la “périodicité” reste constante et on pourrait estimer cette période via une simple transformée de Fourier.
• Le processus \(AR(2)\) ne présente pas de brusques variations d’amplitude, la courbe est lisse. Par contre l’amplitude varie dans de larges limites et la phase glisse continuellement.
• Quand l’observation présente un tel comportement le modèle AR est préférable à un modèle sinusoïde plus bruit. C’est ce qui a poussé George U. Yule à utiliser en 1927 un processus AR pour modéliser la série chronologique du nombre de taches solaires, plutôt que la méthode du périodogramme de Schuster.

Résumé

Processus \(MA(q)\)

• Filtrage RIF d’un bruit blanc.
• Fonction d’autocovariance à support compact : utile pour modéliser des signaux à mémoire courte.
• DSP qui présente des anti-résonances.

Processus \(AR(p)\)

• Filtrage RII d’un bruit blanc.
• Fonction d’autocovariance infinie mais qui tend vers \(0\) avec le décalage \(k\) dans le cas d’un filtre stable : utile pour modéliser des signaux à mémoire longue, présentant des trajectoires assez lisses.
• DSP qui présente des résonances.

Processus \(ARMA(p, q)\)

• Un processus ARMA est la mise en cascade d’un \(MA(q)\) et d’un \(AR(p)\).
• Il s’agit donc du filtrage d’un bruit blanc par un filtre dont la fonction de transfert est une fraction rationnelle, avec des pôles et des zéros.
• La DSP du processus présente ainsi à la fois des résonances et des anti-résonances, ce qui est très flexible en terme de modélisation, et explique la grande universalité du modèle ARMA.

Application en restauration audio

Restauration audio

• Les « vieux enregistrements » audio présentent souvent des défauts dus à des techniques d’enregistrement obsolètes ou à la détérioration physique du support (disque, bande, …).
• La restauration audio consiste à supprimer ces défauts pour améliorer la qualité sonore.
• Dans la suite on s’intéressera à la restauration de craquements : défauts d’assez grande amplitude et de durée très brève (inférieure à la milliseconde).

Signal original :


Signal restauré :

Modèle de signal bruité

• Modèle de signal bruité : \[ \hspace{2cm} Y(t) = \begin{cases} X(t) &\qquad \text{si absence de craquement ;} \\ X(t) + A \delta(t-t_0)& \qquad \text{si présence d'un craquement d'amplitude $A$ en $t_0$.} \end{cases} \]

• Une étude d’un signal bruité montre que les clics sont rares (~1% des échantillons) et leur durée est courte (1 à 10 échantillons).

Modèle AR de signal utile

• On suppose que \(\{X(t)\}_t\) suit un modèle \(AR(p)\) : \[ \hspace{-5cm} X(t) + \sum_{k=1}^p a_k X(t-k) = W(t), \] avec \(\{W(t)\}_t \sim BB(0, \sigma^2)\).

• Cette modélisation stationnaire n’est valable que sur une courte fenêtre temporelle (quelques dizaines de ms).

• Interprétation en terme de prédiction linéaire : \[ X(t) = \hat{X}(t) + W(t), \qquad \hat{X}(t) = - \sum_{k=1}^p a_p X(t-k), \] où \(\hat{X}(t)\) est la prédiction de \(X(t)\) à partir de ses \(p\) échantillons passés, et \(W(t)\) est l’erreur ou résidu de prédiction.

Détection de craquement

• Dans le domaine de la transformée en \(z\), le modèle AR s’écrit :

  \[ X(z) = \frac{1}{A(z)} W(z), \qquad A(z) = 1 + \sum_{k=1}^p a_k z^{-k}. \]

• Supposons les paramètres du modèle AR connus¹ et appliquons le filtre de fonction de transfert \(A(z)\) au signal bruité \(Y(t)\).

  • En cas d’absence de craquement, on a \(A(z)Y(z) = A(z)X(z) = W(z)\) et donc nous obtenons un bruit blanc en sorti, d’où le nom de filtre de blanchiment.
  • En présence de craquement, nous obtenons un bruit blanc auquel s’ajoute une impulsion correspondant au craquement.

• On pourra donc détecter le craquement en seuillant le résultat de ce filtrage.

1. Vous verrez dans un travail préparatoire de TP comment les estimer.

Reconstruction du signal

• On travaille maintenant avec des réalisations des signaux aléatoires précédents, d’où le changement de notation en lettres minuscules.

• On suppose que l’ensemble des clics ont été détectés.

• Soit \(\ell\) l’indice temporel d’un échantillon en présence de craquement. Le signal reconstruit à cet instant est donné par :

  \[\hat{x}(\ell) = - \sum_{k=1}^p a_p x(\ell-k).\]

• En pratique, le modèle de bruit impulsionnel ne correspond pas bien à ce que nous avons observé sur la figure précédente. Le craquement n’est pas parfaitement localisé en temps, il s’étale sur plusieurs échantillons.

• On considère donc une plage de \(m\) échantillons à reconstruire.

• Soit \(\{x_{\ell}, x_{\ell+1}, ..., x_{\ell+m-1}\}\) l’ensemble des échantillons à reconstruire.

• Par une approche de type moindres carrés, on minimise par rapport à cet ensemble de coefficients l’erreur de prédiction :

  \[ \sum_{t=p+1}^T \left( x(t) - \hat{x}(t) \right)^2 = \sum_{t=p+1}^T \left( x(t) + \sum_{k=1}^p a_p x(t-k) \right)^2, \]

  où \(T\) est la longueur de la trame de signal traitée.

• On obtient une solution en forme close qui s’exprime à partir des échantillons non corrompus et des paramètres du modèle AR.

Algorithme final

• L’algorithme final prend en en entrée un signal dégradé et le divise en trames de quelques dizaines de millisecondes avec recouvrement.
• Sur chaque trame, les paramètres du modèle AR sont estimés et utilisés pour la détection puis pour la reconstruction.
• Le signal est ensuite reconstruit grâce à une procédure d’addition-recouvrement.
• Le processus est itéré si besoin (on suppose un unique craquement par trame).

Pour plus d’informations, voir : L. Oudre, Est-il possible de restaurer automatiquement des signaux audio corrompus par du bruit impulsionnel ?, Colloque GRETSI 2013.