Estimateur cohérent

11.2.3 Réduction de la variance

Un examen sommaire des estimations spectrales pour les séries de jauges de fil dans les figures 5 et 6 révèle une variabilité substantielle entre les fréquences, à tel point qu’il est difficile de discerner la structure globale dans les estimations spectrales sans une bonne quantité d’étude. Tous les estimateurs spectraux directs souffrent de cette instabilité inhérente, qui peut s’expliquer en considérant les propriétés distributives de S^X(d)(f). Premièrement, si f n’est pas trop proche de 0 ou de f(N) et si SW(⋅) satisfait à une condition de régularité légère, alors 2S ^w(d)(f) / Sw(f) = dx22; c’est-à-dire que le rv 2S ^w(d)(f) / Sw(f) est approximativement égal en distribution à un chi carré rv avec 2 degrés de liberté. Si le rétrécissement n’est pas utilisé, f est considéré comme “pas trop proche” de 0 ou f(N) si 1 /(n-p) Δt < f < f(N) -1/(n-p) Δt; si on utilise un effilage, il faut remplacer 1/(n-p)Δt par un terme plus grand, reflétant la largeur accrue du lobe central de la fenêtre spectrale (par example, le terme pour la conicité des données de Hanning est d’environ 2/(n-p)Δt donc f est “pas trop proche” si 2/(n-p)Δt < f < f(N)-2/(n-p)Δt).

Puisqu’un chi carré rv xv2 avec v degrés de liberté a une variance de 2u, nous avons l’approximation V = Sw2(f). Ce résultat est indépendant du nombre de Wt, on a: contrairement à des statistiques telles que la moyenne de l’échantillon des VR gaussiennes indépendantes et identiquement distribuées, la variance de S^W(d)(f) ne diminue pas à 0 lorsque la taille de l’échantillon n-p augmente (sauf dans le cas sans intérêt Sw(f) = 0). Ce résultat explique le caractère saccadé des estimations spectrales directes présentées dans les figures 5 et 6. En terminologie statistique, S^W(d)(f) est un estimateur incohérent de Sw(f).

Nous décrivons maintenant trois approches pour obtenir un estimateur cohérent de Sw(f). Chaque approche est basée sur la combinaison de VR qui, sous des hypothèses appropriées, peuvent être considérées comme des estimateurs non corrélés approximativement par paires de SW(f). En bref, les trois approches consistent à

1

lisser S^W(d)(f) sur toutes les fréquences, ce qui donne ce qu’on appelle un estimateur spectral de fenêtre de décalage;

2

{}} (ou {Wt}) en un certain nombre de segments (dont certains peuvent se chevaucher), calculer une estimation spectrale directe pour chaque segment, puis faire la moyenne de ces estimations ensemble, ce qui donne l’estimateur spectral de segmentaveraging chevauché (WOSA) de Welch;

3

calculez une série d’estimations spectrales directes pour {Wt} en utilisant un ensemble de données orthogonales, puis faites la moyenne de ces estimations ensemble, ce qui donne l’estimateur spectral multitaper de Thomson.

Estimateurs spectraux de fenêtre de décalage Un estimateur spectral de fenêtre de décalage de Sw (⋅) prend la forme

(11,15) S^W(lw)(f) =∫-f(N) f(N) Wm(f-f’) S^W(d)(f’) df’

où Wm (⋅) est une fenêtre de lissage dont les propriétés de lissage sont contrôlées par le paramètre de lissage m. En d’autres termes, l’estimateur S^W(lw)(⋅) est obtenu en convoluant une fenêtre de lissage avec l’estimateur spectral direct S^w(d)(⋅). Une fenêtre de lissage typique a à peu près la même apparence qu’une fenêtre spectrale. Il existe un lobe central dont la largeur peut être ajustée par le paramètre de lissage m: plus ce lobe central est large, plus S^W(w)(⋅) sera lisse. Il peut également y avoir un ensemble de lobes latéraux gênants qui provoquent des fuites de fenêtre de lissage. La présence d’une fuite de fenêtre de lissage est facilement détectée en superposant des tracés de S^W(lw) (⋅) et S^W(d) (⋅) et en recherchant des plages de fréquences où la première ne semble pas être une version lissée de la seconde.

