Laboratorium EEG/AR 1

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).

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⁻¹ otrzymujemy

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

gdzie [math]H(\omega)[/math] jest macierzą przejścia. MVAR jest modelem typu "czarna skrzynka", gdzie na wejściu występują szumy, na wyjściu sygnały, a system jest opisany przez macierz przejścia. Zawiera on informacje o własnościach widmowych sygnałów i związkach między nimi.

Na podstawie macierzy [math]H(\omega)[/math] można obliczyć macierz gęstości widmowej zawierającą widma mocy dla pojedynczych kanałów jak również funkcje wzajemnej gęstości mocy pomiędzy kanałami. Stosując tego typu podejście, w którym wszystkie sygnały generowane przez pewien proces są rozpatrywane jednocześnie, można policzyć z macierzy spektralnej nie tylko koherencje zwykłe pomiędzy dwoma kanałami, ale również koherencje wielorakie opisujące związek danego kanału z pozostałymi i koherencje cząstkowe opisujące bezpośrednie związki między dwoma kanałami po usunięciu wpływu pozostałych kanałów. W przypadku gdy pewien kanał 1 będzie wpływał na kanały 2 i 3, obliczając koherencję zwykłą znajdziemy związek między 2 oraz 3, chociaż nie są one ze sobą bezpośrednio powiązane, natomiast koherencja cząstkowa nie wykaże związku między nimi.

Macierz [math]H(\omega)[/math] jest niesymetryczna, a jej wyrazy pozadiagonalne mają sens przyczynowości Grangera, co oznacza, że uwzględnienie wcześniejszej informacji zawartej w jednym z sygnałów zmniejsza błąd predykcji drugiego sygnału. Opierając się na tej własności zdefiniowano Kierunkową Funkcję Przejścia (DTF, directed transfer function) jako znormalizowany element pozadiagonalny [math]H(\omega)[/math]. DTF opisuje kierunek propagacji i skład widmowy rozchodzących się sygnałów.

Otrzymamy w ten sposób całościowy opis zmian wszystkich sygnałów jednocześnie. Co ciekawe, obliczona na tej podstawie funkcja charakteryzująca zależności między sygnałami [math]s_i[/math] (funkcja przejścia) nie jest symetryczna, w przeciwieństwie do np. korelacji. Dzięki temu może służyć wnioskowaniu nie tylko o sile zależności między poszczególnymi sygnałami składowymi, ale też o kierunku przepływu informacji między nimi. W przybliżeniu odpowiada to informacji, w którym z sygnałów struktury odpowiadające danej częstości pojawiają się wcześniej.

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.

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

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]

Ć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. Aby dostosować je do naszych celów dokonamy na nich następujących operacji:

• zastosujemy montaż do połączonych uszu (kanały A1 i A2);
• zmniejszymy częstość próbkowania z 512 do 128 Hz (stosując 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).

Ć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. 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 sprawdź współczynniki dopasowanych modeli AR w dziedzinie czasu.

Ć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 C₁₂ i C₂₁.

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

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

Ćwiczenie 4

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) + szum1;
   sygnał_w_trzecim_kanale(t) = 0,3 * sygnał_z_pierwszego_kanału(t−2) + 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.

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).

Polecenie

Zaimplementuj funkcję obliczającą koherencję dla pary kanałów. 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)