Laboratorium EEG/CSP

Z Brain-wiki

Laboratorium_EEG/BSS

Prezentacja

slajdy z prezentacji

Ślepa separacja źródeł

Rozważmy N-kanałowy sygnał EEG. Próbkę tego sygnału możemy przedstawić jako punkt w przestrzeni rozpiętej przez osie, z których każda reprezentuje wartość potencjału w jednym kanale. Cały sygnał tworzy w tej przestrzeni chmurę punktów. Rozciągłość tej chmury w danym kierunku mówi nam o wariancji (zmienności) sygnału w tym kierunku.

Taki zbiór punktów wygodniej jest analizować w układzie współrzędnych zgodnym z osiami głównymi macierzy kowariancji. W dalszej części rozważań założymy, że te przestrzenie, w których rozważamy sygnały są przestrzeniami wektorowymi, a pojedyncze próbki wielokanałowego sygnału są wektorami.

Kowariancja.png

Filtry przestrzenne i ślepa separacja źródeł

Sygnał EEG jest superpozycją aktywności elektrycznej wielu źródeł. Jak można estymować aktywność samych źródeł?

Mieszanie.png

Niech:

[math]s(t)[/math] - aktywność niezależnych źródeł,
[math]x(t)[/math] mierzony sygnał
[math]A[/math] macierz przejścia taka, że:
[math]x(t) = A s(t)[/math] (*)
[math]s(t) = A^{-1}x(t) = P x(t)[/math]

Macierz kowariancji dla sygnałów [math]x(t)[/math] estymujemy tak:

[math] C_x = E[x(t)x(t)^T][/math]

Podstawiając (*) mamy:

[math] C_x = E[x x^T] = E[As(As)^T] = A E[s s^T] A^T = A C_s A^T[/math]

Z założenia, że źródła są niezależne wynika, że macierz [math]C_s[/math] jest diagonalna. Przekształcając powyższe równanie możemy zapisać:

[math]A^{-1} C_x (A^T)^{-1} = P C_x P^T = C_s[/math]

Odwzorowanie [math]P = A^{-1}[/math] diagonalizuje macierz [math]C_x[/math].

Powyższe rozumowanie jest słuszne w przypadku gdy mamy do czynienia z sygnałem stacjonarnym, tzn. jego macierz kowariancji jest niezależna od czasu, czyli przez cały czas aktywna jest ta sama konfiguracja źródeł niezależnych. W przypadku gdy tak nie jest to konstrukcję filtra przestrzennego można oprzeć o jednoczesną diagonalizację macierzy kowariancji odpowiadających różnym stanom osoby badanej.

Diagonalizacja.png

Common Spatial Pattern

Koncepcja

Dla ustalenia uwagi możemy myśleć o eksperymencie wywołującym potencjał P300. Mamy w nim dwie sytuacje eksperymentalne. Oznaczmy [math]T[/math] (target) próby, w których pojawił się oczekiwany bodziec, zaś [math]NT[/math] (non-target) gdy pojawił się bodziec standardowy. Chcielibyśmy znaleźć taki montaż (czyli taką kombinację liniową kanałów), który maksymalizuje stosunek mocy (wariancji) sygnałów rejestrowanych w dwóch rożnych warunkach eksperymentalnych.

Formalizm

Metoda ta polega na znalezieniu takiego kierunku [math]w[/math] w przestrzeni sygnałów, że sygnał z warunku [math]T[/math] rzutowany na ten kierunek ma dużą wariancje a sygnał z warunku [math]NT[/math] ma wariancję małą.

Rzutowanie sygnału [math]x(t)[/math] na kierunek [math]w[/math] odbywa się przez policzenie iloczynu skalarnego dla każdej chwili czasu [math]t[/math]:

[math] s_w(t) = w^T x(t)[/math]

Wariancja tego rzutowanego sygnału to:

[math] \mathrm{var}(s_w) = E[s_w s_w^T] = E[ w^T x (w^T x)^T] = w^T E[x x^T] w = w^T C_x w [/math]

Zatem znalezienie właściwego kierunku rzutowania można wyrazić jako szukanie maksimum wyrażenia [math] J(w) [/math] (jest to tzw. iloraz Rayleigha):

