Laboratorium EEG/AR 1

Wielokanałowe modele AR

Model AR opisuje wartość sygnału w chwili [math]t[/math] jako kombinację liniową jego wartości w chwilach poprzednich (oraz szumu). W przypadku wielowymiarowym możemy włączyć do tego opisu wartości wszystkich sygnałów [math]s_i[/math], czyli wektora [math]\vec{s}(t)[/math]. Wielozmienny model AR (MVAR, multivariate autoregressive) można wówczas opisać wzorem:

[math] \vec{s}(t)=\sum_{i=1}^p A(i) \vec{s}(t-i) + \vec{\epsilon}(t) , [/math]

gdzie [math]\vec{\epsilon}(t)[/math] będzie wektorem szumów, zaś [math]A(i)[/math] będą macierzami współczynników modelu. Przechodząc do przestrzeni częstości otrzymamy:

[math] \vec{s}(\omega)=A^{-1}(\omega)\vec{\epsilon}(\omega)=H(\omega)\vec{\epsilon}(\omega), [/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ść

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;
• 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

• 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 3

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 4

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)