Topicos Avançados - Aula 05 - Análise no Domínio da Frequência: de Fourier a Wavelets

PRINCÍPIOS E TÉCNICAS DE ELETROENCEFALOGRAFIA EM NEUROCIÊNCIA

AULA 05 - Análise no Domínio da Frequência: de Fourier a Wavelets

Instituto de Ciência e Tecnologia

Graduação em Engenharia Biomédica

Prof. Dr. Adenauer G. Casali

Laboratório de Neuroengenharia e Computação

casali@unifesp.br

PRINCÍPIOS E TÉCNICAS DE EEG EM NEUROCIÊNCIA

Transformada Discreta de Fourier e FFT
As limitações da FFT: spectral leakage, dilema viés x variância
Métodos de janelamento: sFFT, Welch e Multitaper
Do janelamento para wavelets: Decomposição Tempo-Frequência

Adenauer G. CASALI

AULA 05

Nesta aula, nós veremos...

Transformada de Fourier de Tempo Discreto

X(\omega)=\int_{-\infty}^{\infty} x(t)e^{-j\omega t}dt

|X(\omega)|

\angle X(\omega)

Amplitude

Fase

1. Transformada Discreta de Fourier e a FFT

X(\Omega) = \sum_{n=-\infty}^{\infty}x[n]e^{-j\Omega n}

|X(\Omega)|

\angle X(\Omega)

Amplitude

Fase

\omega = \frac{\Omega}{T}

Frequência física (rad/s)

Frequência normalizada (rad)

Intervalo de amostragem

|\omega|\leq \frac{\pi}{T}

Amostragem

\omega_a = \frac{2\pi}{T}

Frequência de amostragem (rad/s)

Transformada de Fourier

Princípios e Técnicas de EEG

Aula 05

Mas como é possível calcular a Transformada de Fourier no computador se a frequência normalizada é, em geral, uma função contínua?

X(\Omega) = \sum_{n=-\infty}^{\infty}x[n]e^{-j\Omega n}

Discreto

Contínuo

Um sinal de EEG: Transformada de Fourier de tempo discreto é função de variável contínua

Princípios e Técnicas de EEG

Aula 05

1. Transformada Discreta de Fourier e a FFT

Discretizar a Frequência envolve tornar o sinal periódico.

: não periódico

cópias de : sinal periódico!

=período

x[n]

x^{(N)}[n]

: sinal construído repetindo-se a cada pontos

x[n]

Princípios e Técnicas de EEG

Aula 05

1. Transformada Discreta de Fourier e a FFT

Discreto!

x^{(N)}[n]

X^{(N)}(\Omega) = 2\pi\sum_{k=-\infty}^{\infty}a_k\delta(\Omega - k\frac{2\pi}{N})

a_k = \frac{1}{N}\sum_{n=\langle N\rangle}x^{(N)}[n]e^{-jk\frac{2\pi}{N}n}

X^{(N)}[k]= N a_k

=\sum_{n=\langle N\rangle}x^{(N)}[n]e^{-jk\frac{2\pi}{N}n}

= \sum_{n=0}^{N-1}x[n]e^{-jk\frac{2\pi}{N}n}

Qual o espectro desta versão replicada do sinal?

Coeficientes da série de Fourier

Número do Harmônico

Frequência Fundamental

Só frequências múltiplas da fundamental

Princípios e Técnicas de EEG

Aula 05

1. Transformada Discreta de Fourier e a FFT

X(\Omega) = \sum_{n=-\infty}^{\infty}x[n]e^{-j\Omega n}

Discreto

Contínuo

X^{(N)}[k] = \sum_{n=0}^{N-1}x[n]e^{-jk\frac{2\pi}{N}n}

DTFT:

Discreto

"N-point Discrete Fourier Transform" (DFT)

"Discrete-Time Fourier Transform" (DTFT)

DFT:

A DFT é a DTFT de infinitas cópias do sinal original, espaçadas entre si de N pontos

corresponde à frequência normalizada de

k\frac{2\pi}{N}

A escolha de N determina a resolução (e o número) de frequências discretas

(ou à frequência física de )

k\frac{\omega_a}{N}

(-\frac{N}{2} < k\leq \frac{N}{2})

Princípios e Técnicas de EEG

Aula 05

1. Transformada Discreta de Fourier e a FFT

"N-point Discrete Fourier Transform" (DFT)

DFT:

FFT: Fast Fourier Transform - modo otimizado de calcular a DFT

A DFT é periódica (em k) com período N:
Para cada k: temos que calcular N produtos de números complexos e somá-los.
O número total de produtos de números complexos é

N^2