Si l’on a utilisé un filtre de pré-blanchissement AR, on peut alors postcolorer S^W(lw)(⋅) pour obtenir un estimateur de Sx(⋅), à savoir,

SX(pc)(f) = S^W(tw)(f) | 1−∑k = 1pϕke-i2πfkΔt|2

Les propriétés statistiques de S^W (lw) (.) sont traçables en raison du résultat d’échantillon important suivant. Si S^W(d)(⋅) est en fait le périodogramme (c’est−à-dire que nous n’avons pas effilé les valeurs de Wt), l’ensemble des vr S^W(d)(j/(n-p)Δt), j = 1,2,…, j,, sont approximativement non corrélés deux à deux, chaque rv étant proportionnel à un χ22, rv (ici J est le plus grand entier tel que J/(n-p) < 1/2). Si nous avons utilisé l’effilage pour former Sw(d)(⋅), une déclaration similaire est vraie sur un plus petit ensemble de VR définis sur une grille plus grossière de fréquences également espacées — à mesure que le degré d’effilage augmente, le nombre de VR approximativement non corrélés diminue. Sous l’hypothèse que le sdf Sw (⋅) varie lentement d’une fréquence à l’autre (le pré-blanchissement aide à le rendre vrai) et que le lobe central de la fenêtre de lissage est suffisamment petit par rapport aux variations de Sw (⋅), il s’ensuit que S^W(d)(f) en Eq. (11.15) peut être approximé par une combinaison linéaire de rv non corrélées χ22. Un argument standard de “degrés de liberté équivalents” peut ensuite être utilisé pour approximer la distribution de S^W(lw)(f). (voir Eq. (11.17) plus tard).

Il existe deux façons pratiques de calculer S^W (lw) (⋅). La première façon est de discrétiser l’égalisation. (11.15), donnant un estimateur proportionnel à une convolution de la forme ΣkWm(f−fk’) SW(d)(fk’), où les valeurs offk’ sont un ensemble de fréquences également espacées. La deuxième façon est de rappeler que “la convolution dans un domaine de Fourier équivaut à la multiplication dans l’autre” pour réécrire Eq. (11,15) comme

(11,16) S ^W(lw)(f) =ττ = -(n−p−1) n−p−1wt, mC ^τ.w(e) e-l2n / τΔι

où C^ τ.W(d), est l’estimateur acvs donné en Eq. (11.9) correspondant à S^W(d) (.), et {wt.m} est une fenêtre de décalage (ceci peut être considéré comme la transformée de Fourier inverse de la fenêtre de lissage Wm(⋅)). En fait, parce que S ^ W(d) (.) est un polynôme trigonométrique, toutes les circonvolutions discrètes de la forme ΣkWm(f-fk’) S^W(d)(fk’) peuvent également être calculées via Eq. (11.16) avec un choix approprié de valeurs wt, m (pour plus de détails, voir Rubrique 6.7). Nos deux méthodes pratiques de calcul S^W(l, w) (.) donnent ainsi des estimateurs équivalents. À moins que la convolution discrète ne soit suffisamment courte, Eq. (11.16) est plus rapide en calcul à utiliser.

La théorie statistique suggère que, sous des hypothèses raisonnables

(11,17) vS ^W(lw)(f) Sw(f) = dxv2

à une bonne approximation, où v est appelé les degrés de liberté équivalents pour S^W (lw)(f) et est donné par v = 2 (n−p) BwΔt /Ch . Ici Bw est une mesure de la largeur de bande de la fenêtre de lissage Wm (⋅) a)n d peut être calculée via BW = 1/Δt∑T =-(n−p−1) n-p-1wT, m2;;d’autre part, Ch ne dépend que de la conicité appliquée aux valeurs de Wtand peut être calculée via Ch =(n−p)∑1 = p + 1nht4Note que, si nous ne rétrécissons pas explicitement, alors ht = 1 / n−pand donc Ch > 1; pour une conicité de données typique, l’inégalité de Cauchy nous indique que Ch > 1 (par exemple, Ch≈1,94 pour la conicité de données Hanning). Les degrés de liberté équivalents pour S ^ W (lw) (f) augmentent donc à mesure que nous augmentons la bande passante de la fenêtre de lissage et diminuent à mesure que nous augmentons le degré de réduction. L’équation (11.17) nous dit que E≈SW(f) et que V≈SW2(f)/v, donc l’augmentation de v diminue V.

