™️ 🧓 ⚜️ Análisis de wavelet. Parte 2 👨‍🚀 🏣 ☝🏽

Introduccion

Esta publicación discute el análisis wavelet de series de tiempo. La idea principal de la transformada wavelet corresponde a los detalles de muchas series de tiempo que demuestran la evolución en el tiempo de sus características principales: el valor promedio, la varianza, los períodos, las amplitudes y las fases de los componentes armónicos. La gran mayoría de los procesos estudiados en diversos campos del conocimiento tienen las características anteriores.

El propósito de esta publicación es describir el método de transformación de wavelet continua de series de tiempo usando la biblioteca PyWavelets.

Un poco de historia

Geofísico D. Morlaix a finales de los años 70 del siglo XX. Encontré el problema de analizar señales de sensores sísmicos que contenían un componente de alta frecuencia (actividad sísmica) durante un corto período de tiempo y componentes de baja frecuencia (estado de calma de la corteza terrestre) durante un largo período. La transformación Window Fourier le permite analizar el componente de alta frecuencia o el componente de baja frecuencia, pero no ambos componentes a la vez.

Por lo tanto, se propuso un método de análisis en el que el ancho de la función de ventana para frecuencias bajas aumentó, y para frecuencias altas disminuyó. La nueva transformación de ventana se obtuvo como resultado del estiramiento (compresión) y el desplazamiento de tiempo de una función generadora (llamada función de escala - función de escala, escala). D. Morlaix llamó a esta función generadora la wavelet.

Wavelet D. Morlaix

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

Transformación continua de wavelet (CWT)

Transformada wavelet de un solo nivel:

coefs, frecuencias = pywt.cwt (datos, escalas, wavelet)

donde:

datos: (matriz_) señal de entrada;

escalas (array_like) : - La escala wavelet a usar. Puede usar la relación f = scale2frequency (scale, wavelet) / samples_period y determinar qué frecuencia física es f. Aquí, f en hertz cuando sample_period en segundos;

wavelet (objeto o nombre Wavelet ) : - Wavelet para usar;

sample_period (float): - Período de muestreo para frecuencias de salida (parámetro opcional). Los valores calculados para los coeficientes son independientes del período de muestreo (es decir, las escalas no se escalan por muestreo de período);

coefs (array_like) : - Wavelet continua: transformación de la señal de entrada para una escala y wavelet determinadas;

frecuencias (tipo matriz) : - Si la unidad del período de muestreo es la segunda y está configurada, las frecuencias se muestran en hercios. De lo contrario, se supone un período de muestreo de 1.

Nota: El tamaño de las matrices de coeficientes depende de la longitud de la matriz de entrada y las longitudes de las escalas dadas.

Obtenga una lista de nombres de wavelet compatibles con 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 ']

Para el análisis wavelet, la elección de la wavelet madre es una de las tareas clave. Por lo tanto, antes de cada selección de wavelet, simplemente es necesario usar el siguiente programa que le permite determinar las propiedades dinámicas de la wavelet:

Listado 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()

A continuación, consideraremos las funciones wavelet y sus propiedades utilizando el programa anterior:

Sombrero mexicano 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)$

Wavelet "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)$

Onda Morlet compleja "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}$

donde B es el ancho de banda; C es la frecuencia central.

En lo sucesivo (sin indicación especial) los valores de B, C se establecen con un punto flotante.

Ondas gaussianas "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}}$

donde: P - se inserta un número entero de 1 a 8 en el nombre de la wavelet, por ejemplo, "gaus8"; C- constante de normalización incorporada.

Ondas Gaussianas integradas "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}$

donde: P - se inserta un número entero de 1 a 8 en el nombre de la wavelet, por ejemplo, "gaga8" y corresponde al orden de la derivada de la función wavelet; C- constante de normalización incorporada.

Shannon wavelets "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}$

donde: B es el ancho de banda; C es la frecuencia central.

CWT en PyWavelets se aplica a datos discretos: convoluciones con muestras de la integral wavelet. Si la escala es demasiado pequeña, obtenemos un filtro discreto con un muestreo inadecuado, lo que conduce a un suavizado, como se muestra en el gráfico de la wavelet 'cmor1.5-1.0'.

En la columna izquierda, los gráficos muestran filtros discretos utilizados en la convolución para varias escalas. La columna derecha es el espectro de potencia de Fourier correspondiente de cada filtro. Con las escalas 1 y 2 para el gráfico 'cmor1.5-1.0' se puede ver que el suavizado ocurre debido a la violación de la restricción Nyquist. El fenómeno indicado puede ocurrir en otras wavelets al elegir C y B, por lo tanto, al trabajar con wavelets, debe usar el programa - Listado 1.

Para relacionar una escala dada con una frecuencia de señal específica, se debe conocer el período de muestreo de la señal. Usando la función pywt.scale2frequency () , puede convertir la lista de escalas a las frecuencias correspondientes. La elección correcta de las escalas depende de la wavelet seleccionada, por lo que se debe utilizar pywt.scale2frequency () para tener una idea del rango de frecuencia correspondiente de la señal.

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

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

Wavelet - Análisis de series temporales usando PyWavelets CWT

El conjunto de datos el-Nino es un conjunto de datos de series de tiempo utilizado para rastrear El Niño y contiene mediciones trimestrales de la temperatura de la superficie del mar desde 1871 hasta 1997. Para comprender el poder de un gráfico de escala, visualicémoslo para un conjunto de datos el-Nino junto con los datos de la serie de tiempo original y su transformación de Fourier.

Para el análisis wavelet de series de tiempo, se deben realizar los siguientes pasos:

1. Seleccione la wavelet madre: Seleccione la wavelet Morlet compleja "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}$

Ancho de banda - B y frecuencia central - De lo cual determinaremos en el siguiente paso.

2. Determine la frecuencia central tomando dt = 0.25 para nuestra serie de tiempo:

 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. Encontramos la transformada de Fourier de la wavelet madre cmor1.0-0.5 (utilizamos el listado 1):
La wavelet continua se evaluará en todo el rango [-8.0, 8.0]

4. Seleccione los datos para la serie de tiempo:

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

Los datos son mediciones trimestrales de la temperatura de la superficie del mar desde 1871 hasta 1997. Para estos datos: t0 = 1871 dt = 0.25

5. Construimos una serie de tiempo con una señal y su promedio móvil en un gráfico:

Listado 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. Llevamos a cabo la transformada de Fourier y la modalidad de espectro de la serie temporal:

Listado 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. Determinamos las escalas: escalas = arange (1, 128); niveles = [2 ** - 4, 2 ** - 3, 2 ** - 2, 2 ** - 1, 2 ** 0, 2 ** 1, 2 ** 2, 2 ** 3].

8. Construya un gráfico de escala:

Listado 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)

El gráfico de escala muestra que la mayor parte del espectro de potencia se concentra en un período de 2 a 8 años, esto corresponde a una frecuencia de 0,125 a 0,5 Hz. En el espectro de Fourier, la moda también se concentra alrededor de estos valores de frecuencia. Sin embargo, la transformada wavelet también nos da información temporal, pero la transformada de Fourier no.

Por ejemplo, en el diagrama de escala, vemos que antes de 1920 había muchas fluctuaciones, mientras que entre 1960 y 1990 no había tantas. También podemos ver que con el tiempo hay un cambio de períodos más cortos a más largos.

Usando el módulo de escalagrama

Scaleogram es una herramienta conveniente para analizar datos 1D con transformación de wavelet continua basada en la biblioteca PyWavelets. El módulo está diseñado para proporcionar una herramienta confiable para el análisis rápido de datos.

Scaleogram tiene las siguientes características:

firma de llamada simple para principiantes
ejes legibles e integración pura matplotlib
muchas opciones para zoom, filtro de espectro, implementación de barra de colores, etc.
soporte para frecuencia y unidades de frecuencia de acuerdo con la marca
velocidad, utiliza el algoritmo N * log (N) para transformaciones
portabilidad: probado con python2.7 y python3.7
mensajes de error detallados y documentación con ejemplos

Sin embargo, en los ejemplos de uso, el acceso a los datos no se registra correctamente:

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

Aparece la siguiente advertencia:

nino3 = pd.read_table (nino3)
que conduce a una advertencia:
Advertencia (del módulo de advertencias):
Archivo "C: /Users/User/Desktop/wavelet/pr.py", línea 7
nino3 = pd.read_table (nino3)
FutureWarning: read_table está en desuso, use read_csv en su lugar, pasando sep = '\ t'

El acceso a los datos debe escribirse así:

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

Para mostrar el uso del escalograma y comparar los resultados con los resultados del ejemplo anterior, utilizamos los mismos datos, en mediciones trimestrales de la temperatura de la superficie del mar desde 1871 hasta 1997. Para estos datos: t0 = 1871 dt = 0.25

Listado 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 comparamos el gráfico de escala con el escaneo resultante, entonces, con la excepción de la paleta de colores, son idénticos y, por lo tanto, muestran la misma dinámica en las series de tiempo.
El listado 4 es más simple que el listado 3 y tiene un mejor rendimiento.

Cuando el espectro de frecuencia y el espectro de potencia según Fourier no permiten obtener información adicional sobre la dinámica de las series de tiempo, entonces es preferible el uso de un escalograma.

Conclusiones

Se proporciona un ejemplo de análisis wavelet de una serie de tiempo mediante la biblioteca PyWavelets.

Análisis de wavelet. Parte 2

Introduccion

Transformación continua de wavelet (CWT)

Wavelet - Análisis de series temporales usando PyWavelets CWT

Usando el módulo de escalagrama

Conclusiones

More articles: