Laboratorium EEG/CSP

Z Brain-wiki

Laboratorium_EEG/BSS

Plan zajęć

Materiał realizowany jest w ciągu 3 tygodni (6 spotkań po ~2,5 h).

Sesja Temat Kluczowe sekcje na tej stronie
1 Wykład BSS + CSP; ćwiczenie symulacyjne Ślepa separacja źródeł, Common Spatial Pattern, Ćwiczenie symulacyjne
2 Dane P300: wczytanie, cięcie epok, wizualizacja ERP, implementacja CSP Analiza wstępna, Analiza CSP
3 Mapki topograficzne; separacja cech Analiza CSP (cd.), Wybór i separacja cech
4 cosSinCSP jako uogólnienie CSP; dane SSVEP Filtry przestrzenne dla SSVEP
5 Wykład ICA; komponenty alfa ICA jako filtr przestrzenny, Komponenty alfa
6 Artefakty (ICLabel/MARA); synteza CSP/cosSinCSP/ICA Identyfikacja artefaktów

Sesja 1: Ślepa separacja źródeł i CSP — wprowadzenie

Czas
~60 min wykład + ~90 min ćwiczenie symulacyjne

Zajęcia wprowadzają do problemu filtracji przestrzennej sygnałów EEG. Punktem wyjścia jest ogólny problem ślepej separacji źródeł (BSS), z którego CSP wyłania się jako szczególny przypadek dla dwóch warunków eksperymentalnych. Filtry przestrzenne omawiane na tych zajęciach stanowią fundament klasycznych metod BCI (P300, SSVEP, motor imagery) oraz są bezpośrednim odpowiednikiem warstw przestrzennych w sieciach neuronowych takich jak ShallowConvNet i EEGNet — do tej analogii wrócimy w sesji 6.

Slajdy do wykładu (PDF)

Zakres wykładu
  • model generatywny EEG: [math]x(t) = As(t)[/math], macierz kowariancji, idea diagonalizacji
  • BSS jako ogólny problem separacji źródeł
  • CSP jako szczególny przypadek: dwa warunki eksperymentalne, iloraz Rayleigha, uogólnione zagadnienie własne
  • interpretacja filtrów i topografii źródeł
Ćwiczenie symulacyjne

Zob. Ćwiczenie symulacyjne w sekcji CSP poniżej.


Ćwiczenie symulacyjne

Zadania do samodzielnej realizacji

Poniższy kod dostarcza gotowej infrastruktury: generacji sygnałów, macierzy mieszającej i wizualizacji. Zadaniem jest samodzielna implementacja kluczowych kroków analizy CSP.

Zadanie 1 (obowiązkowe)
Obliczenie macierzy kowariancji

Napisz funkcję srednia_kowariancja(X, ind), która dla danej tablicy sygnałów X (wymiary: powtórzenia × kanały × próbki) i wektora indeksów czasowych ind zwraca średnią macierz kowariancji uśrednioną po powtórzeniach. Dla każdego powtórzenia zastosuj funkcję cov i znormalizuj wynik przez jego ślad (trace). Sprawdź wynik porównując z macierzami R_B i R_E obliczonymi w dostarczonym kodzie.

Zadanie 2 (obowiązkowe)
Implementacja CSP

Korzystając z funkcji eig oraz macierzy kowariancji obliczonych w zadaniu 1, wyznacz macierz filtrów przestrzennych W i odpowiadające im wartości własne. Odtwórz sygnały źródłowe dla wszystkich powtórzeń. Narysuj chmury punktów (amplituda źródła 1 vs amplituda źródła 2) osobno dla warunków baseline i ERP — przed transformacją CSP i po niej. Co zmieniło się w kształcie chmur?

Zadanie 3 (obowiązkowe)
Interpretacja wartości własnych

Wartości własne znajdują się na przekątnej macierzy Lambda. Który filtr najbardziej różnicuje warunki baseline i ERP? Jaki jest związek między wartością własną a separacją widoczną na wykresach punktowych?

Zadanie 4 (dla chętnych)
Wpływ macierzy mieszającej

W dostarczonym kodzie macierz P (filtr przestrzenny) wynosi: 

P = [1 2; 1.5 1.3]

Zmień ją kolejno na macierz bliską jednostkowej (np. P = [1 0.1; 0.1 1]) oraz na macierz z dużymi elementami pozadiagonalnymi (np. P = [1 5; 4 1]). Jak zmienia się kształt chmury punktów przed transformacją CSP? Czy CSP poprawnie oddziela źródła w obu przypadkach?

Zadanie 5 (dla chętnych)
Wpływ szumu

W kodzie szum jest wyłączony: n = 0*randn(size(tmp)). Usuń mnożnik 0* i zbadaj jak poziom szumu wpływa na separację źródeł. Przy jakim stosunku sygnału do szumu CSP przestaje skutecznie rozdzielać warunki?

Pytanie do dyskusji

Dostarczony kod używa eig(R_E, R_B). Alternatywą jest eig(R_E, (R_E+R_B)/2). Czym różnią się te dwie wersje? W jakich sytuacjach eksperymentalnych druga wersja mogłaby być bardziej uzasadniona?

Szkielet kodu do uzupełnienia

Poniższy szkielet zawiera gotową infrastrukturę symulacji. Uzupełnij brakujące fragmenty oznaczone komentarzem % TUTAJ zgodnie z zadaniami 1–3. Pełny kod referencyjny dostępny jest poniżej — zajrzyj do niego dopiero po samodzielnej próbie. 

%% Parametry symulacji — nie modyfikuj tej sekcji
Fs = 100;
T = 1;
t = 0:1/Fs:T-1/Fs;
N_rep = 100;
N_chan = 2;
s = zeros(N_rep, N_chan, length(t));
X = zeros(N_rep, N_chan, length(t));

P = [1 2; 1.5 1.3];   % filtr przestrzenny (prawdziwy)
A = P^(-1);            % topografia źródeł

for r = 1:N_rep
    s1 = sin(2*pi*11*t + pi/2) + 0.02*randn(size(t));
    s2 = exp(-((t-0.8)/0.05).^2) + 0.01*randn(size(t));
    s(r,1,:) = s1;
    s(r,2,:) = s2;
    tmp = squeeze(s(r,:,:));
    n = 0*randn(size(tmp));   % Zadanie 5: zmień 0* na 1* aby włączyć szum
    X(r,:,:) = A*tmp + n;
end

baseline_ind = find(t < 0.5);
ERP_ind      = find(t >= 0.5);

%% Zadanie 1: Oblicz średnie macierze kowariancji
% Dla każdego powtórzenia r:
%   - wytnij fragment sygnału X(r,:,baseline_ind) i analogicznie ERP
%   - oblicz macierz kowariancji funkcją cov() — pamiętaj o transpozycji
%   - znormalizuj przez ślad (funkcja trace())
%   - dodaj do sumy
% Na końcu podziel przez N_rep

R_B = zeros(N_chan, N_chan);
R_E = zeros(N_chan, N_chan);

for r = 1:N_rep
    B = squeeze(X(r,:,baseline_ind));
    % TUTAJ: oblicz kowariancję B, znormalizuj, dodaj do R_B

    E = squeeze(X(r,:,ERP_ind));
    % TUTAJ: oblicz kowariancję E, znormalizuj, dodaj do R_E
end

% TUTAJ: podziel R_B i R_E przez N_rep

%% Zadanie 2: Rozwiąż uogólnione zagadnienie własne
% Użyj funkcji eig(A,B) gdzie A=R_E, B=R_B
% W wyniku otrzymasz macierz wektorów własnych W (w kolumnach)
% i macierz Lambda z wartościami własnymi na przekątnej

% TUTAJ: [W, Lambda] = ...

disp('Wartości własne (przekątna Lambda):')
disp(diag(Lambda))

%% Zadanie 2 cd.: Odtwórz sygnały źródłowe
% Dla każdego powtórzenia r zastosuj filtr: S(r,:,:) = W' * X(r,:,:)

S = zeros(N_rep, N_chan, length(t));

for r = 1:N_rep
    % TUTAJ: S(r,:,:) = ...
end

%% Zadanie 3: Wizualizacja — chmury punktów przed i po CSP
x_base_k1 = X(:,1,baseline_ind);
x_base_k2 = X(:,2,baseline_ind);
x_ERP_k1  = X(:,1,ERP_ind);
x_ERP_k2  = X(:,2,ERP_ind);