[math]J(w) = \frac{w^T C_T w}{w^T C_{NT} w} [/math]

Ekstremum tego ilorazu można znaleźć poprzez policzenie gradientu [math]J(w)[/math] i przyrównanie go do zera:

[math] \nabla J(w) = \frac{ 2 C_{T} w \left(w^T C_{NT} w\right) -2C_{NT} w \left(w^T C_{T} w \right)}{\left(w^T C_{NT} w\right)^2} = \frac{ 2}{w^T C_{NT} w}\left[ C_{T}w -\frac{w^T C_{T} w}{w^T C_{NT} w} C_{NT} w \right][/math]

Przyrównując to wyrażenie do zera dostajemy do rozwiązania tzw. uogólnione zagadnienie własne:

[math] C_{T}w =\frac{w^T C_{T} w}{w^T C_{NT} w} C_{NT} w [/math]

We wzorze tym liczba [math]\lambda=\frac{w^T C_{T} w}{w^T C_{NT} w}[/math] spełniająca to równanie jest uogólnioną wartością własną, wtedy [math]w[/math] jest uogólnionym wektorem własnym odpowiadającym tej wartości.

Aby znaleźć [math] \lambda[/math] i [math]w[/math] możemy wykorzystać w Matlabie funkcję eig. Funkcja ta rozwiązuje (również) uogólnione zagadnienia własne postaci AwBw dostarczając w wyniku macierz wektorów własnych (w kolumnach) oraz macierz zawierającą na przekątnej odpowiadające im wartości własne.


Ćwiczenie symulacyjne

Zastosowanie filtra CSP do detekcji potencjału P300

Analiza wstępna

Poszczególne etapy analizy proszę kodować w osobnych funkcjach. Funkcje te powinny być wywoływane z nadrzędnego skryptu, który powinien umożliwic wykoanie całości analiz.

  • Wczytać dane kalibracyjne do Matlaba i pociąć je na realizacje typu T — „target” (związane z wystąpieniami litery „B”) i NT — „non-target” (pozostałe litery) o długości −200 do +800 ms wokół triggerów. Dla każdej realizacji odjąć trend liniowy.
  • Sygnał zmontować wzgl. „połączonych uszu” i wyświetlić średnie przebiegi dla warunku T i NT w układzie topograficznym — wykorzystać w tym celu funkcję plottopo z pakietu Eeglab.

Poniżej zaprezentowany jest przykładowy skrypt do cięcia danych wokół znaczników. Działa on z plikami zawartymi w archiwum:

Plik:KalibracjaP300.tar.gz

Korzysta z funkcji pomocniczych dostępnych w dystrybucji obci w katalogu

/usr/share/openbci/analysis/matlab_obci_signal_processing

Openbci można pobrać z https://github.com/BrainTech/openbci

ZADANIE: Analiza CSP

Link do Read menager [1]

  • Wykonać analizę CSP wzmacniającą potencjał P300.
  • Zaprezentować średnią ze wszystkich kanałów źródłowych z warunku target (jeden kolor) i non-target (inny kolor) w subplotach ułożonych w prostokątnej siatce. Zaobserwować dla którego kanału średnie różnią się najbardziej. Czy jest związek tego kanału z wartościami własnymi?
  • Dla kanału najbardziej różnicującego wykonać mapki topograficzne (do wykonania tych mapek wykorzystać funkcję topoplot z pakietu eeglab) wektorów odpowiadających:
    • filtrowi przestrzennemu
    • rzutu topograficznego źródła na elektrody.

Filtry przestrzenne dla SSEP

ICA jako filtr przestrzenny

Definicja

Independent Component Analysis (ICA) jest metodą statystycznej analizy sygnałów, która dokonuje dekompozycji wielokanałowych zapisów na składowe niezależne w sensie statystycznym. Dwie składowe s1 i s2 są niezależne jeżeli wiedza o wartości s1 nie daje żadnych informacji o możliwych wartościach s2. ICA może być wyrażona przez prosty model generatywny:

x = Ds
gdzie x = {x1, x2, ..., xn} jest zmierzonym n kanałowym sygnałem, D jest macierzą mieszającą zaś s = {s1, s2, ..., sn} jest aktywnością n źródeł. Podstawowym założeniem dotyczącym s jest to, że si są statystycznie niezależne. Aby wyestymować model musimy też założyć, że składowe mają niegaussowskie rozkłady wartości (Hyvärinen, 2000).

Dodatkowo model ten zakłada następujące fakty:

  1. Sygnał jest liniową mieszaniną aktywności źródeł
  2. Sygnały pochodzące z każdego ze źródeł są niezależne od pozostałych
  3. Źródła oraz proces ich mieszania są stacjonarne, tzn, ich momenty statystyczne nie zależą od czasu
  4. Energie (wariancje) źródeł nie mogą być wyznaczone jednoznacznie. Dzieje się tak ponieważ pomnożenie amplitudy i-tego źródła może być uzyskane poprzez przemnożenie albo si albo przez przemnożenie i-tej kolumny macierzy D. Naturalnym rozwiązaniem tej niejednoznaczności jest wprowadzenie konwencji, że komponenty są normowane tak aby miały wariancję 1: E[si] = 1.
  5. Kolejność komponentów jest dowolna. Bo jeśli w ten sam sposób zmienimy kolejność komponentów w s i kolumn w D to dostaniemy dokładnie ten sam sygnał x.

Głównym wyzwaniem w analizie ICA jest estymacja macierzy mieszającej D. Gdy jest ona znana to komponenty mogą być wyliczone w następujący sposób:

s = D−1x

Estymacja

Znalezienie niezależnych komponentów może być rozważane w świetle Centralnego Twierdzenia Granicznego jako poszukiwanie komponentów o możliwie nie gaussowskim rozkładzie. Aby zrozumieć to podejście prześledźmy heurystykę zaproponowaną przez (Hyvärinen, 2000). Dla prostoty załóżmy, że poszukiwane źródła niezależne mają identyczne rozkłady. Zdefiniujmy

y = wTx.

Zauważmy, że jeśli

wT

jest jedną z kolumn macierzy

D−1,

to y jest jednym z poszukiwanych komponentów. Zamieniając zmienne

z = DTw

możemy napisać

y = wTx = wTDs = zTs.

Uwidacznia to fakt, że y jest liniową kombinacją składowych si z wagami danymi przez zi.

Z centralnego twierdzenia granicznego wynika, że suma niezależnych zmiennych losowych ma bardziej gaussowski charakter niż każda z tych zmiennych osobno. Liniowa kombinacja staje się najmniej gaussowska gdy z ma tylko jeden element niezerowy. W tym przypadku y jest proporcjonalny do si. Zatem problem estymacji modelu ICA może być sformułowany jako problem znalezienia takiego wektora w, który maksymalizuje niegaussowskość

y = wTx.

Maksymalizacja niegaussowskości y daje jeden niezależny komponent odpowiadający jednemu z 2n maksimów (bo mamy si i −si) w krajobrazie optymalizacyjnym. Aby znaleźć wszystkie niezależne komponenty musimy znaleźć wszystkie maksima. Ponieważ komponenty są nieskorelowane, to poszukiwania kolejnych komponentów można kontynuować w podprzestrzeni ortogonalnej do już znalezionych komponentów.

Obliczenia

Intuicyjna heurystyka poszukiwania najbardziej niegaussowskich składowych może być użyta do wyprowadzenia różnych funkcji kosztu, których optymalizacja daje model ICA, np. kurtoza.

[math]kurt(y) = E\{y^4\} - 3(E{y^2})^2[/math]

Inną miarą gassowskości jest neg-entropia, którą można wyprowadzić z entropii: Entropia jest miarą średniego zdziwienia wynikiem obserwacji zmiennej losowej: [math]H(Y) = - \sum_i P(Y= a_i) \log(P(Y=a_i)) [/math]

Negentropia jest zdefiniowana:

[math] J(y) = H(y_{gauss}) -H(y)[/math] gdzie [math] y_{gauss} [/math] jest gassuwską zmienną losową o takiej samej kowaiancji jak [math] y [/math].

Negentropia jest skomplikowana obliczeniow, więc w praktyce używana jest formuła przybliżona:

[math] J(y) \varpropto [E\{G (y)\} - E\{G(\nu)\}][/math]

