Laboratorium EEG/CSP

Laboratorium_EEG/BSS

Ślepa separacja źródeł

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

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

Filtry przestrzenne i ślepa separacja źródeł

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

Niech:

[math]s(t)[/math] - aktywność niezależnych źródeł,

[math]x(t)[/math] mierzony sygnał

[math]A[/math] macierz przejścia taka, że:

[math]x(t) = A s(t)[/math] (*)

[math]s(t) = A^{-1}x(t) = P x(t)[/math]

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

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

Podstawiając (*) mamy:

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

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

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

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

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

Common Spatial Pattern

Koncepcja

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

Formalizm

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

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

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

Wariancja tego rzutowanego sygnału to:

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

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

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

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

[math] \nabla J(w) = \frac{ 1 C_{T} w+w^T C_{T} 1}{w^T C_{NT} w}-\frac{w^T C_{T} w\left( 1 C_{NT} w+w^T C_{NT} 1\right)}{\left(w^T C_{NT} w\right)^2}[/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:

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

Liczba [math] \lambda = \frac{w^T C_{T} w}{w^T C_{NT} w} C_{NT} w [/math] jest uogólnioną wartością własną, zaś [math]w[/math] jest uogólnionym wektorem własnym odpowiadającym tej wartości.

Aby znaleźć [math] \lambda[/math] i [math]w[/math] wystarczy rozwiązać zagadnienie własne. W matlabie możemy w tym celu wykorzystać funkcję eig

Ćwiczenie symulacyjne

% symulowany eksperyment składa się z sinusoidy udającej alfę spoczynkową i
% gausa udającego potencjał wywołany
% źródła te są symulowane niezależnie a potem mieszane przez macierz L
% symulujemy źródła
% s1 - symuluje alfę
% s2 - symuluje "potencjał wywołany" (ERP)

%ustawiamy parametry do symulacji sygnałów
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));

% filtr przestrzenny - z takimi wagami trzeba wziąść kanały EEG aby odzyskać sygnały źródłowe
P = [1 2
    1.5 1.3]; 


% topografie - z takimi wagami źródła dokładają się do poszczególnych elektrod
A = P^(-1); 


for r =1:N_rep % tworzymy kolejne realizacje "eksperymentu"
    s1 = sin(2*pi*11*t +pi/2+ 0*2*pi*rand())+ 0.02*randn(size(t));  % źródło alfa
    s2 = exp(-((t-0.8)/0.05).^2)+ 0.01*randn(size(t));              % źródło ERP
       
    s(r,1,:) = s1;
    s(r,2,:) = s2;
    tmp = squeeze(s(r,:,:));
    n = 0*randn(size(tmp));
    X(r,:,:) = A*tmp +n; % rzutujemy sygnały źródłowe na elektrody s -> x
    
end

% wycinamy warunki 
% baseline_ind   to indeksy pierwszej połowy każdego powtórzenia "baseline"
% ERP_ind        to indeksy drugiej połowy każdego powtórzenia zawierająca "ERP"
baseline_ind = find(t<0.5);
ERP_ind = find(t>=0.5);

x_baseline_kan_1 = X(:,1,baseline_ind);
x_baseline_kan_2 = X(:,2,baseline_ind);

x_ERP_kan_1 = X(:,1,ERP_ind);
x_ERP_kan_2 = X(:,2,ERP_ind);



% liczymy średnie macierze kowariancji:
R_E = zeros(N_chan,N_chan);
R_B = zeros(N_chan,N_chan);
for r =1:N_rep
    B = squeeze(X(r,:,baseline_ind)); 
    R_B = R_B + B*B' ;
    
    E = squeeze(X(r,:,ERP_ind));
    R_E = R_E + E*E' ;
end

R_B = R_B/N_rep;
R_E = R_E/N_rep;

% rozwiązujemy uogólnione zagadnienie własne
[W,Lambda]=eig(R_E,R_B); % możliwa jest też otymalizacja wzg. średniej macierzy kowariancji (R_B+R_A)/2);

% odzyskujemy sygnały źródeł
for r =1:N_rep
    S(r,:,:) = W'*squeeze(X(r,:,:));
end

% pobieramy wycinki odpowiadające obu częściom eksperymentu z estymowanych
% źródeł
s_baseline_estymowany_kan1 = squeeze(  S(:,1,baseline_ind));
s_baseline_estymowany_kan2 = squeeze(  S(:,2,baseline_ind));

s_ERP_estymowany_kan1 = squeeze(S(:,1,ERP_ind));
s_ERP_estymowany_kan2 = squeeze(S(:,2,ERP_ind));

%%%%%%%%%%%%%% Iloustracje %%%%%%%%%%%%%%%%%%%%%%%
% ilustracja sygnałów mierzonych
figure(1);clf
subplot(2,2,1);
    plot(t(baseline_ind),(squeeze(X(:,1,baseline_ind)))','b'); hold on
    plot(t(ERP_ind),(squeeze(  X(:,1,ERP_ind)))','r'); hold off
    xlabel('elektroda 1')
    title('ilustracja sutuacji pomiarowej -\newline znane są potencjały na elektrodach w dwóch warunkach eksperymentalnych')
subplot(2,2,3);
    plot(t(baseline_ind),(squeeze(X(:,1,baseline_ind)))','b'); hold on
    plot(t(ERP_ind),(squeeze(  X(:,2,ERP_ind)))','r'); hold off
    xlabel('elektroda 2')
subplot(1,2,2)
    plot(x_baseline_kan_1(:),x_baseline_kan_2(:),'b.');
    hold on
    plot(x_ERP_kan_1(:),x_ERP_kan_2(:),'r.');
    xlim([-2,2])
    ylim([-2,2])
    axis equal
    
    % wektor własny odpowiadajacy największej wartości własnej jest
    % kierunkiem najbardziej różnicującym warunki eksperymentalne
    disp('wartości własne znajdują się na diagonali macierzy Lambda')
    disp(Lambda)
    % rysujemy wersory jednostkowe w kierunkach wektorów własnych
    w1 = W(:,1); 
    w1 = w1/norm(w1);
    w2 = W(:,2); 
    w2 = w2/norm(w2);
    line([0, w1(1) ],[0,w1(2)],'Color',[0,0.3,0])
    text(w1(1),w1(2),'wektor własny1')
    line([0, w2(1) ],[0,w2(2)],'Color',[1,0,1])
    text(w2(1),w2(2),'wektor własny2')
    
    
    xlabel('Amplituda na elektrodzie 1')
    ylabel('Amplituda na elektrodzie 2')   
    legend('baseline','ERP')

% Ilustracja estymowanych źródeł
figure(2);clf
subplot(2,2,1);
    plot(t(baseline_ind),(squeeze(S(:,1,baseline_ind)))','b');hold on
    plot(t(ERP_ind),(squeeze(S(:,1,ERP_ind)))','r');hold off
    xlabel('estymowane zrodlo  1')
    title('ilustracja estymacji -\newline estymowane są potencjały żródeł w dwóch warunkach eksperymentalnych')    
subplot(2,2,3);
    plot(t(baseline_ind),(squeeze(S(:,2,baseline_ind)))','b');hold on
    plot(t(ERP_ind),(squeeze(S(:,2,ERP_ind)))','r');hold off
    xlabel('estymowane zrodlo  2');
subplot(1,2,2)
    plot(s_baseline_estymowany_kan1(:),s_baseline_estymowany_kan2(:),'b.');
    hold on
    plot(s_ERP_estymowany_kan1(:),s_ERP_estymowany_kan2(:),'r.');   
    
    xlabel('Amplituda estym. źródła 1')
    ylabel('Amplituda estym. źródła 2')   
    legend('baseline','ERP')

Zastosowanie filtra CSP do detekcji potencjału P300

Eksperyment

Przygotowanie do badania:

• założyć czepek z elektrodami w systemie 10-20
• elektrody referencyjne: M1 i M2
• elektroda GND w pozycji AFz

Przygotowanie scenariuszy obci

• w terminalu uruchomić obci srv
• w terminalu uruchomić obci_gui --preset brain2013
• w interfejsie GUI zapisujemy scenariusze do własnego katalogu np "P300"
  • "P-Brain2013 Signal (with ID)" jako np. "Sygnal"
  • "P-Brain2013 Calibration p300 " jako "kalibracjaP300"
  • "P-Brain 2013 p300" jako Labirynt
• w przeglądarce plików otwórz katalog P300.Powinien on zawierać pliki: Sygnal.ini, kalibracjaP300.ini, Labirynt.ini oraz katalogi Sygnal_configs, kalibracjaP300_configs, Labirynt_configs.
• edytujemy parametry peerów.
  • z katalogu "P300/Sygnal_configs" kopiujemy plik "amplifier.ini" do katalogu "P300" jako "global_amplifier.ini" . To pozwoli nam zmieniać ustawienia wzmacniacza dla wszystkich scenariuszy jednocześnie.
  • w naszych scenariuszach zapisanych w plikach Sygnal.ini, kalibracjaP300.ini, Labirynt.ini podmieniamy ścieżkę config = w sekcji [peers.amplifier] tak aby

pokazywała na ten skopiowany plik "global_amplifier.ini"

• • w tym pliku "global_amplifier.ini" podmieniamy linijki (tak aby były to listy faktycznie wykorzystywanych kanałów) z:
    • active_channels
    • channel_names
  • dodajemy też linijkę: sampling_rate = 256
• wchodzimy po kolei do katalogów: Sygnal_configs, kalibracjaP300_configs, Labirynt_configs i odnajdujemy plik switch_backup.ini. W tym pliku ustawiamy parametr new_scenario na pusty. To spowoduje, że scenariusze te nie będą uruchamiać kolejnych scenariuszy po zakończeniu działania.

Przeprowadzenie badania:

• obejrzeć sygnal w svarogu. (uruchamiamy go z terminala poleceniem svarog). poprawić ewentualnie źle kontaktujące elektrody.
• w interfejsie GUI zakończyc scenariusz " "
• w interfejsie GUI wybrać scenariusz " " Scenariusz ten przeprowadzi kalibrację interfejsu BCI- dane z tej kalibracji będziemy potem wykorzystywać do dalszych analiz.
• po przeprowadzeniu kalibracji uruchomiony zostanie automatycznie scenariusz przejścia przez labirynt
• danych z kalibracji potrzebować będziemy kilka zestawów. W zależności od chęci prosze powtórzyć 3-4 krotnie scenariusz Kalibracja+labirynt, lub samą kalibrację.

Analiza wstępna

• 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ół trigerów. Dla kazdej realizacji odjąć trend liniowy.
• Sygnał zmontować wzg. "średnich uszu" i wyświetlić średnie przebiegi dla warunku T i NT w układzie topograficznym - wykorzystać w tym celu funkcję topoplot z pakietu eeglab

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

Plik:KalibracjaP300.tar.gz

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

obci/analysis/matlab_obci_signal_processing

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

% ustalamy nzawy plików z danymi
nazwaPliku = 'p_6301423_calibration_p300.obci';
nameOfXMLFile = strcat(nazwaPliku,'.xml');
nameOfTagFile = strcat(nazwaPliku,'.tag'); %tagi = znaczniki zdarzeń
namesOfDataFiles = strcat(nazwaPliku,'.raw');

% inicjujemy obiekt rm
rm = ReadManager(nameOfXMLFile,namesOfDataFiles,nameOfTagFile);

% obieramy przydatne parametry i znaczniki
numberOfChannels  = rm.get_param('number_of_channels');
namesOfChannels   = rm.get_param('channels_names');
samplingFrequency = rm.get_param('sampling_frequency');
tagsStruct        = rm.get_tags();

% tworzenie list znaczników Target i NonTarget
numberOfStruct = length(tagsStruct);
targetTimeStamps = [];
NonTargetTimeStamps = [];
for structNumber = 1:numberOfStruct % iterujemy się przez tagi
    if(strcmp(tagsStruct(structNumber).name,'blink')) % szukamy tagów o nazwie 'blink'
        index = tagsStruct(structNumber).children.index; % tu jest numer pola stymulacji
        target= tagsStruct(structNumber).children.target;% tu jest numer pola na którym wyświetlany jest target
        if index == target % warunek na to, że mamy do czynienia z tagiem target
            targetTimeStamps = [targetTimeStamps tagsStruct(structNumber).start_timestamp]; %dodajemy timeStamp do listy targetów
        else
            NonTargetTimeStamps = [NonTargetTimeStamps tagsStruct(structNumber).start_timestamp];%dodajemy timeStamp do listy non-targetów
        end
        
    end
end



% pobieramy próbki
samples = double(rm.get_samples()); % konwersja na double jest potrzebna żeby dobrze funkcjonowało filtrowanie
samples=samples(1:8,:); % odrzucamy kanały, które nie mają EEG
numberOfChannels =8;

% filtrujemy dolnoprzepustowo aby odrzucić artefakty sieci i część
% artefaktów mięśniowych
[b,a] = cheby2(6,80,25 /(samplingFrequency/2),'low');
for ch = 1:numberOfChannels
    samples(ch,:)=filtfilt(b,a,samples(ch,:));
end

% montujemy dane do wspólnej średniej (common average)
M = -ones(8,8)/8;
M=M+eye(8,8)*9/8;
samples = 0.0715*M*samples;

% wycinamy dane wokół znaczników
PRE = -0.2; % czas przed tagiem w sek.
POST = 0.8; % czas po tagu w sek.
wycinek = floor(PRE*samplingFrequency:POST*samplingFrequency); % tablica ze "standardowymi' indeksami do cięcia

% pobieramy targety
TargetSignal = zeros(length(targetTimeStamps),numberOfChannels, length(wycinek)); % tablica na sygnały target
for trialNumber = 1:length(targetTimeStamps)
    trigerOnset = floor(targetTimeStamps(trialNumber)*samplingFrequency);
    tenWycinek = wycinek + trigerOnset;
    if tenWycinek(1)>0 && tenWycinek(end)<=size(samples,2) % test czy wycinek który chcemy pobrać nie wystaje poza dostępny sygnał
        tmpSignal = samples(:,tenWycinek);
        tmpSignal = detrend(tmpSignal')'; % usuwanie liniowego trendu - przy krótkich wycinkach działa lepiej niż filtrowanie górnoprzepustowe
        TargetSignal(trialNumber, :,:) = tmpSignal;
    end
end

% pobieramy non-targety
NonTargetSignal = zeros(length(NonTargetTimeStamps),numberOfChannels, length(wycinek));% tablica na sygnały non-target
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
%
% dla ilustracji podglądamy średnie po powtórzeniach ze wszystkich targetó∑
% i non-targetów
plot(squeeze(mean(TargetSignal,1))','r');
hold on
plot(squeeze(mean(NonTargetSignal,1))','b')

Analiza CSP

• Korzystając z danych kalibracyjnych wykonac 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.
• Do wykonania tych mapek wykorzystać funkcję topoplot z pakietu eeglab
• Zbadać powtarzalność topografii pomiędzy plkiami konfiguracyjnymi.

Filtry przestrzenne dla SSEP

Teoria

Proszę zapoznać się z koncepcją filtra przestrzennego dla SSVEP zaprezentowaną tu: http://www.eurasip.org/Proceedings/Eusipco/Eusipco2009/contents/papers/1569193209.pdf

Eksperyment ASSR

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

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

Rejestracja sygnału

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

Analiza

JZ: zmieniłbym analizę na czas-częstość i zrobił porównanie montażu usznego do filtra G.G. Moliny

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

1. Z sygnału wycinamy fragmenty od 0 do 5 sek. dla wszystkich elektrod położone nad korą słuchową.
2. Dla każdej realizacji obliczamy widma metodą Welcha.
3. Otrzymane zespolone widma uśredniamy po realizacjach.
4. Sprawdzamy czy w uśrednionym widmie występuję maksimum w częstości modulacji tj. 40 Hz.

Wyznaczenie przebiegu czasowego ERD i ERS

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

Transformacja Hjortha

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