👔 👨🏻‍🚒 👩🏼‍⚕️ Análise wavelet. Parte 2 👩🏾‍🎓 🥉 👴🏿

1. Introdução

Esta publicação discute a análise de séries temporais de wavelet. A idéia principal da transformada wavelet corresponde às especificidades de muitas séries temporais que demonstram a evolução no tempo de suas principais características - valor médio, variância, períodos, amplitudes e fases dos componentes harmônicos. A grande maioria dos processos estudados em vários campos do conhecimento tem as características acima.

O objetivo desta publicação é descrever o método de transformação de wavelet contínua de séries temporais usando a biblioteca PyWavelets.

Um pouco de história

Geofísico D. Morlaix no final dos anos 70 do século XX. Encontrei o problema de analisar sinais de sensores sísmicos que continham um componente de alta frequência (atividade sísmica) por um curto período de tempo e componentes de baixa frequência (estado calmo da crosta terrestre) por um longo período. A transformação de Fourier da janela permite analisar o componente de alta frequência ou o componente de baixa frequência, mas não os dois componentes ao mesmo tempo.

Portanto, foi proposto um método de análise em que a largura da função da janela aumentava para baixas frequências e diminuía para altas frequências. A transformação da nova janela foi obtida como resultado do alongamento (compressão) e do deslocamento de tempo de uma função geradora (a chamada função de escala - função de escala, escalamento). Essa função geradora foi chamada de wavelet por D. Morlaix.

Wavelet D. Morlaix

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

Transformação Contínua de Wavelet (CWT)

Transformada wavelet de nível único:

coefs, frequências = pywt.cwt (dados, escalas, wavelet)

onde:

data: (array_like) - sinal de entrada;

scale (array_like) : - A escala de wavelet a ser usada. Você pode usar a razão f = scale2frequency (scale, wavelet) / sampling_period e determinar qual frequência física f é. Aqui, f em hertz quando sampling_period em segundos;

wavelet (objeto ou nome da Wavelet ) : - Wavelet para usar;

sampling_period (float): - Período de amostragem para frequências de saída (parâmetro opcional). Os valores calculados para coefs são independentes do período de amostragem (ou seja, as escalas não são dimensionadas pela amostragem de período);

coefs (array_like) : - wavelet contínuo - transformação do sinal de entrada para uma determinada escala e wavelet;

frequencies (array_like) : - Se a unidade do período de amostragem for o segundo e estiver definida, as frequências serão exibidas em hertz. Caso contrário, é assumido um período de amostragem 1.

Nota: O tamanho das matrizes de coeficiente depende do comprimento da matriz de entrada e dos comprimentos das escalas fornecidas.

Obtenha uma lista dos nomes de wavelets compatíveis com o cwt disponíveis:

 >>> 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 a análise de wavelets, a escolha da wavelet mãe é uma das principais tarefas. Portanto, antes de cada seleção de wavelet, é simplesmente necessário usar o seguinte programa que permite determinar as propriedades dinâmicas da wavelet:

Listagem 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 seguir, consideraremos as funções da wavelet e suas propriedades usando o programa acima:

Wavelet mexh de chapéu mexicano:

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

Wavelet complexo de Morlet "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}$

onde B é a largura de banda; C é a frequência central.

A seguir (sem indicação especial) os valores de B, C são definidos com um ponto flutuante.

Wavelets 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}}$

onde: P - um número inteiro de 1 a 8 é inserido no nome da wavelet, por exemplo "gaus8"; C - constante de normalização incorporada.

Wavelets Gaussianos Integrados “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}$

onde: P - um número inteiro de 1 a 8 é inserido no nome da wavelet, por exemplo, "gaus8" e corresponde à ordem da derivada da função wavelet; C - constante de normalização incorporada.

Wavelets de 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}$

onde: B é a largura de banda; C é a frequência central.

O CWT no PyWavelets se aplica a dados discretos - convoluções com amostras da integral da wavelet. Se a escala é muito pequena, obtemos um filtro discreto com amostragem inadequada, o que leva à suavização, conforme mostrado no gráfico da wavelet 'cmor1.5-1.0'.

Na coluna da esquerda, os gráficos mostram filtros discretos usados na convolução para várias escalas. A coluna da direita é o espectro de potência Fourier correspondente de cada filtro. Nas escalas 1 e 2 para o gráfico 'cmor1.5-1.0', pode-se observar que a suavização ocorre devido à violação da restrição de Nyquist. O fenômeno indicado pode ocorrer em outras wavelets ao escolher C e B, portanto, ao trabalhar com wavelets, você deve usar o programa - Listagem 1.

Para relacionar uma determinada escala a uma frequência de sinal específica, o período de amostragem do sinal deve ser conhecido. Usando a função pywt.scale2frequency () , você pode converter a lista de escalas nas frequências correspondentes. A escolha correta das escalas depende da wavelet selecionada; portanto, pywt.scale2frequency () deve ser usado para ter uma idéia da faixa de frequência correspondente do sinal.

 >>> 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álise de séries temporais usando CWT PyWavelets

O conjunto de dados el-Nino é um conjunto de dados de séries temporais usado para rastrear o El Nino e contém medições trimestrais de temperatura da superfície do mar de 1871 a 1997. Para entender o poder de um gráfico de escala, vamos visualizá-lo para um conjunto de dados el-Nino, juntamente com os dados originais da série temporal e sua transformação de Fourier.

Para a análise wavelet de séries temporais, as seguintes etapas devem ser executadas:

1. Selecione a wavelet mãe: Selecione a complexa wavelet Morlet "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}$

Largura de banda - B e frequência central - A partir da qual determinaremos na próxima etapa.

2. Determine a frequência central tomando dt = 0,25 para nossa série temporal:

 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,666666667]

3. Encontramos a transformação de Fourier da mãe wavelet cmor1.0-0.5 (usamos a listagem 1):
A wavelet contínua será avaliada em todo o intervalo [-8,0, 8,0]

4. Selecione os dados para a série temporal:

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

Os dados são medições trimestrais da temperatura da superfície do mar de 1871 a 1997. Para esses dados: t0 = 1871 dt = 0,25

5. Construímos uma série temporal com um sinal e sua média móvel em um gráfico:

Listagem 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. Realizamos a modalidade de transformação e espectro de Fourier a partir da série temporal:

Listagem 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 as escalas: escalas = arange (1, 128); níveis = [2 ** - 4, 2 ** - 3, 2 ** - 2, 2 ** - 1, 2 ** 0, 2 ** 1, 2 ** 2, 2 ** 3].

8. Crie um gráfico de escala:

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

O gráfico de escala mostra que a maior parte da potência do espectro está concentrada durante um período de 2 a 8 anos, o que corresponde a uma frequência de 0,125 a 0,5 Hz. No espectro de Fourier, a moda também se concentra em torno desses valores de frequência. No entanto, a transformação wavelet também nos fornece informações temporárias, mas a transformação Fourier não.

Por exemplo, no diagrama de escala, vemos que antes de 1920 havia muitas flutuações, enquanto entre 1960 e 1990 não havia tantas. Também podemos ver que, com o tempo, há uma mudança de períodos mais curtos para períodos mais longos.

Usando o módulo scaleogram

O Scaleogram é uma ferramenta conveniente para analisar dados 1D com transformação de wavelet contínua com base na biblioteca PyWavelets. O módulo foi projetado para fornecer uma ferramenta confiável para análise rápida de dados.

O Scaleogram possui os seguintes recursos:

assinatura simples para iniciantes
eixos legíveis e pura integração com matplotlib
muitas opções para zoom, filtro de espectro, implementação da barra de cores, etc ...
suporte para unidades de frequência e frequência de acordo com a marcação
velocidade, usa o algoritmo N * log (N) para transformações
portabilidade: testado com python2.7 e python3.7
mensagens de erro detalhadas e documentação com exemplos

No entanto, nos exemplos de uso, o acesso aos dados não é registrado corretamente:

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

O seguinte aviso aparece:

nino3 = pd.read_table (nino3)
levando a um aviso:
Aviso (do módulo de avisos):
Arquivo “C: /Users/User/Desktop/wavelet/pr.py”, linha 7
nino3 = pd.read_table (nino3)
FutureWarning: read_table está obsoleta, use read_csv em vez disso, passando sep = '\ t'

O acesso aos dados deve ser escrito assim:

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

Para mostrar o uso do escalograma e comparar os resultados com os resultados do exemplo anterior, usamos os mesmos dados - em medições trimestrais da temperatura da superfície do mar de 1871 a 1997. Para esses dados: t0 = 1871 dt = 0,25

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

Se compararmos o gráfico de escala com a digitalização resultante, com exceção da paleta de cores, eles são idênticos e, portanto, mostram a mesma dinâmica na série temporal.
A Listagem 4 é mais simples que a Listagem 3 e tem melhor desempenho.

Quando o espectro de frequências e o espectro de potência de acordo com Fourier não permitem obter informações adicionais sobre a dinâmica das séries temporais, é preferível o uso de um escalograma.

Conclusões

Um exemplo de análise wavelet de uma série temporal por meio da biblioteca PyWavelets é dado.

Análise wavelet. Parte 2

1. Introdução

Transformação Contínua de Wavelet (CWT)

Wavelet - Análise de séries temporais usando CWT PyWavelets

Usando o módulo scaleogram

Conclusões

More articles: