⛳️ 👩🏿‍🏭 💩 Analyse des ondelettes. 2e partie 🖍️ 😎 🤐

Présentation

Cette publication traite de l'analyse en ondelettes des séries chronologiques. L'idée principale de la transformée en ondelettes correspond aux spécificités de nombreuses séries chronologiques qui démontrent l'évolution dans le temps de leurs principales caractéristiques - la valeur moyenne, la variance, les périodes, les amplitudes et les phases des composantes harmoniques. La grande majorité des processus étudiés dans divers domaines de la connaissance présentent les caractéristiques ci-dessus.

Le but de cette publication est de décrire la méthode de transformation continue en ondelettes de séries chronologiques à l'aide de la bibliothèque PyWavelets.

Un peu d'histoire

Géophysicien D. Morlaix à la fin des années 70 du XXe siècle. J'ai rencontré le problème de l'analyse des signaux des capteurs sismiques qui contenaient une composante haute fréquence (activité sismique) pendant une courte période et des composantes basse fréquence (état calme de la croûte terrestre) pendant une longue période. La transformée de Fourier de fenêtre vous permet d'analyser la composante haute fréquence ou la composante basse fréquence, mais pas les deux à la fois.

Par conséquent, une méthode d'analyse a été proposée dans laquelle la largeur de la fonction de fenêtre pour les basses fréquences augmentait et pour les hautes fréquences, elle diminuait. La nouvelle transformation de fenêtre a été obtenue à la suite de l'extension (compression) et du décalage temporel d'une fonction génératrice (dite fonction de mise à l'échelle - fonction de mise à l'échelle, scalet). Cette fonction génératrice a été appelée l'ondelette par D. Morlaix.

Wavelet D. Morlaix

from pylab import* import scaleogram as scg axes = scg.plot_wav('cmor1-1.5', figsize=(14,3)) show()

Transformation en ondelettes continue (CWT)

Transformation en ondelettes à un niveau:

coefs, fréquences = pywt.cwt (données, échelles, ondelettes)

où:

données: (array_like) -Signal d'entrée ;

échelles (array_like) : - L'échelle d'ondelettes à utiliser. Vous pouvez utiliser le rapport f = scale2frequency (scale, wavelet) / d'échantillonnage_période et déterminer quelle fréquence physique f est. Ici, f en hertz lors de la période d'échantillonnage en secondes;

ondelette (objet ou nom ondelette) : - ondelette à utiliser;

échantillonnage_période (float): - Période d'échantillonnage pour les fréquences de sortie (paramètre facultatif). Les valeurs calculées pour les coefs sont indépendantes de la période d'échantillonnage (c'est-à-dire que les échelles ne sont pas mises à l'échelle par échantillonnage de période.);

coefs (array_like) : - Ondelette continue - transformation du signal d'entrée pour une échelle et une ondelette données;

fréquences (array_like) : - Si l'unité de la période d'échantillonnage est la seconde et est définie, les fréquences sont affichées en hertz. Sinon, une période d'échantillonnage de 1 est supposée.

Remarque: La taille des tableaux de coefficients dépend de la longueur du tableau d'entrée et des longueurs des échelles données.

Obtenez une liste des noms d'ondelettes compatibles cwt disponibles:

 >>> import pywt >>> wavlist = pywt.wavelist(kind='continuous') >>> wavlist

['cgau1', 'cgau2', 'cgau3', 'cgau4', 'cgau5', 'cgau6', 'cgau7', 'cgau8', 'cmor', 'fbsp', 'gaus1', 'gaus2', ' gaus3 ',' gaus4 ',' gaus5 ',' gaus6 ',' gaus7 ',' gaus8 ',' mexh ',' morl ',' shan ']

Pour l'analyse en ondelettes, le choix de l'ondelette mère est l'une des tâches clés. Par conséquent, avant chaque sélection d'ondelettes, il est simplement nécessaire d'utiliser le programme suivant qui vous permet de déterminer les propriétés dynamiques de l'ondelette:

Listing 1

 import numpy as np import pywt import matplotlib.pyplot as plt vav='cmor1.5-1.0' wav = pywt.ContinuousWavelet(vav) #  ,      print("       [{}, {}]".format( wav.lower_bound, wav.upper_bound)) width = wav.upper_bound - wav.lower_bound scales = [1, 2, 3, 4, 10, 15] max_len = int(np.max(scales)*width + 1) t = np.arange(max_len) fig, axes = plt.subplots(len(scales), 2, figsize=(12, 6)) for n, scale in enumerate(scales): #       cwt int_psi, x = pywt.integrate_wavelet(wav, precision=10) step = x[1] - x[0] j = np.floor( np.arange(scale * width + 1) / (scale * step)) if np.max(j) >= np.size(int_psi): j = np.delete(j, np.where((j >= np.size(int_psi)))[0]) j = j.astype(np.int) # normalize int_psi     int_psi /= np.abs(int_psi).max() #     filt = int_psi[j][::-1] # CWT          #          . nt = len(filt) t = np.linspace(-nt//2, nt//2, nt) axes[n, 0].plot(t, filt.real, t, filt.imag) axes[n, 0].set_xlim([-max_len//2, max_len//2]) axes[n, 0].set_ylim([-1, 1]) axes[n, 0].text(50, 0.35, 'scale = {}'.format(scale)) f = np.linspace(-np.pi, np.pi, max_len) filt_fft = np.fft.fftshift(np.fft.fft(filt, n=max_len)) filt_fft /= np.abs(filt_fft).max() axes[n, 1].plot(f, np.abs(filt_fft)**2) axes[n, 1].set_xlim([-np.pi, np.pi]) axes[n, 1].set_ylim([0, 1]) axes[n, 1].set_xticks([-np.pi, 0, np.pi]) axes[n, 1].set_xticklabels([r'$-\pi$', '0', r'$\pi$']) axes[n, 1].grid(True, axis='x') axes[n, 1].text(np.pi/2, 0.5, 'scale = {}'.format(scale)) axes[n, 0].set_xlabel('time (samples)') axes[n, 1].set_xlabel('frequency (radians)') axes[0, 0].legend(['real', 'imaginary'], loc='upper left') axes[0, 1].legend(['Power'], loc='upper left') axes[0, 0].set_title('filter: wavelet - %s'%vav) axes[0, 1].set_title(r'|FFT(filter)|$^2$') plt.show()

Ensuite, nous considérerons les fonctions d'ondelettes et leurs propriétés en utilisant le programme ci-dessus:

Chapeau mexicain wavelet mexh:

v a r p h i l e f t (t r i g h t) = f r a c 2 s q r t 3 s q r t [4] p i c d o t e x p^{- f r a c t^{2} 2} c d o t l e f t (1 - t^{2} r i g h t)

$\ varphi \ left (t \ right) = \ frac {2} {\ sqrt {3} \ sqrt [4] {\ pi}} \ cdot exp ^ {- \ frac {t ^ {2}} {2} } \ cdot \ left (1-t ^ {2} \ right)$

Ondelette "morl" Morlaix:

v a r p h i l e f t (t r i g h t) = e x p^{- f r a c t^{2} 2} c d o t c o s (5 t)

$\ varphi \ left (t \ right) = exp ^ {- \ frac {t ^ {2}} {2}} \ cdot cos (5t)$

Ondelette de Morlet complexe "cmorB-C"

v a r p h i l e f t (t r i g h t) = f r a c 1 s q r t p i B e x p^{- f r a c t^{2} B} c d o t e x p^{j 2 p i C t}

$\ varphi \ left (t \ right) = \ frac {1} {\ sqrt {\ pi} B} exp ^ {- \ frac {t ^ {2}} {B}} \ cdot exp ^ {j2 \ pi C t}$

où B est la bande passante; C est la fréquence centrale.

Ci-après (sans indication particulière) les valeurs de B, C sont fixées avec une virgule flottante.

Ondelettes gaussiennes «gausP»

v a r p h i l e f t (t r i g h t) = C c d o t e x p^{- t^{2}}

$\ varphi \ left (t \ right) = C \ cdot exp ^ {- t ^ {2}}$

où: P - un entier de 1 à 8 est inséré dans le nom de l'ondelette, par exemple, "gaus8"; C- constante de normalisation intégrée.

Ondelettes gaussiennes intégrées «cgauP»

v a r p h i l e f t (t r i g h t) = C c d o t e x p^{- t^{2}} c d o t e x p^{- j t}

$\ varphi \ left (t \ right) = C \ cdot exp ^ {- t ^ {2}} \ cdot exp ^ {- jt}$

où: P - un entier de 1 à 8 est inséré dans le nom de l'ondelette, par exemple, "gaus8" et correspond à l'ordre de la dérivée de la fonction d'ondelette; C- constante de normalisation intégrée.

Ondelettes Shannon "shanB-C"

v a r p h i l e f t (t r i g h t) = s q r t B c d o t f r a c s i n (p i B t) p i B T c d o t e x p^{j 2 p i C t}

