Laboratorium EEG/EEGLAB: Różnice pomiędzy wersjami
Linia 173: | Linia 173: | ||
Jak już zdecydujemy się na konkretne wartości to ... | Jak już zdecydujemy się na konkretne wartości to ... | ||
− | + | === Wracamy na chwilę do <tt>run_me.m</tt> === | |
====Automatyczne usowanie artefaktow==== | ====Automatyczne usowanie artefaktow==== |
Wersja z 12:06, 14 mar 2018
Laboratorium_EEG/EEGLAB
Zajęcia bazują na podręczniku praktycznym użytkowania pakietu EEGLAB dostępnym na stronach:
- rozdziały 1-8 z części I: http://sccn.ucsd.edu/wiki/EEGLAB#The_EEGLAB_Tutorial_Outline
- rozdziały z części: http://sccn.ucsd.edu/wiki/Chapter_02:_Writing_EEGLAB_Scripts
Spis treści
Analiza potencjałów wywołanych
Wprowadzenie
W tej części zajęć naszym celem jest przećwiczenie pracy z danymi z eksperymentu w którym badamy potencjały wywołane na grupie osób. Dane pochodzą z eksperymentu opisanego w artykule: Effects of Valence and Origin of Emotions in Word Processing Evidenced by Event Related Potential Correlates in a Lexical Decision Task
Skupimy się na badaniu różnic pomiędzy przetwarzaniem słów i pseudo-słów.
Przebieg eksperymentu
Uczestnicy zostali poinformowani o celu eksperymentu i charakterze pomiaru EEG. Zachęciliśmy ich, aby utrzymywali wygodną postawę i kontrolowali mruganie. Protokół zapewniał 3-s przerwy na normalne mruganie co 10 prób, jak również dwie dłuższe przerwy, których czas trwania jest kontrolowany przez uczestnika, dla odpoczynku i poprawy postawy. Długie przerwy występowały co 270 prób.
Zadaniem było odczytanie bodźców, które pojawiły się na środku ekranu, i zaklasyfikowanie ich jako słów lub pseudo-słów przez naciśnięcie oznaczonych klawiszy na klawiaturze. Rejestrowano typ i opóźnienie odpowiedzi. Pojedynczy blok eksperymentalny zawierał 135 słów i 135 pseudo-słów; ten blok powtórzono trzy razy. Słowa i pseudo słowa były wyświetlane w losowej kolejności we wszystkich blokach. Próby przebiegały w następujący sposób: (1) punkt fiksacji wyświetlany przez 500 ms; (2) bodziec wyświetlany, dopóki uczestnik nie odpowie; (3) pusty ekran wyświetlany przez losowo zmienny odstęp czasu między 1000 a 1100 ms.
Przygotujmy katalogi do analizy. Przydadzą się następujące:
- DANE
- SKRYPTY
- RYSUNKI
Dane dostępne są tutaj.
Proszę je ściągnąć i wypakować do katalogu DANE. Powinno pojawić się 66 plików, po 2 na każdą osobę. Pliki *.set to pliki matlabowe zawierające struktury EEG, pliki *fdt zawierają właściwe sygnały.
SKRYPTY
Etapy analizy, które musimy wykonać to:
- FILTROWANIE i CIECIE sygnałów
- DOBÓR PROGÓW DLA ARTEFAKTOW na podstawie oglądania sygnału i prób z różnymi progami do autoamtu (abnormal values / abnorlmal trends)
- Automatyczne usuwanie artefaktów
- Przygotowanie danych do analiz wg. przynależności do czynnika grupującego
- Wybór obszarów zainteresowania (grup elektrod do analizy)
- Analiza wariancji
- Zilustrowanie wyników analizy
Aby nasze analizy były łatwo powtarzalne przygotujemy główny skrypt np: run_me.m, do którego będziemy dodawać stopniowo skrypty i funkcje kolejnych etapów analizy. Dla czytelności kodu i możliwości uruchamiania osobno poszczególnych etapów posłużymy się technologią komórek w Matlabie (komórkę rozpoczynają dwa procenty '%%');
Zacznijmy zatem pisać run_me.m. Ustawiamy zmienne ze ścieżkami, wybieramy, które dane wchodzą do analizy (część została odrzucona ze względu na nadmiar artefaktów), częstości odcięcia filtrów dolno- i górno-przepustowego, czasy wokół momentu wyświetlenia bodźca, informacja, które kanały zawierają EEG, zakres czasu do korekty poziomu "0":
%% USTAWIENIA:
directoryIn = '../DANE/'; % lokalizacja folderu z surowymi danymi
directoryOut = '../DANE/'; % lokalizacja do zapisu przygotowanych danych
addName = '_epoch'; % dodatek do nazwy pliku przy zapisywaniu przygotowanych danych
subjects = [1:6,8:11 ,15:24,26:37]; % lista badanych do wczytania (z poprawnymi danymi)
LowFreqCut = 0.1; % granica odciecia filtru gornoprzepustowego
HighFreqCut = 30; % granica odciecia filtru dolnoprzepustowego
startT = -0.3; % start do wyciecia epok [s]
stopT = 1.0; % stop do wyciecia epok [s]
eegChan = 1:19; % numery kanalow EEG
baseline = [-200, 0]; % baseline [od,do] [ms] do odjecia
%% FILTROWANIE i CIECIE
eeglab
filtruj_i_tnij(EEG , ALLEEG , CURRENTSET, directoryIn, directoryOut, addName, subjects, LowFreqCut, HighFreqCut, startT, stopT, eegChan, baseline)
filtruj_i_tnij.m
Musimy teraz napisać funkcję filtruj_i_tnij. Musi ona kolejno:
- wczytać nasze zbiory danych,
- stworzyć montaż kanałów do średniej ze wszystkich kanałów, tzw. common average reference (CAR),
- przefiltrować je,
- kanał po kanale,
- wyciąć interesujące nas epoki,
- wykonać dla nich korektę poziomu "0",
- zapisać tak przetworzone dane do nowych "*.set"-ów.
Szkielet tej funkcji wygląda następująco, proszę uzupełnić brakujące fragmenty:
function filtruj_i_tnij(EEG, ALLEEG, CURRENTSET, directoryIn, directoryOut, addName, subjects, LowFreqCut, HighFreqCut, startT, stopT, eegChan, baseline)
% zapamiętujemy stan wejściowy zmiennych globalnych eeglaba
EEG0 = EEG;
ALLEEG0 = ALLEEG;
CURRENTSET0 = CURRENTSET;
events = {'Pseudo','ANeg','ANeu','APos','ONeg','ONeu','OPos','RNeg','RNeu','RPos'}; % typ eventow do wyciecia
for sub = subjects
% tworzymy nazwę setu dla osoby z numerem 'sub'
if sub < 10
nameOfSubj = ['B0' num2str(sub)];
else
nameOfSubj = ['B' num2str(sub)];
end
% wczytujemy set z danymi
EEG = pop_loadset('filename',[nameOfSubj '.set'],'filepath',directoryIn);
% uaktualniamy zestaw danych globalnych dla eeglaba:
[ALLEEG, EEG, CURRENTSET] = eeg_store( ALLEEG, EEG, 0 );
EEG = eeg_checkset( EEG);
% zostawiamy tylko potrzebne kanaly EEG , te które są wymienione w zmiennej eegChan:
EEG.data = EEG.data(....); % wycianmy dane
EEG.nbchan = .... ; % uaktualniamy liczbę kanałów
EEG.chanlocs = EEG.chanlocs(1,...); % obcianamy lstę położeń kanałów
% referencja - odjecie sredniej ze wszystkich kanalow
EEG = pop_reref(.... );
EEG = eeg_checkset(EEG);
% filtracja danych - zastosujemy filtr Butterwartha 2-go rzędu
[b1,a1] = butter(2, .... , 'high');
[b2,a2] = butter(2, ... , 'low');
[b3,a3] = butter(2, 2*[49.5 50.5]/EEG.srate, 'stop'); % filtr pasmowozaporowy "notch" na sieć
for channel = 1:size(EEG.data,1)
EEG.data(channel,:) = ... % górnoprzepustowy
EEG.data(channel,:) = ... % dolno
EEG.data(channel,:) = ... % notch
end
EEG = eeg_checkset( EEG );
% ekstrakcja epok:
% zauważmy, że podajemy w tej funkcji listę eventów wokół których mają być robione wycinki i okresy wycinków
% przy okazji nadajemy setowi nazwę z dodanym stringiem addName:
EEG = pop_epoch( EEG, events, [startT, stopT], 'newname', [nameOfSubj addName], 'epochinfo', 'yes');
[ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, 1,'gui','off');
EEG = eeg_checkset( EEG );
% usuniecie baseline
EEG = pop_rmbase(...);
EEG = eeg_checkset( EEG );
% zapisanie setu
EEG = pop_saveset( EEG, 'filename',[nameOfSubj addName],'filepath',directoryOut);
disp(['filtrowanie i cięcie ----------- dane ' nameOfSubj ' : GOTOWE ----------------'])
end
% przywracamy stan wejściowy zmiennych globalnych eeglaba
EEG0 = EEG;
EEG = EEG0;
ALLEEG = ALLEEG0;
CURRENTSET = CURRENTSET0;
end
DOBÓR PROGOW DLA ARTEFAKTÓW
Usuwanie artefaktów stanowi jeden z bardziej żmudnych ale niezbędnych etapów przygotowania danych do dalszej analizy. Jeśli mamy dostatecznie dużo danych to najbardziej bezpiecznym podejściem jest usunięcie tych epok (spośród wcześniej wyciętych), które mają artefakty. Do najczęstszych i najbardziej szkodliwych w analizie potencjałów wywołanych należą:
- mrugnięcia
- płynięcie sygnału
- mięśnie
Zwykle usuwanie artefaktów wykonujemy dwuetapowo:
- Dla ustalonych progów na wykrycie artefaktu oznaczamy automatycznie złe fragmenty
- Dokonujemy wizualnego przeglądu tak oznaczonych danych i wykonujemy ewentualne korekty, w miejscach nietypowych zaburzeń.
Tu dla oszczędzenia czasu zrobimy tylko procedurę uproszczoną:
- wczytamy jeden set, np.: B01_epoch.set
- popracujemy chwilę z narzędziem: Tools -> Reject data epochs -> Reject data (all methods)
- link do tutoriala EEGlab na ten temat: https://sccn.ucsd.edu/wiki/Chapter_01:_Rejecting_Artifacts#Rejecting_artifacts_in_epoched_data
Naszym celem jest takie dobranie progów
- abnormalValue: wartosc [uV] powyżej/poniżej(na minusie) której sygnał uznawany jest za artefakt
- abnormalTrend: maksymalne nachylenie trendu [uV/epokę], R squared limit (0-1)]
Jak już zdecydujemy się na konkretne wartości to ...
Wracamy na chwilę do run_me.m
Automatyczne usowanie artefaktow
%% Automatyczne usowanie artefaktow
abnormalValue = 60; % wartosc [uV] powyzej/ponizej(na minusie) ktorej oznaczane sa artefakty
abnormalTrend = [50, 0.3]; % [max slope [uV/epoch], R squared limit (0-1)]
markArtefacts(EEG , ALLEEG , CURRENTSET, directoryOut, addName, subjects, abnormalValue, abnormalTrend)
%% przygotowanie danych
data = prepareData(EEG , ALLEEG , CURRENTSET, directoryOut, addName, subjects);
% data - macierz z danymi o rozmiarze [2(1.Words|2.Pseudo) x subjects x kanaly x probki]
%% GFP
timeMarks = [];
plotGFP(data,startT,stopT,EEG.srate,timeMarks,EEG.chanlocs)
timeMarks = [70,120,225,290,375,670];
plotGFP(data,startT,stopT,EEG.srate,timeMarks,EEG.chanlocs)
%% ANOVA
ROI.channels = [[18,13];[15,10];[19,17];[9,4];[11,6]];
ROI.labels = {'LF','CF','RF','LP','RP'};
dataANOVA = prepareDataForANOVA(data,timeMarks,ROI,startT,stopT,EEG.srate);
% dataANOVA - [5 okien czasowych x obserwacje x |1.srednia amplituda 2. subject 3.slowo(0)/pseudo(1) 4.ROI|]
for window = 1:(length(timeMarks)-1)
dANOVA = squeeze(dataANOVA(window,:,:));
stats = anova(dANOVA(:,1),dANOVA(:,2),dANOVA(:,3),dANOVA(:,4),{'type','ROI'});
stats2 = PostHocStats(stats,dANOVA,ROI);
disp(['------------- Time window: ' num2str(timeMarks(window)) '-' num2str(timeMarks(window+1)) '-------------'])
disp(stats)
disp(stats2)
end
%% rysowanie ERP w ROI
plotERP_ROI(data,startT,stopT,EEG.srate,timeMarks,ROI)
Laboratorium_EEG/EEGLAB