X^{(N)}[k] = \sum_{n=0}^{N-1}x[n]e^{-jk\frac{2\pi}{N}n}

X^{(N)}[k+N] = X^{(N)}[k]

Princípios e Técnicas de EEG

Aula 05

1. Transformada Discreta de Fourier e a FFT

2. Limitações da FFT

"Spectral Leakage"

Mas esta não é a única limitação!

Princípios e Técnicas de EEG

Aula 05

Considere o seguinte exemplo:

y[n] - 0.75y[n-1] +0.5y[n-2] = r[n]

Processo gerado por um filtro IIR de segunda ordem:

Ruído branco (espectro = 1)

Espectro de y é a reposta em frequência do filtro:

Princípios e Técnicas de EEG

Aula 05

2. Limitações da FFT

Considere o seguinte exemplo:

y[n] - 0.75y[n-1] +0.5y[n-2] = r[n]

Processo gerado por um filtro IIR de segunda ordem:

Amostra com 1000 pontos do sinal

FFT

Espectro verdadeiro

Espectro ruidoso!

Princípios e Técnicas de EEG

Aula 05

FFT (e DFT) são estimadores inconsistentes da densidade espectral de potência!

2. Limitações da FFT

Este é um problema de viés x variância!

Princípios e Técnicas de EEG

Aula 05

Pouca resolução espectral

(trecho curto)

Ruído amostral (inconsistência estatística)

"trade-off"

3. Métodos de Janelamento

Estratégia para reduzir o ruído na estima do espectro:

"janelar" o sinal e mediar os espectros

short-FFT (sFFT): média entre os espectros construídos com a FFT de pedaços do sinal

Mas cortar pedaços produz "spectral leakage"!

Usar janelas suaves nas bordas!

Método de Welch:

Princípios e Técnicas de EEG

Aula 05

Este é um problema de viés x variância!

Princípios e Técnicas de EEG

Aula 05

Pouca resolução espectral

(trecho curto)

Ruído amostral (inconsistência estatística)

"trade-off"

Escolher uma janela específica é assumir um compromisso específico entre viés e variância

QUE TAL USAR VÁRIAS?

3. Métodos de Janelamento

Fonte: Prerau et al., 2017

Método "Multitaper": usar uma base ortogonal de janelas

Funções ortogonais!

(DPSS: "Discrete Prolate Spheroidal Sequences")

sFFT

Welch

Multitaper

Princípios e Técnicas de EEG

Aula 05

3. Métodos de Janelamento

Matematicamente:

Princípios e Técnicas de EEG

Aula 05

P(k) = \frac{1}{L}\sum_{l=1}^{L}\Bigl|\sum_{n=0}^{N-1}v_{l}[n]x[n]e^{-jk\frac{2\pi}{N} n}\Bigr|^2

"Tapers"

Pergunta possível: qual sequência discreta de janelas possui a maior fração possível de sua energia espectral concentrada dentro de uma banda limitada? (isto é, quais janelas possuem espectros mais estreitos possíveis?)

Todo sinal finito sofre de vazamento espectral devido ao espectro largo da janela que produz o truncamento do sinal no tempo.

V(k)

Transformada de Fourier da janela

\lambda = \frac{\sum_{k=k_1}^{k_2}|V(k)|^2}{\sum_{k=0}^{N}|V(k)|^2}

Fração da energia espectral que está na banda de largura desejada

MAXIMIZAR

3. Métodos de Janelamento

Slepian (1978):

Princípios e Técnicas de EEG

Aula 05

Se N for o N-point empregado e se W for meia largura da banda desejada, existem aproximadamente

L\approx 2NW

janelas (tapers) com espectro quase perfeitamente concentrado.

Estas soluções formam as chamadas DPSS (Discrete Prolate Spheroidal Sequences)

Propriedade importante das DPSS: ORTOGONALIDADE

\sum_{n=0}^{N-1}v_{l}[n]v_{s}[n] = 0 \:\:(s\neq l)

(Dentro da banda os espectros dos tapers também são aproximadamente ortogonais).

3. Métodos de Janelamento

Parâmetro fundamental:

Princípios e Técnicas de EEG

Aula 05

L = 2NW

NW pequeno = poucas janelas, alta resolução espectral (mais variância)

NW grande = muitas janelas, menor variância (baixa resolução espectral)

Na prática: ao controlar o parâmetro W, o multitaper dá mais controle da suavização do espectro (que é feita em uma banda de largura 2W). Isso reduz variância, mas pode misturar picos próximos!

3. Métodos de Janelamento

Princípios e Técnicas de EEG

Aula 05

3. Métodos de Janelamento

