Pracownia EEG/Potencjały wywołane

Pracownia EEG / Potencjały wywołane

Wstęp

Potencjały wywołane

Wstęp

Potencjały wywołane EEG (ang. evoked potentials, EP) są śladami odpowiedzi mózgu na bodźce. W przypadku MEG poprawnie mówi się o polach wywołanych (ang. evoked fields, EF). Zwykle reakcja ta jest mała (wyjątek stanowią kompleksy K) i w zapisie pojedynczej realizacji reakcji na bodziec najczęściej niewidoczna wśród czynności pochodzącej od wielu innych procesów zachodzących w tym samym czasie w mózgu. Ich wyodrębnienie z tak zwanego tła EEG/MEG, czyli manifestacji elektrycznej innych, trwających w tym samym czasie w mózgu procesów (aktywny prąd niezależny), wymaga w dniu dzisiejszym zapisu odpowiedzi na szereg powtórzeń tego samego bodźca (rys. 1). Pierwsze potencjały wywołane pokazał pod koniec lat 40-tych Dawson (Dawson 1947) wykonując superpozycję pojedynczych realizacji na kliszy fotograficznej. Z czasem skomplikowane urządzenia analogowe zastąpione zostały przez cyfrowe uśrednianie kolejnych fragmentów EEG, zsynchronizowanych według momentu wystąpienia bodźca, łatwo realizowane za pomocą komputera.

Uśrednianie potencjałów wywołanych. a) ciągły zapis EEG z wyróżnionymi momentami wystąpienia bodźca (w tym przypadku słuchowego) b) kilkadziesiąt kolejnych odcinków, wyciętych jako sekunda EEG od momentu wystąpienia kolejnych bodźców, ustawione jeden pod drugim c) uśredniony potencjał wywołany — widać m. in. załamek P300 ok. 300 milisekund po bodźcu

Techniki te opierają się na założeniu, że zawarta w EEG odpowiedź mózgu na każdy z kolejnych bodźców jest niezmienna, a EEG odzwierciedlające pozostałe procesy traktowane jest jak nieskorelowany z nią proces stochastyczny. Zależnie od rodzaju potencjałów wywołanych, założenia te są mniej lub bardziej nieuzasadnione; podważa je choćby powszechnie znany efekt habituacji, polegający na osłabieniu późnych potencjałów wywołanych kolejnymi powtórzeniami bodźca.

Istota potencjałów wywołanych jest przedmiotem otwartej dyskusji i dziesiątków prac, dotykających od lat podstawowych w tej dziedzinie pytań: czym jest potencjał wywołany, który widzimy w uśrednionym przebiegu? Czy naprawdę wynika z deterministycznie powtarzanej, jednakowej odpowiedzi pojawiającej się po każdym powtórzeniu bodźca niezależnie od „tła”, czy może wynika z reorganizacji faz tego właśnie „tła” EEG, czy może kombinacji tych dwóch efektów? Jest to wspaniałe pole dla zastosowań zaawansowanych metod modelowania i analizy sygnałów. W ostatnich latach powstają dziesiątki prac na ten temat, a od czasu do czasu również krytyczne artykuły wykazujące, że dotychczasowe odkrycia są raczej artefaktami stosowanych metod a nie wynikają z własności analizowanych danych (Yeung et al., 2004). Jak widać jest tu wciąż bardzo wiele do zrobienia, gdyż zrozumienie tego mechanizmu stoi na drodze do „świętego Graala”, którym w tej dziedzinie jest wyodrębnienie pojedynczych potencjałów wywołanych.

Niezależnie od tego, w neurofizjologii klinicznej nazwą potencjał wywołany określa się krzywą widoczną po uśrednieniu odpowiedzi na kilkanaście do kilku tysięcy bodźców (Szelenberger 2000). Gromadzona od dziesięcioleci wiedza o behawioralnych i klinicznych korelatach potencjałów wywołanych opiera się na rozpoznawaniu w przebiegach uśrednionych tak zwanych załamków, czyli przejściowych wzrostów lub spadków potencjału (na przykład na rys. 1 widać załamek P300). Nazwy załamków składają się zwykle z litery „P” (od ang. positive), jeśli wychylenie jest dodatnie, lub „N” (od ang. negative), jeśli wychylenie jest ujemne, oraz liczby. Uwaga:

• Nie ma niestety jednej ustalonej konwencji prezentacji potencjałów koniecznie trzeba więc zwracać uwagę na orientację wykresu, a wykonując wykres samemu - zamieszczać informację o orientacji kierunku dodatniego,gdyż kierunek dodatni może być wykresie skierowany do góry lub do dołu.
• Liczba następująca po N lub P określa przybliżoną liczbę milisekund od wystąpienia bodźca, czyli tak zwaną latencję. Nie można jednak ścisłego związku między tą liczbą a liczbą milisekund np. potencjał P300 w zależności od szczegółów paradygmatu doświadczalnego może wystąpić znacznie później. np. w eksperymentach językowych zdarza się że P300 występuje później niż N400 :-).

Rys. 2 przedstawia schematycznie najważniejsze załamki rozpoznawane w potencjałach słuchowych, czyli wywoływanych bodźcem dźwiękowym. Dla potencjałów wzrokowych i somatosensorycznych istnieją podobne klasyfikacje.

Orientacyjny schemat załamków rozpoznawanych w słuchowych potencjach wywołanych, na podstawie (Szelenberger 2000). Skala czasu logarytmiczna. Najszybsze (do 12 ms) składowe egzogenne to potencjały pniowe (BAEP), oznaczane jako fale I-VII. Litery „P” i „N” oznaczają dodatnie i ujemne wychylenia związane z dalszymi załamkami

Najszybsze (czyli o najmniejszej latencji) składowe to potencjały egzogenne, odzwierciedlające wstępne fazy przekazu informacji. W potencjałach słuchowych są to potencjały pnia (ang. brainstem auditory evoked potentials, BAEP), składające się z siedmiu fal występujących pomiędzy 1 a 12 ms od bodźca. Są one generowane np. w nerwie słuchowym czy pniu mózgu, i wykorzystywane w klinicznej diagnostyce integralności dróg słuchowych. Podobnie jak inne wczesne potencjały o latencjach poniżej 100 ms, są praktycznie niezależne od stanu uwagi.

Po potencjałach o średniej latencji (między 12 a 50 ms) zaczynają się późne potencjały słuchowe (od 100 ms od podania bodźca), odzwierciedlające bardziej złożone reakcje na bodziec. Potencjały ok. 100 ms obserwowane w okolicach sensorycznych odzwierciedlają odpowiedź tejże kory na bodziec. Ich amplituda zależy od parametrów fizycznych bodźca oraz od uwagi, np. amplituda załamka N100 wzrasta w stanie skupienia uwagi czy pobudzenia emocjonalnego. Pierwszy i najpopularniejszy całkowicie endogenny załamek to P300, pojawiający się po rozpoznaniu bodźca oczekiwanego, na którym skupiamy uwagę: na przykład, jeśli w serii pojawiających się losowo liter mamy zliczać wystąpienia jednej z nich, to na średniej odpowiedzi na pojawienie się tej litery pojawi się załamek P300. Załamek ten występuje również w odpowiedzi na bodźce wyczekiwane, niespodziewane i ważne dla osoby (np. własne imię).

Uśrednianie w dziedzinie czasu

Jak już było wspomniane powyżej podstawową techniką stosowaną do analizy potencjałów wywołanych jest uśrednianie wielu realizacji odpowiedzi na bodziec wyrównanych względem jakiegoś charakterystycznego zdarzenia, np. momentu podania bodźca albo momentu behawioralnej reakcji na bodziec (np. wciśnięcie przycisku). Podejście to bazuje na trzech założeniach:

• Reakcja jest czasowo związana z bodźcem. Tzn. w kolejnych realizacjach występuje ona zawsze z tym samym opóźnieniem.
• Reakcja na bodziec skutkuje zawsze pojawieniem się w sygnale składowej o stałym kształcie (mówimy o stałej morfologii).
• Spontaniczna czynność EEG — ta która nie dotyczy przetwarzania interesującego nas bodźca — jest niezależnym, stacjonarnym szumem o średniej zero.

Zgodnie z powyższymi założeniami mierzony w i-tej realizacji sygnał można wyrazić jako:

[math]x_i(t) = s(t) + n_i(t)[/math]

Uśrednianie po N realizacjach daje:

[math]\bar x (t) = \frac{1}{N} \sum_{i=1}^N x_i(t) = \frac{1}{N} \left(N s(t) + \sum_{i=1}^N n_i(t) \right) [/math]