[math]\nu [/math] jest zmienną losową ze standardowego rozkładu normalnego , a G są pewnymi niekwadratowymi funkcjami.

W algorytmie FastICA extremum negentropii jest znajdowane w procedurze bazującej na optymalizacji Netwona. (szczegóły np.: sekcja 6 w https://www.cs.helsinki.fi/u/ahyvarin/papers/NN00new.pdf)


Procedura wykorzystywana w eeglabie („runica”, Makeig 1996) dąży do minimalizacji informacji wzajemnej. Oba podejścia są w przybliżeniu równoważne (Hyvärinen, 2000), chociaż owo przybliżenie dla sygnałów elektrofizjologicznych nie zostało to jeszcze w pełni wyeksplorowane. Dla sygnałów o niskiej wymiarowości i spełniających dokładnie założenia ICA wszystkie powszechnie wykorzystywane algorytmy dają niemal identyczne wyniki.

Bardzo ważna uwaga
ogólną zasadą jest, że jeśli estymujemy N stabilnych komponentów (z N-kanałowych danych) to musimy dysponować kN2 punktami danych w każdym kanale, gdzie N2 jest liczbą elementów w macierzy D, którą ICA próbuje wyestymować, k jest liczbą całkowitą. Nie ma dobrych oszacowań teoretycznych na wielkość k, z praktycznych obserwacji wynika, że rośnie ona z liczbą kanałów.

Możliwe zastosowania

Najczęściej ICA jest stosowana jako narzędzie do:

  • usuwania artefaktów z sygnałów EEG (ruchy oczu i mięśnie)
  • wydobywania składowych do dalszej analizy (Onton, 2006)
  • jako analiza wstępna do lokalizacji źródeł (Grau, 2007).
  • ICA jest także stosowana w analize sygnałów EKG i EMG.

Bibliografia

Bazowa praca:

Nieco prościej opisana wersja z przykładami:

  • Hyvärinen, A. and Oja, E. (2000). Independent component analysis: Algorithms and applications. Neural Networks, 13(4-5):411–430.

https://www.cs.helsinki.fi/u/ahyvarin/papers/NN00new.pdf

ZADANIE: Wydobywanie interesujących komponentów

Dane do tej części ćwiczeń proszę pobrać i rozpakować w swoim katalogu: http://www.fuw.edu.pl/~jarekz/LabEEG/Dane_do_ICA_alfa.tar.gz

Pochodzą one z eksperymentu, w którym osoba badana siedziała z zamkniętymi oczami słuchając nagrania czytanego spokojnym głosem. Metadane opisujące sygnał znajdują się w pliku Miro.xml, zaś lokalizacje elektrod w pliku Miro-10-20-Cap.locs.

Proszę:

  • wczytać dane do eeglaba
  • wyedytować lokalizację elektrod
  • usunąć kanały nie zawierające EEG
  • zmienić referencje na średnią z kanałów A1 i A2
  • przefiltrować filtrem FIR górnoprzepustowym z częstością odcięcia 0,5 Hz
  • obejrzeć wstępnie przygotowane dane
  • policzyć ICA na całym sygnale
  • obejrzeć właściwości otrzymanych komponentów
    • Czy są wśród nich takie, które zawierają znaczny udział rytmu alfa?
    • Jaka jest ich topografia?
  • usunąć wszystkie komponenty nie zawierające alfy
  • odtworzyć z tych komponentów sygnał na elektrodach
  • wykonać dekompozycję ICA kilkukrotnie (co najmniej 3) i porównać wyniki
    • Czy uzyskiwane komponenty są powtarzalne?
    • Swoje wyniki porównać też z sąsiednimi grupami.

ZADANIE: Identyfikacja artefaktów

Proszę pobrać dane:

Pochodzą one z eksperymentu w którym osoba badana czytała słowa o różnych właściwościach wzbudzania emocji.

  • wczytaj je do eeglaba
  • wczytaj lokalizację kanałów z pliku Arousal-10-20-Cap.locs
  • obejrzyj przebiegi czasowe
  • odrzuć kanał z diodą (21) i z GSR (20)
  • zrób dekompozycję ICA
  • obejrzyj topografię komponentów
  • zidentyfikuj komponenty odpowiadające mruganiu i aktywności mięśniowej.
UWAGA
Aktualnie do wykrywania komponentów artefaktowych warto posłużyć się wtyczkami do eeglaba dostępnymi przez stronę:

https://sccn.ucsd.edu/eeglab/plugin_uploader/plugin_list_all.php

  • ICLabel
  • MARA

Eksperyment ASSR

W eksprymencie wykorzystujemy układ do generacji potencjałów słuchowych stanu ustalonego (ASSR). Wejście układu ASSR typu mini-jack wkładamy w wyjście słuchawkowe w laptopie. Drugie wejście układu ASSR wkładamy do wyjścia triggera we wzmacniaczu. Uruchamiamy plik dźwiękowy MM40tr.wav. Można go znalezc w: http://www.fuw.edu.pl/~suffa/LabEEG/MM40tr.wav

Stymulacja dźwiękowa składa sie z fali nośnej o częstości 400 Hz modulowanej z częstością 40 Hz. Plik dźwiękowy zawiera 5 sekund ciszy i 5 sekund stymulacji, powtórzone 40 razy.

Rejestracja sygnału

  1. Zakładamy czepek i elektrody w systemie 10-10, dbamy o to by opory pomiędzy elektrodami były poniżej 5 kΩ i różnice pomiędzy oporami różnych elektrod nie przekraczały 20%.
  2. Oklejamy kwadrat 3×3 elektrod na korze słuchowej z lewej strony (elektrody FT7, FC5, FC3, T7, C5, T3, TP7, CP5, CP3), 3×3 elektrod na korze słuchowej z prawej strony (elektrody FT8, FC6, FC4, T8, C6, T4, TP8, CP6, CP4), elektrody Fz, Cz, Pz i Oz, elektrody referencyjne A1 i A2. W sumie powinno być 24 elektrody.
  3. Elektrodę GND mocujemy na pozycji AFz.
  4. Sygnał rejestrujemy z częstością 2048 Hz.
  5. Do rejestracji stosujemy scenariusz 'ASSR' w interfejsie obci_gui.

Analiza

Początek stymulacji dźwiękowej oznaczymy jako 0. Poniższą analizę zastosuj dla sygnałów w referencji do uśrednionych odprowadzeń usznych A1 i A2. Wyznaczenie pasma częstości odpowiedzi ASSR

  1. Z sygnału wycinamy fragmenty od 0 do 5 sek. dla wszystkich elektrod położone nad korą słuchową (odcinki podczas stymulacji oraz bez niej).
  2. Dla każdej realizacji (odpowiedniego typu) obliczamy widma metodą Welcha.
  3. Otrzymane zespolone widma uśredniamy po realizacjach.
  4. Sprawdzamy czy w uśrednionym widmie występuję maksimum w częstości modulacji tj. 40 Hz i czy jest różnica między odcinkami ze stymulacją i bez niej.

Wyznaczenie przebiegu czasowego ERD i ERS

  1. Zaprojektuj filtry pasmowo przepustowe (Czebyszewa 2 rodzaju) zgodne z wyznaczonym pasmem. Zbadaj funkcje przenoszenia i odpowiedzi impulsowej.
  2. Powycinaj sygnały od −5 do +10 sekund (wszystkie kanały). Przefiltruj każdą realizację.
  3. Oblicz moc chwilową za pomocą transformaty Hilberta (kwadrat modułu transformaty Hilberta).
  4. Uśrednij moc chwilową po realizacjach.
  5. Oblicz względną zmianę mocy chwilowej względem czasu −4 do −2 s. W ten sposób otrzymasz przebieg ERD i ERS w czasie.
  6. Wykreśl ERD i ERS w układzie topograficznym. (Rozmieść subploty tak, aby z w przybliżeniu odpowiadały pozycjom elektrod).

Transformacja Hjortha

Transformacja Hjortha jest przybliżeniem numerycznym transformacji Laplace'a, czyli drugiej pochodnej przestrzennej. Obliczamy ją jako różnicę potencjału pomiędzy daną elektrodą i średnią z czterech sąsiednich elektrod. Przelicz potencjały z elektrod, w których występuję odpowiedź ASSR na montaż Hjortha i powtórz analizę opisaną powyżej.