Laboratorium EEG/CSP

Z Brain-wiki

Laboratorium_EEG/BSS

Spis treści

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?

Wykreśl także wizualizację przebiegów czasowych sygnałów X i sygnałów S uzyskanych po rozmieszaniu.

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);

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

%% Zadanie 1: Oblicz średnie macierze kowariancji
% korzystając z funkcji srednia_kowariancja(X, ind)

% 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







%% 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')
%% Wizualizacja przebiegów czasowych X i S
% TUTAJ:

%% 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 średniej ze wszystkich elektrod
  • 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

Dane: https://www.fuw.edu.pl/~jarekz/LabEEG/KalibracjaP300.tar.gz

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 wektorów odpowiadających:
    • filtrowi przestrzennemu
    • rzutu topograficznego źródła na elektrody.

Szkielet kodu do analizy wstępnej

%% KROK 1: Wczytanie danych
% ReadManager to pomocnicza klasa do odczytu plików w formacie OpenBCI.
% Przyjmuje trzy pliki: XML z metadanymi, RAW z próbkami i TAG ze znacznikami zdarzeń.
nazwaPliku = 'p_6301423_calibration_p300.obci';
nameOfXMLFile  = strcat(nazwaPliku, '.xml');
nameOfTagFile  = strcat(nazwaPliku, '.tag');
namesOfDataFiles = strcat(nazwaPliku, '.raw');

rm = ReadManager(nameOfXMLFile, namesOfDataFiles, nameOfTagFile);

% Pobieramy podstawowe parametry nagrania
numberOfChannels  = rm.get_param('number_of_channels');
namesOfChannels   = rm.get_param('channels_names');
samplingFrequency = rm.get_param('sampling_frequency');

% Pobieramy listę wszystkich znaczników zdarzeń z nagrania
tagsStruct = rm.get_tags();

%% KROK 2: Identyfikacja znaczników Target i Non-Target
% Każdy znacznik typu 'blink' odpowiada jednemu mignięciu w eksperymencie P300.
% Pole 'index' mówi który symbol mignął, pole 'target' który symbol był oczekiwany.
% Jeśli index == target to było to mignięcie docelowe (Target), w przeciwnym razie
% Non-Target.
targetTimeStamps    = [];
NonTargetTimeStamps = [];

for structNumber = 1:length(tagsStruct)
    if strcmp(tagsStruct(structNumber).name, 'blink')
        index  = tagsStruct(structNumber).children.index;
        target = tagsStruct(structNumber).children.target;
        if index == target
            targetTimeStamps    = [targetTimeStamps    tagsStruct(structNumber).start_timestamp];
        else
            NonTargetTimeStamps = [NonTargetTimeStamps tagsStruct(structNumber).start_timestamp];
        end
    end
end

%% KROK 3: Wczytanie i filtrowanie sygnału
samples = double(rm.get_samples());
samples = samples(1:8,:);   % odrzucamy kanały bez EEG
numberOfChannels = 8;

% Filtr dolnoprzepustowy Czebyszewa 2. rodzaju — odcięcie 25 Hz
[b, a] = cheby2(6, 80, 25/(samplingFrequency/2), 'low');
for ch = 1:numberOfChannels
    samples(ch,:) = filtfilt(b, a, samples(ch,:));
end

%% KROK 4: Montaż względem średniej (common average reference)
% Odejmujemy od każdego kanału średnią ze wszystkich kanałów w danej chwili.
% Jest to prostsze i bardziej ogólne niż montaż względem połączonych uszu.
M = -ones(numberOfChannels, numberOfChannels) / numberOfChannels;
M = M + eye(numberOfChannels) * (numberOfChannels+1) / numberOfChannels;
samples = 0.0715 * M * samples;

%% KROK 5: Cięcie epok
PRE  = -0.2;
POST =  0.8;
wycinek = floor(PRE*samplingFrequency : POST*samplingFrequency);

TargetSignal    = zeros(length(targetTimeStamps),    numberOfChannels, length(wycinek));
NonTargetSignal = zeros(length(NonTargetTimeStamps), numberOfChannels, length(wycinek));