Wartość oczekiwana średniego sygnału wynosi:

[math]\mathrm{E}\left[ \bar x(t) \right] = s(t) [/math]

gdyż dla szumu o średniej zero mamy:

[math]\mathrm{E}\left[\frac{1}{N} \sum_{i=1}^N n_i(t)\right] = 0[/math].

Wariancja [math]\bar x(t)[/math] wynosi:

[math]\sigma^2_{\bar x(t)} = \mathrm{E} \left[ \left( \frac{1}{N} \sum_{i=1}^N n_i(t)\right)^2 \right] \approx \frac{1}{N} \sigma^2_{n(t)}[/math]

gdyż [math]s(t)[/math] jest deterministyczne.

Dla potencjałów pojawiających się w pierwszych kilkudziesięciu milisekundach po bodźcu model ten można uznać za poprawny, gdyż wykazują one zależność głównie od parametrów fizycznych bodźca. Dla późniejszych składowych staje się on coraz bardziej wątpliwy ze względu na występujące korelacje ze stanem (np. uwagi) badanego.

Dla późniejszych składowych można by postulować bardziej ogólną wersję modelu potencjału wywołanego:

[math]x_i(t) = s_i(t) + n_i(t)[/math]

Ten model uwidacznia, że do opisu potencjału wywołanego potrzebna jest nie tylko średnia, ale i wyższe momenty rozkładu (np. drugi moment centralny czyli wariancja).

Wariancja w przypadku szumu skorelowanego

W poprzednim paragrafie pokazaliśmy, że w przypadku niezależnego szumu wariancja uśrednionego sygnału (w danej chwili czasu) maleje jak [math]\frac{1}{N}[/math]. Jednak w przypadku gdy kolejne próbki szumu są skorelowane (np. występuje silna aktywność rytmiczna alfa) sytuacja się komplikuje. Intuicyjnie łatwo możemy to sobie wyobrazić w granicznym przypadku. Załóżmy, że tło dla potencjałów o stałej morfologii stanowi sinusoida. Jeśli będziemy podawać bodźce w odstępach będących wielokrotnością okresu tej sinusoidy to w uśrednianych fragmentach sygnału owa sinusoida jest tak samo powtarzalna jak badany potencjał i uśrednianie nie prowadzi do poprawy stosunku amplitudy sygnału do amplitudy tła.

Stopień zależności pomiędzy próbkami można zmierzyć przy pomocy funkcji autokorelacji. Jeśli funkcję autokorelacji przybliżymy przez:

[math]R_{xx}(\tau) = \sigma^2 \exp(-\beta |\tau|) \cos(2 \pi f_0 \tau)[/math]

gdzie [math]\sigma^2[/math] jest wariancją szumu, [math]f_0[/math] jest średnią częstością aktywności rytmicznej, [math]\beta/ \pi[/math] jest szerokością pasma tej czynności, zaś [math]\tau[/math] opóźnieniem, to wariancję potencjału wywołanego można wyrazić jako:

[math]\sigma^2_{\bar x(t)} = \frac{\sigma^2}{N} \left[ \frac{1-\exp(-2 \beta T)}{1 - 2 \exp(- \beta T) \cos(2 \pi f_0 T) +\exp(-2 \beta T)} \right] [/math]

gdzie [math]T[/math] jest odstępem między bodźcami. Widać stąd, że obecność czynności rytmicznej w tle wpływa na stosunek sygnału do szumu. Z powyższego równania widać, że stosunek wariancji dąży do [math]\frac{\sigma^2}{N}[/math] gdy [math]\beta T[/math] staje się duże.

Periodyczna stymulacja może ponadto prowadzić do wzbudzenia czynności rytmicznej o częstości równej częstości pobudzania bądź jej harmonicznej.

Wariancja latencji

Jedną z form zmienności potencjału wywołanego w pojedynczych realizacjach jest wariancja latencji.

Czerwona linia: Uśrednianie 100 potencjałów wywołanych o latncjach z rozkładu normalnego N(300,30 ), Niebieskie linie: przykładowe pojedyncze realizacje wzięte do średniej

Zmienność latencji prowadzi do:

• zmniejszenia amplitudy [math]\bar x[/math]
• wzrostu rozciągłości czasowej załamka

O tych konsekwencjach zmienności latencji należy w szczególności pamiętać porównując średnie potencjały uzyskane w różnych warunkach eksperymentalnych. Gdy zauważymy, że średnie potencjały różnią się amplitudą należy rozważyć możliwość, że w istocie warunki eksperymentalne wpływają na różnice w wariancji latencji.

Naturalne jest, że zmienność latencji jest bardziej widoczna w przypadku załamków trwających krótko niż w przypadku załamków trwających długo. Wariancja latencji rożnie wraz z jej wartością, stąd wśród późnych składowych potencjałów wywołanych nie obserwuje się komponentów krótkotrwałych.

Habituacja

Zmienność w odpowiedziach na kolejno występujące po sobie bodźce można podzielić na komponent stochastyczny i deterministyczny. Komponent deterministyczny odzwierciedla zmianę czułości układu nerwowego na kolejne bodźce. Nowy bodziec często wywołuje reakcje pobudzenia widoczną w wielu parametrach takich jak przewodnictwo skóry, rytm serca i reaktywność kory mózgowej. Powszechnie obserwowanym zjawiskiem jest spadek reaktywności na powtarzające się bodźce zwany habituacją. Niekiedy obserwuje się także zjawisko przeciwne — wzrost reaktywności na bodziec — zwane sensytyzacją.

Klasycznym sposobem na badanie tych zjawisk w potencjałach wywołanych jest uśrednianie w podzbiorach lub uśrednianie blokowe. Wymaga ono zmiany paradygmatu rejestracji danych. Załóżmy, że interesującym nas bodźcem jest błysk światła. Chcemy zaobserwować ewentualne zmiany amplitudy i latencji w ciągu 100 potencjałów będących reakcją na błysk powtarzający się co 500 ms.

Uśrednianie w podzbiorach

W metodzie tej musimy założyć, że proces habituacji jest powolny, tzn. zmiany pomiędzy kolejnymi reakcjami są bardzo małe. W takim przypadku możemy podzielić długą serię bodźców na ciąg pod-serii i uśrednienie przeprowadzić dla każdej pod-serii osobno.

Ilustracja uśredniania w podzbiorach

Uśrednianie blokowe

W tym celu eksperyment wykonujemy w wielu seriach. Jedna seria składa się ze 100 błysków powtarzanych co 500 ms. Po serii następuje przerwa trwająca np. 5 s (długość przerwy wyraża nasze oszacowanie czasu potrzebnego na zanik habituacji). Po przerwie seria jest powtarzana. Powtórzeń serii wykonujemy N. Następnie uśredniamy fragmenty sygnałów odpowiadające sobie w kolejnych seriach.

W tym paradygmacie zakładamy, że każda seria rejestrowana jest w tych samych warunkach. Zaniedbujemy efekty związane ze zmęczeniem i zmianami w np. stopniu koncentracji badanego.

Ilustracja uśredniania w blokach

Załamki w krzywej potencjału wywołanego vs ukryte komponenty

Komponenta potencjału wywołanego to bardzo ważne ale jednocześnie źle określone pojęcie. Krzywa otrzymana w wyniku uśredniania potencjałów wywołanych przedstawia serię załamków dodatnich i ujemnych, ale trzeba być świadomym, że ten przebieg napięcia rejestrowany przez elektrody przedstawia tak naprawdę sumę potencjałów pochodzących od wielu dość niezależnie działających źródeł w mózgu. Przyczynek pochodzący od jednego funkcjonalnego źródła będziemy określać jako komponentę.

Komentarza wymaga też pojęcie źródła funkcjonalnego. Czynność elektryczną synchronicznie działających neuronów piramidalnych w określonym fragmencie kory mózgowej z pewnej odległości można przybliżyć potencjałem dipola prądowego. Taki aktywny fragment kory stanowi źródło dające przyczynek do mierzonego na powierzchni głowy potencjału. W najprostszej sytuacji taki przyczynek może być obserwowany jako komponenta. Mózg jest jednak obiektem niezwykle złożonym i gęsto połączonym. Wiele obszarów korowych może być funkcjonalnie połączonych i aktywowanych synchronicznie w celu realizacji określonego zadania — o takim funkcjonalnym zespole też możemy myśleć jako o pojedynczym źródle, choć może ono mieć złożoną strukturę przestrzenną. Często zatem myślimy o komponentach ukrytych potencjału wywołanego jako o aktywności takich funkcjonalnych zespołów.

Po tym wstępie widać już chyba, że do zrozumienia co dzieje się w czasie reakcji na bodziec, którą obserwujemy w potencjale wywołanym, bardziej pożądane byłoby badanie komponent ukrytych niż poszczególnych załamków krzywej potencjału wywołanego, która stanowi ważoną sumę owych komponentów. Zanim przejdziemy do omówienia technik służących „wyłuskiwaniu” komponent prześledźmy kilka symulacji i wnioski z nich płynące.

Załamki i komponenty nie są tym samym.

Tę samą krzywą można otrzymać na nieskończenie wiele sposobów przez sumowanie rożnych kombinacji domniemanych komponentów. Przykładowy rysunek pokazuje, jak tę samą krzywą można złożyć z dwóch rożnych zestawów komponentów (kod wytwarzający przykładową symulację pokazany jest poniżej).

Panel a) przedstawia hipotetyczną krzywą ERP, panel b) jedna z możliwych dekompozycji sygnału w panelu a), panel c) inna możliwa dekompozycja sygnału z panelu a)

# -*- coding: utf-8 -*-
from __future__ import division
import numpy as np
from matplotlib.pyplot import plot, show ,subplot, figure, title, xlabel, ylabel, ylim

Fs = 1
t = np.arange(0,700)

C1 = np.exp(-(t-100)**2/(2*30**2))
C2 = np.exp(-(t-200)**2/(2*50**2))
C3 = np.exp(-(t-350)**2/(2*100**2))
ERP = 0.5*C1-C2+C3

C2p = np.concatenate((C2[0:200], np.ones(30), C2[200:700-30] ))
C3p = ERP - 0.5*C1 + C2p
ERP_v2 = 0.5*C1-C2p+C3p

subplot(3,1,1)
plot(t,ERP, t[::10],ERP_v2[::10],'g.')
ylabel('a)',rotation='horizontal')
ylim([-1, 1.5])

subplot(3,1,2)
plot(t,0.5*C1,t,-C2,t,C3)
ylabel('b)',rotation='horizontal')
ylim([-1, 1.5])

subplot(3,1,3)
plot(t,0.5*C1,t,-C2p,t,C3p)
ylabel('c)',rotation='horizontal')
ylim([-1, 1.5])

show()

Zmiana jednej komponenty może powodować zmiany w amplitudzie i latencji więcej niż jednego załamka.

Poniższy rysunek przedstawia sytuację gdy zmiana amplitudy jednego komponentu może znacząco wpływać na amplitudę i latencję innych załamków.

Panel a) przedstawia dwie hipotetyczne krzywe ERP, panel b) panel b i c komponenty z których powstały owe krzywe. Zmiana amplitudy komponentu oznaczonego kolorem powoduje zmianę amplitudy oznaczonych załamków ale także zmianę latencji wcześniejszego załamka.

Poniżej kod odtwarzający powyższy przykład:

# -*- coding: utf-8 -*-
from __future__ import division
import numpy as np
from matplotlib.pyplot import plot, show, subplot, figure, title, xlabel, ylabel, ylim

Fs = 1
t = np.arange(0,700)

C1 = np.exp(-(t-100)**2/(2*30**2))
C2 = np.exp(-(t-200)**2/(2*50**2))
C3 = np.exp(-(t-350)**2/(2*100**2))
ERP1 = 0.5*C1-C2+C3
ERP2 = 0.5*C1-C2+1.5*C3


subplot(3,1,1)
plot(t,ERP1,'b', t,ERP2,'g')
plot(t[np.argmin(ERP1)],np.min(ERP1),'r.' )
plot(t[np.argmax(ERP1)],np.max(ERP1),'m.' )
plot(t[np.argmin(ERP2)],np.min(ERP2),'r.' )
plot(t[np.argmax(ERP2)],np.max(ERP2),'m.' )

ylabel('a)',rotation='horizontal')
ylim([-1, 1.6])

subplot(3,1,2)
plot(t,0.5*C1,'k',t,-C2,'k',t,C3,'b')
ylabel('b)',rotation='horizontal')
ylim([-1, 1.6])

subplot(3,1,3)
plot(t,0.5*C1,'k',t,-C2,'k',t,1.5*C3,'g')
ylabel('c)',rotation='horizontal')
ylim([-1, 1.6])