L’approximation en Eq. (11.17) peut être utilisé pour construire un intervalle de confiance pour SW(f) de la manière suivante.Soit nv(α) le point de pourcentage α×100% de la distribution xv2 ; c’est-à-dire P = α.L’intervalle de confiance A100 (1 − 2α) % pour Sw(f) est approximativement donné par

(11.18)

Les points de pourcentage ην (α) sont tabulés dans de nombreux manuels ou peuvent être calculés à l’aide d’un algorithme donné par Best et Roberts

L’intervalle de confiance de (11,18) est gênant dans la mesure où sa longueur est proportionnelle à S ^ W (lw) (f). D’autre part, l’intervalle de confiance correspondant pour 10.log10(Sw(f)) (c’est-à-dire, SW(f) sur une échelle de décibels) est juste

qui a une largeur indépendante de S^W(lw)(.). C’est la raison d’être pour tracer les estimations de fds sur une échelle de décibels (ou logarithmique).

Un nombre déconcertant de fenêtres de décalage différentes a été discuté dans la littérature (voir). Nous ne donnons ici qu’un seul exemple, la célèbre fenêtre fag de Parzen (Parzen):

wt.m = 1 -6τ∼2 +6| τ∼|3, | τ|≤m/22 (1- τ∼) 3, m/2 <| τ/≤m0, / τ/ > m

où m est considéré comme un entier positif et τ = τ / m. Cette fenêtre de décalage est facile à calculer et a des lobes latéraux dont l’enveloppe se désintègre en f-4, de sorte que la fuite de fenêtre de lissage est rarement un problème. Pour une bonne approximation, la bande passante de la fenêtre de lissage pour la fenêtre de décalage de Parzen est donnée par Bw = 1,85/(mΔt). À mesure que m augmente, la bande passante de la fenêtre de lissage diminue et l’estimateur de fenêtre de décalage résultant devient moins lisse. Les degrés de liberté équivalents associés sont donnés approximativement par v = 3,71 (n-p)/(mCh). La fenêtre de décalage de Parzen pour m = 32 et sa fenêtre de lissage associée sont représentées sur la figure 7.

Fig.7. Fenêtre de décalage de Parzen (a) et la fenêtre de lissage correspondante (b) pour m = 32. La bande passante de la fenêtre de lissage Estbw = 0,058.

À titre d’exemple, la figure 8(a) montre un estimateur de fenêtre de décalage postcolorée pour les données de jauge d’onde de fil (la courbe solide), ainsi que l’estimateur spectral direct postcoloré correspondant (les points, ceux-ci représentent la même estimation que celle montrée à la figure 6 (b)). La fenêtre de décalage de Parzen a été utilisée ici avec une valeur de m = 237 pour le paramètre de fenêtre de lissage (le degré de liberté équivalent correspondant v est de 64). Cette valeur a été choisie après quelques expérimentations et semble produire un estimateur de fenêtre de décalage qui capture toutes les caractéristiques spectrales importantes indiquées par l’estimateur spectral direct pour les fréquences comprises entre 0,4 et 4,0 Hz (notez cependant que cet estimateur efface assez mal le pic entre 0,0 et 0,4 Hz). Nous avons également tracé un entrecroisement dont la hauteur verticale représente la longueur d’un intervalle de confiance de 95% pour 10 ⋅ log10(SX(f)) (basé sur l’estimateur de fenêtre de décalage postcolorée) et dont la largeur horizontale représente la largeur de bande de fenêtre de lissage BW

Fig.8. Estimation spectrale de fenêtre de décalage de Parzen post-couleur – courbe solide sur le graphique (a) – et estimation spectrale WOSA — courbe solide sur (b) – pour les séries chronologiques de jauge d’onde de fil. Le paramètre de fenêtre de lissage pour la fenêtre de décalage Parzen était m = 237, ce qui donne v = 64 degrés de liberté équivalents. L’estimation spectrale WOSA a été formée à l’aide d’un cône de données Hanning sur des blocs avec 256 points de données, avec des blocs adjacents se chevauchant de 50%. Les degrés de liberté équivalents pour cette estimation sont v = 59.

