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

Z Brain-wiki
 
(Nie pokazano 79 wersji utworzonych przez 2 użytkowników)
Linia 1: Linia 1:
−
=Funkcja kowariancji i korelacji=
+
[[Laboratorium_EEG]]/MVAR
  
−
==Wstęp==
+
=Wielokanałowe modele AR=
  
−
W celu scharakteryzowania zależności wzajemnej dwóch sygnałów losowych, stosuje się funkcję kowariancji, zdefiniowaną w następujący sposób:
+
[[Model autoregresyjny (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):
  
−
<equation id="uid98">
 
 
<math>
 
<math>
−
\gamma _{xy} (\tau ) = \mathrm{cov}(x(t),y(t-\tau ))=\mathrm{E}[(x(t)-\mu _x)(y(t-\tau )-\mu _y)]
+
\let\hdots \dots
 +
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>
 
</math>
−
</equation>
 
  
−
gdzie:
+
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ść
  
−
<equation id="uid99">
 
 
<math>
 
<math>
−
\begin{array}{l}
+
A(f)X(f)=E(f)
−
\mu _x = \mathrm{E}[x(t)]\\
 
−
\mu _y = \mathrm{E}[y(t)]\\ \end{array}
 
 
</math>
 
</math>
−
</equation>
 
  
−
W przypadku sygnałów ciągłych estymację tę można zapisać w poniższy sposób:
+
gdzie, podobnie jak poprzednio, odpowiednie wielkości stają się wektorami ''k''-wierszowymi i macierzami rozmiaru ''k''&times;''k'':
  
−
<equation id="uid100">
 
 
<math>
 
<math>
−
\gamma _{xy} (\tau ) = \frac{1}{T}\int _0^{T}(x(t)-\mu_x)(y(t-\tau)-\mu_y)dt
+
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>
 
</math>
−
</equation>
 
  
−
natomiast dla sygnałów dyskretnych jako:
+
Mnożąc powyższe równanie lewostronnie przez macierz ''A''<sup>&minus;1</sup> otrzymujemy
  
−
<equation id="uid101">
 
 
<math>
 
<math>
−
\gamma _{xy}(k) = \frac{1}{N-1}\sum _{i=0}^{N-k}(x(i+k)-x_s)(y(i)-y_s)
+
X(f)=A^{-1}(f)E(f)=H(f)E(f)
 
</math>
 
</math>
−
</equation>
 
  
−
W odróżnieniu od funkcji autokowariancji, funkcja kowariancji nie musi mieć maksimum dla przesunięcia <math>\tau =0</math>. Ponadto posiada ona następującą cechę:
+
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'':
  
−
<equation id="uid102">
 
 
<math>
 
<math>
−
\gamma _{xy}(-\tau ) = \gamma _{yx}(\tau )
+
X(t)=\sum_{j=1}^\infty{A(j)X(t-j)}+E(t)
 
</math>
 
</math>
−
</equation>
 
  
−
Funkcję kowariancji można znormalizować:
+
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''
  
−
<equation id="uid103">
 
 
<math>
 
<math>
−
\rho (k) = \frac{\mathrm{E}[(x(t)-\mu _x)(y(t-\tau )-\mu _y)]}{\sqrt{\mathrm{E}[(x(t)-\mu _x)^2]\mathrm{E}[(y(t)-\mu _y)^2]}} = \frac{\gamma _{xy}}{\sigma_x\sigma_y}
+
X(t)=\sum_{j=1}^\infty{A(j)X(t-j)}+\sum_{j=1}^\infty{B(j)Y(t-j)}+N(t)
 
</math>
 
</math>
−
</equation>
 
−
Otrzymaną funkcję nazywamy funkcją korelacji.
 
−
Jednym z zastosowań funkcji korelacji jest wyznaczanie czasu przejścia sygnału przez dany układ liniowy. Funkcja korelacji pomiędzy sygnałem na wejściu układu i sygnałem na jego wyjściu osiągnie wartość maksymalną dla przesunięcia <math>\tau </math> równego czasowi, jaki potrzebował sygnał na pokonanie danego układu. Niestety, taka metoda wyznaczania opóźnienia obarczona jest pewną wadą &mdash; w przypadku gdy prędkość sygnału bądź jego droga zależą od częstości, wtedy na wykresie funkcji korelacji nie uzyskamy wyraźnego maksimum.
 
  
 +
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==
  
−
==Zadanie : Funkcja kowariancji i korelacji==
+
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).
−
Zaimplementuj funkcję obliczającą funkcję kowariancji dla różnych sygnałów ''x'' i ''y'' (równanie 13) skorzystaj przy tym z własności opisanej równaniem (14).  
 
−
Przykładowe wywołanie:
 
−
<source lang = python>
 
−
a = np.array([1,2,3])
 
−
b = np.array([-1,-2,-3])
 
  
−
print koreluj(a,b,2)
+
Załóżmy, że dla układu dwukanałowego dopasowaliśmy model AR rzędu 1 otrzymując
−
</source>
 
−
powinno dać w wyniku:
 
−
[ 0.5 0.  -1.  0.  0.5]
 
  
−
<!--
+
<math>\left( \begin{array}{c} X_1(t)\\X_2(t) \end{array}\right)=
−
{{hidden begin|title=Przykładowe rozwiązanie:}}
+
\left( \begin{array}{cccc} a_{11} & a_{12}\\
−
<source lang = python>
+
0 & a_{22} \end{array}\right)
−
import numpy as np
+
\left( \begin{array}{c} X_1(t-1)\\X_2(t-1) \end{array}\right)
−
import pylab as py
+
+\left( \begin{array}{c} E_1(t)\\E_2(t) \end{array}\right)
 +
</math>
  
−
def koreluj(x,y,max_tau):
+
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 ''a''<sub>12</sub> oznacza, że poprzednie próbki sygnału ''X''<sub>2</sub> mają wpływ na bieżącą próbkę sygnału ''X''<sub>1</sub>. Natomiast
−
    x = x - np.mean(x)
+
współczynnik ''a''<sub>21</sub> wynosi 0, widzimy więc, że poprzednie próbki sygnału ''X''<sub>1</sub> nie wpływają na bieżącą próbkę sygnału ''X''<sub>2</sub>. Wnioskujemy więc, że w naszym układzie występuje wpływ ''X''<sub>2</sub>&rarr;''X''<sub>1</sub>, a nie ma wpływu ''X''<sub>1</sub>&rarr;''X''<sub>2</sub>. 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.
−
    y = y - np.mean(y)
 
−
    cor = np.zeros(2*max_tau+1)
 
−
    cor[max_tau] = np.sum(x[:]*y[:])
 
−
    for i in range(1,max_tau+1):
 
−
        cor[max_tau+i] = np.sum(x[i:]*y[:-i])
 
−
        cor[max_tau-i] = np.sum(y[i:]*x[:-i])
 
−
    N= len(x)
 
−
    cor = cor /(N-1)
 
−
    return cor
 
  
−
a = np.array([1,2,3])
+
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:
−
b = np.array([-1,-2,-3])
 
−
for i in range(3):
 
−
    print koreluj(a,b,i)
 
−
</source>
 
−
{{hidden end}} -->
 
−
Z danych zarejestrowanych w trakcie czuwania z zamkniętymi oczami wybierz sygnały z następujących kanałów: Fp1, P3, Pz, P4, Fp2, O1, O2.
 
  
−
<ol>
+
<math>
 +
X(f)=A^{-1}(f)E(f)=H(f)E(f)
 +
</math>
  
−
<li>
+
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''.  
−
Dla każdego kanału oblicz funkcję autokorelacji, zaś  dla każdej pary kanałów oblicz funkcję korelacji wzajemnej. Wyniki zaprezentuj w formie kwadratowej macierzy wykresów (za pomocą funkcji subplot, tak jak na przykładowym rys. (rys. <xr id="uid9"> %i</xr>)). Na przekątnej macierzy narysuj funkcję autokorelacji odpowiednich kanałów, poza przekątną &mdash; funkcję korelacji wzajemnej. Wskaż kanały, które są najbardziej skorelowane ze sobą. Czy możliwe jest wyznaczenie opóźnienia sygnału pomiędzy tymi kanałami?
 
  
 +
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ą.
  
−
<li>
+
Wersja znormalizowana
−
Powtórz punkt 1, tym razem jednak funkcję autokorelacji i korelacji wzajemnej oblicz na sygnałach przefiltrowanych filtrem wąskopasmowym w paśmie alfa charakterystycznym dla badanej osoby. ([[%C4%86wiczenia_7#Funkcje_do_projektowania_filtr.C3.B3w_IIR_dost.C4.99pne_w_module_scipy.signal|przypomnienie konstrukcji filtrów]])
+
* <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>
  
−
<li>
+
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].
−
Oszacuj istotność statystyczną zależności między parami kanałów. Twoją hipotezą zerową jest brak istotnej korelacji pomiędzy sygnałami zarejestrowanymi przez dwie różne elektrody EEG. Hipoteza alternatywna to występowanie zależności pomiędzy tymi sygnałami. Podanie estymatorów wariancji funkcji korelacji jest bardzo trudne, dlatego jednym ze sposobów oszacowania progu powyżej którego wartość funkcji korelacji można byłoby uznać za istotną statystycznie, jest zastosowanie metody ''bootstrap''. Teoretycznie, funkcja korelacji policzona dla dwóch rzeczywistych, nieskorelowanych sygnałów, powinna wynosić 0 dla każdego przesunięcia <math>\tau</math>. Tak jest jednak w przypadku sygnałów nieskończonych; w analizie sygnałów takowych nie spotkamy.
 
  
−
Dokonując losowej zamiany kolejności próbek, możemy doprowadzić do wytworzenia sygnałów zależnych losowo, które jednak ze względu na skończony czas trwania, dadzą niezerową funkcję korelacji. Poziom losowych fluktuacji tej funkcji oszacujemy wykonując następujące kroki:
+
Dla obu wersji funkcji DTF wartości bliskie zeru mówią o braku związku między danymi sygnałami.
−
<ol type="A">
 
−
<li> Losowa zamiana kolejności próbek w analizowanych sygnałach. Jeżeli pomiędzy dwoma sygnałami istnieją jakieś zależności, losowa zamiana próbek doprowadzi do zniszczenia tych związków. W ten sposób uzyskujemy sygnały, które teoretycznie są nieskorelowane.
 
−
<li> Obliczenie funkcji  korelacji wzajemnej dla sygnałów policzonych w punkcie A.
 
−
<li> Powtórzenie kroków A i B wiele (np. 1000) razy.
 
−
<li> Oszacowanie 95 % przedziału ufności dla wartości średniej funkcji korelacji wzajemnej dla danego przesunięcia <math>\tau</math> korzystając z otrzymanego w kroku C empirycznego rozkładu wartości tych funkcji dla sygnałów niezależnych. 
 
−
<li> Powtórzenie kroków A-D dla kolejnych przesunięć <math>\tau</math>.
 
−
<li> Sprawdzenie, dla których przesunięć <math>\tau </math> funkcje autokorelacji i korelacji obliczone dla oryginalnych sygnałów uzyskały wartości wyższe niż wartości progowe oszacowane dla sygnałów o losowych zależnościach.
 
−
</ol>
 
  
−
Procedura opisana powyżej ma jednak pewną wadę. Staramy się w niej oszacować poziom przypadkowych korelacji pomiędzy dwoma sygnałami dla kolejnych przesunięć <math>\tau </math>, co jest niczym innym jak wielokrotnym powtórzeniem pewnego testu. Obserwowanie korelacji dla wielu par kanałów równocześnie również prowadzi do zwiększenia szansy na zaobserwowanie ekstremalnie dużych fluktuacji.
+
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.
−
Występuje tu zatem ''problem wielokrotnych porównań''.
 
−
Przypominamy, iż może to doprowadzić do przypadkowego uznania wyników jako &bdquo;istotnych&rdquo; statystycznie. Np. jeśli pojedynczy test wykonujemy na poziomie istotności 5% to dopuszczamy odrzucenie w 1 przypadku na 20 hipotezy zerowej pomimo, iż jest ona prawdziwa. Z drugiej jednak strony, jeśli powtórzymy wykonywany test 20 razy, to oczekujemy uzyskania 1 przypadku, w którym poziom <math>p</math> będzie mniejszy od 5% co jest przesłanką za odrzuceniem hipotezy zerowej.  
 
  
−
W przypadku wykonywania serii testów należałoby więc zastosować odpowiednie poprawki, np. [http://www.bmj.com/content/310/6973/170.full korektę Bonferroniego] czy [http://en.wikipedia.org/wiki/False_discovery_rate false discovery rate (FDR)]. Innym rozwiązaniem w analizowanym przez nas problemie jest zastosowanie tzw. statystyk wartości ekstremalnych, które prowadzą do następujących zmian w procedurze (nie działa dla funkcji autokorelacji ze względu na jej normalizację do 1 dla zerowego przesunięcia):
+
=Ćwiczenia=
  
−
<ol type="A">
+
==Wstęp do ćwiczeń==
−
<li> Losowa zmiana kolejności próbek w analizowanych sygnałach (we wszystkich analizowanych kanałach). Jeżeli pomiędzy dwoma sygnałami istnieją jakieś zależności, losowa zamiana próbek doprowadzi do zniszczenia tych związków. W ten sposób uzyskujemy sygnały, które teoretycznie są nieskorelowane.
 
−
<li> Obliczenie funkcji korelacji dla sygnałów otrzymanych w punkcie A.
 
−
<li>    Zapamiętanie maksymalnej wartości bezwzględnej funkcji korelacji z punktu B (maksimum bierzemy po wszystkich przesunięciach i po wszystkich parach kanałów).
 
−
<li> Powtórzenie kroków A-C 1000 razy. Uzyskamy w ten sposób rozkład maksymalnych wartości funkcji korelacji możliwych do zaobserwowania dla sygnałów niezależnych.
 
−
<li>    Wyznaczenie 95 centyla rozkładu wartości maksymalnych.
 
−
<li> Nałożenie na rysunki funkcji korelacji uzyskane w Zadaniu 2 poziomych linii symbolizujących poziom zależności dwóch sygnałów o losowych zależnościach i sprawdzenie, dla których przesunięć <math>\tau </math> wartości funkcji korelacji przekraczają estymowane progi istotności statystycznej.
 
−
</ol>
 
  
−
</ol>
+
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).
  
−
[[Plik:Korelacje_wzajemne.png|700px|center|thumb|<figure id="uid9" />Przykład wyniku analizy korelacji wzajemnych dla sygnału niefiltrowanego z naniesionymi granicami możliwych fluktuacji.]]
+
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]).
  
−
==Wzajemna gęstość widmowa sygnałów i koherencja==
+
Przez &bdquo;szum&rdquo; rozumieć będziemy wektor wartości losowych zwracanych przez funkcję <tt>randn</tt>.
  
−
====Zadanie 4: Wzajemna gęstość widmowa sygnałów i koherencja====
+
==Ćwiczenie 1==
−
=====Wstęp=====
 
−
Podobnie jak w przypadku twierdzenia Chinczyna dla pojedynczego sygnału, możliwe jest policzenie transformaty Fouriera funkcji kowariancji. Uzyskana w ten sposób wielkość nazywa się funkcją wzajemnej gęstości mocy widmowej sygnału:
 
  
−
<equation id="uid122">
+
Z zestawu danych do obliczania ICA (poprzedni rozdział) wybierz jeden kanał EEG, zawierający wyraźną czynność alfa.
−
<math>
+
<br><br>
−
S_{xy}(f) = \int _{-\infty }^{\infty }\gamma_{xy}(\tau )e^{-2\pi i f \tau}d\tau </math>
+
Przytnij wybrany odcinek do długości 2000 próbek. Podziel wybrany odcinek danych przez ich odchylenie standardowe.
−
</equation>
+
<br><br>UWAGA:
 +
* Sygnał po wycięciu z oryginalnego nagrania musimy jeszcze przefiltrować filtrem górnoprzepustowym (o granicy odcięcia np. 1 Hz) aby z sygnału usunąć składową stałą i inne efekty „płynięcia” o niskich częstościach. Możemy zastosować np. filtr Butterwortha, przyzwyczajajmy się również do filtrowania funkcją <tt>filtfilt</tt>.
 +
* „szum” w opisie konstrukcji drugiego kanału wytwarzamy funkcją <tt>randn</tt>.
 +
* Po skonstruowaniu kanału 2 również dzielimy go przez jego wariancję aby oba sygnały miały jednakową wariancję (jednostkową).
 +
 
 +
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*0.15
 +
 
 +
* Zestaw 2
 +
  Kanał 1 to nasz wybrany kanał EEG
 +
  Kanał 2 = szum
  
−
W celu dalszego omówienia własności funkcji wzajemnej mocy widmowej sygnałów funkcję tę zapiszemy w postaci:
+
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.
  
−
<equation id="uid123">
+
==Ćwiczenie 2==
−
<math>
 
−
\begin{array}{l}
 
−
S_{xy}(f) = |S_{xy}(f)|e^{i\phi _{xy}(f)}\\
 
−
\\
 
−
\phi _{xy} = \arg(S_{xy})
 
−
\end{array} </math>
 
−
</equation>
 
−
<!-- \mathrm{arc\,tg}\left[\frac{\mathrm{Im}(S_{xy}(f))}{\mathrm{Re}(S_{xy}(f))}\right]-->
 
  
−
Wartość bezwzględna funkcji wzajemnej gęstości mocy widmowej osiąga największą wartość dla '''częstości''', w których sygnały <math>x(t)</math> i <math>y(t)</math> są ze sobą skorelowane. Funkcja wzajemnej mocy widmowej sygnałów pozbawiona jest zatem wady, która charakteryzowała funkcję korelacji, to jest problemu z wyznaczeniem czasu transmisji sygnału, w przypadku gdy czas ten zależał od częstości. Przy pomocy funkcji wzajemnej mocy widmowej, czas ten można oszacować przy pomocy fazy tej funkcji &mdash; <math>\phi _{xy}(f)</math>. Jeśli funkcja wzajemnej mocy widmowej została wyznaczona pomiędzy sygnałami na wejściu i wyjściu układu liniowego, to faza ta reprezentuje przesunięcie fazowe sygnału przy przejściu przez układ. Czas tego przejścia można oszacować za pomocą następującej wyrażenia:
+
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 ''m'' = 1...''p''.
 +
Dla każdego rzędu ''p'' oblicz dla kanałów 1 i 2 oraz 2 i 1 miarę DC [3]:
  
−
<equation id="uid124">
 
 
<math>
 
<math>
−
\tau = \frac{\phi _{xy}(f)}{2\pi f}
+
\mathrm{DC}_{ij}^2=\sum_{m=1}^p{a_{ij}^2(m)}
 
</math>
 
</math>
−
</equation>
 
  
−
Podobnie jak w przypadku funkcji autokorelacji i korelacji wzajemnej, funkcję wzajemnej gęstości mocy widmowej można znormalizować:
+
==Ć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ą &pi;/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 <tt>(2,1000)</tt>).
 +
 
 +
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ą 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:
  