show()

Instrukcje do ćwiczeń

Zadanie 1: uśrednianie symulowanych sygnałów

w tym ćwiczeniu chcemy pokazać efekty:

1. uśredniania sygnałów związanych fazowo z bodźcem i nieskorelowanych z bodźcem / wpływ liczby realizacji,
2. uśredniania sygnałów z fluktuacją latencji,
3. wpływ filtrów o różnych funkcjach odpowiedzi impulsowej na morfologię uśrednionego potencjału.

Sygnały testowe

Proszę zaimplementować funkcje generujące następujące przebiegi czasowe (oznaczenia: F_s — częstość próbkowania, T — czas trwania w sekundach, analogicznie do tego):

• sinusoida o zadanej częstości f i fazie [math]\phi[/math]: y = sin(f, phi, Fs, T),
• funkcja Gaussa o zadanym położeniu i szerokości y = g_1(t0, sigma, Fs, T)

[math]g_1(t) = \exp\left(-\frac{(t-t_0)^2}{2 \sigma^2}\right)[/math],

• pochodną funkcji Gaussa o zadanym położeniu i szerokości y = g_2(t0, sigma, Fs, T),
• połowę funkcji Gaussa o zadanym położeniu i szerokości y = g_3(t0, sigma, Fs, T),

[math]g_{3}(t) = \left\{ \begin{array}{lll} 0 & \mathrm{dla}& t\lt t_0 \\ \exp\left(-\frac{(t-t_0)^2}{2 \sigma^2}\right)& \mathrm{dla}& t \ge t_0 \end{array} \right. [/math].

W dalszych zadaniach tej sekcji przyjmiemy:

• T = 1 s,
• F_s = 128 Hz,
• t₀ = 0,3 s,
• σ = 0,02 s.

Uśrednianie w modelu addytywnym: sygnał + szum

Dla funkcji [math]g_1[/math], [math]g_2[/math] i [math]g_3[/math] wykonaj N symulacji sygnału opisanego daną funkcją z dodanym szumem gaussowskim (numpy.random.randn). Zaobserwuj zmianę wariancji uśrednianego sygnału wraz z N.

Uśrednianie w modelu addytywnym: sygnał + szum + niezależna czynność rytmiczna

Powtórz poprzednie symulacje dokładając w każdej realizacji sinusoidę o losowo wybranej fazie. Symulacje wykonaj dla różnych częstości sinusoidy, takich których połówka okresu zbliżona jest do szerokości funkcji Gaussa, oraz znacząco od niej różnych.

Uśrednianie w modelu z resetem fazy

Zbuduj sygnał składający się z dwóch fragmentów sinusoidy o częstości f:

• pierwszy fragment ma losową fazę,
• drugi fragment ma fazę równą 0.

Wykonaj N realizacji takiego sygnału, przy czym każda realizacja ma częstość losowaną z rozkładu Gaussa o średniej 10 i pewnej wariancji (czyli [math]f \in \mathrm{N}(10,\sigma)[/math]). Uśrednij otrzymane realizacje. Zaobserwuj zależność średniej od liczby uśrednianych realizacji i od wariancji częstości.

Uśrednianie potencjału z fluktuacją latencji

Wygeneruj sygnały składające się z funkcji [math]g_1[/math], [math]g_2[/math] albo [math]g_3[/math] z fluktuującą latencją, tzn. z wartością [math]t_0[/math] pochodzącą z jakiegoś rozkładu prawdopodobieństwa. U nas dla ustalenia uwagi niech będzie to rozkład Gaussa. Wykonaj zestaw symulacji obrazujący zależność od liczby realizacji i od wielkości fluktuacji, mierzonej przez wariancję rozkładu. Czy jest związek pomiędzy wariancją fluktuacji [math]t_0[/math] a szerokością załamków?

Wpływ filtrów

1. Zaprojektuj następujące filtry (przyjmij częstość próbkowania 128 Hz):
   • dolnoprzepustowy: Butterwortha z częstością odcięcia: 30, 40, 50 Hz, i rzędach od 1 do 5 (funkcja butter w module scipy.signal),
   • górnoprzepustowy: Butterwortha z częstością odcięcia: 0,1, 0,5, 2 i 5 Hz, i rzędach od 1 do 5,
   • notch: Czybyszewa II rodzaju z pasmem zaporowym ustawionym na 50 Hz, i rzędach od 1 do 3 (funkcja cheby2),
2. Zbadaj funkcję odpowiedzi impulsowej i funkcję przenoszenia tych filtrów (w skali dB). Zaobserwuj jak długa jest odpowiedź impulsowa i jak daleko od częstości odcięcia zaburzone jest pasmo przenoszenia.
3. Zastosuj powyższe filtry do funkcji [math]g_1[/math], [math]g_2[/math] albo [math]g_3[/math]. Filtrowanie przeprowadź w jedną (lfilter) oraz w dwie strony (filtfilt). Zaobserwuj związek między zniekształceniem wprowadzanym przez filtr a funkcją odpowiedzi impulsowej filtru.
4. W ćwiczeniu tym posłuż się symulowanymi sygnałami uzyskanymi z modelu addytywnego: „sygnał+szum”. Dla funkcji [math]g_1[/math], [math]g_2[/math] albo [math]g_3[/math] wykonaj N = 60 realizacji sygnału opisanego daną funkcją z dodanym szumem gaussowskim (numpy.random.randn). Każdą z realizacji przefiltruj, a następnie uśrednij uzyskane sygnały.

Zadanie 2: rejestracje i analiza potencjałów wzrokowych

Eksperyment 1. Obserwacja bodźca wzrokowego

Proszę zapoznać się z typowym wyglądem załamków w przypadku potencjałów wzrokowych, jaki jest przedstawiony na rysunku (cały artykuł jest dostępny tu).

Rejestracja

• Paradygmat: z odległości około 40 cm badany obserwuje wyświetalny na ekranie biały kwadrat. Bodźce trwają 300 ms i powtarzają się co 0,5 sekundy. Bodziec poprzedza punkt fiksacji pokazywany przez 300 ms. Natężenie bodźca ma przyjąć jedną z 4 wartości (manipulujemy jego przeźroczystością: 1, 0,7, 0,3 i 0,1, gdzie 1 oznacza brak przeźroczystości). Każda z intensywności bodźca ma pojawić się 100 razy podczas badania, bodźce podawane mają być w sposób losowy. Zadanie poprzedza krótka instrukcja, a kończy podziękowanie.

• Pomiar
  • Zakładamy czepek i elektrody w systemie 10-20, dbamy o to by opory pomiędzy elektrodami były poniżej 10 kΩ i różnice pomiędzy oporami różnych elektrod nie przekraczały 20%.
    
  • Elektrodę GND mocujemy w otworze czepka odpowiadającym pozycji AFz (na linii centralnej 10% w stronę nosa od elektrody Fz)
  • Montujemy także elektrody M1 i M2 i FCz.

Analiza

1. W sygnale wyszukujemy wzrokowo odcinki z artefaktami i oznaczamy je tagami w SVAROGu.
2. Filtrujemy sygnał w paśmie 1-20 Hz.
3. Z ciągłego zapisu wycinamy fragmenty od −300 ms do +800 ms (zero oznacza moment wystąpienia bodźca). W czasie wycinania badamy wycinane fragmenty pod względem zakresu amplitudy (wszystkie kanały powinny mieć amplitudę w granicach ±75 μV) i zawartości artefaktów od mrugania (w tym celu trzeba oglądając sygnał w SVAROG-u ustalić próg na amplitudę sygnału w kanałach Fp1 i Fp2).
4. W przypadku, gdy dane są zanieczyszczone wieloma mrugnięciami (zwłaszcza tylnych okolicach), należy skorzystać z metody ICA w SVAROGU, informacje tutaj: [1])
5. Porównujemy zestawy danych otrzymane po wzrokowym i automatycznym (metodą opisaną powyżej) usunięciu artefaktów.
6. Uśredniamy wycięte sygnały.
7. Rysujemy średni potencjał we wszystkich kanałach.
8. Pierwszy sposób analizy: Identyfikujemy załamki P1 i N2. Mierzymy amplitudę P1 jako różnicę między ekstremalną wartością załamka P1 a ekstremalną wartością załamka N2.
9. Drugi sposób analizy: Identyfikujemy załamki P1 i N2. Mierzymy amplitudę P1 w stosunku do zera.
10. Porównujemy uzyskane amplitudy w 4 warunkach za pomocą odpowiedniego testu statystycznego (porównanie warunków ze skrajnymi wartościami amplitudy) — zastosuj test permutacyjny (opisany poniżej) lub test t-Studenta/nieparametryczny odpowiednik).

