Kodu okruchy - mmagnuski/AnalizaDanychEEG01 GitHub Wiki
Co to są funkcje?
Funkcje to procedury, które przyjmują pewne argumenty wejściowe (input), coś z nimi robią i coś nam zwracają (output). Funkcje w matlabie używamy pisząc w command window coś w stylu:
output = funkcja(input);
Przykładem funkcji jest średnia (ang. mean). Funkcji mean
używamy na liście wartości liczbowych (wektorze) aby otrzymać średnią liczb. Lista liczb jest więc argumentem wejściowym, średnia jest natomiast tym, co funkcja mean
zwraca. Spróbujcie:
wartosci = [1, 3, 8];
srednia = mean(wartosci);
Bardzo przydatną funkcją, z której będziemy dziś korzystać jest funkcja plot
. Pozwala ona wyświetlić listę wartości (wektor) jako linię na wykresie. Najpierw stwórzmy wektor (listę wartości liczbowych):
wektor = [3, 8, 1, 2, 5, 6, 3, 9];
Następnie skorzystajmy z funkcji plot
aby wyświetlić wartości wektora jako linię:
plot(wektor);
Funkcja plot
obsługuje wiele różnych opcji (chociażby kolor, czy grubość linii) - pełną dokumentację (tej i dowolnej innej funkcji) możemy zobaczyć pisząc:
help plot
co wypisuje nam tekst dokumentacji w command window. Możemy też napisać:
doc plot
co otworzy dokumnetację w odzielnym oknie jako html.
W następnych krokach użyjemy funkcji mean
oraz plot
aby policzyć i narysować potencjał wywołany dla wybranej elektrody.
Wywoływanie ERPów
Interfejs graficzny EEGlab'a jest świetny do preprocessingu i szybkiego eksplorowania danych, nie jest jednak dobry do analizy. Aby przeprowadzić porządną analizę potrzebujemy zajrzeć pod maskę interfejsu graficznego i skorzystać z komend tekstowych.
Event Related Potential (ERP) to dumna nazwa dla dosyć prostej analizy - uśrednienia powtórzeń tak aby otrzymać dla każdej elektrody zmiany w czasie średniej wartości napięcia. Nauczymy się liczyć ERPy na dwa sposoby: trudniejszy (za pomocą funkcji mean
) oraz prostszy (za pomocą funkcji erpy_daj
, którą napisałem dla was :rat: ).
mean
ERPy funkcją Wasze dane eeg dla aktywnego datasetu są przechowywane przez EEGlaba w zmiennej EEG
. EEG
przechowuje dużo różnych informacji takich jak próbkowanie sygnału (EEG.srate
) czy liczba kanałów (EEG.nbchan
). Sam sygnał jest przechowywany w EEG.data
jako trójwymiarowa macierz (elektrody na czas na epoki), możemy więc policzyć średnią po trzecim wymiarze macierzy (tzn. uśrednić trzeci wymiar) i uzyskać tym samym dla każdej elektrody ERPa (macierz elektrody na czas):
erpy = mean(EEG.data, 3);
Ok, fajnie, brzmi profesjonalnie (i pewnie niezbyt zrozumiale), ale co jeżeli chcemy wyświetlić teraz sygnał dla konkretnej elektrody? Wtedy trzeba zaadresować macierz erpy
, którą otrzymaliśmy w poprzednim kroku. Musimy wtedy wiedzieć, która jest w kolejności wszystkich elektrod ta interesująca nas elektroda. W systemie EGI jest to dosyć proste (pod warunkiem, że mamy w sygnale wszystkie elektrody), elektroda E5
jest piąta a E42
czterdziesta druga. Dla elektrody E5
wyświetlamy ERPa tak:
plot(erpy(5,:));
Gdybyśmy jednak chcieli policzyć ERPa dla konkretnego warunku, na konkretnej elektrodzie - sprawa robi się jeszcze trudniejsza... Tu przychodzi nam z pomocą erpy_daj
!
erpy_daj
Po prostu erpy_daj
ułatwia życie studenta neurokognitywistyki (innych kierunków też!) pozwalając policzyć erpa wygodnie i rach ciach. Np. aby otrzymać ERPa dla elektrody E23
, dla warunku face_0
piszemy (pamiętaj, że erpy_daj
musi być dodane do ścieżki matlaba!):
erp = erpy_daj(EEG, 'E23', 'face_0');
WOW, jakie to proste i fajne! Od teraz ERPy są moją pasją! :hamster: :hammer: :hamburger:
Następnie erpy wyświetlamy tak:
plot(erp);
Zauważcie, że na osi x mamy numer próbki czasowej, a nie czas. Funkcja plot
nie ma pojęcia o czasie naszych kolejnych próbek! :scream: Możemy jej jednak informacje o wartościach na osi x podać. Wygląda to wtedy tak: plot(x,y)
, gdzie x
to informacja o pozycji na osi x, a y
to pozycja na osi y. W naszym wypadku oś x to czas, oś y natomiast - średnie napięcie. Informacja o czasie jest przechowywana w EEGlabie w EEG.times
, piszemy więc:
plot(EEG.times, erp);
Mamy więc opisany na osi x czas, ale zechcieliśmy na jednym wykresie namalować dwa różne erpy - jeden dla warunku face_0
a drugi dla warunku face_180
- oba dla elektrody E23
. Najpierw policzmy erpy:
erp_0 = erpy_daj(EEG, 'E23', 'face_0');
erp_180 = erpy_daj(EEG, 'E23', 'face_180');
Ok, teraz wyplotujmy je na jednym wykresie, erp dla warunku twarze prosto (0) - na zielono; dla warunku do góry nogami (180) - na czerwono. Najpierw otwórzmy nowe okienko i poprośmy matlaba aby nie nadpisywał tzn. pozwalał na dorysowywanie kolejnych linii (erpów w naszym wypadku) bez usuwania poprzednich:
figure
hold on
Aby wyplotować dane w kolorze używamy komendy plot(y, c)
lub plot(x, y, c)
gdzie c
to tekstowy skrót koloru: 'r'
- czerwony, 'g'
- zielony (itp.). A więc możemy zrobić tak:
plot(EEG.times, erp_0, 'g');
plot(EEG.times, erp_180, 'r');
Czad, nie? :smile:
Widmo nie takie straszne
- co i jak w interfejsie:
plot
:arrow_right:channel spectra and maps
eegh
prawdę Ci powie: funkcjapop_spectopo
- co nam zwraca
pop_spectopo
- rysujemy widma, kredki każdy ma
kodem będzie to tak:
Liczymy widmo przed i po:
widmo_pre = pop_spectopo(EEG, 1, [-500 0], 'EEG', 'plot', 'off');
[widmo_post, freq] = pop_spectopo(EEG, 1, [0 500], 'EEG', 'plot', 'off');
Plotujemy:
elec = 30; % tutaj numer elektrody
figure; hold on;
plot (freq, widmo_pre(elec, :), 'g');
plot (freq, widmo_post(elec, :),'r');
legend ('przed bodzcem', 'po bodzcu');
Jeżeli chcecie zrobić wykres widma tylko dla jakiegoś zakresu częstotliwości, przyda się:
% tutaj ustawiamy zakres częstotliwości:
rng = [0, 75];
% szukamy próbek widma najbliższych zakresowi, jakiego chcemy:
cnt = @(x) x(1) : x(2);
freq_rng = cnt( arrayfun(@(x) find(min(abs(freq - x)) == abs(freq - x)), rng) );
% plotujemy widmo tak:
plot(freq(freq_rng), widmo(elec, freq_rng), 'g');