3. De la théorie à l'expérimentation

La cyclostationnarité suppose l'existence de plusieurs réalisations d'un même processus. Or, en pratique, les signaux (issus d'une carte d'acquisition par exemple) correspondent à une seule réalisation de longueur finie. On va donc être conduit tout comme dans le cas stationnaire à parler d'ergodicité et d'estimation.

Le but de ce chapitre est dans un premier temps de généraliser la notion d'ergodicité aux processus cyclostationnaires et d'expliquer comment estimer les différentes grandeurs caractéristiques à l'ordre $1$ et $2$. Dans un deuxième temps, un exemple de diagnostic sur les signaux d'engrenages illustre l'exploitation de la cyclostationnarité à l'ordre $1$.

3.1 Cycloergodisme

Lorsque l'on ne dispose que d'une seule réalisation, il existe deux solutions pour exploiter la cyclostationnarité.

• Il est possible d'utiliser le formalisme Fraction of Time Probability construit par Gardner [#!gardner:cyclostationarity!#]. Dans ce formalisme, les signaux ne sont plus des réalisations d'un processus aléatoire mais une série temporelle (ou série angulaire selon nos notations). Ce formalisme conduit notamment à redéfinir la notion de fonction de répartition. Il traite le problème partiellement puisqu'il supprime la notion de réalisation au profit d'une série temporelle.

• Il est possible de caractériser nos signaux à partir d'une seule réalisation en étendant la notion d'ergodicité aux processus cyclostationnaires. On parle alors de cycloergodicité [#!boyles:cycloergodic!#].

Dans les deux cas, la réalisation devra contenir un grand nombre de cycles. Nous allons ici définir et utiliser l'approche cycloergodique, plus classique que celle de Gardner.

Comme la notion d'ergodisme fait appel à la notion de moyenne temporelle, la notion de cycloergodisme fera appel à la notion de moyenne cyclique ou synchrone.

3.1.1 Moyenne cyclique (synchrone)

Soit $\pstoch{x}{\theta}{Z}$ un processus stochastique cyclostationnaire complexe ou réel ayant comme période cyclique $\Theta$, contenant un nombre de cycles important, et, $x\left(\theta\right)$ une réalisation $\omega$ particulière.

On découpe le signal en $K$ blocs consécutifs de longueur $\Theta$. La moyenne synchrone est alors donnée par :

$$\left<x\left(\theta\right)\right>_\Theta^K = \frac{1}{K} \s... ...^{K-1}x\left[mod\left(\theta + k \Theta,K\Theta\right)\right]$$ (3.1)

Figure: Calcul de la moyenne synchrone

Où $mod\left(a,b\right)$ est le reste de la division entière de $a$ par $b$. La fonction $mod$ permet de définir la moyenne synchrone pour tout $\theta$. Cette opération est schématisée sur la figure .

Puisque la moyenne synchrone extrait une composante périodique, qui admet une série de Fourier, elle peut être également écrite lorsque le nombre de blocs est infini :

$$\left<x\left(\theta\right)\right>_{\Theta}=\sum_{k=-\infty}... .../2}^{W/2}x\left(u\right) e^{-2 \pi j k u/\Theta} du\right\}$$ (3.2)

Où le terme entre crochet représente la transformée de Fourier pour les fréquences $\frac{1}{\Theta}$ et ses harmoniques. La composante périodique est ensuite générée par la série de Fourier. La moyenne synchrone apparaît donc comme un filtre en peigne sélectionnant uniquement la fréquence $1/\Theta$ et ses multiples.

3.1.2 Définition du cycloergodisme (au sens fort)

3.1.2.1 Cas où il existe $1$ seul cycle

Soit un processus stochastique cyclique $\pstoch{x}{\theta}{R}$ de période cyclique $\Theta$, la $\omega^{\text{\tiny ième}}$ réalisation de ce processus $x_{\omega}\left(\theta\right)$, et une fonction certaine $g$. $X\left(\Theta\right)$ est cycloergodique au sens fort à la période cyclique $\Theta$ si la moyenne synchrone (cyclique) $\left<g\left[x_{\omega}\left(\theta\right)\right]\right>_\Theta^K$ converge presque sûrement^3.1 vers un signal certain quand $K$ tend vers l'infini.

Si le processus $X\left(\theta\right)$ est de plus cyclostationnaire à la période cyclique $\Theta$, ce signal certain correspond au moment d'ordre $n$ du processus $X\left(\theta\right)$ pour $g\left(x\right)=x^n$.

En pratique, la cycloergodicité sera supposée (ou établie d'après un modèle). En effet, la vérifier nécessite de posséder différentes réalisations.

3.1.2.2 Cas où il existe plusieurs cycles

Dans le cas polycyclostationnaire, le signal n'a pas de période mais plusieurs ``périodes'' incommensurables entre elles. Chacune des périodes produit un peigne de Dirac dans le domaine fréquentiel. Puisque ces périodes sont incommensurables, tous les peignes sont distincts. Dès lors, la moyenne synchrone agit comme un filtre en peigne et, permet de sélectionner uniquement la période $\Theta$ prise comme référence.

Afin de définir la cycloergodicité dans ce contexte, nous allons définir l'extracteur des composantes poly-périodiques associées aux K périodes cycliques $\left\{\Theta_k\right\}$ :

$$\left<x\left(\theta\right)\right>_{\left\{\Theta_k\right\}}=\sum_{k=1}^{K} \left<x\left(\theta\right)\right>_{\Theta_k}$$ (3.3)

Nous pouvons alors définir le cycloergodisme pour plusieurs fréquences cycliques :

Soit un processus stochastique poly-cyclique $\pstoch{x}{\theta}{R}$ de périodes cycliques $\left\{\Theta_k\right\}$, $x_{\omega}\left(\theta\right)$ $\omega^{\text{\tiny ième}}$ réalisation de ce processus, et une fonction certaine $g$. $X\left(\Theta\right)$ est poly-cycloergodique au sens fort aux périodes cycliques $\left\{\Theta_k\right\}$ si $\left<g\left[x_{\omega}\left(\theta\right)\right]\right>_{\left\{\Theta_k\right\}}$ converge presque sûrement vers un signal certain.

L'introduction du cycloergodisme va nous permettre de développer les estimateurs basés sur une réalisation $x\left(\theta\right)$ d'un processus cycloergodique $X\left(\theta\right)$ afin de caractériser la cyclostationnarité aux ordres $1$ et $2$.

Les estimations, basées sur une seule réalisation, peuvent également être réalisées en supposant l'unicité du processus générateur pour chacun des cycles considérés. Sous cette hypothèse, il est possible d'assimiler un cycle à une réalisation c'est à dire de substituer une moyenne d'ensemble par une moyenne de cycle.

3.2 Caractérisation de la cyclostationnarité aux ordres $1$ et $2$

3.2.1 Caractérisation à l'ordre 1

Pour caractériser la cyclostationnarité à l'ordre $1$, nous allons chercher à estimer la composante (presque) périodique d'un signal (presque) cyclostationnaire.

3.2.1.1 Moyenne synchrone

Le premier outil destiné à extraire cette composante périodique est la moyenne synchrone [#!braun:extraction!#] qui nous a servi à introduire la cycloergodicité dans le paragraphe précédent. La moyenne synchrone, bien qu'en apparence simple à calculer présente quelques pièges :

1. la période doit être connue de manière exacte. Dans le cas contraire, il est possible d'utiliser des outils tels que le spectre, l'analyse cepstrale (présentée au chapitre ), ou la corrélation pour mettre en évidence cette période, et la calculer avec une bonne précision. Si elle n'est pas connue exactement, l'estimateur de la moyenne synchrone tendra vers 0 (sauf si cette période correspond par hasard à une autre période présente dans le signal).
2. la période doit correspondre à un nombre entier de points. Cette condition n'est pas nécessairement vérifiée même si l'on travaille dans le domaine angulaire. Par exemple, dans le cas d'un engrenage, si la période de l'arbre moteur correspond à un nombre entier, la période en sortie du réducteur, multipliée par le rapport de réduction ($2/3$ par exemple) ne sera pas entière.

   Dans, un tel cas où les fréquences d'entrée et de sortie sont commensurables (rapport $p/q$, $p$ et $q$ entiers), il suffit tout simplement de sur-échantillonner d'un facteur $p$ et de le sous-échantillonner d'un facteur $q$ afin d'obtenir un nombre de points entier par tour de l'arbre de sortie (les outils de conversion de fréquence d'échantillonnage sont décrits dans [#!crochiere:multirate!#]). Il est également possible d'utiliser une interpolation pour effectuer le ré-échantillonnage.

3. Le signal utilisé pour le calcul de la moyenne synchrone doit être tronqué afin d'avoir un nombre fini de périodes.

Ces conditions tendent à favoriser l'échantillonnage angulaire qui fournit un nombre de points constant par tour de l'arbre de référence (suppression de l'effet des fluctuations de vitesses). Des méthodes permettant de transformer des signaux à variable générique temporelle en signaux à variable générique angulaire ont été proposées au chapitre . En angulaire, les signaux sont ``synchronisés'' par rapport à la période cyclique.

Il est possible d'écrire la moyenne synchrone sous une forme similaire à l'équation () pour un signal de longueur $L$, multiple de la taille d'un bloc de $\Theta$ échantillons :

$$\left<x\left(n_\theta\right)\right>_\Theta^K=\sum_{k=0}^{\T... ...1}{L} \sum_{m=0}^{L-1} x\left(m\right) e^{-2\pi j k m/\Theta}$$ (3.4)

A la différence du cas continu, le nombre de points fini, conduira à employer un filtre en peigne avec des ``dents élargies''. On ne sélectionnera donc plus aussi précisément une fréquence donnée. La réponse fréquentielle du filtre en peigne équivalent est donnée par [#!braun:extraction2!#]:

$$H\left(f\right)=\frac{1}{K} \frac{\sin\left(\pi f \Theta K\... ...eft(\pi f \Theta\right)} e^{-j \pi f \Theta \left(K-1\right)}$$ (3.5)

Cette réponse est une fonction de Dirichlet qui a un maximum à la fréquence $1/\Theta$ ainsi qu'à ses harmoniques (voir figure ). La largeur des dents est inversement proportionnelle au nombre de cycles $K$. Dès lors, plus le nombre de cycles est grand, plus le filtre sera sélectif. Quand le nombre de cycles tend vers l'infini, ce filtre tend vers le filtre en peigne de l'équation ().

Figure: Module du filtre équivalent à la moyenne synchrone

3.2.1.2 Seuil pour la moyenne synchrone :

La moyenne synchrone étant calculée sur un nombre fini de points, il n'est pas toujours facile de pouvoir différencier une composante ``synchrone'' qui ressortira grâce au moyennage constructif, d'une composante asynchrone subissant un moyennage destructif. Nous avons alors proposé dans [#!bonnardot:aide!#] d'utiliser un critère afin de détecter les valeurs résultant d'un moyennage constructif. Ce dernier est intéressant lorsque la composante périodique a une amplitude faible vis à vis du signal ou que le nombre de cycles utilisés pour le moyennage est faible. Lorsque le nombre de cycles est important (quelques centaines), le bruit est réduit fortement, alors un examen visuel du signal est souvent suffisant.

Ce critère est basé sur la définition d'un seuil. Toute valeur au dessus de ce seuil sera considérée comme une composante synchrone (avec un risque de se tromper donné).

La moyenne synchrone dépend de l'angle, donc le seuil en dépendra également. Pour définir un tel seuil, nous modélisons le signal comme la somme d'une composante déterministe $m\left(\theta\right)$ de période $\Theta$ et d'une composante résiduelle $r\left(\theta\right)$ incluant le bruit ainsi que les contributions cyclostationnaires aux ordres supérieurs. Le signal devra avoir une composante continue nulle, c'est-à-dire $\frac{1}{K\Theta}\sum_{\theta=0}^{K \Theta}x\left(\theta\right)=0$. Cette condition n'est pas restrictive puisque l'on s'intéresse uniquement à la partie synchrone. On pourra également si besoin est, rajouter la moyenne classique au seuil et au signal après traitement.

$$x\left(\theta\right)=m\left(\theta\right)+r\left(\theta\right)$$ (3.6)

Où la moyenne synchrone calculée sur $K$ cycles est $\hat{m}\left(\theta\right)=\left<x\left(\theta\right)\right>_\Theta^K$.

Soit $a_{\theta_1}\left(k\right)=x\left(\theta_1+k \Theta\right)$, le signal obtenu en gardant les échantillons distants de $\Theta$ à partir de $\Theta_1$. Nous supposerons pour établir le seuil que pour un angle $\theta_1$ donné, $a_{\theta_1}\left(k\right)$ suit une loi normale de moyenne $m\left(\theta_1\right)$ et d'écart type $\sigma_r\left(\theta_1\right)$. Après moyennage synchrone, $\hat{m}\left(\theta_1\right)=\left<a_{\theta_1}\left(k\right)\right>^K$ suivra donc une loi normale de moyenne $m\left(\theta_1\right)$ et d'écart type $\sigma_r\left(\theta_1\right)/\sqrt{K}$.

L'écart type de $x\left(\theta\right)$ notée $\hat{\sigma}_b\left(\theta\right)$ sera estimée sur $K$ cycles par :

$$\hat{\sigma}_b\left(\theta\right)=\sqrt{\frac{K}{K-1}\left<... ...a\right)-\hat{m}\left(\theta\right)\right]^2\right>_\Theta^K}$$ (3.7)

Nous allons ensuite tester si $m\left(\theta\right)\neq 0$, c'est-à-dire si la moyenne synchrone existe pour un angle donné. Étant donné que l'on utilise un estimateur pour l'écart type, la quantité $\frac{\hat{m}\left(\theta\right)}{\hat{\sigma}_r\left(\theta\right)/\sqrt{K}}$ suivra une loi de Student $t_\beta$ à $K-1$ degrés de liberté.

Pour valider l'hypothèse $x\left(\theta\right)$ cyclostationnaire à l'ordre $1$ avec un risque de se tromper de $\alpha$, on exigera :

$$\left\vert\frac{\hat{m}\left(\theta\right)}{\hat{\sigma}_b\left(\theta\right)/\sqrt{K}}\right\vert\ge t_{1-\alpha}$$ (3.8)

On ne peut rien conclure si le signal est en dessous du seuil : il s'agit soit de bruit, soit d'une moyenne synchrone trop faible (auquel cas il faut augmenter le nombre de cycles).

Cette approche peut être discutable en présence de défauts puisque le kurtosis (c'est-à-dire la non gaussianité) de la série $a_{\theta_1}\left(\theta\right)$ a été utilisée pour bâtir un indicateur de défaut [#!raad:contributions!#]. Néanmoins, elle reste pertinente avant et lors de l'apparition du défaut, car une fois le défaut établi, le seuil n'est plus utile puisque ce défaut devient clairement identifiable dans la moyenne synchrone.

La figure illustre l'utilisation du seuil. La figure (a) montre un cycle d'un signal accéléromètrique issu d'un engrenage auquel un bruit gaussien a été ajouté afin de bien mettre en évidence le rôle du seuil.

La figure (b) montre la moyenne synchrone et le seuil associé, calculés sur $K=2500$ tours. Étant donné le nombre de cycles, on obtient un seuil très bas. Ce signal pourra donc être considéré comme étant la composante périodique. Il sera comparé au signal (c) qui est une moyenne synchrone sur $K=30$ tours. Le seuil qui est ici beaucoup plus élevé, va indiquer une zone correspondant à une moyenne synchrone nulle avec un risque de se tromper de $\alpha=5\;\%$. En conséquence, toute valeur à l'intérieur de la zone pourra être nulle et ne sera pas interprétable (ceci est facilement vérifiable en comparant les courbes (b) et (c) aux alentours du $350^\text{\tiny ème}$ échantillon).

La courbe (d) suggère une autre utilisation du seuil : la détection de cyclostationnarité. Pour cette courbe, une période de 502 points a été utilisée au lieu de 512 points (1 tour). Dès lors, on a un signal qui n'est plus cyclostationnaire à l'ordre 1 pour la période cyclique 502. Il est donc naturellement en dessous du seuil. En tirant parti de cette remarque, il est alors possible en comptabilisant le pourcentage de signal au-dessus du seuil pour différentes périodes de détecter le ou les cycles présents dans le signal. Cette utilisation est suggérée à titre indicatif puisqu'il existe de meilleurs outils pour trouver les fréquences cycliques (cepstre, corrélation spectrale, ...).

Figure: Application du seuil à la moyenne synchrone

3.2.1.3 Signaux presques-cyclostationnaires :

Il convient de distinguer deux cas:

• les fréquences sont incommensurables : dans ce cas il ``suffira'' de calculer les moyennes synchrones respectives suivant les différentes périodes. La moyenne est sans biais car elle est destructrice pour toutes les autres périodes. De tels signaux d'apparence simple sont la cause de nombreuses difficultés :

  o

  Si le nombre d'échantillons par cycle est constant pour une source périodique, il ne sera pas constant pour l'autre, car, les deux périodes sont incommensurables. Il sera donc nécessaire de resynchroniser au préalable ces signaux sur la période étudiée avant de calculer la moyenne synchrone. Le chapitre présente un exemple de signaux possèdant deux fréquences cycliques différentes.

  o

  Si ces fréquences sont proches, la séparation est rendue plus ardue lorsque le nombre de cycles disponibles est faible. En effet, le filtre en peigne associé à la moyenne synchrone a des dents assez ``larges'' et ne pourra donc pas facilement séparer deux fréquences proches.
  

• les fréquences sont commensurables : deux approches sont alors possibles.
  1. Considérer les signaux comme cyclostationnaires selon une période commune. Cette approche, bien que légitime conduit à une perte d'information puisque l'on n'exploite plus les périodes respectives mais une seule période commune. Sur les engrenages industriels, cette période commune est en général la plus longue possible afin d'éviter l'excitation des résonnances du système et la dégradation rapide de la boîte à vitesse.
  2. Une autre approche, produisant des résultats non optimaux consiste à étudier chacune des sous périodes indépendamment. Malheureusement, les périodes étant dans un rapport commensurable la moyenne est biaisée car elle n'est pas totalement destructrice pour les autres périodes même si la taille du signal est infinie. L'efficacité de la moyenne synchrone sera réduite.

     Pour illustrer ce biais nous avons construit sur la figure un signal composé de la somme d'une composante périodique $p_1\left(n\right)$ de période $\Theta_1=400$ échantillons et d'une autre composante périodique $p_2\left(n\right)$ de période $\Theta_2=300$ échantillons. Le maximum de chacune des composantes est d'amplitude unitaire. La période commune du signal est alors $\Theta_{com}=1200$ échantillons. La figure montre la moyenne synchrone calculée relativement à la période $\Theta_1$.

     La moyenne synchrone calculée sur une durée infinie de période commune est formée du motif récurrent de $p_1\left(n\right)$ (flèche notée C) que l'on souhaitait isoler mais également de $4$ motifs récurents issus de $p_2\left(n\right)$ dont l'amplitude a été divisée par quatre (flèches notées D). En effet, dans un cycle commun, il y a $3$ motifs issus de $p_1\left(n\right)$ qui subissent un moyennage constructif et $4$ motifs issus de $p_2\left(n\right)$ qui subissent un moyennage destructif (ce qui explique leur amplitude de 0.25).

     Figure: Moyennage synchrone de signaux de périodes commensurables

     

3.2.2 Caractérisation à l'ordre 2

A l'ordre $2$, il est possible d'utiliser $4$ représentations pour caractériser la $\CSP{2}$. Chaque représentation est liée soit à un repère angulaire, soit à un repère fréquentiel associé $f_\theta$ soit à un repère mixte. Pour le repère angulaire nous utiliserons l'angle $\theta$ ainsi que le retard $\tau_\theta$ pour la corrélation. Bien que classiquement le domaine fréquentiel soit gradué en $Hz$, l'utilisation d'une telle graduation pour des signaux dépendant de l'angle nécessite quelques approximations. Dans le domaine angulaire, le nombre de points par tour $ppt$ est constant. Pour deux fréquences de rotation distinctes, la quantité $1/ppt$ restera constante. Dès lors, on ne pourra plus appeler cette quantité une fréquence. La quantité $\ordre=1/ppt$ sera alors appelé un ordre [#!bruel:order!#], l'ordre $n$ correspond à $n$ fois la fréquence de rotation instantanée. L'échelle fréquentielle en $Hz$ étant en général plus familière, il sera toutefois possible de l'utiliser si les variations de vitesse sont faibles autour d'une vitesse de rotation moyenne constante, nous aurons alors $f=f_\theta$ (de l'ordre de quelques pour cent ou pour mille autour de la fréquence de rotation moyenne). Il est bien sûr nécessaire de mesurer au préalable la vitesse de rotation moyenne.

3.2.2.1 Densité Spectrale de corrélation

Pour définir un spectre, il est nécessaire de calculer la transformée de Fourier du signal. Malheureusement, la transformation de Fourier définie par ${\TF{x}\left(f\right)=\int_{\mathbb{R}}x\left(\theta\right) e^{-2\pi j f \theta} d\theta}$ n'est pas convergente, en moyenne quadratique lorsque $x\left(\theta\right)$ est un signal aléatoire. Il convient alors d'utiliser la décomposition de Cramér présentée dans [#!cramer:mathematical!#,#!blanc:theorie!#] :

$$x\left(\theta\right) = \int_{\mathbb{R}}e^{2 \pi j f \theta} d\TF{X}\left(f\right)$$ (3.9)

Où $d\TF{X}\left(f\right)$ est l'incrément spectral à la fréquence $f$. Pour $x\left(\theta\right)$ centré (approche cumulant), on a :

$$\esp{d\TF{X}\left(f_1\right)d\TF{X^{*}}\left(f_2\right)} = d^2 \rho_{x\left(2\right)}\left(f_1,f_2\right)$$ (3.10)

où $d^2 \rho_{x\left(2\right)}\left(f_1,f_2\right)$ est une mesure de puissance qui quantifie l'interaction entre les deux canaux $f_1$ et $f_2$ telle que :

$$\cumns{x}{2}{\theta_1,\theta_2}=\int_{\mathbb{R}^2} e^{2 \p... ...ta_2\right)} d^2 \rho_{x\left(2\right)}\left(f_1,f_2\right)$$ (3.11)

La mesure $d^2 \rho_{x\left(2\right)}\left(f_1,f_2\right)$ est une différentielle au second ordre d'une fonction certaine $\rho_{x\left(2\right)}\left(f_1,f_2\right)$ homogène à une distribution spectrale de la fonction d'autocovariance $\cumns{x}{2}{\theta_1,\theta_2}$. Par suite, on admettra qu'il existe une densité de corrélation spectrale (DCS) $\wv{x}\left(f_1,f_2\right)$ telle que :

$$d^2 \rho_{x\left(2\right)}\left(f_1,f_2\right) = \wv{x}\left(f_1,f_2\right) df_1 df_2$$ (3.12)

La DCS dépend alors de deux fréquences et s'identifie à un spectre lorsque ces dernières sont identiques : $\wv{x}\left(f,f\right) d^2f=\esp{d\TF{X}\left(f\right)d\TF{X^{*}}\left(f\right)}$, et, dans le cas général, à l'inter-spectre entre deux versions décalées dans le domaine fréquentiel de $\Delta_f=f_2-f_1$ : $\wv{x}\left(f_1,f_1+\Delta_f\right) df_1 df_2=\esp{d\TF{X}\left(f_1\right)d\TF{X^{*}}\left(f_1+\Delta_f\right)}$.

Dans le cas cyclostationnaire, la périodicité du cumulant d'ordre $2$ : $\cumns{x}{2}{\theta_1,\theta_2}=\cumns{x}{2}{\theta_1+\Theta,\theta_2+\Theta}$ appliquée aux équations () et () impliquera :

0

$$\int_{\mathbb{R}^2} e^{2 \pi j \left(f_1 \theta_1 - f_2 \theta_2\right)} d^2 \rho_{x\left(2\right)}\left(f_1,f_2\right) = \int_{\mathbb{R}^2} e^{2 \pi j \left(f_1 \theta_1 - f_2 \theta_2\right)} \wv{x}\left(f_1,f_2\right) df_1 df_2$$ (3.13)

$$= \int_{\mathbb{R}^2} e^{2 \pi j \left(f_1 \theta_1 - f_2 \theta_2\... ...t)} e^{2\pi j \Theta \left(f_1-f_2\right)} \wv{x}\left(f_1,f_2\right) df_1 df_2$$ (3.14)

Soit :

$$d^2 \rho_{x\left(2\right)}\left(f_1,f_2\right) = \left\{ \... ...,\quad k\in\mathbb{Z} \\ 0 & ailleurs \end{array} \right.$$ (3.15)

La mesure de puissance sera donc nulle partout sauf sur l'ensemble des droites d'équation $f_1=f_2+\frac{k}{\Theta}$ (cf. figure ). Il apparaît alors intéressant de faire le changement de repère :

$$\left\{ \begin{array}{l} f=\frac{f_1+f_2}{2} \\ \alpha=f_1-f_2 \end{array} \right.$$ (3.16)

Ce changement de repère permet de faire apparaître la fréquence cyclique $\alpha$. Cette fréquence cyclique prend tout son sens en réécrivant l'équation dans ce nouveau repère :

$$\wv{x}\left(f,\alpha\right) df d\alpha=\esp{dX\left(f+\frac{\alpha}{2}\right) dX^*\left(f-\frac{\alpha}{2}\right)}$$ (3.17)

La DCS apparaît donc comme le produit de deux versions décalées en fréquence d'une quantité $\alpha$ des incréments spectraux à la fréquence $f$. Une valeur non nulle de la DCS caractérisera donc un lien entre les fréquences $f-\frac{\alpha}{2}$ et $f+\frac{\alpha}{2}$.

On peut extraire la partie intéressante de cette mesure de puissance en définissant la densité spectrale cyclique (DSC) qui correspond aux valeurs non nulles de la DCS :

$$\left\{ \begin{array}{lcl} \speccorr{x}^{\alpha}\left(f\r... ...lpha}{\text{\tiny$\frac{k}{\Theta}$}} \end{array} \right.$$ (3.18)

La densité de corrélation est alors continue par rapport à la variable $f$ et discrète par rapport à la fréquence cyclique $\alpha$. Ce résultat est une propriété remarquable des signaux cyclostationnaires.

L'équation () est la forme symétrique de la DCS. Il existe également une forme assymétrique définie par :

$$\wv{x}\left(f,\alpha\right) df d\alpha=\esp{dX\left(f+\alpha\right) dX^*\left(f\right)}$$ (3.19)

L'estimation de la corrélation spectrale se fera en deux étapes comme l'indique la figure :

1. Pour une fréquence cyclique donnée $\alpha$, on effectue deux décalages fréquentiels du signal apodisé en le multipliant respectivement par $e^{+i \pi \alpha \theta}$ et $e^{-i \pi \alpha \theta}$. On obtient alors deux signaux.
2. On calcule ensuite l'interspectre de ces deux signaux. L'interspectre obtenu correspond alors à une coupe de la corrélation spectrale à la fréquence cyclique $\alpha$ choisie.

Comme la corrélation spectrale est basée sur le calcul d'un interspectre, on utilisera les mêmes techniques d'estimation : le périodogramme moyenné et le périodogramme lissé.

La figure illustre ces deux techniques.

La méthode du périodogramme moyenné est basée sur le découpage du signal en $K$ blocs de $M$ points. La corrélation spectrale est tout d'abord calculée sur chacun de ces blocs. Les corrélations spectrales ainsi obtenues sont ensuite moyennées.

La méthode du périodogramme lissé est basée sur le lissage d'une corrélation spectrale calculée sur le signal complet. On sous-échantillonne ensuite cette image selon l'axe fréquentiel.

Figure: Estimation de la corrélation spectrale

D'après les calculs de complexité de [#!bouillaut:approches!#], la méthode du périodogramme moyenné devient plus avantageuse à partir de $4$ moyennes, de plus, elle nécessite moins de mémoire.

Lors de l'utilisation du périodogramme moyenné, il faudra prendre garde à la méthode d'évaluation de l'interspectre. En effet, dans les logiciels comme Matlab, la transformée de Fourier est calculée en considérant que chacun des blocs débute à un instant nul ($n=0$). Dès lors, pour la $k^\text{\tiny ième}$ tranche, on introduira un déphasage de $e^{2 \pi j \alpha k M}$ [#!capdessus:aide!#] qu'il suffira de compenser pour estimer correctement la corrélation spectrale.

Le périodogramme moyenné conduit à faire un compromis entre le biais et la résolution fréquentielle. En effet, le signal étant de taille finie, si l'on augmente le nombre de blocs on diminue le biais, mais, un bloc de taille réduite conduit à une résolution fréquentielle plus faible.

Un dilemme biais variance apparaît pour l'estimateur basé sur le périodogramme lissé lors du choix de la fenêtre de lissage : une fenêtre de lissage de durée importante permet de réduire le biais mais la variance dépend de la forme de cette fenêtre.

Signaux presque-cyclostationnaires

Dans le cas presque-cyclostationnaire, le cumulant d'ordre $2$ peut être approximé uniformément par un polynôme trigonométrique de la forme :

$$\cumns{x}{2}{\theta,\tau} \approx \sum_{n=1}^{N} c_n\left(\tau\right) e^{2\pi j \nu_n \theta}$$ (3.20)

Cette formulation conduira à une mesure de puissance nulle partout sauf sur l'ensemble des droites d'équation $f_1=f_2+\nu_n$. La figure montre ce support pour un signal constitué de deux pseudo-périodes cycliques ainsi que de leurs harmoniques. Dès lors, les valeurs de la DSC non nulles correspondront non plus à une fréquence cyclique et à ces harmoniques mais à l'ensemble des pseudo-périodes $\left\{\nu_n\right\}$ (incluant les harmoniques).

Figure: Support de la DCS

Parmi les applications de la DSC, on pourra citer le degré de cyclostationnarité introduit par Gardner et repris dans [#!zivanovic:degree!#] :

$$D_{CS}=\frac{\sum_{\alpha\neq 0} \int_{\mathbb{R}} \left\ve... ...eft\vert S^0_{x\left(2\right)}\left(f\right)\right\vert^2 df}$$ (3.21)

Ce degré de cyclostationnarité peut être interprété comme la distance avec le processus stationnaire le plus proche [#!gardner:characterization!#]. Une étude et une généralisation des mesures de cyclostationnarité aux ordres supérieurs pourront être trouvés dans [#!raad:contributions!#].

3.2.2.2 Autocovariance

L'autocovariance correspond au cumulant d'ordre 2. Puisqu'elle est périodique dans le contexte cyclostationnaire: $\cumns{X}{2}{\theta_1,\theta_2}=\cumns{X}{2}{\theta_1+\Theta,\theta_2+\Theta}$, il est classique de faire le changement de repère:

$$\left\{ \begin{array}{l} \theta=\frac{\theta_1+\theta_2... ...\\ \tau=\frac{\theta_2-\theta_1}{2} \end{array} \right.$$ (3.22)

On utilise également la variance synchrone définie par $\cumns{X}{2}{\theta,0}$ qui correspond à la puissance instantanée du système.

Cas cyclostationnaire

Dans le cas cyclostationnaire on a $\cumns{X}{2}{\theta,\tau}=\cumns{X}{2}{\theta+m \Theta,\tau+n \Theta}$ quels que soient les entiers $m$ et $n$. L'autocovariance admet donc un développement en série de Fourier dont les coefficients sont donnés par :

$$\cum{X}{2}^\alpha\left(\tau\right)=\frac{1}{\Theta} \int_0^... ...text{ avec $\alpha=\frac{k}{\Theta}, \quad k\in\mathbb{Z}$}$$ (3.23)

Où les fréquences cycliques $\alpha$ sont multiples de la fréquence de base du système : $\alpha=\frac{k}{\Theta}$, k entier.

La fonction $\cum{X}{2}^\alpha\left(\tau\right)$ est alors appelée fonction d'autocovariance cyclique. Elle est reliée à la DSC par :

$$\speccorr{x}^{\alpha}\left(f\right) = \int_\mathbb{R} \cum{X}{2}^\alpha\left(\tau\right) e^{-2 \pi j f \tau} d\tau$$ (3.24)

Cas presque-cyclostationnaire

Dans le cas presque-cyclostationnaire la variance peut être approximée par un polynôme trigonométrique qui admet une décomposition de Fourier-Bohr (équation ). Ce polynôme fournit directement la valeur de la variance cyclique : $\cum{X}{2}^\alpha\left(\tau\right)=c_{\alpha}\left(\tau\right)$.

A la différence du cas cyclostationnaire où l'ensemble des valeurs $\alpha$ correspondait au fondamental et à ses harmoniques, elles correspondent ici à un ensemble de fréquences données.

La variance cyclique sera alors reliée à la DSC par :

$$\speccorr{x}^{\alpha}\left(f\right) \approx \int_\mathbb{R} \cum{X}{2}^\alpha\left(\tau\right) e^{-2 \pi j f \tau} d\tau$$ (3.25)

L'égalité ne peut pas être utilisée ici car le polynôme trigonométrique approche seulement la variance.

3.2.2.3 Spectre de Wigner-Ville

Nous traitons ici des signaux stochastiques. Or, la distribution de Wigner-Ville $\wv{x}\left(\theta,f\right)$ a été originellement conçue pour des signaux déterministes d'énergie finie par :

$$\wv{x}\left(\theta,f\right)=\int_{\mathbb{R}} x^* \left(\th... ...x\left(\theta+\frac{\tau}{2}\right) e^{-2 \pi j f \tau} d\tau$$ (3.26)

Pour pouvoir manipuler des signaux aléatoires, nous allons utiliser le spectre de Wigner-Ville [#!martin:time!#,#!antoni:apports!#] $\swv{X}\left(\theta,f\right)$ qui est défini comme l'espérance mathématique de la distribution de Wigner-Ville du processus aléatoire centré. Ainsi, pour un processus aléatoire $X\left(\theta\right)$ on aura :

$$\swv{X}\left(\theta,f\right)=\esp{\wv{X}\left(\theta,f\right)}$$ (3.27)

Sous réserve de commutativité entre l'espérance et l'intégration, il est possible de lier l'autocovariance et le spectre de Wigner-Ville par :

$$\swv{X}\left(\theta,f\right) = \int_\mathbb{R} \cumns{X}{2}{\theta,\tau} e^{2 \pi j f \tau} d\tau$$ (3.28)

Lorsque le processus stochastique $X\left(\theta\right)$ est (pseudo)-cyclostationnaire, le spectre de Wigner-Ville présente une structure (pseudo)-périodique suivant la variable angulaire $\theta$.

En combinant les équations () et (), on peut écrire :

$$\swv{X}\left(\theta,f\right) = \int_\mathbb{R} e^{-2 \pi j f \tau} \sum_{\alpha\in \left\{\nu_k\right\}} \cum{X}{2}^\alpha\left(\tau\right) e^{-2 \pi j \alpha \theta} d\tau$$ (3.29)

$$= \sum_{\alpha \in \left\{\nu_k\right\}} \speccorr{x}^{\alpha}\left(f\right) e^{-2 \pi j \alpha \theta}$$ (3.30)

Dans le cas cyclostationnaire, on aura $\left\{\nu_k\right\}=\frac{k}{\Theta}$. Dans le cas presque-cyclostationnaire, les fréquences $\nu_k$ ne seront pas forcément multiples entre elles, de plus, il conviendra de remplacer l'égalité par une approximation.

Cette équation établit le lien entre le spectre de Wigner-Ville et la DSC.

La figure résume les relations entre les différentes représentations à l'ordre $2$.

Figure: Relations entre les différentes représentations à l'ordre 2

3.2.3 Alternative dans le cas quasi-cyclostationnaire

Lorsque les signaux sont quasi-cyclostationnaires (il existe des périodes incommensurables entre elles), il est possible d'utiliser l'opérateur d'extraction de composante périodique $\left<x\left(\theta\right)\right>_{\Theta}$ défini par l'équation ().

Comme indiqué précédemment, cet opérateur apparaît comme un filtre en peigne sélectionnant uniquement les composantes de période $\Theta$.

On peut alors généraliser cet opérateur à l'ordre 2 par :

$$\left<\cum{X}{2}\left(\theta_1,\theta_2\right)\right>_\Thet... ...)-\left<X\left(\theta_2\right)\right>_\Theta\right]^*}\right>$$ (3.31)

Nous avons proposé ici deux méthodes pour traiter le cas presque cyclostationnaire : soit de prendre en considération et obtenir des statistiques presque-cyclostationnaires, soit extraire chaque période (sous réserve d'existence c'est-à-dire de quasi-cyclostationnarité) et travailler alors sur des signaux cyclostationnaires. Le choix de la méthode dépendra alors du type de signal et des traitements que l'on veut faire. Par exemple, si chaque période a une signification physique, la dernière méthode deviendra intéressante.

3.3 Application au diagnostic

Ce paragraphe a pour but de montrer un exemple d'exploitation de la cyclostationnarité pour le diagnostic de défaut d'engrenages. Les signaux accéléromètriques ont été acquis au centre d'études et recherches EDF à Chatou sur un banc d'essai de fatigue. Ces signaux ont notamment été exploités dans [#!fontanive:surveillance!#] et [#!elbadaoui:contribution!#]. Nous proposons ici une autre analyse basée sur la cyclostationnarité.

Il s'agit d'un banc d'essai composé d'une roue de $56$ dents et d'une autre de $15$ dents. La fréquence de rotation de la roue menante est de $12,5\;Hz$ et la fréquence d'échantillonnage de $6,4\;kHz$. Afin d'étudier l'apparition d'un défaut, un essai de fatigue a été réalisé. Tout au long de cet essai 15 acquisitions ont été faites. Un défaut de type écaillage a été détecté à partir de la dixième mesure. Le défaut s'est aggravé progressivement jusqu'à la rupture d'une dent lors du quinzième relevé.

Dans les précédentes publications, le défaut a été détecté lors du $10^\text{\tiny ième}$ essai. Par ailleurs, les signaux ont été qualifiés de ``fortement bruités''. Nous allons, en exploitant la cyclostationnarité retrouver ces mêmes résultats. Pour cela, nous avons préalablement re-échantillonné les signaux dans le domaine angulaire en estimant la phase instantanée par démodulation autour de deux fois la fréquence d'engrènement.

La figure montre les $15$ acquisitions pendant une durée correspondant à la période commune (intervalle au bout duquel l'engrenage revient dans la même configuration (mêmes dents en contact). La période commune $\Theta_{com}$ représente $15$ tours de la roue de $56$ dents ou $56$ tours de la roue de $15$ dents (15 et 56 n'ayant pas de multiple commun). Chacun des signaux a été normalisé par rapport à leur valeur crête-à-crête avant affichage afin de faciliter leur comparaison. Le défaut apparaît nettement le $11^\text{\tiny ième}$ jour alors qu'il avait été détecté par un expert le $10^\text{\tiny ième}$ jour.

Figure: Premiers points du signal

En toute rigueur les signaux d'engrenage sont cyclostationnaires, c'est-à-dire qu'il existe une période cyclique. Cette période correspond au cycle commun de l'engrenage. Nous avons préféré les qualifier de polycyclostationnaires puisque physiquement cette période est composée de ``sous'' périodes correspondant aux périodes de rotation des roues qui sont dans un rapport rationnel, il serait alors dommage de ne pas exploiter ces informations. Nous allons exploiter la moyenne synchrone sur la période commune $\Theta_{com}$ afin d'examiner l'apparition au dixième jour.

Si on considère le signal cyclostationnaire à l'ordre $2$, on le modélise sous la forme :

$$x\left(\theta\right)=p\left(\theta\right)+r\left(\theta\right)$$ (3.32)

Où $p\left(\theta\right)$ est un signal périodique de période $\Theta_{com}$ que l'on estimera à l'aide d'une moyenne synchrone. Le signal ne contient malheureusement que $7$ blocs de taille $\Theta_{com}$. La figure montre l'estimation de cette période commune, le résidu est présenté sur la figure . On peut constater que la période commune est très proche du signal et que le résidu n'apporte que peu d'information nouvelle.

Figure: Période commune $p\left(\theta\right)$

Figure: Signal résiduel $r\left(\theta\right)$

Devant l'insuccès de l'approche cyclostationnaire, nous allons utiliser une approche polycyclostationnaire. Pour cela le signal sera modélisé par :

$$x\left(\theta\right)=p_1\left(\theta\right)+p_2\left(\theta\right)+p_{1,2}\left(\theta\right)+r\left(\theta\right)$$ (3.33)

Où le signal est composé de deux signaux périodiques $p_1\left(\theta\right)$ de période $\Theta_1$ et $p_2\left(\theta\right)$ de période $\Theta_2$ dont la période correspond respectivement à la $1^{\text{\tiny {ère}}}$ et à la $2^{\text{\tiny {ème}}}$ roue ; de $p_{1,2}\left(\theta\right)$ de période $\Theta_{com}$ correspondant à l'interaction des $2$ roues ; d'un signal résiduel $r\left(\theta\right)$ non périodique. Nous supposerons que les périodes $\Theta_1$ et $\Theta_2$ ne sont pas multiples entre elles mais dans un rapport rationnel. Dans le cas contraire, cette approche n'aurait aucun intérêt.

La composante précédemment estimée, $p\left(\theta\right)$ correspond à la somme des trois composantes périodiques. L'estimation de $p_1\left(\theta\right)$ se fera en calculant la moyenne synchrone $\left<x\left(\theta\right)\right>_{\Theta_1}$.

Malheureusement, l'estimateur de $p_1\left(\theta\right)$ est biaisé. En effet, comme $p_{1,2}\left(\theta\right)$ est périodique d'une part et que $\Theta_1$ est sous-multiple de cette période, une moyenne synchrone calculée sur un nombre infini de blocs fera toujours apparaître un nombre fini $\alpha=\frac{\Theta_{com}}{\Theta_1}$ de composantes décalées de $p_{1,2}\left(\theta\right)$ lors du moyennage :

$$\lim_{N\rightarrow\infty} \left<x\left(\theta\right)\right>... ...ta+k \Theta_1\right) + r\left(\theta+k \Theta_1\right)\right]$$ (3.34)

Où :

• $\lim_{N\rightarrow\infty} \frac{1}{N\alpha}\sum_{k=0}^{N\alpha-1} p_1\left(\theta+k \Theta_1\right)$ correspondra, par définition à $p_1\left(\theta\right)$,
• $\lim_{N\rightarrow\infty} \frac{1}{N\alpha}\sum_{k=0}^{N\alpha-1} p_2\left(\theta+k \Theta_1\right)=0$ puisque $p_2\left(\theta\right)$ est périodique de période $\Theta_2$ qui n'est pas multiple de $\Theta_1$,
• $\lim_{N\rightarrow\infty} \frac{1}{N\alpha}\sum_{k=0}^{N\alpha-1} p_{1,2}\left(... ...ht)=\frac{1}{\alpha}\sum_{k=0}^{\alpha-1} p_{1,2}\left(\theta+k \Theta_1\right)$ puisque $p_{1,2}\left(\theta\right)$ est périodique de période $\Theta_{com}$ et que $\alpha$ périodes sont contenues dans une période $\Theta_{com}$,
• $\lim_{N\rightarrow\infty} \frac{1}{N\alpha}\sum_{k=0}^{N\alpha-1} r\left(\theta\right)=0$ par définition puisque $r\left(\theta\right)$ n'est pas périodique de période $\Theta_1$.

Soit :

$$\lim_{N\rightarrow\infty} \left<x\left(\theta\right)\right>... ...2}\left(\theta+k \Theta_1\right) \neq p_1\left(\theta\right)$$ (3.35)

Néanmoins, si les nombres de dents sont élevés et qu'ils sont premiers entre eux, le nombre $\alpha$ sera également élevé. Il est tout de même important de noter que l'estimateur de $p_1\left(\theta\right)$ est quoiqu'il en soit biaisé.

Pour l'engrenage de EDF, nous allons essayer d'évaluer la moyenne associée à l'arbre 1 :

$$p_1\left(\theta\right)+p_{1,2}\left(\theta\right)$$ (3.36)

Pour cela, nous avons évalué la moyenne synchrone suivant la période commune, pour estimer la contribution :

$$p_1\left(\theta\right)+p_2\left(\theta\right)+p_{1,2}\left(\theta\right)$$ (3.37)

Puis nous avons retranché $p_2\left(\theta\right)$ estimé par moyennage synchrone sur une période $\Theta_2$ ($\alpha=56$). Le résultat est présenté sur la figure .

Figure: Signal après traitement

Bien que l'estimation soit biaisée $\left(\alpha=56\right)$, le défaut est nettement visible et au $10^\text{\tiny ième}$ jour (il se manifeste par l'apparition de pics dans le signal).

Ce signal a été qualifié de fortement bruité, nous avons établi ici que cette source de bruit provenait de la contribution de la deuxième roue. Cette composante, fortement énergétique lors des premières acquisitions masque les premiers signes de l'apparition du défaut. Lorsque le défaut progresse, l'énergie associée à la roue 1 porteuse du défaut devient plus importante et le fait ainsi ressortir sur le signal non traité.

De nombreuses autres approches ou méthodes complémentaires pourront être envisagées. Nous nous sommes néanmoins limités à cette dernière puisqu'elle est très simple et qu'elle illustre bien le problème du calcul de la moyenne synchrone dans le cas d'engrenages.

3.4 Conclusion

Nous avons défini dans ce chapitre la notion de cycloergodicité qui nous permet d'utiliser une seule réalisation de durée suffisamment longue. En pratique, nous supposons implicitement cette propriété vérifiée. Nous avons également défini les outils permettant de caractériser la cyclostationnarité et la quasi-cyclostationnarité aux ordres $1$ et $2$. Enfin, un exemple montre les limites de l'approche cyclostationnaire et l'intérêt de l'approche polycyclostationnaire pour les signaux d'engrenages.

Le chapitre suivant va illustrer l'utilisation de la cyclostationnarité dans le cadre de signaux pseudo-quasicyclostationnaires sur les signaux de roulement.