EEG/Pracownia EEG/SSVEP 2

Koncepcja drgania uogólnionego. Transformata Hilberta

Wstęp

Sygnałem najczęściej występującym w przyrodzie oraz najczęściej stosowanym w technice jest sygnał harmoniczny o postaci:

[math]x(t) = A \sin(\omega t + \varphi_0)[/math] gdzie:
[math]t[/math] — chwila czasu
[math]A[/math] — amplituda sygnału,
ω — częstość sygnału,
φ[math]_0[/math] — faza początkowa sygnału.

Okazuje się, że istnieje szeroka klasa sygnałów rzeczywistych, które można przedstawić w postaci tzw. drgania uogólnionego: [math]x(t) = A(t) \sin(\Omega(t)\cdot t)[/math]

gdzie:
[math]t[/math] — chwila czasu
[math]A(t)[/math] — amplituda chwilowa sygnału x(t) (jego obwiednia),
Ω(t) — częstość chwilowa sygnału x(t),

Do klasy sygnałów, które mogą być reprezentowane w postaci drgania uogólnionego, należą m.in. wszystkie sygnały o ograniczonej energii i ograniczonej mocy średniej przedziałami ciągłe i bez składowej stałej (tzw. sygnały przestrzeni [math]L^2[/math]). W celu przedstawienia sygnału [math]x(t)[/math] jako drgania uogólnionego należy wpierw wyznaczyć jego sygnał analityczny z(t), który zdefiniowany jest w następujący sposób:

[math]z(t) = x(t) + ix_H(t)\;[/math]

gdzie:
[math]i=\sqrt{-1}[/math]
[math]x_H(t)[/math] — transformata Hilberta sygnału [math]x(t)[/math].

Transformatę Hilberta [math]x_H(t)[/math] sygnału [math]x(t)[/math] i transformatę do niej odwrotną definiujemy jak poniżej: [math] x_H(t) = \frac{1}{\pi}\int_{-\infty}^{\infty} \frac{x(\tau)}{t - \tau}d\tau [/math]
[math] x(t) = -\frac{1}{\pi}\int_{-\infty}^{\infty} \frac{x_H(\tau)}{t - \tau}d\tau [/math]

Uwaga praktyczna: 
do wyznaczania sygnału analitycznego korzysta się z jego następującej własności:
Widmo sygnału analitycznego odpowiadającego sygnałowi rzeczywistemu jest zerowe 
dla ujemnych częstości, zaś dla dodatnich częstości ma podwojoną amplitudę:
[math]x_a = F^{-1}\left(F(x)\cdot 2U\right)[/math]
gdzie F — transformacja Fouriera, a U funkcja schodkowa.
Metoda ta zaimplementowana jest w funkcji scipy.signal.hilbert

Jak można zauważyć, sygnał analityczny jest funkcją zespoloną, w związku z czym można go przedstawić w postaci: [math] z(t) = \left|z(t)\right|e^{i\varphi(t)} [/math] gdzie (patrz równanie (3)):
[math] \begin{array}{l} \left|z(t)\right| = \sqrt{x^2(t) + x^2_H(t)} \\ \\ \varphi(t) = \mathrm{arc\,tg}(\frac{x_H(t)}{x(t)}) \end{array} [/math]

Wielkości te służą do wyznaczania chwilowej fazy φ (wzór powyżej), chwilowej amplitudy A (obwiedni) oraz chwilowej częstości Ω sygnału [math]x(t)[/math]:

[math] \begin{array}{l} A(t)=\left|z(t)\right| \\ \\ \Omega(t) = \frac{d\varphi(t)}{dt} \end{array} [/math]

co umożliwia przedstawienie sygnału x(t) w postaci drgania uogólnionego: [math]x(t) = A(t) \sin(\Omega(t)\cdot t)[/math]

Porównując powyższy wzór ze wzorem na funkcję harmoniczną: [math]x(t) = A \sin(\omega t + \varphi_0)[/math]

widzimy, że sygnały nieharmoniczne charakteryzują się zmienną w czasie amplitudą i częstością oraz nie mają określonej fazy początkowej. Ten ostatni parametr jednak również może być wyznaczony, pod warunkiem iż określimy go względem pewnej stałej w czasie częstości ω₀(t): [math]\varphi(t)=\omega_0\cdot t + \varphi_0(t)[/math] gdzie:
[math]\varphi_0(t)[/math] — faza chwilowa sygnału [math]x(t)[/math].

Faza chwilowa jest zatem zdefiniowana jednoznacznie, ale tylko względem określonej częstości [math]\omega_0[/math]. W przypadku gdy analizujemy sygnały szerokopasmowe, wyznaczenie fazy chwilowej jest możliwe po przefiltrowaniu sygnału filtrem pasmowo-przepustowym.

Ćwiczenia symulacyjne

Proszę zaimplementować funkcje generujące następujące przebiegi czasowe (oznaczenia: F_s — częstość próbkowania, T — czas trwania w sekundach, analogicznie do tego):

• sinusoida o zadanej częstości f i fazie [math]\phi[/math]: y = sin(f, phi, Fs, T),
• funkcja Gabora o zadanym położeniu szerokości i częstości y = gabor(t0, sigma, omega, Fs, T):

[math]g_1(t) = \exp\left(-\frac{(t-t_0)^2}{2 \sigma^2}\right)\cdot \cos\left(\omega(t-t_0) \right)[/math].

Ćwiczenie 1

• Wygeneruj 2-sekundowy odcinek sygnału Gabora o pozycji w czasie t₀ = 1 s, częstości 16 Hz, skali 0,1 s i częstości próbkowania 128 Hz.
• Wyznacz i narysuj amplitudę chwilową sygnału Gabora.

Ćwiczenie 2

• Wyznacz i narysuj amplitudę chwilową sygnału sinusoidalnego, o częstości 16 Hz i częstości próbkowania 128 Hz.
• Wyznacz i narysuj fazę chwilową sygnału sinusoidalnego, korzystając ze związku:

[math]\varphi(t)=\omega_0\cdot t + \varphi_0(t)[/math]

Ćwiczenie 3

Wygeneruj dwa sygnały sinusoidalne o tej samej częstości 32 Hz i częstości próbkowania 128 Hz, ale różnych fazach początkowych. Pierwszy sygnał powinien mieć fazę początkową równą 0, drugi sygnał sinusoidalny powinien mieć fazę początkową równą π/4. Za pomocą transformaty Hilberta wyznacz różnicę faz symulowanych sygnałów.

Ćwiczenie 4

To ćwiczenie jest ku przestrodze. Wytwórz sygnał będący sumą dwóch sinusoid: jednej o częstości 30 i drugiej o częstości 32 Hz. Wykreśl przebieg sygnału i jego amplitudy chwilowej.

Ćwiczenia na danych pomiarowych

Ćwiczenie 1

W zebranych sygnałach SSVEP wybierz zapisy dla trzech różnych częstości stymulacji (po jednym dla każdej częstości). Do analizy wybierz trzy kanały EEG, dla których sygnał SSVEP jest a) bardzo wyraźny; b) widoczny, ale słabszy; c) w zasadzie niewidoczny. Do analizy wybierz fragmenty od 2 sekund przed rozpoczęciem stymulacji do 2 sekund po jej zakończeniu.

Dla każdej wybranej częstości stymulacji wybrane kanały EEG przefiltruj filtrem wąskopasmowym przepuszczającym częstości skupione wokół tej częstości stymulacji. Do przefiltrowanych sygnałów zastosuj transformację Hilberta, wyznacz amplitudę i częstość chwilową. Wyznacz fazę chwilową dla ω₀ równej częstości stymulacji. Wypróbuj dwa sposoby filtrowania: „w jedną stronę” (filter) i „w obie strony” (filtfilt).

Wyrysuj przefiltrowane sygnały wraz z wyliczoną amplitudą chwilową. Do rysunku dodaj wykres sygnału triggera aby widać było początek i koniec stymulacji. Narysuj też wykres zależności częstości chwilowej i fazy chwilowej od czasu.