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

Z Brain-wiki
Linia 69: Linia 69:
 
==== <tt>filtruj_i_tnij.m</tt>====
 
==== <tt>filtruj_i_tnij.m</tt>====
 
Musimy teraz napisać funkcję <tt>filtruj_i_tnij</tt>. Musi ona kolejno wczytać nasze zbiory danych, stworzyć montaż kanałów, przefiltrować je, kanał po kanale, wyciąć interesujące nas epoki, wykonać dla nich korektę poziomu "0", i zapisać tak przetworzone dane do nowych "*.set"-ów. Szkilet tej funkcji wygląda następująco, proszę uzupełnić brakujące fragmenty:
 
Musimy teraz napisać funkcję <tt>filtruj_i_tnij</tt>. Musi ona kolejno wczytać nasze zbiory danych, stworzyć montaż kanałów, przefiltrować je, kanał po kanale, wyciąć interesujące nas epoki, wykonać dla nich korektę poziomu "0", i zapisać tak przetworzone dane do nowych "*.set"-ów. Szkilet tej funkcji wygląda następująco, proszę uzupełnić brakujące fragmenty:
 +
 
<source lang = matlab>
 
<source lang = matlab>
   
+
 
 +
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
  
 
</source>
 
</source>
  
 +
==== Wracamy na chwilę do <tt>run_me.m</tt> ====
  
 
<source lang = matlab>
 
<source lang = matlab>
Linia 84: Linia 154:
 
%Tools -> Reject data epochs -> Reject data (all methods)
 
%Tools -> Reject data epochs -> Reject data (all methods)
 
%link do tutoriala EEGlab: https://sccn.ucsd.edu/wiki/Chapter_01:_Rejecting_Artifacts#Rejecting_artifacts_in_epoched_data
 
%link do tutoriala EEGlab: https://sccn.ucsd.edu/wiki/Chapter_01:_Rejecting_Artifacts#Rejecting_artifacts_in_epoched_data
 
+
</source>
 +
====Automatyczne usowanie artefaktow====
 +
<source lang = matlab>
 
%% Automatyczne usowanie artefaktow
 
%% Automatyczne usowanie artefaktow
  

Wersja z 11:50, 14 mar 2018

Laboratorium_EEG/EEGLAB

Zajęcia bazują na podręczniku praktycznym użytkowania pakietu EEGLAB dostępnym na stronach:

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, przefiltrować je, kanał po kanale, wyciąć interesujące nas epoki, wykonać dla nich korektę poziomu "0", i zapisać tak przetworzone dane do nowych "*.set"-ów. Szkilet 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

Wracamy na chwilę do run_me.m

%% DOBOR PROGOW DLA ARTEFAKTOW na podstawie ogladanie sygnalu i prob z roznymi progami do autoamtu (abnormal values / abnorlmal trends)

EEG = pop_loadset('filename','B01_epoch.set','filepath',directoryOut);
[ALLEEG, EEG, CURRENTSET] = eeg_store(ALLEEG, EEG, 0);
EEG = eeg_checkset( EEG);
eeglab redraw
%Tools -> Reject data epochs -> Reject data (all methods)
%link do tutoriala EEGlab: https://sccn.ucsd.edu/wiki/Chapter_01:_Rejecting_Artifacts#Rejecting_artifacts_in_epoched_data

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