$\ varphi \ left (t \ right) = \ sqrt {B} \ cdot \ frac {sin (\ pi B t)} {\ pi B T} \ cdot exp ^ {j 2 pi C t}$

où: B est la bande passante; C est la fréquence centrale.

CWT dans PyWavelets s'applique aux données discrètes - convolutions avec des échantillons de l'intégrale d'ondelettes. Si l'échelle est trop petite, alors nous obtenons un filtre discret avec un échantillonnage inadéquat, conduisant à un lissage, comme le montre le graphique pour l'ondelette 'cmor1.5-1.0'.

Dans la colonne de gauche, les graphiques montrent des filtres discrets utilisés dans la convolution pour différentes échelles. La colonne de droite est le spectre de puissance de Fourier correspondant de chaque filtre. Avec les échelles 1 et 2 pour le graphique «cmor1,5-1,0», on peut voir que le lissage se produit en raison de la violation de la contrainte de Nyquist. Le phénomène indiqué peut se produire dans d'autres ondelettes lors du choix de C et B, par conséquent, lorsque vous travaillez avec des ondelettes, vous devez utiliser le programme - Listing 1.

Afin de relier une échelle donnée à une fréquence de signal spécifique, la période d'échantillonnage du signal doit être connue. En utilisant la fonction pywt.scale2frequency () , vous pouvez convertir la liste des échelles aux fréquences correspondantes. Le choix correct des échelles dépend de l'ondelette sélectionnée, donc pywt.scale2frequency () doit être utilisé pour se faire une idée de la plage de fréquence correspondante du signal.

 >>> import pywt >>> dt = 0.01 # 100 Hz sampling >>> frequencies = pywt.scale2frequency('cmor1.5-1.0', [1, 2, 3, 4]) / dt >>> frequencies

tableau ([100., 50., 33.33333333, 25.])

Ondelettes - Analyse de séries chronologiques à l'aide de CWT PyWavelets

L'ensemble de données el-Nino est un ensemble de données de séries chronologiques utilisé pour suivre El Nino et contient des mesures trimestrielles de la température de la surface de la mer de 1871 à 1997. Pour comprendre la puissance d'un graphique à l'échelle, visualisons-le pour un ensemble de données el-Nino avec les données de série temporelle d'origine et sa transformée de Fourier.

Pour l'analyse en ondelettes de séries chronologiques, les étapes suivantes doivent être effectuées:

1. Sélectionnez l'ondelette mère: Sélectionnez l'ondelette Morlet complexe "cmorB-C":

v a r p h i l e f t (t r i g h t) = f r a c 1 s q r t p i B e x p^{- f r a c t^{2} B} c d o t e x p^{j 2 p i C t}

$\ varphi \ left (t \ right) = \ frac {1} {\ sqrt {\ pi} B} exp ^ {- \ frac {t ^ {2}} {B}} \ cdot exp ^ {j2 \ pi C t}$

Bande passante - B et fréquence centrale - à partir de laquelle nous déterminerons à l'étape suivante.

2. Déterminez la fréquence centrale en prenant dt = 0,25 pour notre série chronologique:

 import pywt dt = 0.25 # 4 Hz sampling scale=range(1,4) frequencies = pywt.scale2frequency('cmor1.0-0.5', scale) / dt print(frequencies)

[2. 1. 0.66666667]

3. On trouve la transformée de Fourier de l'ondelette mère cmor1.0-0.5 (on utilise le listing 1):
Les ondelettes continues seront évaluées sur toute la plage [-8,0, 8,0]

4. Sélectionnez les données pour la série chronologique:

 pd.read_csv("http://paos.colorado.edu/research/wavelets/wave_idl/sst_nino3.dat", sep = "|")

Les données sont des mesures trimestrielles de la température de la surface de la mer de 1871 à 1997. Pour ces données: t0 = 1871 dt = 0,25

5. Nous construisons une série chronologique avec un signal et sa moyenne mobile sur un graphique:

