out 072013
 
piwavelet

piwaveletUm grupo da National Oceanography Centre (NOC) da Universidade de Southampton  desenvolve um pacote para análise de coerência Wavelet (Ondeletas) em Matlab (http://noc.ac.uk/using-science/crosswavelet-wavelet-coherence). Como eu queria usar esse pacote, mas não tinha o menor interesse em usar o Matlab, eu acabei desenvolvendo uma API que permite usar as saídas desse pacote em um programa escrito em python. A página do meu projeto pode ser acessada em http://duducosmos.github.io/PIWavelet/. Basicamente o que fiz foi rodar o pacote original em Matlab no octave 3.6, os resultados obtidos são passados para um array python e assim vai.

 

 

 

A ideia aqui é fazer um breve tutorial sobre Wavelets e, dado dois conjuntos de sinal, como medir o grau de relação entre os dois.

Tradicionalmente o que se faz é, dado, por exemplo, dois conjuntos de dados, buscar qual a correlação linear entre eles. Mas fica a questão, e se a relação entre os dados não for linear? E se essa relação muda com o tempo?

Uma ferramente recente para resolver tais questões, em particular para análise de séries temporais, é o uso de Wavelets e da técnica de coerência, que mede o grau de relação entre os sinais ao longo do tempo.

 

Wavelet (Ondeleta)

 

A Wavelet, ou ondeleta em português, é uma função com média zero e que é definida em frequência e tempo. Pode se caracterizar a ondeleta pelo modo como ela se localiza no espaço de tempo e  de frequência. Aqui vale ressaltar que de acordo com  o princípio de incerteza de Heisenberg que sempre existirá uma incerteza intrínseca em o quão pequeno pode ser a banda de tempo e frequência.

 

Uma ondeleta particular é a de Morlet, a qual é definida como:

\psi_{0}(\eta)=\pi^{-1/4}e^{i\omega_{0}\eta}e^{-\frac{1}{2}\eta^{2}} (1)

 

sendo \omega_{0} a frequencia adimensional e \etao tempo adimensional.

 

Transformação contínua de ondeleta

 

A ideia da transformação contínua de ondeleta (Continuous Wavelet Transform - CWT) está em aplicar a ondeleta como um filtro passa bandas em uma série temporal. A ondeleta sofre estreitamento no tempo variando sua escala (s), tal que \eta=s.t, sendo que sua normalização tem unidade de energia.  A  CWT de uma função f é definida pela transformação integral:

\Psi(a,b) = \int_{-\infty}^{\infty}f(u)\psi^{*}_{a,b}(u)du, para a data-recalc-dims=0" />. (2)

 

Sendo:

\psi^{*}_{a,b}(u) = \frac{1}{a}\psi\left(\frac{u-b}{a}\right) (3)

representa uma família de funções wavelets, chamada de wavelet mãe. O parâmetro a refere-se a escala, b é a translação ou parâmetro de localização da wavelet mãe e \psi^{*}_{a,b}(u) é o complexo conjugado de   \psi_{a,b}(u).

Exemplo e uso do código:

No exemplo abaixo é gerado um sinal aleatório. Esse dado é normalizado e em seguida se gera o gráfico da CTW desse sinal.

[python]
import numpy as np
import matplotlib.pyplot as plt
from piwavelet import piwavelet
y1 = np.random.rand(50) #Generation of the Random Signal
y1 = (y1-y1.mean())/y1.std() #Normalization of the Signal 1
wave, scales, freqs, coi, fft, fftfreqs=piwavelet.cwt(y1) # If you want to know the individual properties.'
piwavelet.plotWavelet(y1,title='test',label='Random Signal',units='days')
[/python]

 

Na figura abaixo está a imagem resultante.

 

Decifrando o código:

Primeiramente são importados os módulos a serem usados

[python]
import numpy as np
import matplotlib.pyplot as plt
from piwavelet import piwavelet

[/python]

Aqui é gerado uma série de dados aleatória

[python]
y1 = np.random.rand(50) #Generation of the Random Signal

[/python]

 

O dado é normalizado de forma que é subtraída a média do mesmo e dividido pelo desvio padrão

[python]
y1 = (y1-y1.mean())/y1.std() #Normalization of the Signal 1

[/python]

 

Finalmente são obtidos a CTW, a escala, a frequência, o cone de influência, a transformada de Fourrier (TF) do sinal e a frequências da TF.

[python]
wave, scales, freqs, coi, fft, fftfreqs=piwavelet.cwt(y1) # If you want to know the individual properties.'
piwavelet.plotWavelet(y1,title='test',label='Random Signal',units='days')
[/python]

Coerência Wavelet :

O objetivo da coerência wavelet é o de identificar em que banda de frequência e em que intervalo de tempo duas séries temporais estão relacionadas. Para esse tipo de análise é preciso realizar uma suavização do espectro cruzado antes de calcular a coerência, pois de outra forma o resultado sempre será igual a 1. [3]

 

Dada duas séries X e Y, com as tranformadas wavelets W^{X}_{n}(s) e W^{Y}_{n}(s), sendo n o índice de tempo e s a escala, o espectro cruzado wavelet é dado por:

W^{XY}_{n}(s) = W^{x}_{n}(s)W^{Y*}_{n}(s) (4)

sendo que o simbolo * denota o complexo conjugado.

A coerência wavelet quadrada é definida como o valor absoluto quadrado o espectro cruzado de wavelets, normalizado pelos espectros de potência suavizados das wavelets [3]:

R^{2}_{n}(s) = \frac{|S(s^{-1}W^{XY}_{n}(s))|^{2}}{S(|s^{-1}W^{X}_{n}(s)|^{2})S(|s^{-1}W^{Y}_{n}(s)|^{2})} (5)

sendo que S(W) denota suavização no tempo e na escala. O fator s^{-1} é usado para converter o espectro em densidade de energia.

 A função de suavização vai depender da escala usada e da wavelet mãe considerada [3].

A função de suavização é definida como [5]:

S(W) = S_{scale}(S_{time}(W(s,t))) (6)

sendo que S_{scale} denota suavização ao logo do eixo da escala e S_{time} no do tempo.

 Para o caso da Wavelet de Morlet, tem-se [5]:

S_{time}(W)|_{s} = \left(W(t,s) c_{1}e^{-t^{2}/2s^{2}}\right)|_{s} (7)

S_{scale}(W)|_{t} = \left(W(t,s) c^{2}\Pi(0.6s)\right)|_{t} (8)

sendo que c_{1} e c_{2} constantes de normalização e \Pi é a função retângulo. O fator 0.6 é empiricamente determidado pelo comprimento de decorrelação de escala para a wavelet de Morelet [3].

Apresentação

Abaixo está uma apresentação sobre esse post.

Referências:

[1] Mallat, Stephane G. (1999). A wavelet tour of signal processing

[2] Addison, Paul S. The illustrated wavelet transform handbook

[3] Torrence, Christopher and Compo, Gilbert P. (1998). A Practical Guide to Wavelet Analysis

[4] Grinsted, A., Moore, J.C., Jevrejeva, S., 2004, Nonlin. Processes Geophys., 11, 561-566, doi:10.5194/npg-11-561-2004

[5] Jevrejeva, S., Moore, J.C., Grinsted, A., 2003, J. Geophys. Res., 108(D21), 4677, doi:10.1029/2003JD003417

[6] Torrence, C., Webster, P. ,1999, J.Clim., 12, 2679-2690

 Leave a Reply

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong>

(required)

(required)

%d blogueiros gostam disto: