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

Z Brain-wiki
 
(Nie pokazano 29 pośrednich wersji utworzonych przez tego samego użytkownika)
Linia 18: Linia 18:
 
</math>
 
</math>
  
Uzyskujemy w ten sposób tzw. model wielokanałowy (wielozmienny, ang. ''multichannel'', ''multivariate'', MVAR).
+
Uzyskujemy w ten sposób tzw. model wielokanałowy (wielozmienny, ang. ''multichannel'', ''multivariate'', MVAR). Należy odróżniać taki model od wersji wielowymiarowej, w której opisywany proces ''X'' zależy od wielu zmiennych, np. w wersji dwuwymiarowej szukalibyśmy ''X''(''x'', ''y''). Modelami wielowymiarowymi nie będziemy się tu zajmować.
  
 
Podobnie jak w przypadku jednokanałowym możemy przetransformować model do przestrzeni częstości uzyskując zależność
 
Podobnie jak w przypadku jednokanałowym możemy przetransformować model do przestrzeni częstości uzyskując zależność
Linia 26: Linia 26:
 
</math>
 
</math>
  
gdzie, podobnie jak poprzednio, odpowiednie wielkości stają się wektorami ''k''-wierszowymi i macierzami rozmiaru ''k&times;''k'':
+
gdzie, podobnie jak poprzednio, odpowiednie wielkości stają się wektorami ''k''-wierszowymi i macierzami rozmiaru ''k''&times;''k'':
  
 
<math>
 
<math>
Linia 92: Linia 92:
 
Macierz ''H'' produkuje z transformat Fouriera sygnałów szumowych (posiadających z definicji widmo płaskie, niezależne od częstości) widmo badanych sygnałów. Zawiera więc ona potrzebną nam informację o zależnościach miedzy naszymi sygnałami w dziedzinie częstości. Jest to macierz niesymetryczna, jej element ''H<sub>ij</sub>''(''f'') mówi o intensywności wpływu sygnału ''X<sub>j</sub>'' na sygnał ''X<sub>i</sub>'' w częstości ''f''.  
 
Macierz ''H'' produkuje z transformat Fouriera sygnałów szumowych (posiadających z definicji widmo płaskie, niezależne od częstości) widmo badanych sygnałów. Zawiera więc ona potrzebną nam informację o zależnościach miedzy naszymi sygnałami w dziedzinie częstości. Jest to macierz niesymetryczna, jej element ''H<sub>ij</sub>''(''f'') mówi o intensywności wpływu sygnału ''X<sub>j</sub>'' na sygnał ''X<sub>i</sub>'' w częstości ''f''.  
  
Z elementów macierzy ''H'' zbudowana jest kierunkowa funkcja przejścia (ang. directed transfer function, DTF). Istnieje kilka wariantów tej funkcji, poniżej prezentowane są wzory opisujące wersję znormalizowaną i nieznormalizowaną.
+
Z elementów macierzy ''H'' zbudowana jest kierunkowa funkcja przejścia (ang. ''directed transfer function'', DTF). Istnieje kilka wariantów tej funkcji, poniżej prezentowane są wzory opisujące wersję znormalizowaną i nieznormalizowaną.
  
 
Wersja znormalizowana
 
Wersja znormalizowana
Linia 103: Linia 103:
  
 
Dla obu wersji funkcji DTF wartości bliskie zeru mówią o braku związku między danymi sygnałami.
 
Dla obu wersji funkcji DTF wartości bliskie zeru mówią o braku związku między danymi sygnałami.
 +
 +
Zauważmy jeszcze, że model AR traktuje wszystkie kanały modelowanego układu na raz (jednocześnie). Jest to sytuacja zupełnie inna niż w przypadku badania każdej pary kanałów osobno.
  
 
=Ćwiczenia=
 
=Ćwiczenia=
Linia 108: Linia 110:
 
==Wstęp do ćwiczeń==
 
==Wstęp do ćwiczeń==
  
Do ćwiczeń w tym rozdziale używać będziemy zestawu danych, które służyły w poprzednim rozdziale do wyznaczania komponentów ICA. Aby dostosować je do naszych celów dokonamy na nich następujących operacji:
+
Do ćwiczeń w tym rozdziale używać będziemy zestawu danych, które służyły w poprzednim rozdziale do wyznaczania komponentów ICA (http://www.fuw.edu.pl/~jarekz/LabEEG/Dane_do_ICA_alfa.tar.gz). Aby dostosować je do naszych celów dokonamy na nich następujących operacji:
 +
* wybierzemy kanały EEG;
 
* zastosujemy montaż do połączonych uszu (kanały A1 i A2);
 
* zastosujemy montaż do połączonych uszu (kanały A1 i A2);
 
* zmniejszymy częstość próbkowania z 512 do 128 Hz (stosując antyaliasingowy filtr dolnoprzepustowy funkcją filtfilt i biorąc co czwartą próbkę filtrowanego sygnału);
 
* zmniejszymy częstość próbkowania z 512 do 128 Hz (stosując antyaliasingowy filtr dolnoprzepustowy funkcją filtfilt i biorąc co czwartą próbkę filtrowanego sygnału);
* przefiltrujemy sygnał górnoprzepustowo z granicą odcięcia 1 Hz (stosując funkcję filtfilt).
+
* przefiltrujemy sygnał (po zmianie) górnoprzepustowo z granicą odcięcia 1 Hz (stosując funkcję filtfilt).
 +
 
 +
Aby dopasować współczynniki modelu MVAR do posiadanych danych użyjemy funkcji <tt>licz_wsp_AR.m</tt> zawartej w pakiecie: [http://www.fuw.edu.pl/~moira/prezentacja/MVAR_Matlab.zip MVAR_Matlab.zip] &mdash; należy go ściągnąć i wypakować do katalogu roboczego. Funkcja ta zwraca współczynniki ''A'' modelu MVAR oraz macierz wariancji szumu ''V''. Aby uzyskać macierz przejścia modelu ''H'' należy zaimplementować samemu procedurę transformacji współczynników do dziedziny częstości (patrz: [http://brain.fuw.edu.pl/edu/index.php/Pracownia_EEG/AR_1#Kilka_s.C5.82.C3.B3w_o_transformacji_Z tu]).
 +
 
 +
Przez &bdquo;szum&rdquo; rozumieć będziemy wektor wartości losowych zwracanych przez funkcję <tt>randn</tt>.
  
 
==Ćwiczenie 1==
 
==Ćwiczenie 1==
  
Z zestawu danych do obliczania ICA (poprzedni rozdział) wybierz jeden kanał EEG, zawierający wyraźną czynność alfa. Przytnij wybrany odcinek do długości 2000 próbek.
+
Z zestawu danych do obliczania ICA (poprzedni rozdział) wybierz jeden kanał EEG, zawierający wyraźną czynność alfa.
 +
<br><br>
 +
Przytnij wybrany odcinek do długości 2000 próbek. Podziel wybrany odcinek danych przez ich odchylenie standardowe.
 +
 
 
Wygeneruj dwa zestawy danych:
 
Wygeneruj dwa zestawy danych:
 
* Zestaw 1
 
* Zestaw 1
Linia 129: Linia 139:
 
==Ćwiczenie 2==
 
==Ćwiczenie 2==
  
Dla obu zestawów danych z ćwiczenia 1 sprawdź współczynniki dopasowanych modeli AR w dziedzinie czasu.
+
Dla obu zestawów danych z ćwiczenia 1 dopasuj dwukanałowe modele AR dla rzędów od 1 do 10. Sprawdź, ile wynoszą współczynniki ''a''<sub>12</sub> i ''a''<sub>21</sub> dopasowanych modeli AR w dziedzinie czasu dla ''j'' = 1...''p''.
 +
Dla każdego rzędu ''p'' oblicz dla kanałów 1 i 2 oraz 2 i 1 miarę DC [3]:
 +
 
 +
<math>
 +
\mathrm{DC}_{ij}^2=\sum_{m=1}^p{a_{ij}^2(m)}
 +
</math>
  
 
==Ćwiczenie 3==
 
==Ćwiczenie 3==
Linia 140: Linia 155:
 
Ustal optymalny rząd modelu AR (tym razem dwukanałowego) i oblicz macierz gęstości widmowej mocy oraz koherencji między tymi sygnałami. Narysuj moduł i fazę koherencji ''C''<sub>12</sub> i ''C''<sub>21</sub>.
 
Ustal optymalny rząd modelu AR (tym razem dwukanałowego) i oblicz macierz gęstości widmowej mocy oraz koherencji między tymi sygnałami. Narysuj moduł i fazę koherencji ''C''<sub>12</sub> i ''C''<sub>21</sub>.
  
Dla tego zestawu kanałów oblicz i narysuj normalizowaną i nienormalizowaną fukcję DTF.
+
Dla tego zestawu kanałów oblicz i narysuj normalizowaną i nienormalizowaną funkcję DTF.
  
 
Zmień fazę początkową drugiego sygnału. Jak zmienia się funkcja koherencji? Co dzieje się z funkcją DTF?
 
Zmień fazę początkową drugiego sygnału. Jak zmienia się funkcja koherencji? Co dzieje się z funkcją DTF?
  
 
==Ćwiczenie 4==
 
==Ćwiczenie 4==
 +
Ćwiczenie to opisuje problem tzw. wspólnego źródła. W takiej sytuacji prawdziwa analiza wielokanałowa wykazuje swoją przewagę nad badaniem zależności parami (dla każdej pary kanałów osobno).
  
 
Wygeneruj układ trzech sygnałów w następujący sposób:
 
Wygeneruj układ trzech sygnałów w następujący sposób:
  
 
     jako pierwszego kanału użyj sygnału z ćwiczenia 1;
 
     jako pierwszego kanału użyj sygnału z ćwiczenia 1;
     sygnał_w_drugim_kanale(t) = 0,4 * sygnał_z_pierwszego_kanału(t−1) + szum1;
+
     sygnał_w_drugim_kanale(t) = 0,4 * sygnał_z_pierwszego_kanału(t−1) + 0,4 * szum1;
     sygnał_w_trzecim_kanale(t) = 0,3 * sygnał_z_pierwszego_kanału(t−2) + szum2.
+
     sygnał_w_trzecim_kanale(t) = 0,3 * sygnał_z_pierwszego_kanału(t−2) + 0,4 * szum2.
  
 
Oblicz macierz koherencji zwyczajnych dla tego układu i na ich podstawie wyznacz zależności między kanałami. Powtórz to samo dla koherencji cząstkowych.
 
Oblicz macierz koherencji zwyczajnych dla tego układu i na ich podstawie wyznacz zależności między kanałami. Powtórz to samo dla koherencji cząstkowych.
  
 
Oblicz dla tego zestawu danych funkcje DTF.
 
Oblicz dla tego zestawu danych funkcje DTF.
 +
 +
Oblicz funkcję DTF tylko dla kanałów 2 i 3 (model dwukanałowy). Czym różni się wynik od odpowiedniej funkcji dla układu wszystkich trzech kanałów?
  
 
Wyniki wszystkich obliczeń przedstaw na rysunkach.
 
Wyniki wszystkich obliczeń przedstaw na rysunkach.
Linia 160: Linia 178:
 
==Ćwiczenie 5==
 
==Ćwiczenie 5==
  
Oblicz funkcje DTF dla wszystkich kanałów EEG z przygotowanego zestawu danych do ICA (dla pełnej długości w czasie każdego kanału).
+
Oblicz funkcje DTF dla wszystkich kanałów EEG z przygotowanego zestawu danych do ICA (dla pełnej długości w czasie każdego kanału). Zaobserwuj przepływy w paśmie częstości alfa.
 +
 
 +
==Ćwiczenie 6==
 +
 
 +
Oblicz funkcje DTF dla sygnałów zebranych podczas eksperymentu z ruchu palcem (Pracownia EEG: https://brain.fuw.edu.pl/edu/index.php/Pracownia_EEG/ERDS#Cwiczenia). Do analizy wybierz kanały położone nad korą ruchową rąk lewej i prawej: C3, Cz, C4, F3, F4. Oblicz wyniki dla poszczególnych realizacji eksperymentu, a następnie je uśrednij po realizacjach. Porównaj przepływy dla okresu referencyjnego (&minus;4..&minus;2 s) oraz dla okresu po ruchu (+0,5..+2,5 s).
 +
Porównaj również wyniki dla normalizowanej i nienormalizowanej wersji funkcji DTF.
 +
 
 +
Podziel cały wybrany odcinek danych eksperymentalnych (&minus;5..+5 s) na okna o długości 1 sekundy nakładające się na siebie z zakładką 0,5 s. Policz dla każdego okna nienormalizowane funkcje DTF (uśredniając po realizacjach). Narysuj otrzymane wyniki w postaci macierzy map czas-częstość obrazujących zmiany intensywności przepływów w czasie.
 +
 
 +
=Literatura=
 +
* Kamiński M, Blinowska KJ. A new method of the description of the information flow in brain structures. Biol Cybern 1991; 65:203–10.
 +
* Granger CWJ. Investigating causal relations in by econometric models and crossspectral methods. Econometrica 1969; 37:424–38.
 +
* Kamiński M, Ding M, Truccolo W, Bressler S. Evaluating causal relations in neural systems: Granger causality, directed transfer function and statistical assessment of significance. Biol Cybern 2001; 85:145–57.

Aktualna wersja na dzień 21:09, 23 kwi 2019

Laboratorium_EEG/MVAR

Wielokanałowe modele AR

Model AR opisuje wartość sygnału w chwili czasu t jako kombinację liniową jego wartości w chwilach poprzednich oraz szumu.

[math]X(t)=\sum_{j=1}^p {A(j)X(t-j)}+E(t) [/math]

Wzór powyższy może posłużyć do jednoczesnego opisu wielu sygnałów tworzących układ wielokanałowy. Mamy wtedy (dla k kanałów):

[math]X(t)=\left( \begin{array}{c} X_1(t)\\X_2(t)\\ \vdots\\X_k(t) \end{array}\right),\ E(t)=\left( \begin{array}{c} E_1(t)\\E_2(t)\\ \vdots\\E_k(t) \end{array}\right),\ A(j)= \left( \begin{array}{cccc} A_{11}(j) & A_{12}(j) & \hdots & A_{1k}(j)\\ A_{21}(j) & A_{22}(j) & \hdots & A_{2k}(j)\\ \vdots & \vdots & \ddots & \vdots \\A_{k1}(j) & A_{k2}(j) & \hdots & A_{kk}(j) \end{array}\right) [/math]

Uzyskujemy w ten sposób tzw. model wielokanałowy (wielozmienny, ang. multichannel, multivariate, MVAR). Należy odróżniać taki model od wersji wielowymiarowej, w której opisywany proces X zależy od wielu zmiennych, np. w wersji dwuwymiarowej szukalibyśmy X(x, y). Modelami wielowymiarowymi nie będziemy się tu zajmować.

Podobnie jak w przypadku jednokanałowym możemy przetransformować model do przestrzeni częstości uzyskując zależność

[math] A(f)X(f)=E(f) [/math]

gdzie, podobnie jak poprzednio, odpowiednie wielkości stają się wektorami k-wierszowymi i macierzami rozmiaru k×k:

[math] X(f)=\left( \begin{array}{c} X_1(f)\\X_2(f)\\ \vdots\\X_k(f) \end{array}\right),\ E(f)=\left( \begin{array}{c} E_1(f)\\E_2(f)\\ \vdots\\E_k(f) \end{array}\right),\ A(f)=\left( \begin{array}{cccc} A_{11}(f) & A_{12}(f) & \hdots & A_{1k}(f)\\ A_{21}(f) & A_{22}(f) & \hdots & A_{2k}(f)\\ \vdots & \vdots & \ddots & \vdots \\A_{k1}(f) & A_{k2}(f) & \hdots & A_{kk}(f) \end{array}\right) [/math]

Mnożąc powyższe równanie lewostronnie przez macierz A−1 otrzymujemy

[math] X(f)=A^{-1}(f)E(f)=H(f)E(f) [/math]

Macierz H nazywana jest macierzą przejścia (ang. transfer matrix) modelu.

Przyczynowość

Aby móc efektywnie opisywać zależności przyczynowe między sygnałami musimy w matematyczny sposób opisać, w jaki sposób tego typu zależności rozumiemy. W przypadku sygnałów biomedycznych będziemy musieli uwzględnić ich stochastyczny charakter. Zakładamy również, że na wartość danych w chwili obecnej mogą wpływać jedynie czynniki z chwil poprzednich (skutek nie może wyprzedzać przyczyny).

Przyczynowość Grangera

Szeroko stosowaną definicją przyczynowości między sygnałami jest definicja podana przez Grangera (Clive Granger, ekonomista i matematyk amerykański). Bazuje ona na przewidywalności szeregów czasowych. Wyobraźmy sobie, że próbujemy przewidzieć wartość sygnału (szeregu czasowego) X w chwili czasu t używając wartości tego samego sygnału zmierzonych w chwilach poprzednich, wziętych z pewnymi współczynnikami A:

[math] X(t)=\sum_{j=1}^\infty{A(j)X(t-j)}+E(t) [/math]

W powyższym wzorze wielkość E(t) jest uzyskanym przez nas błędem dopasowania.

Spróbujmy teraz przewidzieć wartość X(t), ale wzbogacając formułę po prawej stronie poprzednimi wartości innego sygnału Y

[math] X(t)=\sum_{j=1}^\infty{A(j)X(t-j)}+\sum_{j=1}^\infty{B(j)Y(t-j)}+N(t) [/math]

Uzyskujemy nowy błąd dopasowania N. Zasada sformułowana przez Grangera mówi, że jeśli po dodaniu poprzednich próbek drugiego sygnału nasze przewidywanie sygnału pierwszego się poprawiło, to kanał drugi (Y) jest przyczynowym dla pierwszego (X). Matematycznie jakość predykcji badamy porównując wariancje szumów N i E:

  • jeśli var(N) < var(E) to Y jest przyczynowy dla X.

Tak rozumianą relację przyczynowości między sygnałami nazywamy przyczynowością Grangera (lub przyczynowością w sensie Grangera, ang. Granger causality).

Zauważmy, że predykcja nie pogorszy się (pierwotne wyrażenie z poprzednimi wartościami X wciąż we wzorze jest), może albo pozostać bez zmian (jeśli Y nie ma nic wspólnego z X) albo się poprawić.

Funkcja DTF

W procesie dopasowywania współczynników modelu AR staramy się je tak dobrać, aby jak najlepiej opisać nasz sygnał, czyli aby jak najmniejsza część wariancji sygnału zostawała dla składowej losowej E(t). Możemy więc interpretować tę czynność jak w przypadku przewidywania wartości sygnału X na podstawie jego poprzednich wartości, a składową losową traktować jak błąd predykcji. Tak więc na podstawie badania współczynników dopasowanego modelu możemy wnioskować o zależnościach przyczynowych w naszych danych (wielokanałowych).

Załóżmy, że dla układu dwukanałowego dopasowaliśmy model AR rzędu 1 otrzymując

[math]\left( \begin{array}{c} X_1(t)\\X_2(t) \end{array}\right)= \left( \begin{array}{cccc} a_{11} & a_{12}\\ 0 & a_{22} \end{array}\right) \left( \begin{array}{c} X_1(t-1)\\X_2(t-1) \end{array}\right) +\left( \begin{array}{c} E_1(t)\\E_2(t) \end{array}\right) [/math]

Analiza współczynników dopasowanego modelu również dostarcza nam informacji o zależnościach między sygnałami w analizowanym zestawie. Obecność niezerowego współczynnika a12 oznacza, że poprzednie próbki sygnału X2 mają wpływ na bieżącą próbkę sygnału X1. Natomiast współczynnik a21 wynosi 0, widzimy więc, że poprzednie próbki sygnału X1 nie wpływają na bieżącą próbkę sygnału X2. Wnioskujemy więc, że w naszym układzie występuje wpływ X2X1, a nie ma wpływu X1X2. Widzimy również, że każdy kierunek możliwego wpływu jest niezależny od drugiego i mogą wystąpić każdy osobno lub oba równocześnie, zależnie od charakteru danych. Tego typu zależności nie bylibyśmy w stanie wykryć posługując się np. analizą korelacji.

Informacja zawarta we współczynnikach modelu A mówi nam o zależnościach w dziedzinie czasu. Aby móc coś powiedzieć o zależnościach w dziedzinie częstości wykorzystamy macierz przejścia H modelu AR wyrażonego w dziedzinie częstości:

[math] X(f)=A^{-1}(f)E(f)=H(f)E(f) [/math]

Macierz H produkuje z transformat Fouriera sygnałów szumowych (posiadających z definicji widmo płaskie, niezależne od częstości) widmo badanych sygnałów. Zawiera więc ona potrzebną nam informację o zależnościach miedzy naszymi sygnałami w dziedzinie częstości. Jest to macierz niesymetryczna, jej element Hij(f) mówi o intensywności wpływu sygnału Xj na sygnał Xi w częstości f.

Z elementów macierzy H zbudowana jest kierunkowa funkcja przejścia (ang. directed transfer function, DTF). Istnieje kilka wariantów tej funkcji, poniżej prezentowane są wzory opisujące wersję znormalizowaną i nieznormalizowaną.

Wersja znormalizowana

  • [math]\mathrm{DTF}_{ij}(f)=\mathrm{DTF}_{j\rightarrow i}(f)=\frac{\left| H_{ij}(f) \right|^2}{\sum_{m=1}^k{\left| H_{im}(f) \right|^2} }[/math]

Wersja nieznormalizowana

  • [math]\mathrm{NDTF}_{ij}(f)=\mathrm{NDTF}_{j\rightarrow i}(f)=\left| H_{ij}(f) \right|^2[/math]

Wersja nieznormalizowana to po prostu kwadrat modułu odpowiedniego elementu macierzy H. Jego wartość jest wprost proporcjonalna do intensywności związku między sygnałami; zależy ona od wariancji danych i może być dowolnie duża. Wersja znormalizowana opisuje stosunek wpływu z kanału o indeksie j do kanału o indeksie i (w częstości f) w stosunku do wszystkich wpływów do kanału o indeksie i (w tej częstości). W ten sposób wartość ta jest znormalizowana do przedziału [0, 1].

Dla obu wersji funkcji DTF wartości bliskie zeru mówią o braku związku między danymi sygnałami.

Zauważmy jeszcze, że model AR traktuje wszystkie kanały modelowanego układu na raz (jednocześnie). Jest to sytuacja zupełnie inna niż w przypadku badania każdej pary kanałów osobno.

Ćwiczenia

Wstęp do ćwiczeń

Do ćwiczeń w tym rozdziale używać będziemy zestawu danych, które służyły w poprzednim rozdziale do wyznaczania komponentów ICA (http://www.fuw.edu.pl/~jarekz/LabEEG/Dane_do_ICA_alfa.tar.gz). Aby dostosować je do naszych celów dokonamy na nich następujących operacji:

  • wybierzemy kanały EEG;
  • zastosujemy montaż do połączonych uszu (kanały A1 i A2);
  • zmniejszymy częstość próbkowania z 512 do 128 Hz (stosując antyaliasingowy filtr dolnoprzepustowy funkcją filtfilt i biorąc co czwartą próbkę filtrowanego sygnału);
  • przefiltrujemy sygnał (po zmianie) górnoprzepustowo z granicą odcięcia 1 Hz (stosując funkcję filtfilt).

Aby dopasować współczynniki modelu MVAR do posiadanych danych użyjemy funkcji licz_wsp_AR.m zawartej w pakiecie: MVAR_Matlab.zip — należy go ściągnąć i wypakować do katalogu roboczego. Funkcja ta zwraca współczynniki A modelu MVAR oraz macierz wariancji szumu V. Aby uzyskać macierz przejścia modelu H należy zaimplementować samemu procedurę transformacji współczynników do dziedziny częstości (patrz: tu).

Przez „szum” rozumieć będziemy wektor wartości losowych zwracanych przez funkcję randn.

Ćwiczenie 1

Z zestawu danych do obliczania ICA (poprzedni rozdział) wybierz jeden kanał EEG, zawierający wyraźną czynność alfa.

Przytnij wybrany odcinek do długości 2000 próbek. Podziel wybrany odcinek danych przez ich odchylenie standardowe.

Wygeneruj dwa zestawy danych:

  • Zestaw 1
   Kanał 1 to nasz wybrany kanał EEG
   Kanał 2 = (kanał 1 opóźniony o 1 próbkę)*0,6 + szum
  • Zestaw 2
  Kanał 1 to nasz wybrany kanał EEG
  Kanał 2 = szum

Dla obu zestawów danych sprawdź stosując metodę przyczynowości Grangera, który sygnał możemy uznać za przyczynowy dla drugiego sygnału. W tym celu w każdym zestawie dopasuj kolejno jednokanałowe modele AR oraz model dwukanałowy i porównaj otrzymane wariancje szumu.

Ćwiczenie 2

Dla obu zestawów danych z ćwiczenia 1 dopasuj dwukanałowe modele AR dla rzędów od 1 do 10. Sprawdź, ile wynoszą współczynniki a12 i a21 dopasowanych modeli AR w dziedzinie czasu dla j = 1...p. Dla każdego rzędu p oblicz dla kanałów 1 i 2 oraz 2 i 1 miarę DC [3]:

[math] \mathrm{DC}_{ij}^2=\sum_{m=1}^p{a_{ij}^2(m)} [/math]

Ćwiczenie 3

  • Wygeneruj dwa sygnały sinusoidalne o długości 1000 próbek każdy, o tej samej częstości 32 Hz i częstości próbkowania 128 Hz, ale różnych fazach początkowych.
  • Pierwszy sygnał powinien mieć fazę początkową równą 0, drugi sygnał sinusoidalny powinien mieć fazę początkową równą π/4.
  • Do drugiego z sygnałów dodaj małą (o amplitudzie ok 0,2 amplitudy sinusoidy) składową losową (czyli dodatkowy niezależny szum biały).
  • Z tak otrzymanych sygnałów utwórz jeden sygnał dwukanałowy (macierz o rozmiarze (2,1000)).

Ustal optymalny rząd modelu AR (tym razem dwukanałowego) i oblicz macierz gęstości widmowej mocy oraz koherencji między tymi sygnałami. Narysuj moduł i fazę koherencji C12 i C21.

Dla tego zestawu kanałów oblicz i narysuj normalizowaną i nienormalizowaną funkcję DTF.

Zmień fazę początkową drugiego sygnału. Jak zmienia się funkcja koherencji? Co dzieje się z funkcją DTF?

Ćwiczenie 4

Ćwiczenie to opisuje problem tzw. wspólnego źródła. W takiej sytuacji prawdziwa analiza wielokanałowa wykazuje swoją przewagę nad badaniem zależności parami (dla każdej pary kanałów osobno).

Wygeneruj układ trzech sygnałów w następujący sposób:

   jako pierwszego kanału użyj sygnału z ćwiczenia 1;
   sygnał_w_drugim_kanale(t) = 0,4 * sygnał_z_pierwszego_kanału(t−1) + 0,4 * szum1;
   sygnał_w_trzecim_kanale(t) = 0,3 * sygnał_z_pierwszego_kanału(t−2) + 0,4 * szum2.

Oblicz macierz koherencji zwyczajnych dla tego układu i na ich podstawie wyznacz zależności między kanałami. Powtórz to samo dla koherencji cząstkowych.

Oblicz dla tego zestawu danych funkcje DTF.

Oblicz funkcję DTF tylko dla kanałów 2 i 3 (model dwukanałowy). Czym różni się wynik od odpowiedniej funkcji dla układu wszystkich trzech kanałów?

Wyniki wszystkich obliczeń przedstaw na rysunkach.

Ćwiczenie 5

Oblicz funkcje DTF dla wszystkich kanałów EEG z przygotowanego zestawu danych do ICA (dla pełnej długości w czasie każdego kanału). Zaobserwuj przepływy w paśmie częstości alfa.

Ćwiczenie 6

Oblicz funkcje DTF dla sygnałów zebranych podczas eksperymentu z ruchu palcem (Pracownia EEG: https://brain.fuw.edu.pl/edu/index.php/Pracownia_EEG/ERDS#Cwiczenia). Do analizy wybierz kanały położone nad korą ruchową rąk lewej i prawej: C3, Cz, C4, F3, F4. Oblicz wyniki dla poszczególnych realizacji eksperymentu, a następnie je uśrednij po realizacjach. Porównaj przepływy dla okresu referencyjnego (−4..−2 s) oraz dla okresu po ruchu (+0,5..+2,5 s). Porównaj również wyniki dla normalizowanej i nienormalizowanej wersji funkcji DTF.

Podziel cały wybrany odcinek danych eksperymentalnych (−5..+5 s) na okna o długości 1 sekundy nakładające się na siebie z zakładką 0,5 s. Policz dla każdego okna nienormalizowane funkcje DTF (uśredniając po realizacjach). Narysuj otrzymane wyniki w postaci macierzy map czas-częstość obrazujących zmiany intensywności przepływów w czasie.

Literatura

  • Kamiński M, Blinowska KJ. A new method of the description of the information flow in brain structures. Biol Cybern 1991; 65:203–10.
  • Granger CWJ. Investigating causal relations in by econometric models and crossspectral methods. Econometrica 1969; 37:424–38.
  • Kamiński M, Ding M, Truccolo W, Bressler S. Evaluating causal relations in neural systems: Granger causality, directed transfer function and statistical assessment of significance. Biol Cybern 2001; 85:145–57.