Laboratorium EEG/CSP: Różnice pomiędzy wersjami

Z Brain-wiki
 
(Nie pokazano 47 pośrednich wersji utworzonych przez tego samego użytkownika)
Linia 4: Linia 4:
 
[https://brain.fuw.edu.pl/edu/images/2/2f/BSS.pdf slajdy z prezentacji]
 
[https://brain.fuw.edu.pl/edu/images/2/2f/BSS.pdf slajdy z prezentacji]
  
=Ślepa separacja źródeł=
+
=Ślepa separacja źródeł (BSS)=
 +
{{hidden begin|title=Wstęp teoretyczny do BSS}}
 
Rozważmy ''N''-kanałowy sygnał EEG.
 
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.  
 
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.  
Linia 51: Linia 52:
 
: <math>J(w) = \frac{w^T C_T w}{w^T C_{NT} w}  </math>
 
: <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:
 
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}</math>
+
:<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>
ponieważ macierze kowariancji są symetryczne
 
::<math>\nabla J(w) =  \frac{  1}{w^T  C_{NT} w}\left[    C_{T} w+ C_{T}w  -\frac{w^T  C_{T} w}{w^T  C_{NT} w} \left(  C_{NT} w+ C_{NT}w \right) \right]</math>
 
::<math>= \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:
 
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>
 
:<math>      C_{T}w  =\frac{w^T  C_{T} w}{w^T  C_{NT} w}  C_{NT} w  </math>
Linia 60: Linia 58:
  
 
Aby znaleźć <math> \lambda</math> i <math>w</math> możemy wykorzystać w Matlabie funkcję <tt>eig</tt>. Funkcja ta rozwiązuje (również) uogólnione zagadnienia własne postaci ''Aw''=&lambda;''Bw'' 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.
 
Aby znaleźć <math> \lambda</math> i <math>w</math> możemy wykorzystać w Matlabie funkcję <tt>eig</tt>. Funkcja ta rozwiązuje (również) uogólnione zagadnienia własne postaci ''Aw''=&lambda;''Bw'' 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.
 +
{{hidden end}}
  
 
<!--
 
<!--
Linia 93: Linia 92:
 
-->
 
-->
  
===Ćwiczenie symulacyjne ===
+
====Ćwiczenie symulacyjne ====
 
{{hidden begin|title=kod przykładowy}}
 
{{hidden begin|title=kod przykładowy}}
 
<source lang  = matlab>
 
<source lang  = matlab>
Linia 188: Linia 187:
 
     title('ilustracja sytuacji pomiarowej -\newline znane są potencjały na elektrodach w dwóch warunkach eksperymentalnych')
 
     title('ilustracja sytuacji pomiarowej -\newline znane są potencjały na elektrodach w dwóch warunkach eksperymentalnych')
 
subplot(2,2,3);
 
subplot(2,2,3);
     plot(t(baseline_ind),(squeeze(X(:,1,baseline_ind)))','b'); hold on
+
     plot(t(baseline_ind),(squeeze(X(:,2,baseline_ind)))','b'); hold on
 
     plot(t(ERP_ind),(squeeze(  X(:,2,ERP_ind)))','r'); hold off
 
     plot(t(ERP_ind),(squeeze(  X(:,2,ERP_ind)))','r'); hold off
 
     xlabel('elektroda 2')
 
     xlabel('elektroda 2')
Linia 241: Linia 240:
  
 
==Zastosowanie filtra CSP do detekcji potencjału P300==
 
==Zastosowanie filtra CSP do detekcji potencjału P300==
 +
{{hidden begin|title=Eksperyment}}
 
===Eksperyment===
 
===Eksperyment===
 
* Proszę zapoznać się z instrukcją: http://laboratorium-eeg.braintech.pl/rozdz10.html
 
* Proszę zapoznać się z instrukcją: http://laboratorium-eeg.braintech.pl/rozdz10.html
Linia 287: Linia 287:
 
# Danych z kalibracji potrzebować będziemy kilka zestawów.  Proszę powtórzyć kilkukrotnie scenariusz &bdquo;kalibracjaP300&rdquo;. Przed każdym uruchomieniem trzeba zmienić string w pliku <tt>file_id_name</tt> np. na <tt>test???</tt> gdzie <tt>???</tt> oznacza kolejne numery.
 
# Danych z kalibracji potrzebować będziemy kilka zestawów.  Proszę powtórzyć kilkukrotnie scenariusz &bdquo;kalibracjaP300&rdquo;. Przed każdym uruchomieniem trzeba zmienić string w pliku <tt>file_id_name</tt> np. na <tt>test???</tt> gdzie <tt>???</tt> oznacza kolejne numery.
 
-->
 
-->
 
+
{{hidden end}}
 
===Analiza wstępna===
 
===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.
 
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 &mdash; &bdquo;target&rdquo; (związane z wystąpieniami litery &bdquo;B&rdquo;) i NT &mdash; &bdquo;non-target&rdquo; (pozostałe litery) o długości &minus;200 do +800 ms wokół triggerów. Dla każdej realizacji odjąć trend liniowy.
 
* Wczytać dane kalibracyjne do Matlaba i pociąć je na realizacje typu T &mdash; &bdquo;target&rdquo; (związane z wystąpieniami litery &bdquo;B&rdquo;) i NT &mdash; &bdquo;non-target&rdquo; (pozostałe litery) o długości &minus;200 do +800 ms wokół triggerów. Dla każdej realizacji odjąć trend liniowy.
* Sygnał zmontować wzgl. &bdquo;połączonych uszu&rdquo; i wyświetlić średnie przebiegi dla warunku T i NT w układzie topograficznym &mdash; wykorzystać w tym celu funkcję <tt>plottopo</tt> z pakietu Eeglab.
+
* Sygnał zmontować wzgl. &bdquo;połączonych uszu&rdquo; i wyświetlić średnie przebiegi dla warunku T i NT w układzie topograficznym &mdash; można wykorzystać w tym celu poniższy fragment kodu.
 +
{{hidden begin|title= fragment kodu do rysowania dwóch ERPów -- zamiast funkcji plottopo}}
 +
<source lang = matlab>
 +
%% stworzenie osi
 +
      % zakładam, że dane do rysowania są w dwóch strukturach EEG1 i EEG2 (struktury eeglab),
 +
      % te dane są podzielone na realizacje i mają kształt (kanały , czas, epoki)
 +
      % i że mamy już załadowane położenia elektrod we współrzędnych
 +
      % biegunowych
 +
      w = 0.125;
 +
      h = 0.125;
 +
      sc = 0.8;
 +
      figure()
 +
    for chanNum = 1:length(EEG1.chanlocs)
 +
        r =EEG1.chanlocs(chanNum).radius;
 +
        theta = EEG1.chanlocs(chanNum).theta;
 +
        x = r* sin(theta/180*pi)*sc+0.46;
 +
        y = r* cos(theta/180*pi)*sc+0.46;
 +
        ax(chanNum) = axes('Position', [x, y, w, h]);
 +
    end
 +
%%uśredniam EEG po powtórzeniach aby otrzymać ERP
 +
  ERP1 =squeeze( mean( EEG1.data,3)) ;
 +
  ERP2 =squeeze( mean( EEG2.data,3)) ;
 +
%% rysowanie uśrednionych potencjałów z dwóch struktur EEG
 +
for chanNum = 1:length(EEG1.chanlocs)
 +
    hold(ax(chanNum),'on');
 +
    plot(ax(chanNum), EEG1.times, ERP1(chanNum ,:));
 +
    plot(ax(chanNum), EEG2.times, ERP2(chanNum ,:));
 +
    title(ax(chanNum),EEG1.chanlocs(chanNum).labels );
 +
    hold(ax(chanNum),'off');
 +
end
 +
</source>
 +
{{hidden end}}
  
Poniżej zaprezentowany jest przykładowy skrypt do cięcia danych wokół znaczników. Działa on z plikami zawartymi w archiwum:
+
<small>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]]  
 
: [[Plik:KalibracjaP300.tar.gz]]  
 
Korzysta z funkcji pomocniczych dostępnych w dystrybucji obci w katalogu  
 
Korzysta z funkcji pomocniczych dostępnych w dystrybucji obci w katalogu  
 
: /usr/share/openbci/analysis/matlab_obci_signal_processing
 
: /usr/share/openbci/analysis/matlab_obci_signal_processing
Openbci można pobrać z https://github.com/BrainTech/openbci
+
Openbci można pobrać z https://github.com/BrainTech/openbci</small>
  
 
{{hidden begin|title=przykładowy skrypt}}
 
{{hidden begin|title=przykładowy skrypt}}
Linia 390: Linia 421:
 
{{hidden end}}
 
{{hidden end}}
  
===Analiza CSP===
+
===ZADANIE: Analiza CSP===
Link do Read menager [https://drive.google.com/open?id=0BzwQ_Lscn8yDS3RXNWdBbkxEQ2c]
 
  
* Korzystając z danych kalibracyjnych wykonać analizę CSP wzmacniającą potencjał P300.
+
Przegląd badań o P300: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2715154/
 +
 
 +
Link do Read menager [https://drive.google.com/file/d/1OgOduK5Zn7GgNl5XdCyLWXHJ7wTJIWcC/view?usp=share_link]
 +
<!---
 +
[https://drive.google.com/open?id=0BzwQ_Lscn8yDS3RXNWdBbkxEQ2c]
 +
--->
 +
* 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?
 
* 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?
  
Linia 399: Linia 435:
 
** filtrowi przestrzennemu  
 
** filtrowi przestrzennemu  
 
** rzutu topograficznego źródła na elektrody.
 
** rzutu topograficznego źródła na elektrody.
 +
<!---
 +
* Zbadać powtarzalność topografii pomiędzy plikami konfiguracyjnymi.
 +
--->
  
* Zbadać powtarzalność topografii pomiędzy plikami konfiguracyjnymi.
+
{{hidden begin|title=Wybór i separacja cech}}
  
 
===Wybór i separacja cech===
 
===Wybór i separacja cech===
Linia 410: Linia 449:
 
**  i sumujemy je w wybranym zakresie czasu.
 
**  i sumujemy je w wybranym zakresie czasu.
 
* Zaobserwuj jak zmienia się separacja w grupach target i non-target.
 
* Zaobserwuj jak zmienia się separacja w grupach target i non-target.
 +
{{hidden end}}
  
 
==Filtry przestrzenne dla SSEP ==
 
==Filtry przestrzenne dla SSEP ==
 
+
{{hidden begin|title=Filtry przestrzenne dla SSEP}}
 
=== Teoria===
 
=== Teoria===
 
Ciekawa koncepcja filtra przestrzennego dla SSVEP zaprezentowana jest  tu: http://www.eurasip.org/Proceedings/Eusipco/Eusipco2009/contents/papers/1569193209.pdf
 
Ciekawa koncepcja filtra przestrzennego dla SSVEP zaprezentowana jest  tu: http://www.eurasip.org/Proceedings/Eusipco/Eusipco2009/contents/papers/1569193209.pdf
Linia 483: Linia 523:
 
     PP=0;
 
     PP=0;
 
     for rep = 1:numberOfTrials
 
     for rep = 1:numberOfTrials
         x = signal(rep,i,:);
+
         x = squeeze(signal(rep,i,:));
 
         [Pxx,ff] = pwelch(x, X.sampling, 1, X.sampling, X.sampling);
 
         [Pxx,ff] = pwelch(x, X.sampling, 1, X.sampling, X.sampling);
 
         PP =PP + Pxx;
 
         PP =PP + Pxx;
Linia 495: Linia 535:
 
     PP=0;
 
     PP=0;
 
     for rep = 1:numberOfTrials
 
     for rep = 1:numberOfTrials
         s = S(rep,i,:);
+
         s = squeeze(S(rep,i,:));
 
         [Pss,ff]=pwelch(s, X.sampling, 1, X.sampling, X.sampling);
 
         [Pss,ff]=pwelch(s, X.sampling, 1, X.sampling, X.sampling);
 
         PP =PP + Pss;
 
         PP =PP + Pss;
Linia 506: Linia 546:
  
  
 +
====ZADANIE: Analiza danych z eksperymentu własnego ====
 +
# Przefiltruj sygnały EEG w paśmie 1-45 Hz za pomocą procedury <tt>filtfilt</tt>.
 +
# Na podstawie sygnału trigger oraz danych zapisanych w pliku wyodrębnij sygnały EEG zarejestrowane w trakcie stymulacji z odpowiednimi częstościami.
 +
# Uśrednij sygnały odpowiadające stymulacji tą samą częstością.
 +
# Obejrzyj uśrednione sygnały. (Zaprezentuj je).
 +
 +
#* Sposób I:
 +
#**Dla każdej realizacji wyestymuj przy pomocy metody Welcha widmo mocy sygnału EEG.
 +
#**Dla każdej częstości stymulacji wyznacz poziom tła na podstawie widm pochodzących ze stymulacji innymi częstościami. Np. dla stymulacji częstością 10 Hz poziom tła można wyznaczyć jako 95 centyl ze zbioru wartości widma w częstości 10 Hz dla stymulacji częstościami {4, 7, 13, 16, 20, 25, 30, 35, 40} Hz.
 +
#**Dla każdej częstości stymulacji wyznacz miarę odpowiedzi SSVEP (amplitudę widma odpowiadającą częstości stymulacji powyżej poziomy tła dla pojedynczej próby).
 +
#**Zaprezentuj widma otrzymane przy stymulacjach różnymi częstościami wraz z korytarzem odpowwiadajacym 95% przedziałowi ufności wyznaczonemu dla widm z pojedynczych realizacji.
 +
#**Sporządź wykres odpowiedzi SSVEP od częstości z zaznaczeniem przedziałów ufności i poziomu tła.
 +
#* Sposób II:
 +
#**Wyestumuj filtr CSP-SSVEP wspólny dla wszystkich częstości stymulacji
 +
#**Powtórz kroki ze sposobu I dla uzyskanych komponentów.
 +
#**Zaprezentuj filtry przestrzenne i topografie dla poszczególnych komponentów
 +
#**Dla przypomnienia:
 +
#***Filtr to zestaw współczynników z jakimi należy zsumować sygnały z poszczególnych kanałów EEG aby dostać komponenty odpowiadające hipotetycznym źródłom nieskorelowanym. Filtr można zilustrować na głowie, przypisując poszczególnym pozycjom elektrod wagi równe współrzędnym wektora w (kolumna macierzy W).
 +
#***Topografia źródła to zestaw współczynników z jakimi docierają one do poszczególnych kanałów EEG. Topografia zawarta jest w wierszach macierzy odwrotnej do W.
 +
#**Można rysunki wykonać jako macierze 5x5 i w pozycji elektrody kolorujemy proporcjonalnie do współczynnika.
 +
 +
===SSVEP-BCI===
 +
W zajęciach tych przydadzą nam się informacje z:
 +
https://brain.fuw.edu.pl/edu/index.php/Laboratorium_EEG/Wprowadzenie_do_syg_online
 +
====Wstęp====
 +
Naszym celem będzie stworzenie prostego interfejsu wykorzystującego zjawisko SSVEP.
 +
Po dwóch stronach monitora zamocujemy diody. Każda będzie migać ze swoją ustaloną częstością (np. 16 i 22 Hz) - warto zadbać aby nie były to częstości powiązane ze sobą harmonicznie, a z drugiej strony aby, biorąc pod uwagę wyniki poprzedniego zadania, dawały dobra odpowiedź SSVEP.
 +
 +
Eksperyment będzie miał dwie cześci: sesję kalibracje i sesję on-line. Pomiędzy tymi sesjami będziemy uczyć kalsyfikator. Na podstawie części kalibracyjnej ustalimy jakie parametry przetwarznego on-line sygnału świadczą o patrzeniu się na diodę z lewej a jakie na tą z prawej strony ekranu.
 +
 +
Potem w sesji online będziemy porównywać (za pomocą predykcji klasyfikatora) rejestrowany sygnał z tymi warościami kalibracyjnymi i na tej podstawie system będzie zwracał informację o wyborze lewej lub prawej diody. To powinno umożliwić prostą komunikację na zasadzie pytanie i odpowiedź TAK/NIE.
 +
 +
==== Sesja kalibracycjna====
 +
* Zakładamy czepek
 +
* Częstość próbkowania ustawiamy na 256Hz
 +
* Na podstawie kodu: https://brain.fuw.edu.pl/edu/index.php/Laboratorium_EEG/Wprowadzenie_do_syg_online#.C4.86wiczenie:_Wykorzystanie_pomiaru_EMG_do_sterowania_on-line
 +
tworzymy programik, który zamiast pętli while tworzy odpowiednie pętle for aby:
 +
* pięciokrotnie zarejestrować sekwencję trzech warunków eksperymentalnych:
 +
** patrz 5s na diodę z lewej (warunek l)
 +
** patrz 5s na środek ekranu (warunek s)
 +
** patrz 5s na diodę z prawej (warunek p)
 +
* w czasie tego patrzenia:
 +
** odbieramy próbki w pakietach o długości 0.5s
 +
** każdy pakiet filtrujem (technika filtrowania on-line) w pasmach dookoła wybranych dwóch częstości (PASMO_LEWE / PASMO_PRAWE)
 +
** przefiltrowany pakiet przeliczamy na RMS
 +
** zbieramy sześć zbiorów wyników -  zestawy RMSów związane z każdym z waunków kalibracyjnych (LEWY/SPOCZYNEK/PRAWY) dla PASMO_LEWE i PASMO_PRAWE.
 +
* Normalizujemy RMSy. Dla pasma lewego obliczamy (RMS_(l/p/s) - np.mean(RMS_s)) / np.sdt(RMS_s) i analogicznie dla pasma prawego. Małe indeksy oznaczaają tu warunki. Zapamiętujemy współczynniki normalizacyjne.
 +
** Ogladamy rozkłady uzyskanych wielkości (znormalizowanych RMSów).
 +
** Uczymy klasyfikator np. regresję logistyczną (https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html) do rozpoznawania, z którym warunkiem patrzenia mamy do czynienia.
 +
 +
==== Sesja online====
 +
* wczytujemy wyuczony model i kalibracyjne wartości np.mean(RMS_s)) i np.sdt(RMS_s) dla pasm lewego i prawego
 +
* Odbieramy próbki w pakietach o długości 0.5s
 +
** każdy pakiet filtrujem (technika filtrowania on-line) w pasmach dookoła wybranych dwóch częstości (PASMO_LEWE / PASMO_PRAWE)
 +
** przefiltrowany pakiet przeliczamy na RMS
 +
** Normalizujemy RMSy. Dla pasma lewego obliczamy (RMS_(l/p/s) - np.mean(RMS_s)) / np.sdt(RMS_s) i analogicznie dla pasma prawego. Małe indeksy oznaczaają tu warunki.
 +
** robimy predykcję klasyfikatora, z którym warunkiem patrzenia mamy do czynienia. Wyświetlamy na ekranie komunikat.
 +
* Testujemy czy powyższy schemat analizy pozwala na komunikację.
 +
{{hidden end}}
 
<!--
 
<!--
 
===Eksperyment ASSR===
 
===Eksperyment ASSR===
Linia 543: Linia 642:
  
 
==ICA jako filtr przestrzenny==
 
==ICA jako filtr przestrzenny==
 +
{{hidden begin|title=Wstęp teoretyczny do ICA}}
 
===Definicja ===
 
===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.
 
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.
Linia 595: Linia 695:
  
 
===Obliczenia===
 
===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. Procedura wykorzystywana w eeglabie (&bdquo;runica&rdquo;, 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.
+
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 (&bdquo;runica&rdquo;, 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.
 
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.
  
Linia 608: Linia 731:
  
 
===Bibliografia===
 
===Bibliografia===
 +
Bazowa praca:
 +
* A. Hyvärinen. Fast and Robust Fixed-Point Algorithms for Independent Component Analysis. IEEE Transactions on Neural Networks 10(3):626-634, 1999 http://www.cs.helsinki.fi/u/ahyvarin/papers/TNN99_reprint.pdf
 +
 +
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
 +
 
* Grau, C., Fuentemilla, L., Marco-Pallars, J. (2007). Functional neural dynamics underlying auditory event-related n1 and n1 suppression response. Neuroimage, 36(6):522–31.
 
* Grau, C., Fuentemilla, L., Marco-Pallars, J. (2007). Functional neural dynamics underlying auditory event-related n1 and n1 suppression response. Neuroimage, 36(6):522–31.
* Hyvärinen, A. and Oja, E. (2000). Independent component analysis: Algorithms and applications. Neural Networks, 13(4-5):411–430.
 
 
* Makeig, S., Bell, A., Jung, T.-P., Sejnowski,T. (1996). Independent component analysis of electroencephalographic data. W: Touretzky, D., Mozer, M., and Hasselmo, M., editors, Advances in Neural Information Processing Systems, volume 8, pages 145–151. MIT Press, Cambridge, MA.
 
* Makeig, S., Bell, A., Jung, T.-P., Sejnowski,T. (1996). Independent component analysis of electroencephalographic data. W: Touretzky, D., Mozer, M., and Hasselmo, M., editors, Advances in Neural Information Processing Systems, volume 8, pages 145–151. MIT Press, Cambridge, MA.
 
* Onton, J., Makeig, S. (2006). Information-based modeling of event-related brain dynamics. Prog Brain Res., 159:99–120.
 
* Onton, J., Makeig, S. (2006). Information-based modeling of event-related brain dynamics. Prog Brain Res., 159:99–120.
Linia 615: Linia 744:
 
* http://sccn.ucsd.edu/~arno/indexica.html
 
* http://sccn.ucsd.edu/~arno/indexica.html
 
* http://cis.legacy.ics.tkk.fi/aapo/papers/IJCNN99_tutorialweb/
 
* http://cis.legacy.ics.tkk.fi/aapo/papers/IJCNN99_tutorialweb/
 +
{{hidden end}}
 +
{{hidden begin|title=Wydobywanie interesujących komponentów}}
  
=== Wydobywanie interesujących komponentów ===
+
=== ZADANIE: Wydobywanie interesujących komponentów ===
  
 
Dane do tej części ćwiczeń proszę pobrać i rozpakować w swoim katalogu:
 
Dane do tej części ćwiczeń proszę pobrać i rozpakować w swoim katalogu:
Linia 639: Linia 770:
 
** Czy uzyskiwane komponenty są powtarzalne?  
 
** Czy uzyskiwane komponenty są powtarzalne?  
 
** Swoje wyniki porównać też z sąsiednimi grupami.
 
** Swoje wyniki porównać też z sąsiednimi grupami.
 +
{{hidden end}}
  
=== Identyfikacja artefaktów ===
+
===ZADANIE: Identyfikacja artefaktów ===
 
Proszę pobrać dane:
 
Proszę pobrać dane:
  
Linia 656: Linia 788:
 
* obejrzyj topografię komponentów
 
* obejrzyj topografię komponentów
 
* zidentyfikuj komponenty odpowiadające mruganiu i aktywności mięśniowej.
 
* 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
 +
 +
====W raporcie: ====
 +
* zaprezentuj fragmenty sygnału zawierającego artefakty oczne i mięśniowe przed i po zastosowaniu czyszczenia poprzez usuwanie komponentów zdominowanych przez artefakty.
 +
* zaprezentuj topografię i przebiegi czasowe komponentów zidentyfikowanych jako artefakty oczne i mięśniowe.
 +
 +
  
 
<!--
 
<!--
Linia 663: Linia 807:
 
-->
 
-->
  
 +
<!---
 
==Eksperyment ASSR==
 
==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
 
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
Linia 694: Linia 839:
 
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.
 
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.
 
Przelicz potencjały z elektrod, w których występuję odpowiedź ASSR na montaż Hjortha i powtórz analizę opisaną powyżej.
 +
--->

Aktualna wersja na dzień 08:47, 14 maj 2024

Laboratorium_EEG/BSS

Prezentacja

slajdy z prezentacji

Ślepa separacja źródeł (BSS)


Ć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 — można wykorzystać w tym celu poniższy fragment kodu.

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

Przegląd badań o P300: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2715154/

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

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

W raporcie:

  • zaprezentuj fragmenty sygnału zawierającego artefakty oczne i mięśniowe przed i po zastosowaniu czyszczenia poprzez usuwanie komponentów zdominowanych przez artefakty.
  • zaprezentuj topografię i przebiegi czasowe komponentów zidentyfikowanych jako artefakty oczne i mięśniowe.