• • Opis testu permutacyjnego

Pytanie: Czy intensywność bodźca wpływa na amplitudę potencjału P1? Jak stwierdzić czy amplitudy załamków P1 istotnie różnią się pomiędzy seriami?

Musimy przeprowadzić test statystyczny. Możemy skonstruować go np. tak: Mamy dwa warunki I i II. Jako statystykę przyjmiemy różnicę pomiędzy amplitudą załamka zmierzoną w serii z dwoma diodami (P2_II) a amplitudą zmierzoną w serii z jedną diodą (P2_I):

S = P2_II − P2_I

Postawmy hipotezę zerową, że nie ma różnicy w amplitudzie pomiędzy warunkiem II a I oraz hipotezę alternatywną, że różnica taka występuje.

H₀: S = 0

H₁: S ≠ 0

Zgodnie z H₀ dla wartości S nie ma różnicy czy realizacja została zarejestrowana w warunkach I czy II, zatem można losowo poprzypisywać numery warunków do realizacji i w ten sposób uzyskać rozkład statystyki S pod warunkiem prawdziwości H₀, czyli obliczyć średnie potencjały dla „pomieszanych” tagów, wyliczyć średnie potencjały dla sztucznie wygenerowaych grup, i wyliczyć różnicę amplitud dla tych średnich, czyli wartość statystyki dla hipotezy H₀.

Zaimplementuj test permutacyjny zgodny z powyższą koncepcją.

Eksperyment 2: rejestracja i analiza potencjału P3

Rejestracja

1. Paradygmat:
   • z odległości około 40 cm badany obserwuje kwadrat, który występuje w dwóch kolorach.
   • Bodźce trwają 300 ms, poprzedzane punktem fiksacji wyświetlanym na 300 ms, przerwa między próbami jest losowana i wynosi od 0.8-1.2 s.
   • Proporcja wystapień bodźców koloru pierwszego i drugiego wynosi 1:5.
   • Zadaniem badanego jest liczenie kwadratów w rzadkim kolorze.
   • Wykonujemy tyle realizacji aby zarejestrować 100 powtórzeń bodźców rzadkich.
2. Zakładamy czepek i elektrody w systemie 10-20, dbamy o to by opory pomiędzy elektrodami były poniżej 5 kΩ i różnice pomiędzy oporami różnych elektrod nie przekraczały 20%.
   • Elektrodę GND mocujemy w otworze czepka odpowiadającym pozycji AFz (na linii centralnej 10% w stronę nosa od elektrody Fz)
   • Montujemy także elektrody M1 i M2, będziemy ich używać w montażu i w czasie analizy jako referencji typu połączone uszy.
   • Częstość próbkowania 512 Hz.

Analiza

1. Filtrujemy sygnał w paśmie 1-20 Hz.
2. Z ciągłego zapisu wycinamy fragmenty od −300 ms do +900 ms (zero oznacza moment wystąpienia bodźca). W czasie wycinania badamy wycinane fragmenty pod względem zakresu amplitudy (wszystkie kanały powinny mieć amplitudę w granicach ±75 μV) i zawartości artefaktów od mrugania (w tym celu trzeba oglądając sygnał w SVAROG-u ustalić próg na amplitudę sygnału w kanałach Fp1 i Fp2).
3. Uśredniamy osobno realizacje dla każdego koloru kwadratu.
4. Nakładamy na siebie na rysunkach oba typy średnich.
5. Wyznaczamy amplitudę P3 i liczymy test statystyczny pomiędzy warunkami.

Zadanie 3. Analiza Heartbeat-Evoked Potential (HEP)

Celem tego zadania jest zaobserwowanie występowania ERP związanego z działaniem serca (szczegóły można doczytać tutaj: [2]. W zadaniu tym wykorzystany dane spoczynkowe zebrane podczas pierwszej rejestracji.

• w pierwszym kroku sygnał montujemy i filtrujemy;
• następnie zaznaczamy fragmenty sygnału, w których widzimy artefakty;
• w drugim robimy detekcję załamków R z EKG — po ich podstawie będziemy cieli sygnał (−200-600 ms po R);
• następnie uśredniamy otrzymane odcinki.