Projecto de Física Computacional (2018-19)

25 de Dezembro de 2018
Projecto de Física Computacional (MEFT/IST 2018-19)
Dep. de Física, IST/U. Lisboa
25 a 31 de Dezembro 2018 (svn commit limit: 18H)

Docentes: Fernando Barão (responsável), Rui Coelho, Miguel Orcinha, Carlos Couto


Instruções muito importantes


Índice

Decomposição empírica de um sinal temporal

Decomposição empírica de um sinal temporal

Introdução

A análise de sinais em tempo é feita habitualmente com recurso a transformadas de Fourier que transformam a base de tempo do sinal, em frequência. O método de Fourier utiliza ondas harmónicas de amplitude constante como funções de base para a decomposição do sinal, assumido como estacionário (média, variância constantes no tempo). No entanto em muitas situações físicas, o sinal é não estacionário e não-linear (amplitude e frequência variáveis no tempo) e portanto a sua análise deve ser feita sem qualquer condição apriorística. Isto é, faz-se uma análise dos dados sem recurso a qualquer modelização.

Como exemplos de aplicação da análise de sinais variáveis com o tempo, poderíamos pensar no clima ao longo dos anos e desta forma analisar o efeito de aquecimento global, a evolução de cotações financeiras, as correntes oceânicas, a quantidade de dióxido de carbono emitida na atmosfera, etc. A decomposição do sinal deverá servir para separar o "ruído" (variações de alta frequência) de tendências existentes nos dados e eventualmente fazer previsões.

Neste projecto vamos desenvolver uma ferramenta de análise de sinais não estacionários, chamada decomposição empírica modal.

A decomposição empírica modal de um sinal (empirical mode decomposition - EMD), também chamada de transformada de Huang, foi inicialmente sugerida por Huang et al. em 1998 (1). Trata-se de uma técnica de análise e decomposição de sinal que trata os sinais como uma sobreposição de oscilações rápidas e lentas (2), sem necessidade de qualquer informação prévia do sinal. Esta técnica, simples, eficiente e iterativa, permite a decomposição do sinal num conjunto de funções temporais chamadas funções modais intrínsicas (intrinsic mode functions - IMF).

Suponhamos então um sinal temporal $s(t)$, discretizado e avaliado em $N$ pontos $t_0, t_1, t_2, \cdots, t_{N-1}$. A sua decomposição temporal em $n$ IMF's $c_j(t)$ , permitirá escrevê-lo como: \begin{equation} s(t) \simeq c_1(t) + c_2(t) + c_3(t) + \cdots = \sum_{j=1}^{n} c_j(t) \end{equation}

A referência (3) constitui uma boa base de partida para uma melhor compreensão do método e dos diferentes detalhes que o envolvem.

Algoritmo EMD para extração das IMF

O método visa decompor um dado sinal de entrada numa soma de vários sinais a determinar (“intrinsic mode functions – IMFs”), de modo iterativo, de tal modo que cada IMF da decomposição capture localmente i.e. em cada instante do tempo, as diferentes escalas de variação temporal do sinal original. A extracção das IMF consegue-se implementando o seguinte algoritmo, considerando o sinal original como possuindo N+1 amostragens temporais, i.e. dado pela função discretizada $s(t_i)$ com $i=0, \cdots, N$:

passo 1: Determinação das funções interpoladoras e função "detalhe" do sinal

Consideremos o sinal inicial $s_k(t) = s(t)$ com k=1 (1ª interação).
Determinar as funções envelope interpoladoras $e_{max}(t)$ e $e_{min}(t)$ que passem respectivamente pelos máximos locais e mínimos locais da função $s_k(t)$. Na literatura científica estão relatados o uso de vários tipos de interpoladores, tendo-se concluído que o uso de interpoladores cubic spline seria o mais indicado.

passos para a definição da função envelope $e_{max}(t)$:

passos para a definição da função envelope $e_{min}(t)$:

Uma nota sobre os pontos fronteira do sinal: poderá ser conveniente fazer o "mirroring" dos últimos extremos da função sinal para ter um bom comportamento na interpoladora.

Construir o "sinal médio" de $m_k(t)$ com $m_k(t) = \left( e_{max}(t) + e_{min} \right) / 2$.

Obter o “detalhe” do sinal $d_k(t)$, dado por $d_k(t_i)= s_k(t_i) – m_k(t)$. A função "detalhe" resulta da subtracção ao sinal, do seu valor médio. Será por isso uma função tendencialmente de "média nula".

Tecnicamente $d_k(t_i)$ seria suposto ser uma função IMF. No entanto, os erros introduzidos pela função interpoladora podem exigir várias iterações até se chegar à função IMF. Este processo de "regularização" dos dados, é conhecido como "sifting process".