Listing 2

 from numpy import * from scipy import * import pandas as pd from pylab import * import pywt def get_ave_values(xvalues, yvalues, n = 5): signal_length = len(xvalues) if signal_length % n == 0: padding_length = 0 else: padding_length = n - signal_length//n % n xarr = array(xvalues) yarr = array(yvalues) xarr.resize(signal_length//n, n) yarr.resize(signal_length//n, n) xarr_reshaped = xarr.reshape((-1,n)) yarr_reshaped = yarr.reshape((-1,n)) x_ave = xarr_reshaped[:,0] y_ave = nanmean(yarr_reshaped, axis=1) return x_ave, y_ave def plot_signal_plus_average(time, signal, average_over = 5): fig, ax = subplots(figsize=(15, 3)) time_ave, signal_ave = get_ave_values(time, signal, average_over) ax.plot(time, signal, label='') ax.plot(time_ave, signal_ave, label = '   (n={})'.format(5)) ax.set_xlim([time[0], time[-1]]) ax.set_ylabel(' ', fontsize=18) ax.set_title(' +   ', fontsize=18) ax.set_xlabel('', fontsize=18) ax.legend() show() df_nino = pd.read_csv("http://paos.colorado.edu/research/wavelets/wave_idl/sst_nino3.dat", sep = "|") N = df_nino.shape[0] t0=1871 dt=0.25 time = arange(0, N) * dt + t0 signal = df_nino.values.squeeze() scales = arange(1, 128) plot_signal_plus_average(time, signal)

6. Nous réalisons la transformée de Fourier et la modalité de spectre à partir de la série chronologique:

Listing 3

 from numpy import * from scipy import * import pandas as pd from pylab import * import pywt def get_ave_values(xvalues, yvalues, n = 5): signal_length = len(xvalues) if signal_length % n == 0: padding_length = 0 else: padding_length = n - signal_length//n % n xarr = array(xvalues) yarr = array(yvalues) xarr.resize(signal_length//n, n) yarr.resize(signal_length//n, n) xarr_reshaped = xarr.reshape((-1,n)) yarr_reshaped = yarr.reshape((-1,n)) x_ave = xarr_reshaped[:,0] y_ave = nanmean(yarr_reshaped, axis=1) return x_ave, y_ave def get_fft_values(y_values, T, N, f_s): f_values = linspace(0.0, 1.0/(2.0*T), N//2) fft_values_ = fft(y_values) fft_values = 2.0/N * abs(fft_values_[0:N//2]) return f_values, fft_values def plot_fft_plus_power(time, signal): dt = time[1] - time[0] N = len(signal) fs = 1/dt fig, ax = subplots(figsize=(15, 3)) variance = std(signal)**2 f_values, fft_values = get_fft_values(signal, dt, N, fs) fft_power = variance * abs(fft_values) ** 2 # FFT power spectrum ax.plot(f_values, fft_values, 'r-', label='FFT ') ax.plot(f_values, fft_power, 'k--', linewidth=1, label=' ') ax.set_xlabel('[ / ]', fontsize=18) ax.set_ylabel('', fontsize=18) ax.legend() show() df_nino = pd.read_csv("http://paos.colorado.edu/research/wavelets/wave_idl/sst_nino3.dat", sep = "|") N = df_nino.shape[0] t0=1871 dt=0.25 time = arange(0, N) * dt + t0 signal = df_nino.values.squeeze() scales = arange(1, 128) plot_fft_plus_power(time, signal)

7. Nous déterminons les échelles: échelles = arange (1, 128); niveaux = [2 ** - 4, 2 ** - 3, 2 ** - 2, 2 ** - 1, 2 ** 0, 2 ** 1, 2 ** 2, 2 ** 3].

8. Créez un graphique à l'échelle:

Listing 4

 from numpy import * import pandas as pd from pylab import * import pywt def plot_wavelet(time, signal, scales, waveletname = 'cmor1.0-0.4', cmap = plt.cm.seismic, title = '-( ) ', ylabel = ' ()', xlabel = ''): dt = time[1] - time[0] [coefficients, frequencies] = pywt.cwt(signal, scales, waveletname, dt) power = (abs(coefficients)) ** 2 period = 1. / frequencies levels = [2**-4 , 2**-3 , 2**-2 , 2**-1 , 2**0 , 2**1 , 2**2 , 2**3] contourlevels = log2(levels) fig, ax = subplots(figsize=(15, 10)) im = ax.contourf(time, log2(period), log2(power), contourlevels, extend='both',cmap=cmap) ax.set_title(title, fontsize=20) ax.set_ylabel(ylabel, fontsize=18) ax.set_xlabel(xlabel, fontsize=18) yticks = 2**arange(np.ceil(log2(period.min())), ceil(log2(period.max()))) ax.set_yticks(log2(yticks)) ax.set_yticklabels(yticks) ax.invert_yaxis() ylim = ax.get_ylim() ax.set_ylim(ylim[0], -1) cbar_ax = fig.add_axes([0.95, 0.5, 0.03, 0.25]) fig.colorbar(im, cax=cbar_ax, orientation="vertical") show() df_nino = pd.read_csv("http://paos.colorado.edu/research/wavelets/wave_idl/sst_nino3.dat", sep = "|") N = df_nino.shape[0] t0=1871 dt=0.25 time = arange(0, N) * dt + t0 signal = df_nino.values.squeeze() scales = arange(1, 128) plot_wavelet(time, signal, scales)