s_base_k1 = squeeze(S(:,1,baseline_ind));
s_base_k2 = squeeze(S(:,2,baseline_ind));
s_ERP_k1  = squeeze(S(:,1,ERP_ind));
s_ERP_k2  = squeeze(S(:,2,ERP_ind));

figure(1); clf
subplot(1,2,1)
    plot(x_base_k1(:), x_base_k2(:), 'b.'); hold on
    plot(x_ERP_k1(:),  x_ERP_k2(:),  'r.');
    axis equal; xlim([-2 2]); ylim([-2 2])
    xlabel('Kanał 1'); ylabel('Kanał 2')
    title('Przed CSP')
    legend('baseline','ERP')

subplot(1,2,2)
    plot(s_base_k1(:), s_base_k2(:), 'b.'); hold on
    plot(s_ERP_k1(:),  s_ERP_k2(:),  'r.');
    axis equal
    xlabel('Źródło 1'); ylabel('Źródło 2')
    title('Po CSP')
    legend('baseline','ERP')

%% Zadanie 4 (dla chętnych): zmień macierz P powyżej i powtórz analizę
%% Zadanie 5 (dla chętnych): włącz szum (zmień 0* na 1* w linii z n=...)


Sesja 2: CSP dla P300 — dane realne

Czas
~150 min (jedno spotkanie)

Pracujemy z danymi kalibracyjnymi z eksperymentu P300. Celem sesji jest przejście pełnego pipeline'u: od surowych danych do odtworzonych źródeł CSP.

Zakres
  • wczytanie danych kalibracyjnych, identyfikacja znaczników T i NT
  • cięcie epok (−200 do +800 ms), usuwanie trendu liniowego
  • montaż względem połączonych uszu
  • wizualizacja ERP w układzie topograficznym
  • obliczenie macierzy kowariancji RT i RNT, normalizacja przez ślad
  • rozwiązanie uogólnionego zagadnienia własnego: [W, Lambda] = eig(R_T, R_NT)
  • odtworzenie sygnałów źródłowych: S = W'*EEG
  • przegląd estymowanych źródeł w dwóch warunkach
Materiały

Zob. Analiza wstępna i Analiza CSP poniżej.

Sesja 3: CSP dla P300 — interpretacja i separacja cech

Czas
~150 min (jedno spotkanie)

Kontynuacja analizy CSP. Celem sesji jest interpretacja wyników w kategoriach fizjologicznych oraz ocena separowalności klas w przestrzeni cech opartej na mocy sygnału.

Zakres
  • mapki topograficzne filtra przestrzennego i topografii źródła (funkcja topoplot z pakietu EEGLAB)
  • związek wartości własnych z siłą różnicowania warunków T i NT
  • wykresy punktowe mocy: oś X — moc składowej z największą wartością własną, oś Y — moc kolejnej składowej; jeden punkt = jedno powtórzenie
  • uśrednianie po 2, 4, 6, 8 i 10 realizacjach — obserwacja poprawy separacji wraz z liczbą uśrednień
  • dyskusja: utożsamianie komponentów CSP z fizjologicznymi źródłami jest ryzykowne; komponenty maksymalizują kontrast między warunkami, niekoniecznie odpowiadają izolowanym generatorom
Materiały

Zob. Analiza CSP i Wybór i separacja cech poniżej.


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żliwić wykonanie 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.

Sesja 4: Filtry przestrzenne dla SSEP — cosSinCSP

Czas
~150 min (jedno spotkanie)

Sesja wprowadza metodę cosSinCSP jako naturalne uogólnienie CSP o wiedzę dziedzinową dotyczącą sygnału SSEP. W klasycznym CSP dwa warunki eksperymentalne definiują dwie macierze kowariancji. W cosSinCSP rolę jednej z macierzy przejmuje projekcja sygnału na przestrzeń sinusów i cosinusów częstości stymulacji i jej harmonicznych. Formalizm ilorazu Rayleigha i uogólnionego zagadnienia własnego pozostaje identyczny.

Dane używane na tych zajęciach pochodzą z eksperymentu SSEP z archiwum — bez sesji rejestracyjnej. Na wykresach widma sygnału przed transformacją CSP dominuje artefakt sieciowy 50 Hz, który praktycznie zasłania odpowiedź na stymulację. Zadaniem jest sprawdzenie czy cosSinCSP potrafi ten sygnał wydobyć.

Zakres
  • idea cosSinCSP: macierz S złożona z sinusów i cosinusów harmonicznych, projekcja sygnału na przestrzeń SSEP i resztę Y
  • implementacja funkcji cosSinCSP w MATLAB
  • analiza danych archiwalnych: widma Welcha dla kanałów oryginalnych i źródeł CSP — porównanie SNR przed i po transformacji
  • prezentacja filtrów przestrzennych i topografii źródeł
Materiały

Zob. Filtry przestrzenne dla SSEP poniżej.

Filtry przestrzenne dla SSEP

Sesja 5: ICA — teoria i wydobywanie interesującego komponentu

Czas
~150 min (jedno spotkanie)

Sesja wprowadza analizę składowych niezależnych (ICA) jako trzecie podejście do problemu ślepej separacji źródeł — obok BSS/CSP i cosSinCSP. Kluczowa różnica: CSP szuka kierunków maksymalizujących kontrast wariancji między warunkami (kryterium drugiego rzędu), ICA szuka składowych niezależnych statystycznie, maksymalizując niegaussowość (kryterium wyższego rzędu). ICA nie wymaga dwóch warunków eksperymentalnych — działa na jednym stacjonarnym sygnale.

Zakres wykładu (~40 min)
  • model generatywny ICA: [math]\mathbf{x} = \mathbf{D}\mathbf{s}[/math], założenie niegaussowości składowych
  • negoentropia i algorytm FastICA
  • porównanie z CSP: kiedy ICA, kiedy CSP
  • ograniczenia: liczba próbek a liczba kanałów ([math]kN^2[/math]), niejednoznaczność kolejności i skali komponentów
Zakres ćwiczenia (~80 min)
  • wczytanie danych do EEGLAB, edycja lokalizacji elektrod
  • filtrowanie górnoprzepustowe (FIR, 0,5 Hz), zmiana referencji
  • obliczenie ICA na całym sygnale
  • identyfikacja komponentów zawierających rytm alfa: topografia, widmo, przebieg czasowy
  • usunięcie komponentów niezawierających alfy, odtworzenie sygnału na elektrodach
  • porównanie wyników między co najmniej trzema uruchomieniami ICA: czy komponenty są powtarzalne?
Materiały

Zob. ICA jako filtr przestrzenny i Wydobywanie interesujących komponentów poniżej.

ICA jako filtr przestrzenny

Sesja 6: ICA — artefakty i synteza metod

Czas
~150 min (jedno spotkanie)

Pierwsza część sesji poświęcona jest praktycznemu zastosowaniu ICA do usuwania artefaktów z danych EEG przy użyciu automatycznych klasyfikatorów komponentów (ICLabel, MARA). Druga część syntetyzuje całość materiału z bloku filtrów przestrzennych: BSS/CSP, cosSinCSP i ICA są trzema odpowiedziami na ten sam ogólny problem separacji źródeł, różniącymi się założeniami i kryterium optymalizacji.

Zakres ćwiczenia — artefakty (~75 min)
  • wczytanie danych Arousal do EEGLAB, usunięcie kanałów nieEEG
  • obliczenie ICA, przegląd topografii komponentów
  • identyfikacja artefaktów ocznych i mięśniowych
 przy użyciu wtyczek ICLabel i MARA
  • porównanie sygnału przed i po usunięciu komponentów artefaktowych
Zakres syntezy (~45 min)
  • tabela porównawcza CSP / cosSinCSP / ICA:
 założenia, dane wejściowe, kryterium optymalizacji, kiedy stosować
  • związek filtrów przestrzennych z architekturami sieci neuronowych:
 warstwa depthwise conv w ShallowConvNet i EEGNet jako
 wyuczony odpowiednik macierzy W z CSP —
 szczegóły w CSP a sieci neuronowe poniżej
  • omówienie raportów, pytania
Materiały

Zob. Identyfikacja artefaktów poniżej.

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.




Źródło: „http://brain.fuw.edu.pl/edu/index.php?title=Laboratorium_EEG/CSP&oldid=11725”