−
<equation id="uid125">
+
    jako pierwszego kanału użyj sygnału z ćwiczenia 1;
−
<math>
+
    sygnał_w_drugim_kanale(t) = 0,4 * sygnał_z_pierwszego_kanału(t−1) + 0,4 * szum1;
−
C_{xy}(f) = \frac{S_{xy}(f)}{\sqrt{S_x(f)S_y(f)}}
+
    sygnał_w_trzecim_kanale(t) = 0,3 * sygnał_z_pierwszego_kanału(t−2) + 0,4 * szum2.
−
</math>
+
 
−
</equation>
+
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==
  
−
Znormalizowaną postać funkcji wzajemnej gęstości mocy widmowej nazywamy funkcją ''koherencji''.  
+
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).  
−
Koherencja jest wielkością zespoloną. Faza koherencji odzwierciedla różnicę faz pomiędzy dwoma sygnałami. Moduł koherencji reprezentuje stopień synchronizacji sygnałów i zawiera się w przedziale od 0.0 do 1.0. Moduł tej funkcji zawiera się w przedziale od 0 do 1. Wartości 0 odpowiada brak synchronizacji pomiędzy sygnałami, zaś wartości 1 pełna synchronizacja dwóch przebiegów czasowych. Należy również zwrócić uwagę na nazewnictwo - często sam moduł koherencji określany jest jako koherencja, w literaturze anglojęzycznej moduł koherencji posiada jednak odrębną nazwę: Magnitude Square Coherence (MSC). Istotny jest również sposób estymacji modułu koherencji, który wyprowadzono w następnym rozdziale, zaś sam estymator reprezentuje wzór (36).
+
Porównaj również wyniki dla normalizowanej i nienormalizowanej wersji funkcji DTF.  
  
−
=====Kilka słów o koherencji=====
+
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.
−
Wzór (<xr id="uid125"/>), definiujący ilościową miarę koherencji, nie uwzględnia stochastycznego charakteru sygnałów. Łatwo zauważyć, że bezpośrednie zastosowanie tego wzoru do obliczenia koherencji dwóch sygnałów o tej samej częstości i różniących się jedynie amplitudą oraz fazą, zawsze da wynik równy 1. Prześledźmy to na następującym przykładzie. <br>
 
−
Dane są dwa sygnały harmoniczne <math>x(t) = A\cos(\Omega t + \phi_x)</math> oraz <math>y(t) = B\cos(\Omega t + \phi_y)</math>.
 
−
Widmo tych sygnałów, wyrażone za pomocą transformaty Fouriera, będzie miało następującą postać: <br>
 
−
<br>
 
−
<math>X(f)=Ae^{-j\phi_x}</math> <br>
 
−
<br>
 
−
<math>Y(f)=Be^{-j\phi_y}</math>, <br>
 
−
<br>
 
−
zaś ich widmo wzajemne: <br>
 
−
<br>
 
−
<math>X(f)\cdot Y^*(f) = A\cdot Be^{-j(\phi_x - \phi_y)}</math>, <br>
 
−
<br>
 
−
gdzie: <math>j=\sqrt{-1}</math>, a * oznacza sprzężenie liczby zespolonej. <br>
 
−
Podstawienie wyrażeń na widmo sygnałów <math>x(t)</math>, <math>y(t)</math> oraz ich widmo wzajemne do wzoru <xr id="id9"/> da koherencję <math>K_{xy}(f) = 1</math> niezależnie od amplitudy sygnałów <math>A</math> i <math>B</math> oraz ich faz <math>\phi_x</math> i <math>\phi_y</math>.
 
−
<br>
 
−
<br>
 
−
W praktyce rzadko jednak mamy do czynienia z sygnałami harmonicznymi. Zwykle mierzone przez nas wielkości mają stochastyczny charakter bądź też ich pomiar jest zaburzany przez różne czynniki.
 
−
Rozważmy teraz najprostszy model pomiaru sygnału, w którym uwzględniono wpływ zakłóceń w postaci białego szumu. Na wejście układu LTI o funkcji impulsowej opisanej wyrażeniem <math>h(t)</math> podamy sygnał <math>x(t)</math> i widmie danym funkcją <math>X(f)</math>. Układ LTI przetworzy sygnał wejściowy na przebieg <math>y(t)</math> o widmie <math>Y(f)</math>. Z uwagi na zaburzenia <math>n(t)</math> o widmie <math>N(f)</math> towarzyszące pomiarowi aparatura nie zarejestruje sygnał <math>y(t)</math> lecz <math>z(t) = y(t) + n(t)</math>. Opisane zależności możemy opisać za pomocą poniższych wzorów:<br>
 
−
<br>
 
−
<math>y(t) = h(t)*x(t)</math> <br>
 
−
<br>
 
−
<math>z(t) = y(t) + n(t)</math> <br>
 
−
<br>
 
−
gdzie: <math>*</math> - operacja splotu. <br>
 
−
Dokonując transformacji powyższych wzorów do dziedziny częstości dostajemy:
 
−
<br>
 
−
<math>Y(f) = H(f)X(f)</math> <br>
 
−
<br>
 
−
<math>Z(f) = Y(f) + N(f)</math> <br>
 
−
<br>
 
−
gdzie: <math>H(f) = \textrm{FFT}\left\{h(t)\right\}</math>. <br>
 
−
<br>
 
−
Wzory te można zapisać w postaci jednej zależności:<br>
 
−
<equation id="LTI_1">
 
−
<math>Z(f) = H(f)X(f) + N(f)</math>
 
−
</equation>
 
−
Załóżmy teraz, że w celu redukcji składowej losowej <math>n(t)</math> wielokrotnie powtarzamy w tych samych warunkach pomiar sygnału <math>z(t)</math>. Za każdym razem na wejściu układu LTI występuje ten sam sygnał <math>x(t)</math>. Układ LTI również przetwarza sygnał wejściowy w ten sam sposób, jednak z uwagi na stochastyczny charakter zakłóceń, otrzymujemy kolejne różniące się do siebie przebiegi <math>z_i(t)</math>. Niech liczbę powtórzeń pomiaru wynosi <math>K</math>. Możemy napisać <math>K</math> równań opisujących relację pomiędzy sygnałem wejściowym, wyjściowym i mierzonym:
 
−
<equation id="LTI_2">
 
−
<math>
 
−
\begin{array}{l}
 
−
Z_1(f) = H(f)X(f) + N_1(f) \\
 
−
\\
 
−
Z_2(f) = H(f)X(f) + N_2(f) \\
 
−
\\
 
−
\vdots \\
 
−
\\
 
−
Z_K(f) = H(f)X(f) + N_K(f) \\
 
−
\end{array}
 
−
</math>
 
−
</equation>
 
−
Przemnóżmy teraz równania (<xr id="LTI_2"/>) obustronnie przez sprzężone widmo sygnału rejestrowanego <math>Z(f)</math>. Dla uproszczenia zapisu operacji dokonamy na jednym, dowolnie wybranym <math>i</math>-tym równaniu:
 
−
<equation id="LTI_3">
 
−
<math>Z_i(f)Z_i^*(f) = \left\{H(f)X(f) + N_i(f)\right\}\cdot Z_i^*(f)</math>
 
−
</equation>
 
−
Na równaniu (<xr id="LTI_3"/>) dokonamy kolejno następujących przekształceń:
 
−
<equation id="equ_1">
 
−
<math>|Z_i(f)|^2 = \left\{H(f)X(f) + N_i(f)\right\}\cdot\left\{H^*(f)X^*(f) + N_i^*(f)\right\}</math>
 
−
</equation>
 
−
<br>
 
−
<equation id="LTI_4">
 
−
<math>|Z_i(f)|^2 = |H(f)|^2|X(f)|^2 + |N_i(f)|^2 + H(f)X(f)N_i^*(f) + N_i(f)H^*(f)X^*(f)</math>
 
−
</equation>
 
−
Dokonajmy teraz uśredniania (<xr id="LTI_4"/>) po kolejnych powtórzeniach pomiaru.
 
−
<equation id="LTI_5">
 
−
<math>\left\langle|Z_i(f)|^2\right\rangle= \left\langle|H(f)|^2|X(f)|^2\right\rangle + \left\langle|N_i(f)|^2\right\rangle + \left\langle H(f)X(f)N_i^*(f)\right\rangle + \left\langle N_i(f)H^*(f)X^*(f)\right\rangle</math>
 
−
</equation>
 
−
Zakładamy, że szum <math>N(f)</math> jest nieskorelowany z sygnałem wejściowym, w związku z czym w wyniku uśredniania dwa ostatnie składniki równania (<xr id="LTI_5"/>) zostaną zredukowane: <math>\left\langle H(f)X(f)N_i^*(f)\right\rangle \approx 0 </math>, <math>\left\langle N_i(f)H^*(f)X^*(f)\right\rangle \approx 0 </math>. Założyliśmy również za każdym razem na wejściu układu liniowego pojawia się ten sam sygnał <math>x(t)</math>, sam układ zaś nie zmienia swoich właściwości, w zwiazku z czym: <math>\left\langle|H(f)|^2|X(f)|^2\right\rangle = |H(f)|^2|X(f)|^2 </math>. Ostatecznie uzyskaliśmy następującą zależność:
 