for trialNumber = 1:length(targetTimeStamps)
    trigerOnset = floor(targetTimeStamps(trialNumber) * samplingFrequency);
    tenWycinek  = wycinek + trigerOnset;
    if tenWycinek(1) > 0 && tenWycinek(end) <= size(samples,2)
        tmpSignal = samples(:, tenWycinek);
        tmpSignal = detrend(tmpSignal')';   % usuwanie trendu liniowego
        TargetSignal(trialNumber,:,:) = tmpSignal;
    end
end

for trialNumber = 1:length(NonTargetTimeStamps)
    trigerOnset = floor(NonTargetTimeStamps(trialNumber) * samplingFrequency);
    tenWycinek  = wycinek + trigerOnset;
    if tenWycinek(1) > 0 && tenWycinek(end) <= size(samples,2)
        tmpSignal = samples(:, tenWycinek);
        tmpSignal = detrend(tmpSignal')';
        NonTargetSignal(trialNumber,:,:) = tmpSignal;
    end
end

% Oś czasu w milisekundach — przyda się do rysowania
times = (wycinek / samplingFrequency) * 1000;

%% KROK 6: Wizualizacja ERP w układzie topograficznym
% Oblicz potencjały wywołane (ERP) jako średnie po powtórzeniach
ERP_T  = squeeze(mean(TargetSignal,    1));   % wymiary: kanały × czas
ERP_NT = squeeze(mean(NonTargetSignal, 1));

% TUTAJ: narysuj ERP_T i ERP_NT dla każdego kanału w subplotach
% ułożonych w siatce odpowiadającej przybliżonym pozycjom elektrod.
% Na każdym subplocie nałóż oba warunki różnymi kolorami.
% Oś X: czas w ms (użyj wektora times), oś Y: amplituda w µV.
% Zadbaj o legendę i tytuły z nazwami kanałów (namesOfChannels).

%% KROK 7: Obliczenie średnich macierzy kowariancji
% Dla każdego powtórzenia oblicz macierz kowariancji funkcją cov(),
% znormalizuj przez jej ślad (trace()), a następnie uśrednij po powtórzeniach.
% Pamiętaj: cov() oczekuje macierzy próbki × kanały, czyli potrzebujesz transpozycji.

R_T  = zeros(numberOfChannels, numberOfChannels);
R_NT = zeros(numberOfChannels, numberOfChannels);

for r = 1:length(targetTimeStamps)
    % TUTAJ: wytnij r-te powtórzenie T, oblicz kowariancję, znormalizuj, dodaj do R_T
end

for r = 1:length(NonTargetTimeStamps)
    % TUTAJ: wytnij r-te powtórzenie NT, oblicz kowariancję, znormalizuj, dodaj do R_NT
end

% TUTAJ: podziel R_T i R_NT przez odpowiednie liczby powtórzeń

%% KROK 8: Rozwiązanie uogólnionego zagadnienia własnego
% Szukamy filtrów przestrzennych W maksymalizujących stosunek mocy
% w warunku T względem NT: eig(R_T, R_NT)
% W kolumnach W znajdują się filtry, na przekątnej Lambda — wartości własne.

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

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

%% KROK 9: Odtworzenie sygnałów źródłowych
% Dla każdego powtórzenia zastosuj filtr przestrzenny:
% s(r,:,:) = W' * x(r,:,:)
% Pamiętaj o wymiarach — squeeze() może być pomocne.

S_T  = zeros(size(TargetSignal));
S_NT = zeros(size(NonTargetSignal));

for r = 1:size(TargetSignal, 1)
    % TUTAJ: S_T(r,:,:) = ...
end

for r = 1:size(NonTargetSignal, 1)
    % TUTAJ: S_NT(r,:,:) = ...
end

%% KROK 10: Wizualizacja estymowanych źródeł
% Narysuj średnie przebiegi źródeł (ERP w przestrzeni CSP) dla obu warunków.
% Ułóż subploty w prostokątnej siatce — jeden subplot na każde źródło.
% Dla którego źródła różnica między T i NT jest największa?
% Czy odpowiada to źródłu z największą wartością własną?

ERP_S_T  = squeeze(mean(S_T,  1));
ERP_S_NT = squeeze(mean(S_NT, 1));

% TUTAJ: narysuj ERP_S_T i ERP_S_NT w subplotach

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
  • 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. Szkielet kodu do analizy wstępnej i Zadanie: Wybór i separacja cech poniżej.



Zadanie: Wybór i separacja cech

  • Przedstaw na rysunkach nałożone na siebie pojedyncze realizacje z warunków target i non-target po rzutowaniu na wektor [math]w[/math] odpowiadający największej i kolejnej wartości własnej.
  • Przedstaw wykresy punktowe takie, że na jednej osi jest moc sygnału (suma kwadratów wartości próbek w wybranym zakresie czasu) odpowiadającego największej wartości własnej, a na drugiej osi kolejnej (mniejszej) wartości własnej; jeden punkt reprezentuje jedno powtórzenie.
  • Wykonaj serię wykresów jak w poprzednim punkcie dla uśrednień sygnałów kolejno po 2, 4, 6, 8 i 10 realizacjach:
    • Liczymy potencjał wywołany dla danej liczby powtórzeń.
    • Następnie podnosimy wartości próbek do kwadratu
    • i sumujemy je w wybranym zakresie czasu.
  • Zaobserwuj jak zmienia się separacja w grupach target i non-target.

Archirtektury sieci do BCI: ConvNet

EEGNet

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

Teoria

Ciekawa koncepcja filtra przestrzennego dla SSEP zaprezentowana jest tu: http://www.eurasip.org/Proceedings/Eusipco/Eusipco2009/contents/papers/1569193209.pdf

Pokrótce można ją rozumieć podobnie do tego co robiliśmy rozważając filtry przestrzenne CSP z tym, że dla SSEP oraz innych potencjałów wywołanych stanu ustalonego możemy skorzystać z dodatkowych informacji dotyczących poszukiwanych źródeł. Wiemy mianowicie, że powinny one oscylować z częstością bodźca i być może jej harmonicznych.

Przyda nam się macierz [math]S[/math] zbudowana tak, że w kolejnych kolumnach znajdują się sinusy i cosinusy kolejnych częstości harmonicznych. Wektory te unormujemy, żeby miały energię równą 1. Innymi słowy macierz [math]S[/math] zbudowana jest z wersorów rozpinających przestrzeń, w której powinien znajdować się sygnał SSEP.

W Matlabie możemy taką macierz zbudować tak: 

% Fs                - częstość próbkowania
% numberOfSamples   - długość sygnału w próbkach
% numberOfHarmonics - liczba harmonicznych, które chcemy włączyć do analizy
t = (0:1:numberOfSamples - 1) / Fs;
S = zeros(numberOfSamples, 2*numberOfHarmonics);

for harmonicNumber = 1:numberOfHarmonics
    c = cos(2*pi*stimulationFrequency*harmonicNumber*t);
    s = sin(2*pi*stimulationFrequency*harmonicNumber*t);
    S(:,(harmonicNumber - 1)*2 + 1) = c/norm(c);
    S(:,(harmonicNumber - 1)*2 + 2) = s/norm(s);
end

Aby w badanym sygnale znaleźć składowe odpowiadające SSEP musimy rzutować sygnał [math]X[/math] (macierz sygnałów kanały × próbki) na przestrzeń rozpiętą przez [math]S[/math]:

[math]A = X S[/math]

Macierz [math]A[/math] zawiera współczynniki będące iloczynami skalarnymi sygnałów i wersorów — mówią one o tym jak dużo sinusa bądź cosinusa o danej częstości jest w pierwotnym sygnale. Komponenty SSEP zawarte w sygnale [math]X[/math] odzyskujemy tak:

[math]\mathrm{SSEP} = A S^T[/math]

Modelujemy rejestrowany sygnał jako:

[math]X = \mathrm{SSEP} + Y[/math]

gdzie:

[math]Y = X - \mathrm{SSEP}[/math]

to wszystkie komponenty sygnału, które nas nie interesują.

Filtr przestrzenny, który chcemy zbudować powinien maksymalizować stosunek wariancji [math]\mathrm{SSEP}[/math] do wariancji [math]Y[/math]. Macierz kowariancji powinna być uśredniona po powtórzeniach, a kowariancja sygnału w każdym powtórzeniu powinna być znormalizowana przez jej ślad. Dalej stosujemy technikę znaną z CSP: maksymalizację ilorazu Rayleigha przez rozwiązanie uogólnionego zagadnienia własnego dla macierzy kowariancji [math]\mathrm{SSEP}[/math] i [math]Y[/math].

Zadanie: implementacja funkcji cosSinCSP

Dane do zadania: Plik:PrzykladoweDaneSSVEP.mat.gz. Pełny przykład z danymi z eksperymentu SSEP: Plik:SSVEP demo csp.tar.gz

Zaimplementuj funkcję cosSinCSP zgodnie z poniższym szkieletem. Prawidłowo zaimplementowana funkcja wraz ze skryptem demonstracyjnym poniżej powinna generować rysunek: podpis grafiki 

function [W, Lambda] = cosSinCSP(signal, stimulationFrequency, numberOfHarmonics, Fs)
% cosSinCSP - filtr przestrzenny dla sygnałów SSEP
%
% Wejście:
%   signal                - dane EEG, wymiary: powtórzenia × kanały × próbki
%   stimulationFrequency  - częstość stymulacji w Hz
%   numberOfHarmonics     - liczba harmonicznych do uwzględnienia
%   Fs                    - częstość próbkowania w Hz
%
% Wyjście:
%   W      - macierz filtrów przestrzennych (kanały × kanały), filtry w kolumnach
%   Lambda - macierz diagonalna wartości własnych

[numberOfTrials, numberOfChannels, numberOfSamples] = size(signal);

%% KROK 1: Zbuduj macierz referencyjną S (sinusy i cosinusy harmonicznych)
% S ma wymiary: próbki × 2*numberOfHarmonics
% Każda kolumna to jeden wersor (unormowany sinus lub cosinus).
% Wróć do sekcji Teoria powyżej i przełóż tamten kod na tę funkcję.

% TUTAJ

%% KROK 2: Oblicz średnie macierze kowariancji C_SSEP i C_Y
% Dla każdego powtórzenia wytnij sygnał x (kanały × próbki),
% rozłóż go na składową SSEP i resztę Y zgodnie z równaniami z sekcji Teoria,
% oblicz macierze kowariancji obu składowych, znormalizuj przez ślad
% i uśrednij po powtórzeniach.

C_SSEP = zeros(numberOfChannels, numberOfChannels);
C_Y    = zeros(numberOfChannels, numberOfChannels);

for trial = 1:numberOfTrials
    x = squeeze(signal(trial, :, :));   % kanały × próbki

    % TUTAJ: oblicz SSEP i Y

    % TUTAJ: oblicz kowariancje, znormalizuj, dodaj do C_SSEP i C_Y
end

% TUTAJ: podziel C_SSEP i C_Y przez numberOfTrials

%% KROK 3: Rozwiąż uogólnione zagadnienie własne
% Znajdź filtry przestrzenne maksymalizujące stosunek mocy SSEP do mocy Y.
% Analogia z CSP: jakie macierze wstawiłeś do eig() tam?

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

end
Wskazówka do interpretacji wyników

Filtr związany z największą wartością własną daje komponent o największym stosunku mocy SSEP do mocy tła. Sprawdź na widmach Welcha czy odpowiada on składowej z wyraźnym pikiem przy częstości stymulacji.

Skrypt demonstracyjny i analiza wyników

% wczytujemy dane
load('PrzykladoweDaneSSVEP.mat');

[numberOfTrials numberOfChannels numberOfSamples] = size(X.data);
namesOfChannels = X.channels;

numberOfHarmonics = 3;
signal = X.data;   % powtórzenie × kanał × próbki

[W, Lambda] = cosSinCSP(signal, X.stimulation, numberOfHarmonics, X.sampling);

%% Odtworzenie sygnałów źródłowych
S = zeros(size(signal));
for powt = 1:size(signal,1)
    S(powt,:,:) = W' * squeeze(signal(powt,:,:));
end

%% Widma Welcha: kanały oryginalne vs źródła CSP
figure('Name', ['Stymulacja: ', num2str(X.stimulation), ' Hz'])
for i = 1:numberOfChannels
    subplot(2, numberOfChannels, i)
    PP = 0;
    for rep = 1:numberOfTrials
        x = squeeze(signal(rep,i,:));
        [Pxx, ff] = pwelch(x, X.sampling, 1, X.sampling, X.sampling);
        PP = PP + Pxx;
    end
    plot(ff(ff<60), PP(ff<60))
    title(namesOfChannels{i})

    subplot(2, numberOfChannels, numberOfChannels+i)
    PP = 0;
    for rep = 1:numberOfTrials
        s = squeeze(S(rep,i,:));
        [Pss, ff] = pwelch(s, X.sampling, 1, X.sampling, X.sampling);
        PP = PP + Pss;
    end
    plot(ff(ff<60), PP(ff<60))
    title(['źródło CSP: ', num2str(i)])
    xlabel(['\lambda = ', num2str(Lambda(i,i), '%.2f')])
end

%% KROK 4: Mapki topograficzne
% Zidentyfikuj komponent z największą wartością własną.
% Narysuj dwie mapki przy użyciu funkcji topoplot() z pakietu EEGLAB:
%
% (a) Filtr przestrzenny: wagi z jakimi kanały EEG składają się na ten komponent.
%     Wskazówka: filtr to kolumna macierzy W.
%
% (b) Topografia źródła: z jakimi wagami komponent dociera do poszczególnych elektrod.
%     Wskazówka: topografia zawarta jest w wierszach macierzy odwrotnej do W —
%     przypomnij sobie analogiczny krok z analizy CSP dla P300.
%
% Czym różnią się te dwie mapki? Którą łatwiej interpretować fizjologicznie?

% TUTAJ: [~, idx] = max(diag(Lambda));
% TUTAJ: topoplot(...)  % filtr
% TUTAJ: topoplot(...)  % topografia

Sesja 5: ICA jako filtr przestrzenny

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?

Teoria

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)2] = 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

w

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)\}]^2[/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 Newtona. (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

Sesja 6: Problem odwrotny w EEG i MEG

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

Sesja wprowadza problem odwrotny jako alternatywne — i komplementarne — podejście do analizy przestrzennej sygnałów EEG i MEG. O ile BSS/CSP i ICA szukają filtrów przestrzennych w przestrzeni elektrod, nie zakładając nic o geometrii mózgu, to metody problemu odwrotnego wyraźnie osadzają źródła w przestrzeni trójwymiarowej, korzystając z fizyki i modelu głowy.

Materiały z naszych zasobów:

Z zasobów zewnętrznych polecam:



Kontrast z BSS: dwa języki opisu tego samego zjawiska

Zarówno BSS/ICA, jak i metody problemu odwrotnego wychodzą od tego samego modelu generatywnego:

[math]\mathbf{x}(t) = \mathbf{A}\,\mathbf{s}(t) + \boldsymbol{\varepsilon}(t)[/math]

gdzie [math]\mathbf{x}(t) \in \mathbb{R}^{N_e}[/math] to potencjały na elektrodach, [math]\mathbf{s}(t) \in \mathbb{R}^{N_s}[/math] to aktywność źródeł, a [math]\mathbf{A}[/math] to macierz mieszająca (ang. lead field lub forward matrix).

Różnica leży w tym, co robimy z macierzą [math]\mathbf{A}[/math]:

BSS / ICA Problem odwrotny
Macierz A nieznana; estymowana z danych obliczona z modelu głowy (forward model)
Liczba źródeł równa liczbie elektrod od 1 dipola do ~10 000 wokseli
Geometria brak założeń przestrzennych źródła zlokalizowane w przestrzeni 3D
Regularyzacja przez kryterium stat. (niezależność, wariancja) przez regularyzację matematyczną (L2, L1, …)
Wynik filtry przestrzenne + komponenty mapa aktywności w przestrzeni mózgu

Kluczowa trudność problemu odwrotnego wynika z jego fundamentalnej niejednoznaczności:

układ równań:

[math]\mathbf{x} = \mathbf{A}\mathbf{s}[/math]

jest niedookreślony — przy [math]N_e \approx 19\text{–}256[/math] elektrodach i [math]N_s \sim 10^4[/math] możliwych lokalizacjach źródeł istnieje nieskończenie wiele rozkładów [math]\mathbf{s}[/math] dokładnie odtwarzających zmierzony [math]\mathbf{x}[/math]. Każda metoda nakłada inne założenie a priori, by wybrać jedno rozwiązanie.

Model do przodu (forward model)

Zanim można rozwiązać problem odwrotny, należy zbudować model do przodu — obliczenie potencjałów na elektrodach wywołanych przez jednostkowe źródło prądowe w danym miejscu mózgu. Wynik zapisuje się jako macierz pola wiodącego (lead field):

[math]\mathbf{A} \in \mathbb{R}^{N_e \times 3N_s}[/math]

gdzie czynnik 3 wynika z trzech ortogonalnych składowych momentu dipolowego.

Modele geometrii głowy
  • sferyczny — analityczny, szybki, niedokładny; przydatny do testów i dydaktyki
  • BEM (Boundary Element Method) — trójwarstwowy model skóry/czaszki/mózgu z siatką trójkątów; standard dla EEG
  • FEM (Finite Element Method) — uwzględnia anizotropię białej istoty mózgowej; wymagany dla MEG i zaawansowanych zastosowań klinicznych
Programy obliczające model do przodu

MNE-Python (mne.make_bem_model), FieldTrip (ft_prepare_headmodel), OpenMEEG.

Fitowanie dipoli (Dipole Fitting)

Idea

Najprostsze założenie: aktywacja EEG pochodzi z jednego (lub kilku) ekwiwalentnych dipoli prądowych. Każdy dipol opisany jest sześcioma parametrami: trzema współrzędnymi położenia [math](x,y,z)[/math] i trzema składowymi momentu dipolowego [math](q_x, q_y, q_z)[/math].

Formalizm

Dla pojedynczego dipola o momencie [math]\mathbf{q}(t)[/math] w lokalizacji [math]\mathbf{r}_0[/math]:

[math]\mathbf{x}(t) = \mathbf{a}(\mathbf{r}_0)\,\mathbf{q}(t) + \boldsymbol{\varepsilon}(t)[/math]

gdzie [math]\mathbf{a}(\mathbf{r}_0) \in \mathbb{R}^{N_e \times 3}[/math] to kolumna macierzy pola wiodącego odpowiadająca lokalizacji [math]\mathbf{r}_0[/math].

Estymacja polega na minimalizacji reszty:

[math]\hat{\mathbf{r}}_0 = \arg\min_{\mathbf{r}} \left\|\mathbf{x}(t) - \mathbf{a}(\mathbf{r})\,\hat{\mathbf{q}}(t)\right\|^2[/math]

Minimalizację prowadzi się zwykle algorytmem Neldera–Meada lub Levenberga–Marquardta po przestrzeni lokalizacji [math]\mathbf{r}[/math].

Zastosowania i ograniczenia

Kiedy stosować
  • gdy źródło jest punktowe lub quasi-punktowe (np. wczesne składowe słuchowych potencjałów wywołanych N20, P25)
  • gdy liczba aktywnych źródeł jest z góry znana i mała (1–3 dipole)
Ograniczenia
  • wrażliwość na punkt startowy — algorytm może utknąć w minimum lokalnym
  • złe dopasowanie gdy aktywacja pochodzi z rozległego obszaru kory
  • wynik silnie zależy od dokładności modelu do przodu

Minimum Norm Estimation (MNE)

Idea

Zamiast zakładać małą liczbę dipoli, rozkładamy aktywność na gęstą siatkę źródeł pokrywającą całą powierzchnię kory i szukamy rozwiązania o minimalnej normie — tj. takiego, które odtwarza dane, angażując jak najmniejszą całkowitą energię źródeł.

Formalizm

Problem optymalizacji z regularyzacją Tichonowa:

[math]\hat{\mathbf{s}} = \arg\min_{\mathbf{s}} \left\|\mathbf{x} - \mathbf{A}\mathbf{s}\right\|^2 + \lambda\left\|\mathbf{s}\right\|^2[/math]

gdzie [math]\lambda \gt 0[/math] to parametr regularyzacji kontrolujący kompromis między wiernością danym a gładkością rozwiązania.

Rozwiązanie analityczne:

[math]\hat{\mathbf{s}} = \mathbf{A}^T\left(\mathbf{A}\mathbf{A}^T + \lambda\mathbf{I}\right)^{-1}\mathbf{x}[/math]

Macierz [math]\mathbf{W}_{\mathrm{MNE}} = \mathbf{A}^T(\mathbf{A}\mathbf{A}^T + \lambda\mathbf{I})^{-1}[/math] jest liniowym filtrem przestrzennym — analogiem macierzy [math]\mathbf{W}[/math] z CSP, lecz skonstruowanym na podstawie modelu do przodu, nie z danych.

Warianty

dSPM (dynamic Statistical Parametric Mapping, Dale 2000)

Normalizuje estymowaną aktywność przez odchylenie standardowe szumu:

[math]F_{\mathrm{dSPM}}(\mathbf{r}, t) = \frac{\hat{s}(\mathbf{r},t)}{\sigma_{\mathrm{noise}}(\mathbf{r})}[/math]

Wynik ma interpretację statystyczną (wartość F lub Z); popularne w neurobiologii poznawczej.

sLORETA (standardized LORETA, Pascual-Marqui 2002)

Normalizuje przez odchylenie standardowe samego estymatora MNE (nie tylko szumu):

[math]F_{\mathrm{sLORETA}}(\mathbf{r}, t) = \frac{\hat{s}(\mathbf{r},t)}{\sqrt{\left[\mathbf{W}_{\mathrm{MNE}}\mathbf{A}\right]_{(\mathbf{r},\mathbf{r})}}}[/math]

Wykazuje zerowe systematyczne błędy lokalizacji w warunkach idealnych.

Wybór [math]\lambda[/math]

Zwykle [math]\lambda = \alpha \cdot \mathrm{tr}(\mathbf{A}\mathbf{A}^T) / \mathrm{tr}(\mathbf{I})[/math] z [math]\alpha \in [0.05, 0.3][/math]; w MNE-Python domyślnie [math]\alpha = 0.1[/math]. Można też stosować kryterium L-krzywej lub walidację krzyżową.

Zastosowania i ograniczenia

Kiedy stosować
  • gdy aktywacja jest rozproszona lub nieznana a priori
  • jako ogólne narzędzie wizualizacji ERP na powierzchni kory
  • przy małej liczbie elektrod
Ograniczenia
  • tendencja do rozmywania rozwiązania — źródła głębokie są systematycznie niedoszacowane
  • brak zakładanej rzadkości (sparsity); wariant L1 (MCE) jest rzadszy, ale trudniejszy obliczeniowo
  • wymaga dobrego modelu do przodu

Beamformery

Idea

Beamformer (dosł. „formowanie wiązki") to adaptacyjny filtr przestrzenny, który dla każdej lokalizacji [math]\mathbf{r}[/math] konstruuje wagi [math]\mathbf{w}(\mathbf{r})[/math] maksymalizujące stosunek mocy sygnału z tej lokalizacji do mocy wszystkich pozostałych źródeł. W odróżnieniu od MNE, beamformer korzysta z macierzy kowariancji danych — jest więc metodą adaptacyjną, zbliżoną duchem do CSP.

LCMV (Linearly Constrained Minimum Variance, Van Veen 1997)

Szukamy filtra minimalizującego wariancję wyjścia przy ograniczeniu zachowującym sygnał z lokalizacji [math]\mathbf{r}[/math]:

[math]\mathbf{w}^*(\mathbf{r}) = \arg\min_{\mathbf{w}} \mathbf{w}^T\mathbf{C}\mathbf{w} \quad \text{p.o.} \quad \mathbf{w}^T\mathbf{a}(\mathbf{r}) = 1[/math]

gdzie [math]\mathbf{C} = E[\mathbf{x}\mathbf{x}^T][/math] to macierz kowariancji danych.

Rozwiązanie analityczne (mnożnik Lagrange'a):

[math]\mathbf{w}(\mathbf{r}) = \frac{\mathbf{C}^{-1}\mathbf{a}(\mathbf{r})}{\mathbf{a}(\mathbf{r})^T\mathbf{C}^{-1}\mathbf{a}(\mathbf{r})}[/math]

Moc w lokalizacji [math]\mathbf{r}[/math]:

[math]P(\mathbf{r}) = \mathbf{w}(\mathbf{r})^T\mathbf{C}\,\mathbf{w}(\mathbf{r}) = \frac{1}{\mathbf{a}(\mathbf{r})^T\mathbf{C}^{-1}\mathbf{a}(\mathbf{r})}[/math]

DICS (Dynamic Imaging of Coherent Sources, Gross 2001)

Wariant beamformera dla analizy częstotliwościowej: zamiast macierzy kowariancji czasowej stosuje cross-spectral density matrix [math]\mathbf{C}(f)[/math] obliczoną dla wybranego pasma częstości. Pozwala mapować koherencję między lokalizacjami na powierzchni kory.

Analogia z CSP

Zauważmy analogię strukturalną:

CSP LCMV Beamformer
Macierz kowariancji [math]\mathbf{R}_T, \mathbf{R}_{NT}[/math] (dwa warunki) [math]\mathbf{C}[/math] (jeden warunek) + [math]\mathbf{a}(\mathbf{r})[/math]
Filtr [math]\mathbf{w}[/math]: maks. stosunek wariancji [math]\mathbf{w}(\mathbf{r})[/math]: min. wariancja przy ograniczeniu
Przestrzeń elektrod przestrzeni mózgu (woksel po wokselu)
Model do przodu nie tak — przez [math]\mathbf{a}(\mathbf{r})[/math]

Kluczowa różnica: CSP optymalizuje globalnie jednym zagadnieniem własnym, beamformer rozwiązuje osobny problem dla każdej lokalizacji.

Zastosowania i ograniczenia

Kiedy stosować
  • analiza ERD/ERS (rytmy ruchowe, alfa, beta) w przestrzeni mózgu
  • lokalizacja źródeł oscylacyjnych w MEG
  • DICS do badania koherencji korowo-korowej
Ograniczenia
  • silnie skorelowane źródła są tłumione (source cancellation) — beamformer zakłada, że źródła są nieskorelowane
  • wymaga dużej liczby próbek do stabilnej estymacji [math]\mathbf{C}^{-1}[/math]
  • przy małej liczbie elektrod macierz [math]\mathbf{C}[/math] może być źle uwarunkowana
Podsumowanie
Inverse problem vs bss.png

Podsumowanie: kiedy stosować którą metodę

Metoda Założenie a priori Dane wejściowe Typowe zastosowanie
Dipole fitting 1–3 punktowe źródła ERP, AEP, SEP lokalizacja wczesnych składowych wywołanych
MNE / dSPM / sLORETA minimalna energia; źródła ciągłe ERP, dowolny sygnał wizualizacja rozproszonej aktywności na korze
LCMV Beamformer źródła nieskorelowane ERD/ERS, dane resting-state lokalizacja źródeł oscylacyjnych w MEG/EEG
DICS j.w. + analiza spektralna dane oscylacyjne mapa koherencji korowo-korowej
CSP / ICA brak modelu geometrycznego dwa warunki / jeden segment klasyfikacja BCI, usuwanie artefaktów
Uwaga praktyczna

Wszystkie metody problemu odwrotnego wymagają przygotowania modelu do przodu. Zaleca się korzystanie z MNE-Python lub FieldTrip — oba pakiety oferują spójne pipeline'y od surowych danych do map korowych i są szeroko stosowane w literaturze, co ułatwia replikację wyników.

Sesja 7: Problem odwrotny — część praktyczna

Czas
~30 min wprowadzenie + ~120 min ćwiczenie

Część praktyczna prowadzi przez pełny pipeline lokalizacji źródeł: od wczytania danych w FieldTrip, przez zbudowanie modelu do przodu, do fitowania dipola i interpretacji wyniku. Sesja jest zaprojektowana tak, aby każdy krok bezpośrednio ilustrował jeden element teorii omówionej w części wykładowej.

Dane

Używamy standardowego zestawu danych FieldTrip (dane somatosensoryczne MEG/EEG, bodziec elektryczny na nadgarstku — składowa N20):

Ten sam zestaw danych jest używany we wszystkich tutorialach FieldTrip wymienionych w tej sesji, co ułatwia porównywanie wyników między krokami.


Instalacja FieldTrip
Wymagania wstępne:

FieldTrip wymaga MATLABa. Przed instalacją sprawdź kompatybilność swojej wersji: MATLAB requirements

Pobranie toolboxa
Pobierz najnowszą wersję ze strony:

Pobierz FieldTrip

Rozpakuj archiwum ZIP w wybranym katalogu, np. ~/matlab/fieldtrip (Linux/macOS).

Ważne
nie zmieniaj wewnętrznej struktury katalogów po rozpakowaniu —

FieldTrip dynamicznie zarządza swoimi podkatalogami.

Konfiguracja ścieżki MATLAB
Pełna instrukcja: Installation and setting up the path

Uruchom poniższe polecenia na początku każdej sesji pracy z FieldTrip (lub wstaw je do pliku startup.m — patrz niżej): 

restoredefaultpath                       % wyczyść aktualną ścieżkę
addpath('~/fieldtrip')                  % podaj tu swoją ścieżkę do FieldTrip
ft_defaults                              % skonfiguruj FieldTrip i minimalne podkatalogi
Uwaga — częsty błąd

Nie używaj addpath(genpath(...)) ani przycisku „Add with Subfolders" w MATLAB GUI. FieldTrip samodzielnie dodaje potrzebne podkatalogi przez ft_defaults i ft_hastoolbox — dodanie wszystkiego ręcznie powoduje konflikty z funkcjami SPM, EEGLAB i wbudowanymi funkcjami MATLABa (w szczególności z funkcjami z katalogu fieldtrip/compat przeznaczonymi tylko dla starych wersji MATLABa).

Trwała konfiguracja przez startup.m

Żeby nie powtarzać powyższych poleceń przy każdym uruchomieniu MATLABa, utwórz plik startup.m w swoim katalogu MATLABa (zwykle Dokumenty/MATLAB/startup.m) z zawartością: 

restoredefaultpath
addpath('C:\fieldtrip')   % zmień na swoją ścieżkę
ft_defaults

MATLAB wykonuje ten plik automatycznie przy każdym uruchomieniu.

Uwaga dla użytkowników EEGLAB
Jeśli korzystasz z EEGLAB i FieldTrip w tej samej sesji,

zawsze ładuj FieldTrip przez restoredefaultpath + addpath + ft_defaults, a nie przez interfejs EEGLAB. Oba toolboxy zawierają funkcje o tych samych nazwach (np. eegplot, topoplot) — kolejność na ścieżce decyduje, która wersja zostanie wywołana.

Weryfikacja instalacji
Sprawdź poprawność instalacji jednym poleceniem:
ft_defaults
ft_version

Powinieneś zobaczyć numer wersji i datę kompilacji, np.: 

FieldTrip version 20250901, path: C:\fieldtrip

Jeśli MATLAB zgłasza błąd Undefined function 'ft_defaults', ścieżka nie została ustawiona poprawnie — wróć do kroku z addpath.

Krok 0: Wejście do FieldTrip

Przed właściwym ćwiczeniem proszę zapoznać się z tutorialem wprowadzającym:

Introduction to the FieldTrip toolbox

Na co zwrócić uwagę przechodząc ten tutorial

Studenci znający EEGLAB natkną się na dwie zasadnicze różnice w filozofii pracy:

  • Struktura cfg zamiast GUI. Każda funkcja FieldTrip przyjmuje strukturę konfiguracyjną cfg i zwraca nową strukturę danych. Nie ma okien dialogowych — cała analiza jest skryptem. Jest to bliższe temu, co robiliśmy w MATLAB-owych sesjach CSP/ICA.
  • Niezmienność danych. FieldTrip nie modyfikuje struktur w miejscu — każde wywołanie funkcji zwraca nową strukturę. Pozwala to łatwo porównywać wyniki różnych parametrów, ale wymaga dbałości o nazewnictwo zmiennych.
Pytanie do dyskusji po lekturze

Jak odpowiadają sobie nawzajem: EEG.data w EEGLAB i struktura data w FieldTrip? Gdzie w strukturach FieldTrip przechowywane są informacje o lokalizacji elektrod?

Krok 1: Model do przodu — BEM dla EEG

Przed fitowaniem dipola musimy obliczyć macierz pola wiodącego [math]\mathbf{A}[/math]. Proszę przejść tutorial:

Construct a BEM headmodel for EEG source analysis

Cel tego kroku

Zrozumieć skąd pochodzi kolumna [math]\mathbf{a}(\mathbf{r})[/math] używana we wszystkich metodach Sesji 7. Każda kolumna [math]\mathbf{a}(\mathbf{r})[/math] mówi: „jeśli w miejscu [math]\mathbf{r}[/math] znajdowałby się dipol o jednostkowym momencie, to jakie potencjały zmierzylibyśmy na każdej z elektrod?"

Kluczowe funkcje FieldTrip w tym kroku
% Segmentacja MRI na trzy tkanki (skóra, czaszka, mózg)
cfg           = [];
cfg.output    = {'brain','skull','scalp'};
segmentedmri  = ft_volumesegment(cfg, mri);

% Budowanie siatki trójkątów BEM
cfg             = [];
cfg.tissue      = {'brain','skull','scalp'};
cfg.numvertices = [3000 2000 1000];
bnd             = ft_prepare_mesh(cfg, segmentedmri);

% Obliczenie modelu do przodu metodą BEM (OpenMEEG)
cfg              = [];
cfg.method       = 'openmeeg';
cfg.conductivity = [0.33 0.0041 0.33];   % mózg, czaszka, skóra [S/m]
headmodel        = ft_prepare_headmodel(cfg, bnd);
Na co zwrócić uwagę
  • Stosunek przewodności czaszka/mózg wynosi tu ~1/80.
 Jak myślisz, jak bardzo błąd w tym parametrze wpłynie na lokalizację źródeł?
  • Funkcja ft_prepare_mesh tworzy trzy zagnieżdżone siatki.
 Jak ma się to do opisu BEM z wykładu?

Krok 2: Fitowanie dipola (ćwiczenie główne)

Czas
~90 min

To jest główne ćwiczenie praktyczne Sesji 7. Proszę przejść tutorial:

Dipole fitting of combined MEG/EEG data

Poniżej znajdziesz szkielet kodu z pytaniami pomocniczymi do samodzielnego wypełnienia — zacznij od tutorialu, a potem wróć do szkieletu i uzupełnij brakujące fragmenty.

Szkielet kodu — fitowanie dipola dla składowej N20

%% KROK A: Wczytanie i preprocessing danych EEG
% Dane: Subject01.zip (patrz wyżej)
% Używamy tylko kanałów EEG, odrzucamy MEG.

cfg                    = [];
cfg.dataset            = 'Subject01.ds';
cfg.trialdef.eventtype = 'backpanel trigger';
cfg.trialdef.eventvalue = 3;           % bodziec na lewym nadgarstku
cfg.trialdef.prestim   = 0.1;
cfg.trialdef.poststim  = 0.3;
cfg                    = ft_definetrial(cfg);

cfg.channel            = 'EEG';
cfg.continuous         = 'yes';
cfg.demean             = 'yes';
cfg.baselinewindow     = [-0.1 0];
data                   = ft_preprocessing(cfg);

%% KROK B: Uśrednianie — potencjał wywołany
% TUTAJ: użyj ft_timelockanalysis aby obliczyć ERP (średnią po powtórzeniach)
% Wynik zapisz jako 'timelock'

% TUTAJ

%% KROK C: Wczytanie modelu do przodu i lokalizacji elektrod
% Zakładamy, że headmodel i elec zostały obliczone w Kroku 1.
% Jeśli nie — można użyć gotowego modelu z paczki danych Subject01.

load('standard_bem.mat');          % headmodel
load('Subject01_elec.mat');        % lokalizacje elektrod (po dopasowaniu do MRI)

%% KROK D: Przygotowanie siatki poszukiwań (grid)
% Definiujemy siatkę regularną 1 cm wewnątrz mózgu.
% ft_prepare_sourcemodel oblicza też macierz pola wiodącego A
% (tzw. leadfield) dla każdego punktu siatki.

cfg                 = [];
cfg.headmodel       = headmodel;
cfg.elec            = elec;
cfg.resolution      = 1;           % odstęp siatki w cm
cfg.unit            = 'cm';
sourcemodel         = ft_prepare_sourcemodel(cfg);

%% KROK E: Fitowanie dipola
% Szukamy jednego dipola najlepiej wyjaśniającego topografię N20.
% Okno czasowe 0.016–0.022 s odpowiada składowej N20.

cfg             = [];
cfg.latency     = [0.016 0.022];   % okno N20 w sekundach
cfg.numdipoles  = 1;
cfg.headmodel   = headmodel;
cfg.elec        = elec;
cfg.sourcemodel = sourcemodel;
cfg.nonlinear   = 'yes';           % optymalizacja nieliniowa po przestrzeni

% TUTAJ: source = ft_dipolefitting(cfg, timelock);

%% KROK F: Wizualizacja wyniku
% Nałożenie lokalizacji dipola na MRI anatomiczne.

load('standard_mri.mat');          % MRI w tym samym układzie współrzędnych

cfg               = [];
cfg.location      = source.dip.pos;   % współrzędne znalezionego dipola [x y z]

% TUTAJ: ft_sourceplot(cfg, mri);

%% KROK G: Ocena dopasowania
% Goodness of fit (GoF) mówi jaka część wariancji danych
% jest wyjaśniona przez znaleziony dipol.
% Wartość > 0.85 uznaje się za dobre dopasowanie.

fprintf('Goodness of fit: %.1f%%\n', source.dip.rv * 100);
% Uwaga: w FieldTrip rv to residual variance (1 - GoF).
% Dobry dipol ma rv < 0.15.

disp('Lokalizacja dipola [x y z] w mm:')
disp(source.dip.pos)

disp('Moment dipolowy [qx qy qz]:')
disp(source.dip.mom)

Zadania do samodzielnej realizacji

Zadanie 1 (obowiązkowe)
Pipeline i lokalizacja

Uzupełnij brakujące fragmenty szkieletu (kroki B i E–F). Gdzie na MRI ląduje znaleziony dipol? Czy odpowiada to oczekiwanej anatomii składowej N20 (pierwotna kora somatosensoryczna S1, okolica bruzdy centralnej)?

Zadanie 2 (obowiązkowe)
Wpływ okna czasowego

Zmień cfg.latency kolejno na:

  • [0.016 0.022] — okno N20
  • [0.025 0.035] — okno P25
  • [0.045 0.060] — okno N45

Czy lokalizacja dipola zmienia się między tymi składowymi? Jak zmienia się GoF? Co to mówi o tym, która składowa jest lepiej opisana modelem jednego dipola?

Zadanie 3 (obowiązkowe)
Dwa dipole

Zmień cfg.numdipoles = 2 i powtórz fitowanie dla okna N20. Czy model dwóch dipoli daje istotnie lepszy GoF niż jeden? Jak interpretujesz dwa znalezione źródła?

Zadanie 4 (dla chętnych)
Wrażliwość na model do przodu

FieldTrip pozwala użyć uproszczonego modelu sferycznego zamiast BEM: ustaw cfg.headmodel na model jednorodnej sfery (cfg.method = 'singlesphere' w ft_prepare_headmodel). O ile zmienia się lokalizacja dipola N20? Czy GoF wzrasta czy maleje?

Pytania do dyskusji
  • Fitowanie dipola wymaga podania okna czasowego — dlaczego?
 Co stałoby się gdybyśmy dopasowywali dipol do całej epoki (−100 do +300 ms)?
  • Algorytm minimalizacji może utknąć w minimum lokalnym.
 Jak cfg.sourcemodel (siatka) pomaga temu zaradzić?
 Jaki jest kompromis między rozdzielczością siatki a czasem obliczeń?
  • Porównaj mapę topograficzną ERP w oknie N20
 z mapą rekonstruowaną z dopasowanego dipola (wywołaj ft_dipolefitting
 z opcją cfg.dipfit.display = 'yes').
 Czy widzisz różnicę? Skąd ona pochodzi?

Krok 3: Dalsze metody — dla zainteresowanych

Poniższe tutoriale stanowią naturalne rozwinięcie ćwiczenia głównego i ilustrują pozostałe metody omówione w teorii Sesji 7. Nie są wymagane na zaliczenie, ale gorąco je polecamy jeśli planujesz używać lokalizacji źródeł w pracy badawczej.

Minimum Norm Estimation

Source reconstruction of event-related fields using minimum-norm estimation

Ten sam sygnał N20 co w ćwiczeniu głównym, ale zamiast jednego dipola otrzymujesz ciągłą mapę aktywności na całej powierzchni kory. Zwróć uwagę jak wybór parametru regularyzacji [math]\lambda[/math] wpływa na rozmycie rozwiązania. Porównaj lokalizację maksimum mapy MNE z lokalizacją dipola z Kroku 2.

Beamformer (LCMV/DICS)

Localizing oscillatory sources in MEG data using a beamformer

Tutorial wykorzystuje ten sam zestaw danych Subject01 i lokalizuje źródła oscylacji gamma wywołanych bodźcem — bezpośrednie rozwinięcie analizy ERD/ERS, którą już znasz. Kluczowa różnica w stosunku do MNE: beamformer korzysta z macierzy kowariancji danych [math]\mathbf{C}[/math], a nie tylko z modelu do przodu. Sprawdź jak wyglądają filtry przestrzenne [math]\mathbf{w}(\mathbf{r})[/math] dla różnych lokalizacji i porównaj je koncepcyjnie z filtrami CSP z Sesji 1–3.




-->

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