y[n] − 2.7607y[n − 1] + 3.8106y[n-2] − 2.6535y[n-3] + 0.9238y[n − 4] = r[n]

Outro exemplo:

Princípios e Técnicas de EEG

Aula 05

3. Métodos de Janelamento

Princípios e Técnicas de EEG

Aula 05

Leakage spectral: energia das frequências baixas vazando para frequências altas

3. Métodos de Janelamento

Princípios e Técnicas de EEG

Aula 05

3. Métodos de Janelamento

No EEG:

Multitaper pode ser mais apropriado para frequências altas e poucas épocas (SNR baixo), pois tem maior "smoothing" na frequência
Welch ou sFFT são mais apropriados para frequências baixas ou mais épocas (SNR alto) pois estes métodos têm mais potencial de separação de picos de frequência
Nenhum dos três métodos é indicado para quando é necessário maior precisão temporal e de fase. Para isso são necessárias análises no domínio Tempo-Frequência.

Princípios e Técnicas de EEG

Aula 05

3. Métodos de Janelamento

4. De Fourier a Wavelets

x(t)

X(\omega)

y(t) = x(t)\ast e^{-j\omega t}

Y(\theta)

Transformada de Fourier:

Y(\theta) = 2\pi X(\theta)\delta(\theta-\omega)

2\pi \delta(\theta-\omega)

= 2\pi X(\omega)\delta(\theta-\omega)

A transformada de Fourier de um sinal em determinada frequência pode ser vista como o resultado da convolução do sinal com uma exponencial complexa de mesma frequência

Considere a seguinte convolução no tempo:

X(\theta)

Tempo

Frequência

Produto

Convolução

Produto

Na frequência:

Princípios e Técnicas de EEG

Aula 05

Filtrando com a convolução:

Transformada de Fourier:

Sinal deve ser estacionário!

O Kernel extrai uma frequência fixa igualmente de todos os tempos

Mas e a ideia de realizar um janelamento?

4. De Fourier a Wavelets

Princípios e Técnicas de EEG

Aula 05

short-time Fast Fourier Transform (sFFT):

Kernel localizado no tempo

Requisito de estacionariedade é mais fraco: sinal deve ser estacionário somente no intervalo do Kernel

Wavelets

4. De Fourier a Wavelets

Princípios e Técnicas de EEG

Aula 05

Wavelets

\psi(t)= g(t)e^{j\omega_0 t}

Frequência da wavelet

"Envelope"

g(t)

e^{j\omega_0 t}

\psi(t)

Source: https://www.youtube.com/watch?v=jnxqHcObNK4

5. Decomposição Tempo-Frequência com Wavelets

Princípios e Técnicas de EEG

Aula 05

Wavelet de Morlet = Senóide x Gaussiana

g(t) = e^{-\bigl(\frac{4{\ln 2}}{\sigma^2}\bigr)t^2}

\frac{g(t_{1})}{g(0)} =\frac{1}{2} = e^{-\bigl(\frac{4{\ln 2}}{\sigma^2}\bigr)t_{1}^2}

-\ln 2 = -\bigl(\frac{4{\ln 2}}{\sigma^2}\bigr) t_{1}^2

t_{1}^2 = \frac{\sigma^2}{4}

t_{1} = \pm \frac{\sigma}{2}

\frac{\sigma}{2}

-\frac{\sigma}{2}

Parâmetro sigma:

\Delta t = \sigma

\psi(t)= g(t)e^{j\omega_0 t}

= intervalo de tempo correspondente a um número k de ciclos oscilando na frequência média da wavelet

\Delta t = \sigma

= \frac{2\pi k}{\omega_0} = kT_0

Parametrização de Morlet

5. Decomposição Tempo-Frequência com Wavelets

Princípios e Técnicas de EEG

Aula 05

Wavelet de Morlet = Senóide x Gaussiana

Fonte: Russel and Han (2016)

T_0=20ms

k=2

5. Decomposição Tempo-Frequência com Wavelets

Princípios e Técnicas de EEG

Aula 05

Princípio da Incerteza

\psi(t)= e^{-\bigl(\frac{4{\ln 2}}{\sigma^2}\bigr)t^2}e^{j\omega_0 t} = e^{j\omega_0 t -\alpha^2t^2}

\Psi(\omega)= \int_{-\infty}^{\infty} e^{j\omega_0 t -\alpha^2 t^2} e^{-j\omega t} dt

= \int_{-\infty}^{\infty} e^{(\frac{(\omega_0-\omega)}{2\alpha} + j\alpha t)^2}e^{-\frac{(\omega_0-\omega)^2 }{4\alpha^2 } }dt