Le graphique à l'échelle montre que la majeure partie de la puissance du spectre est concentrée sur une période de 2 à 8 ans, ce qui correspond à une fréquence de 0,125 à 0,5 Hz. Dans le spectre de Fourier, la mode se concentre également autour de ces valeurs de fréquence. Cependant, la transformée en ondelettes nous donne également des informations temporaires, mais pas la transformée de Fourier.

Par exemple, sur le diagramme à l'échelle, nous voyons qu'avant 1920, il y avait de nombreuses fluctuations, alors qu'entre 1960 et 1990, il n'y en avait pas tant. Nous pouvons également voir qu'au fil du temps, il y a un passage de périodes plus courtes à des périodes plus longues.

Utilisation du module scalogramme

Scaleogram est un outil pratique pour analyser les données 1D avec une transformation en ondelettes continue basée sur la bibliothèque PyWavelets. Le module est conçu pour fournir un outil fiable pour une analyse rapide des données.

Le Scaleogram présente les caractéristiques suivantes:

signature d'appel simple pour débutant
axes lisibles et intégration pure matplotlib
de nombreuses options de zoom, filtre de spectre, implémentation de la barre de couleur, etc ...
prise en charge de la fréquence et des unités de fréquence conformément au marquage
vitesse, utilise l'algorithme N * log (N) pour les transformations
portabilité: testé avec python2.7 et python3.7
messages d'erreur détaillés et documentation avec des exemples

Cependant, dans les exemples d'utilisation, l'accès aux données n'est pas correctement enregistré:

 nino3 = "http://paos.colorado.edu/research/wavelets/wave_idl/sst_nino3.dat"

L'avertissement suivant apparaît:

nino3 = pd.read_table (nino3)
conduisant à un avertissement:
Avertissement (du module des avertissements):
Fichier «C: /Users/User/Desktop/wavelet/pr.py», ligne 7
nino3 = pd.read_table (nino3)
FutureWarning: read_table est déconseillé, utilisez plutôt read_csv, en passant sep = '\ t'

L'accès aux données doit être écrit comme ceci:

 url="http://paos.colorado.edu/research/wavelets/wave_idl/sst_nino3.dat" nino3 = pd.read_csv(url, sep = "|")

Pour montrer l'utilisation du scalogramme et comparer les résultats avec les résultats de l'exemple précédent, nous utilisons les mêmes données - sur les mesures trimestrielles de la température de la surface de la mer de 1871 à 1997. Pour ces données: t0 = 1871 dt = 0,25

Listing 4

 import pandas as pd import pywt from numpy import* import scaleogram as scg from pylab import* url="sst_nino3.dat" nino3 = pd.read_csv(url, sep = "|") data = nino3.values.squeeze() N = data.size t0 = 1871; dt = 0.25 time = t0 + arange(len(data))*dt wavelet = 'cmor1-0.5' ax = scg.cws(time, data, scales=arange(1, 128), wavelet=wavelet, figsize=(14, 3), cmap="jet", cbar=None, ylabel=' ()', xlabel=" []", yscale="log", title='-  \n( )') ticks = ax.set_yticks([2,4,8, 16,32]) ticks = ax.set_yticklabels([2,4,8, 16,32]) show()

Si nous comparons le graphique à l'échelle avec le scan résultant, alors, à l'exception de la palette de couleurs, ils sont identiques et montrent donc la même dynamique dans la série chronologique.
Le listing 4 est plus simple que le listing 3 et offre de meilleures performances.

Lorsque le spectre de fréquence et le spectre de puissance selon Fourier ne permettent pas d'obtenir des informations supplémentaires sur la dynamique de la série temporelle, alors l'utilisation d'un scalogramme est préférée.

Conclusions

Un exemple d'analyse en ondelettes d'une série chronologique au moyen de la bibliothèque PyWavelets est donné.

Analyse des ondelettes. 2e partie

Présentation

Transformation en ondelettes continue (CWT)

Ondelettes - Analyse de séries chronologiques à l'aide de CWT PyWavelets

Utilisation du module scalogramme

Conclusions

More articles: