Laboratorium EEG/CSP
Laboratorium_EEG/BSS
Spis treści
Ś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 ~/.obci/scenarios/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 "~/.obci/scenarios/P300/Sygnal_configs" kopiujemy plik "amplifier.ini" do katalogu "~/.obci/scenarios/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
- w tym pliku "global_amplifier.ini" podmieniamy linijki (tak aby były to listy faktycznie wykorzystywanych kanałów) z:
- 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.
- edytujemy plik "~/.obci/scenarios/P300/kalibracjaP300_configs/clasifier.ini"
- zmieniamy linię
- edytujemy plik "~/.obci/scenarios/P300/kalibracjaP300_configs/clasifier.ini"
- ignore_channels = DriverSaw;AmpSaw;PO7;PO8
- na
- ignore_channels = DriverSaw;AmpSaw;A1;A2
- oraz linię:
- montage_channels = PO7;PO8
- na
- montage_channels = A1;A2
Przeprowadzenie badania:
- Uruchom scenariusz "Sygnał"
- Tworzy on w katalogu domowym plik o nazwie file_id_name.
- Uruchamiamy svaroga z terminala poleceniem svarog. W zakładce sygnały on-line odnajdujemy nazwę naszego scenariusza "Sygnal". Podłączamy się do niego i poprawiamy ewentualnie źle kontaktujące elektrody.
- Jak już jesteśmy zadowoleni z jakości sygnału to zatrzymujemy scenariusz "Sygnal" w obci.
- W pliku file_id_name znajduje się string, który stanowi rdzeń do tworzenia nazw plików, z których korzystają nasze scenariusze. Proszę zmienić ten string np. na: "test1"
- Uruchamiamy scenariusz "kalibracjaP300". Badany będzie oglądał interfejs z trzema literami A B C migającymi w losowej kolejności. Zadaniem jest zliczanie mignięć litery B.
- Po zakończeniu kalibracji uruchamiamy scenariusz Labirynt.
- Danych z kalibracji potrzebować będziemy kilka zestawów. Proszę powtórzyć kilku krotnie scenariusz kalibracjaP300. Przed każdym uruchomieniem trzeba zmienić string w pliku file_id_name np na "test???" gdzie ??? oznacza kolejne numery.
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 każdej 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 zaprezentowany jest przykładowy skrypt do cięcia danych wokół znaczników. Działa on z plikami zawartymi w archiwum:
Korzysta z funkcji pomocniczych dostępnych w dystrybucji obci w katalogu
- /usr/share/openbci/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.
Wybór i separacja cech
- przedstaw na rysunkach nałożone na siebie pojedyncze realizacje target i non-target po rzutowaniu na wetor w odpowiadający największej wartości własnej.
Filtry przestrzenne dla większej ilości warunków
FFDIAG
Analiza ERD/S z użyciem FFDIAG
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
- 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%.
- 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.
- Elektrodę GND mocujemy na pozycji AFz.
- Sygnał rejestrujemy z częstością 2048 Hz.
- 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
- Z sygnału wycinamy fragmenty od 0 do 5 sek. dla wszystkich elektrod położone nad korą słuchową.
- Dla każdej realizacji obliczamy widma metodą Welcha.
- Otrzymane zespolone widma uśredniamy po realizacjach.
- Sprawdzamy czy w uśrednionym widmie występuję maksimum w częstości modulacji tj. 40 Hz.
Wyznaczenie przebiegu czasowego ERD i ERS
- Zaprojektuj filtry pasmowo przepustowe (Chebyszewa 2 rodzaju) zgodne z wyznaczonym pasmem. Zbadaj funkcje przenoszenia i odpowiedzi impulsowej.
- Powycinaj sygnały od -5 do +10 sekund (wszystkie kanały). Przefiltruj każdą realizację.
- Oblicz moc chwilową za pomocą transformaty Hilberta (kwadrat modułu transformaty Hilberta).
- Uśrednij moc chwilową po realizacjach.
- Oblicz względną zmianę mocy chwilowej względem czasu -4 do -2. W ten sposób otrzymasz przebieg ERD i ERS w czasie.
- 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.