=e^{-\frac{(\omega_0-\omega)^2 }{4\alpha^2 } } \int_{-\infty}^{\infty} \frac{e^{(ju)^2}}{\alpha} du

=e^{-\frac{(\omega_0-\omega)^2 }{4\alpha^2 } } \int_{-\infty}^{\infty} \frac{e^{-u^2}}{\alpha} du

\int_{-\infty}^{\infty} e^{-x^2}dx = \sqrt{\pi}

=\frac{\sqrt{\pi}}{\alpha} e^{-\frac{(\omega_0-\omega)^2 }{4\alpha^2 } }

\Psi(\omega)=\frac{\sigma \sqrt{\pi}}{2\sqrt{\ln 2}} e^{-\frac{\sigma^2 }{16 \ln 2 } (\omega_0-\omega)^2 }

= \int_{-\infty}^{\infty} e^{j(\omega_0-\omega) t -\alpha^2 t^2} dt

Na Frequência:

u = -j\frac{(\omega_0-\omega)}{2\alpha} + \alpha t

\alpha^2 = \frac{4\ln 2}{\sigma^2}

5. Decomposição Tempo-Frequência com Wavelets

Princípios e Técnicas de EEG

Aula 05

Princípio da Incerteza

\frac{\Psi(\omega_{1})}{\Psi(\omega_{0})}=\frac{1}{2}=e^{-\frac{\sigma^2 }{16 \ln 2 }(\omega_0-\omega_{1})^2 }

-\ln 2 = -\frac{\sigma^2 }{16 \ln 2 } (\omega_0-\omega_{1})^2

(\omega_0-\omega_{1})^2 = \frac{16 (\ln 2)^2}{\sigma^2 }

(\omega_0-\omega_{1}) = \pm \frac{4 \ln 2}{\sigma}

\Delta \omega

\omega_0 + \frac{4 (\ln 2)}{\sigma }

\omega_0 - \frac{4 (\ln 2)}{\sigma }

\Delta \omega = \frac{8\ln 2}{\sigma}

\Delta t = \sigma

\Delta \omega \Delta t= 8\ln 2

\Delta f \Delta t= \frac{8\ln 2}{2\pi}\approx 0.88

Princípio da Incerteza:

\Psi(\omega)=\frac{\sigma \sqrt{\pi}}{2\sqrt{\ln 2}} e^{-\frac{\sigma^2 }{16 \ln 2 } (\omega_0-\omega)^2 }

5. Decomposição Tempo-Frequência com Wavelets

Princípios e Técnicas de EEG

Aula 05

Source: https://www.youtube.com/watch?v=jnxqHcObNK4

5. Decomposição Tempo-Frequência com Wavelets

Princípios e Técnicas de EEG

Aula 05

Exemplo com ERP

Fonte: Herrman et al. (2014)

5. Decomposição Tempo-Frequência com Wavelets

Princípios e Técnicas de EEG

Aula 05

Transformada Wavelet (formulação moderna):

X(\omega)=\int_{-\infty}^{\infty} x(t)e^{-j\omega t}dt

X_a(b,\tau)=\frac{1}{\sqrt{a}}\int_{-\infty}^{\infty} x(t)\psi^{\ast}_{b}(\frac{t-\tau}{a})dt

\psi_b(x)= \frac{c_{b}}{\sqrt[4]{\pi}}e^{-x^2/2}\bigl(e^{jb x}-e^{-b^2/2})

"Wavelet mãe"

(mother wavelet)

Escala

Centro do tempo

Relacionados ao centro da frequência

Exemplo de wavelet mãe: Morlet

Normalização

Gaussiana

"Carrier"

Termo normalmente insignificante, já que, tipicamente

b>5

5. Decomposição Tempo-Frequência com Wavelets

Princípios e Técnicas de EEG

Aula 05

\psi_b(x)= \frac{c_{b}}{\sqrt[4]{\pi}}e^{-x^2/2}\bigl(e^{jb x}-e^{-b^2/2})

\psi_b(\frac{t-\tau}{a})= \frac{c_{b}}{\sqrt[4]{\pi}}e^{-\frac{(t-\tau)^2}{2a^2}}\bigl(e^{j\frac{b}{a}(t-\tau)}-e^{-b^2/2})

\omega_0 = \frac{b}{a}

a^2=\frac{{\sigma^2}}{8{\ln 2}}

\sigma= \frac{2\pi k}{\omega_0}

b=\frac{{\pi }}{\sqrt{2\ln 2}}k

a=\frac{{\pi k }}{\sqrt{2\ln 2}\omega_0}

f_{min} = \frac{k}{L_{max}}

a=\frac{1}{2\sqrt{2\ln 2} }\frac{k}{f}

a_{max} = \frac{L_{max}}{2\sqrt{2\ln 2}}

Transformada Wavelet de Morlet (formulação moderna):

=\frac{b}{\omega_0}

f =

frequência em Hz

L_{max} =

janela de tempo máxima que podemos usar

k\geq 3

Tipicamente:

5. Decomposição Tempo-Frequência com Wavelets

Princípios e Técnicas de EEG

Aula 05

b = \frac{{\pi }}{\sqrt{2\ln 2}}k \approx 8

Decomposição "Tempo-Frequência":

k = 3 \text{ ciclos}

f_{min} = 6Hz

a_{max} = \frac{1}{2\sqrt{2\ln 2} }\frac{k}{f_{min}} \approx 0.21233

L_{max} = 500ms

X(f,\tau)=\frac{1}{\sqrt{a}}\int_{-\infty}^{\infty} x(t)\psi^{\ast}_{b}(\frac{t-\tau}{a})dt

f = \frac{b}{2\pi a}

f_{max} = 50Hz

a_{min} \approx 0.02546

Exemplo de configuração:

5. Decomposição Tempo-Frequência com Wavelets

Princípios e Técnicas de EEG

Aula 05

Exemplo com potencial evocado real

Frequencia mínima de interesse

Janela disponível

|X(f,\tau)|^2

x(t)

5. Decomposição Tempo-Frequência com Wavelets

Princípios e Técnicas de EEG

Aula 05

Exemplo com potencial evocado real

Muito localizado no tempo, sem resolução na frequência!

|X(f,\tau)|^2

x(t)

5. Decomposição Tempo-Frequência com Wavelets

Princípios e Técnicas de EEG

Aula 05

Exemplo com potencial evocado real

|X(f,\tau)|^2

x(t)

5. Decomposição Tempo-Frequência com Wavelets

Princípios e Técnicas de EEG

Aula 05

Exemplo com potencial evocado real

|X(f,\tau)|^2

x(t)

5. Decomposição Tempo-Frequência com Wavelets

Princípios e Técnicas de EEG

Aula 05

Exemplo com potencial evocado real

|X(f,\tau)|^2

x(t)

5. Decomposição Tempo-Frequência com Wavelets

Princípios e Técnicas de EEG

Aula 05

Exemplo com potencial evocado real

|X(f,\tau)|^2

x(t)

5. Decomposição Tempo-Frequência com Wavelets

Princípios e Técnicas de EEG

Aula 05

Exemplo com potencial evocado real

Espalhando demais no tempo!

|X(f,\tau)|^2

x(t)

5. Decomposição Tempo-Frequência com Wavelets

Princípios e Técnicas de EEG

Aula 05

Exemplo com potencial evocado real

|X(f,\tau)|^2

x(t)

5. Decomposição Tempo-Frequência com Wavelets

Princípios e Técnicas de EEG

Aula 05

Exemplo com potencial evocado real

Baixa frequência sem resolução: uma única wavelet dura 1s

|X(f,\tau)|^2

x(t)

5. Decomposição Tempo-Frequência com Wavelets

Princípios e Técnicas de EEG

Aula 05

Exemplo com EEG espontâneo: dinâmica temporal do alfa

Fonte: Herrman et al. (2014)

5. Decomposição Tempo-Frequência com Wavelets

Princípios e Técnicas de EEG

Aula 05

FFT em Python:

numpy.fft.fft(x,n=N)

Calcula a N-point DFT do sinal x.
Se N for múltiplo de 2, calcula a FFT.
Há outras simetrias que são exploradas para otimizar o cálculo caso N não seja múltiplo de 2.
Retorna a transformada (você precisa saber o que é k para interpretar o resultado!).
Se N for menor que o tamanho do sinal: o sinal é cortado; se N for maior que o tamanho do sinal: o sinal é "zero-padded" (completado com zeros nos intervalos entre as cópias).

X^{(N)}[k]

Princípios e Técnicas de EEG

Aula 05

6. Na prática... (mais sobre isso na aula 18)

Welch em Python:

scipy.signal.welch(x, fs=Fam, window='hann', nperseg=Nwin, noverlap=Nover, nfft=N)

Calcula a N-point DFT de trechos do sinal x e faz a média da potência espectral entre todos os trechos.
Os trechos são janelados com janelas de determinado tipo ("window='hann') tamanho "Nwin" e com tamanho de sobreposição "Nover".

Frequência de amostragem

N-point

Número de amostras da janela

Número de amostras que são sobrepostas