passo 2: Iteração para determinação da IMF (processo "sifting")

Uma função $f(t)$ é considerada uma IMF, caso se considere que possui média zero, o que exige a definição dos seguintes critérios:

Teremos agora que definir o critério de aceitação da IMF. A "sobre-iteração", isto é, a iteração excessiva para a determinação da IMF é um problema aberto (não possui uma solução única e é objecto de publicações). Por isso devemos definir um "razoável" critério de paragem da iteração.

O critério introduzido por Rilling et al. (2), tem o objectivo de garantir globalmente pequenas flutuações na média, permitindo no entanto possíveis flutuações locais. Vamos por isso definir uma função de variação $\sigma(t)$ e estabelecer um critério de aceitação. A função de variação é definida tendo em conta o valor médio da função "detalhe", $m_k(t)$ e a sua amplitude $A_k(t)$, \begin{equation} \sigma(t) = \left| \frac{m_k(t)}{A_k(t)} \right| \end{equation} com:
$m_k(t) = \left( e_{max}(t) + e_{min} \right) / 2$
$A_k(t) = \frac{e_{max}(t) - e_{min}(t)}{2}$

O critério definido por Rilling et al. para paragem do processo de "sifting" é o de exigir que:


resumo dos passos para a determinação da IMF:

passo 3: Caso da iteração sem sucesso para obtenção da IMF

Se o “detalhe” $d_k(t_i)$ não for um sinal com “média nula”, então repete-se o passo 1. Neste caso não se incrementa $k$ e usamos na próxima iteração $s_k(t_i)=d_k(t_i)$.

passo 4: Fim da iteração

A decomposição do sinal cessa quando não for possível encontrar mais IMFs no resíduo $r_k(t)$. Tal acontecerá quando $r_k(t)$ cumprir uma das seguintes condições:

Determinação da frequência e amplitude instantâneas do sinal

Tendo em conta que as Intrinsic Mode Functions que se obtêm no final do algoritmo EMD são basicamente sinais AM/FM (amplitude e frequência variável), o procedimento mais rigoroso para estimar a frequência e amplitude instantânea de cada IMF envolveria a utilização da transformada de Hilbert de cada sinal (IMF), isto é, envolveria transformadas de Fourier. Porém, podemos obter uma evolução aproximada da frequência através da determinação dos instantes de tempo onde um sinal é $0$, e daí fazer a determinação do intervalo de tempo entre dois zeros consecutivos do sinal, que corresponde a metade do período do sinal.

Quanto à amplitude instantânea, podemos estimá-la a partir da semi-diferença entre o envelope dos máximos e mínimos locais do sinal, \begin{equation} A(t) = \frac{e_{max}(t) - e_{min}(t)}{2} \end{equation}

Objectivo e metodologia

Neste trabalho pretende-se primeiramente desenvolver uma biblioteca em C++ para a análise de sinais digitais (sinal definido apenas em determinados instantes de tempo) que implemente um método de análise especialmente útil quando os sinais em estudo se decompõem em várias componentes, cada uma com frequência instantânea variável no tempo – Decomposição Empírica de Sinal (“Empirical Mode Decomposition”). A estas várias componentes, na prática sinais AM/FM, denominam-se “intrinsic mode functions – IMFs”.

A biblioteca pode ser constituída por uma ou mais classes e deve recorrer a código (classes) desenvolvidos pelos alunos ao longo do semestre. Poderão ainda ser utilizados objectos ROOT para display.

Uma vez desenvolvida a biblioteca, deverão analisar os seguintes sinais.

1. Primeiramente, e de forma a validar o método de decomposição EMD implementado, vamos construir a função, \begin{equation} s(t) = A_1 sin(2 \pi f_1 t) + A_2 sin(2 \pi f_2 t) + A_3 sin(2 \pi f_3 t) \end{equation} com:
$A_1=0.5, \quad A_2=0.8, \quad A_3=1.5$
$f_1=1, \quad f_2=3f_1, \quad f_3=5f_1 \quad \text{Hz}$

2. O sol possui um campo magnético cuja orientação vai variando com o tempo de forma cíclica (5).

Um dos indicadores para avaliar a actividade solar é o número de "sun spots" que corresponde a zonas de mais baixa temperatura na superfície do sol (fotoesfera). Para mais detalhes ver aqui.

No ficheiro sunspots.dat encontram-se o número de sun spots medidos em cada mês desde 1749 até aos dias de hoje. Estes dados foram retirados do site http://www.sidc.be/silso/datafiles. Utilize esta contagem de sun spots como sendo o sinal em tempo que pretende decompor.


Bibliografia