Estimateurs spectraux WOSA. Considérons maintenant la deuxième approche commune de la réduction de la variance, à savoir la moyenne des segments chevauchés de Welch (Welch; Carter et ses références). L’idée de base est de diviser une série temporelle en un certain nombre de blocs (i.e., segments), calculer une estimation spectrale directe pour chaque bloc, puis produire l’estimation spectrale WOSA en faisant la moyenne de ces estimations spectrales ensemble. En général, les blocs sont autorisés à se chevaucher, le degré de chevauchement étant déterminé par le degré d’effilage — plus le degré d’effilage est lourd, plus les blocs doivent se chevaucher (Thomson). Ainsi, à l’exception du tout début et de la fin de la série chronologique, les valeurs de données fortement effilées dans un bloc sont légèrement effilées dans un autre bloc, de sorte que nous récupérons intuitivement des “informations” perdues en raison de l’effilage dans un bloc à partir de blocs qui les chevauchent. Parce qu’il peut être implémenté de manière efficace sur le plan informatique (en utilisant l’algorithme de transformée de Fourier rapide) et parce qu’il peut gérer des séries temporelles très longues (ou des séries temporelles avec un spectre variant dans le temps), le schéma d’estimation WOSA est la base de nombreux analyseurs de spectre commerciaux sur le marché.

Pour définir l’estimateur spectral WOSA, soit ns représente une taille de bloc, et soit h1,…, hns soit un cône de données. Nous définissons l’estimateur spectral direct de Sx(f) pour le bloc de valeurs de données contiguës ns commençant à l’indice l comme

S^l, X(d)(f) = Δt |∑t = 1nshtXt−l−1e−l2n / τΔι | 2,1≤l≤n + 1-ns

( il n’y a aucune raison pour laquelle nous ne pouvons pas utiliser une série pré-blanchie {Wt} ici plutôt que Xt, mais le pré-blanchissement est rarement utilisé en conjonction avec WOSA, peut-être parce que le chevauchement de blocs est considéré comme un moyen efficace de compenser les degrés de liberté perdus en raison du rétrécissement). L’estimateur spectral WOSA de SX(f) est défini comme étant

(11.19) S ^X(wosa) (f) = 1nB∑j = 0nB-1S ^js + t. x(d)(f)

où nn est le nombre total de blocs et s est un facteur de décalage entier satisfaisant 0 < s≤ns et s(nB-1) = n-ns (notez que le bloc pour j = 0 utilise des valeurs de données χ1,…, Xns, tandis que le bloc pour j = nB-1 utilise Xn-ns + 1,…, Xe).

Les propriétés statistiques à grand échantillon de S^X(wosa)(f) ressemblent étroitement à celles des estimateurs de fenêtre de décalage. en particulier, nous avons l’approximation que VS ^X(wosa)(f) / Sx(f) = dXv2,, où les degrés de liberté équivalents v sont donnés par

v = 2nB1 + 2∑m = 1nB−1(1−mna) |∑t = 1nshlht + ms|2

( ici ht = 0 par définition pour tous les t > ns). Si nous nous spécialisons dans le cas d’un chevauchement de blocs à 50% (c’est-à-dire s = ns / 2) avec un cône de données de Hanning (une recommandation courante dans la littérature d’ingénierie), cela peut être approximé par la formule simple v≈36nB21 (19nB-1). Ainsi, à mesure que le nombre de blocs nB augmente, les degrés de liberté équivalents augmentent également, donnant un estimateur spectral à variance réduite. À moins que SX (⋅) n’ait un sdf relativement sans caractéristique, nous ne pouvons cependant pas rendre nB arbitrairement petit sans encourir un biais sévère dans les estimateurs spectraux directs individuels principalement en raison de la perte de résolution. (Pour plus de détails sur les résultats ci-dessus, voir la section 6.17.)

La figure 8(b) montre un estimateur spectral WOSA pour les données de jauge d’onde de fil (la courbe solide). Cette série a n = 4096 valeurs de données. Certaines expériences ont indiqué qu’une taille de bloc de ns = 256 et la conicité des données de Hanning sont des choix raisonnables pour estimer le sdf entre 0,4 et 4,0 Hz en utilisant WOSA. Avec un chevauchement de blocs de 50%, le facteur de décalage est s = ns/2 = 128 ; le nombre total de blocs est nB = 1_δ(n−ns) + 1 = 31 ; et v, les degrés de liberté équivalents, est d’environ 59. Les 31 estimations spectrales directes individuelles qui ont été moyennées ensemble pour former l’estimation WOSA sont présentées sous forme de points à la figure 8 (b).

Nous avons également tracé un entrecroisement “bande passante / intervalle de confiance” similaire à celui de la figure 8 (a), mais maintenant la “bande passante” (c’est-à-dire la largeur horizontale) est la distance en fréquence entre des estimations spectrales approximativement non corrélées. La mesure de la bande passante est fonction de la taille du bloc ns et de la conicité des données utilisées dans WOSA. Pour le cône de Hanning, la bande passante est d’environ 1,94/(nsΔt). Les croisements des figures 8(a) et 8(b) sont assez similaires, ce qui indique que les propriétés statistiques de la fenêtre de décalage de Parzen postcolorée et des estimations spectrales WOSA sont comparables: en effet, les estimations réelles concordent étroitement, l’estimation WOSA étant légèrement plus lisse en apparence.

Estimateurs spectraux Multitapers. Une alternative intéressante à l’estimation spectrale de la fenêtre de décalage ou de la WOSA est l’approche multitaper de Thomson. L’estimation spectrale multitaper peut être considérée comme un moyen de produire un estimateur spectral direct avec plus de deux degrés de liberté équivalents (les valeurs typiques sont de 4 à 16). En tant que telle, la méthode multitapeuse est différente dans l’esprit des deux autres estimateurs en ce sens qu’elle ne cherche pas à produire des spectres très lissés. Une augmentation des degrés de liberté de 2 à seulement 10 suffit cependant à réduire la largeur d’un intervalle de confiance de 95% pour le sdf de plus d’un ordre de grandeur et donc à réduire la variabilité de l’estimation spectrale au point où l’œil humain peut facilement discerner la structure globale. Des discussions détaillées sur l’approche multitaper sont données dans et chapitre 7 de. Ici, nous esquissons simplement les idées principales.

L’estimation spectrale multitaper est basée sur l’utilisation d’un ensemble de K cônes de données {ht.k ; t = 1, …, n}, où k va de 0 à K-1. Nous supposons que ces cônes sont orthonormés (c’est-à-dire ∑t = 1nht, jht, k = 1 si j = k et 0 si j≠k). L’estimateur multitapeur le plus simple est défini par

S^X(mt)(f) = 1K∑K = 0K-1S ^K, X(mt)(f) avec ^k, x(mt)(f) Δt /∑t = 1nht, KXte-i2πftΔι|2

( Thomson préconise une pondération adaptative des S^k, X(mt)(f) plutôt que de simplement les faire la moyenne ensemble). Une comparaison de cette définition pour S^k, X(mt)(⋅) avec Eq. (118) montre que S^k, X(mt)(⋅) n’est en fait qu’un estimateur spectral direct, de sorte que l’estimateur multitaper n’est qu’une moyenne d’estimateurs spectraux directs utilisant un ensemble orthonormé de cierges. Dans certaines conditions douces, l’orthonormalité des cierges se traduit dans le domaine fréquentiel par l’indépendance approximative de chaque individu S^k, X(mt)(f); c’est-à-dire S^j.X(mt)(f). L’indépendance approximative implique à son tour que 2KS ^ k, X(mt)(f) / SX(f) = dX22k approximativement, de sorte que les degrés de liberté équivalents pour S^X(mt)(f) sont égaux au double du nombre de contres de données utilisés.

L’astuce clé consiste alors à trouver un ensemble de K séquences orthonormées, dont chacune fait un travail approprié d’effilage. Une approche attrayante consiste à recalculer le problème de concentration qui nous a donné le cône dpss pour une bande passante à résolution fixe de 2W Si nous appelons maintenant ce cône le cône dpss d’ordre zéro et le dénotons par {h,,()}, nous pouvons construire récursivement les K-1 cierges dpss “d’ordre supérieur” restants {ht, k} comme suit. Pour k = 1 , …, K-1, nous définissons la conicité dpss d’ordre k comme l’ensemble de n nombres {ht, k; t = 1,…, n} tels que

1

{ht, k} est orthogonal à chacune des k séquences {ht,()},…, {ht,(k−1)} c’est-à-dire ∑t = 11ht.Jht.k = 0 pour j = 0, …, d-1);

2

{ ht, k} est normalisé de telle sorte quett = 1nht, k2=1;

3

sous réserve des conditions] et 2, la fenêtre spectrale Hk(⋅) correspondant à {ht.k} maximise le rapport de concentration

∫-wwHk(f) df /∫−f(N) f(n) Hk(f) df = λk(n, W)

En termes, sous réserve de la contrainte d’être orthogonale à toutes les coniques dpss d’ordre inférieur, la conicité dpss d’ordre k est ” optimale ” en ce sens restreint que les lobes latéraux de sa fenêtre spectrale sont supprimés autant que possible que mesuré par le rapport de concentration. Les méthodes de calcul des conjectures des données dpss sont discutées au chapitre 8.

Dans une série d’articles, Slepian (et ses références) a largement étudié la nature des dpss. Un fait important qu’il discute est que le rapport de concentration λk(n, W) diminue strictement à mesure que k augmente de telle sorte que λk(n, W) est proche de l’unité pour k < 2NW Δt, après quoi il se rapproche rapidement de 0 avec l’augmentation de k (la valeur 2nWΔt est parfois appelée le nombre de Shannon). Étant donné que λk(n, W) doit être proche de l’unité pour que {ht,k} soit une conicité de données décente, l’estimation spectrale multitaper est limitée à l’utilisation d’au plus – et, en pratique, généralement moins de — 2nWΔt de conicité dpss orthonormale.

Un exemple d’estimation spectrale multitaper est illustré à la figure 9. La colonne de gauche des graphiques montre les diminutions de données dpss d’ordre k pour n = 4096, nW = 4 / Δt et k allant de 0 (graphique du haut) à K-1 = 5 (graphique du bas). Les fines lignes horizontales de chacun de ces graphiques indiquent le niveau zéro, donc, alors que le dpss d’ordre zéro est strictement positif partout (mais assez proche de 0 près de t = 1 et t = n), les cierges d’ordre supérieur supposent des valeurs positives et négatives. Notez également que la conicité d’ordre zéro pèse fortement sur les valeurs des séries chronologiques proches de t = 1 et t = n, mais que ces valeurs reçoivent successivement plus de poids par les coniques d’ordre supérieur (une interprétation du multitapering est que les coniques d’ordre supérieur récupèrent des informations “perdues” lorsqu’une seule conicité de données est utilisée). La courbe solide de la figure 9 (b) montre une estimation spectrale multitaper S^X (mt) (⋅) pour les données de jauge d’onde de fil basées sur ces 6 cierges dpss, tandis que les points montrent les six estimations spectrales directes individuelles S^K.X (mt) (⋅). On notera que le nombre de cierges que nous avons utilisés est inférieur au nombre de Shannon 2nWΔt = 8 et que v, les degrés de liberté équivalents, est ici 2K = 12. L’estimation spectrale multitaper est beaucoup plus hachée en apparence que l’estimation spectrale de fenêtre de décalage de la figure 8 (a) ou l’estimation WOSA de la figure 8 (b), qui ont toutes deux un nombre nettement plus élevé de degrés de liberté équivalents (v = 64 et v = 59, respectivement). Néanmoins, la variabilité de l’estimation spectrale multitapier est suffisamment faible pour que l’œil puisse facilement détecter la structure globale (cf. S^X (mt) (⋅) avec les deux estimations spectrales de la figure 5), et parce qu’elle n’est pas très lissée, l’estimation multitapeuse réussit nettement mieux à capturer la structure spectrale proche de f = 0.

Fig.9. Estimation spectrale multitaper

Sur la base des limites de performance, Bronez [16 soutient que l’estimateur spectral multitapeur a des propriétés statistiques supérieures à WOSA pour les FDS à très hautes plages dynamiques (des recherches supplémentaires sont cependant nécessaires pour vérifier que ces limites se traduisent par un avantage réel dans la pratique). Par rapport au pré-blanchissement, le multitapering est utile dans les situations où les fuites sont un concem, mais il n’est pas pratique de concevoir soigneusement des filtres de pré-blanchissement (cela se produit, par exemple, dans la géophysique d’exploration en raison du volume énorme de séries chronologiques collectées régulièrement). Enfin, nous notons que Thomson et Chave [17 décrivent un schéma attrayant dans lequel le multitapering est utilisé conjointement avec WOSA.

Laisser un commentaire

Votre adresse e-mail ne sera pas publiée.