−
<equation id="LTI_6">
 
−
<math>\left\langle|Z_i(f)|^2\right\rangle= |H(f)|^2|X(f)|^2 + \left\langle|N_i(f)|^2\right\rangle</math>
 
−
</equation>
 
−
Dokonajmy kolejnego przekształcenia równania (<xr id="LTI_2"/>). tym razem przemnożymy obustronnie każde równanie przez sprzężone widmo sygnału wejściowego. W celu uproszczenia zapisu, operację tę wykonamy tylko na jednym dowolnie wybranym <math>i</math>-tym równaniu:
 
−
<equation id="LTI_7">
 
−
<math>Z_i(f)X^*(f) = \left\{H(f)X(f) + N_i(f)\right\}\cdot X^*(f)</math>
 
−
</equation>,
 
−
gdzie: <math>Z_i(f)X^*(f)</math> - to widmo wzajemne sygnałów <math>x(t)</math> i <math>y(t)</math>.
 
−
Proste przekształcenie równania (<xr id="LTI_7"/>) prowadzi do następującego wyrażenia:
 
−
<equation id="LTI_8">
 
−
<math>Z_i(f)X^*(f) = H(f)|X(f)|^2 + N_i(f)X^*(f)</math>
 
−
</equation>
 
−
Uśrednimy teraz równanie (<xr id="LTI_8"/>) po kolejnych realizacjach pomiaru oraz obliczmy moduł uzyskanego wyniku:
 
−
<equation id="LTI_9">
 
−
<math>|\left\langle Z_i(f)X^*(f)\right\rangle| = |H(f)||X(f)|^2 + |\left\langle N_i(f)X^*(f)\right\rangle|</math>
 
−
</equation>
 
−
Brak korelacji pomiędzy szumem <math>n(t)</math> a sygnałem wejściowym <math>x(t)</math> powoduje, że w wyniku uśredniania zostaje zredukowany drugi składnik równania (<xr id="LTI_9"/>): <math>\left\langle N_i(f)X^*(f)\right\rangle \approx 0</math>. Ostatecznie uzyskujemy następującą zależność:
 
−
<equation id="LTI_10">
 
−
<math>|\left\langle Z_i(f)X^*(f)\right\rangle| = |H(f)||X(f)|^2</math>
 
−
</equation>
 
−
która wraz z równaniem  (<xr id="LTI_6"/>) tworzy układ równań opisujących relacje pomiędzy widmami i widmami mocy sygnałów występujących w naszym modelu:
 
−
<equation id="LTI_11">
 
−
<math>
 
−
\left\langle Z_i(f)X^*(f)\right\rangle = |H(f)| |X(f)|^2
 
−
</math>
 
−
</equation>
 
−
::<math>
 
−
\left\langle|Z_i(f)|^2\right\rangle= |H(f)|^2 |X(f)|^2 + \left\langle|N_i(f)|^2\right\rangle
 
−
</math>
 
−
Z pierwszej zależności równania (<xr id="LTI_11"/>) wyznaczmy funkcję przejścia <math>|H(f)|</math>:
 
−
<br>
 
−
<br>
 
−
<math>|H(f)| = \frac{|\left\langle Z_i(f)X^*(f)\right\rangle|}{|X(f)|^2}</math>
 
−
<br>
 
−
<br>
 
−
i podstawy do drugiego równania układu (<xr id="LTI_11"/>). Otrzymujemy:
 
−
<equation id="LTI_12">
 
−
<math>
 
−
\left\langle|Z_i(f)|^2\right\rangle = \left[\frac{|\left\langle Z_i(f)X^*(f)\right\rangle|}{|X(f)|^2}\right]^2 |X(f)|^2 + \left\langle|N_i(f)|^2\right\rangle
 
−
</math>
 
−
</equation>
 
−
Równanie (<xr id="LTI_12"/>) możemy przekształcić do postaci:
 
−
<equation id="LTI_13">
 
−
<math>
 
−
\left\langle|N_i(f)|^2\right\rangle = \left\langle|Z_i(f)|^2\right\rangle - \frac{|\left\langle Z_i(f)X^*(f)\right\rangle|^2}{|X(f)|^2}
 
−
</math>
 
−
</equation>
 
−
a następnie do zależności:
 
−
<equation id="LTI_14">
 
−
<math>
 
−
\left\langle|N_i(f)|^2\right\rangle = \left\langle|Z_i(f)|^2\right\rangle\left[1 - \frac{|\left\langle Z_i(f)X^*(f)\right\rangle|^2}{|X(f)|^2\left\langle|Z_i(f)|^2\right\rangle}\right]
 
−
</math>
 
−
</equation>
 
−
Wyrażenie: <br>
 
−
<equation id="LTI_15">
 
−
<math>
 
−
\mathrm{MSC}_{xz}(f) = \frac{|\left\langle Z_i(f)X^*(f)\right\rangle|^2}{|X(f)|^2\left\langle|Z_i(f)|^2\right\rangle}
 
−
</math>
 
−
</equation>
 
−
nazywana jest '''M'''agnitude '''S'''quare '''C'''oherence pomiędzy sygnałami <math>x(t)</math> i <math>z(t)</math>. W przypadku, gdy wielkość ta jest równa 1 sygnały <math>x(t)</math> i <math>z(t)</math> są w pełni  zsynchronizowane. Wielkość tę uzyskaliśmy dla sygnału na wejściu układu LTI oraz sygnału mierzonego na wyjściu. Funkcję MSC można jednak stosować do dowolnych dwóch sygnałów stochastycznych <math>x(t)</math> i <math>y(t)</math> przy założeniu, że istnieją pomiędzy nimi liniowe zależności:
 
−
<equation id="LTI_16">
 
−
<math>
 
−
\mathrm{MSC}_{xy}(f) = \frac{|\left\langle X_i(f)Y_i^*(f)\right\rangle|^2}{\left\langle|X_i(f)|^2\right\rangle\left\langle|Y_i(f)|^2\right\rangle}
 
−
</math>
 
−
</equation>
 
−
gdzie:
 
−
<math>< ></math> - oznacza wartość średnia,
 
−
<math>X_i(f), Y_i(f) </math> to zespolone widma (policzone np. za pomocą Transformaty Fouriera), wyznaczone odpowiednio dla sygnałów X oraz Y w "i-tej" realizacji eksperymentu lub w "i-tym" oknie czasowym, na który te sygnały zostały podzielone. Wzór (36) reprezentuje estymator wartości bezwzględnej koherencji. Opierając się na podobnym co wyżej rozumowaniu, można wyprowadzić estymator funkcji koherencji, o następującej postaci:
 
−
<equation id="LTI_17">
 
−
<math>
 
−
\mathrm{C}_{xy}(f) = \frac{\left\langle X_i(f)Y_i^*(f)\right\rangle}{(\left\langle|X_i(f)|^2\right\rangle\left\langle|Y_i(f)|^2\right\rangle)^\frac{1}{2}}
 
−
</math>
 
−
</equation>
 
−
Faza koherencji umożliwia nam estymację przesunięcia fazowego pomiędzy sygnałami X i Y, zaś moduł podniesiony do kwadratu funkcji C to MSC.
 
  
−
=====Polecenie 2=====
+
=Literatura=
−
Zaimplementuj funkcję obliczającą koherencję dla pary kanałów.  
+
* Kamiński M, Blinowska KJ. A new method of the description of the information flow in brain structures. Biol Cybern 1991; 65:203–10.
−
<!--Niech argumentami tej funkcji będą dwa wektory zawierające sygnały, zakres częstości, częstość próbkowania. -->
+
* Granger CWJ. Investigating causal relations in by econometric models and crossspectral methods. Econometrica 1969; 37:424–38.
−
Oblicz i narysuj funkcję koherencji dla kolejnych par kanałów (tych samych co w zadaniu 3). Wyniki zaprezentuj w postaci kwadratowej macierzy rysunków. Ponieważ koherencja jest funkcją zespoloną, dobrze jest zaprezentować osobno jej wartość i fazę. Uzyskane wartości bezwzględne koherencje narysuj nad przekątną tej macierzy, a fazę pod przekątną. W celu obliczenia modułu koherencji i jej fazy wykorzystaj wzór 36 (wygenerowane sygnały należy podzielić na pewną liczbę odcinków)
+
* 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ń 08:51, 25 maj 2022

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] \let\hdots \dots 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 X2→X1, a nie ma wpływu X1→X2. 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.

UWAGA:

  • Sygnał po wycięciu z oryginalnego nagrania musimy jeszcze przefiltrować filtrem górnoprzepustowym (o granicy odcięcia np. 1 Hz) aby z sygnału usunąć składową stałą i inne efekty „płynięcia” o niskich częstościach. Możemy zastosować np. filtr Butterwortha, przyzwyczajajmy się również do filtrowania funkcją filtfilt.
  • „szum” w opisie konstrukcji drugiego kanału wytwarzamy funkcją randn.
  • Po skonstruowaniu kanału 2 również dzielimy go przez jego wariancję aby oba sygnały miały jednakową wariancję (jednostkową).

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*0.15
  • 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 m = 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.
Źródło: „http://brain.fuw.edu.pl/edu/index.php?title=Laboratorium_EEG/AR_1&oldid=8852”