Laboratorium EEG/AR 1

Z Brain-wiki

Funkcja kowariancji i korelacji

Wstęp

W celu scharakteryzowania zależności wzajemnej dwóch sygnałów losowych, stosuje się funkcję kowariancji, zdefiniowaną w następujący sposób:

[math] \gamma _{xy} (\tau ) = \mathrm{cov}(x(t),y(t-\tau ))=\mathrm{E}[(x(t)-\mu _x)(y(t-\tau )-\mu _y)] [/math]

gdzie:

[math] \begin{array}{l} \mu _x = \mathrm{E}[x(t)]\\ \mu _y = \mathrm{E}[y(t)]\\ \end{array} [/math]

W przypadku sygnałów ciągłych estymację tę można zapisać w poniższy sposób:

[math] \gamma _{xy} (\tau ) = \frac{1}{T}\int _0^{T}(x(t)-\mu_x)(y(t-\tau)-\mu_y)dt [/math]

natomiast dla sygnałów dyskretnych jako:

[math] \gamma _{xy}(k) = \frac{1}{N-1}\sum _{i=0}^{N-k}(x(i+k)-x_s)(y(i)-y_s) [/math]

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ę:

[math] \gamma _{xy}(-\tau ) = \gamma _{yx}(\tau ) [/math]

Funkcję kowariancji można znormalizować:

[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} [/math]

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ą — 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.


Zadanie : Funkcja kowariancji i korelacji

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:

a = np.array([1,2,3])
b = np.array([-1,-2,-3])

print koreluj(a,b,2)

powinno dać w wyniku:

[ 0.5 0.  -1.   0.   0.5]

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.

  1. 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. %i 1)). Na przekątnej macierzy narysuj funkcję autokorelacji odpowiednich kanałów, poza przekątną — 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?
  2. 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. (przypomnienie konstrukcji filtrów)
  3. 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:
    1. 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.
    2. Obliczenie funkcji korelacji wzajemnej dla sygnałów policzonych w punkcie A.
    3. Powtórzenie kroków A i B wiele (np. 1000) razy.
    4. 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.
    5. Powtórzenie kroków A-D dla kolejnych przesunięć [math]\tau[/math].
    6. 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.

    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. Występuje tu zatem problem wielokrotnych porównań. Przypominamy, iż może to doprowadzić do przypadkowego uznania wyników jako „istotnych” 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. korektę Bonferroniego czy 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):

    1. 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.
    2. Obliczenie funkcji korelacji dla sygnałów otrzymanych w punkcie A.
    3. Zapamiętanie maksymalnej wartości bezwzględnej funkcji korelacji z punktu B (maksimum bierzemy po wszystkich przesunięciach i po wszystkich parach kanałów).
    4. 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.
    5. Wyznaczenie 95 centyla rozkładu wartości maksymalnych.
    6. 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.
Przykład wyniku analizy korelacji wzajemnych dla sygnału niefiltrowanego z naniesionymi granicami możliwych fluktuacji.

Wzajemna gęstość widmowa sygnałów i koherencja

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:

[math] S_{xy}(f) = \int _{-\infty }^{\infty }\gamma_{xy}(\tau )e^{-2\pi i f \tau}d\tau [/math]

W celu dalszego omówienia własności funkcji wzajemnej mocy widmowej sygnałów funkcję tę zapiszemy w postaci:

[math] \begin{array}{l} S_{xy}(f) = |S_{xy}(f)|e^{i\phi _{xy}(f)}\\ \\ \phi _{xy} = \arg(S_{xy}) \end{array} [/math]

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 — [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:

[math] \tau = \frac{\phi _{xy}(f)}{2\pi f} [/math]

Podobnie jak w przypadku funkcji autokorelacji i korelacji wzajemnej, funkcję wzajemnej gęstości mocy widmowej można znormalizować:

[math] C_{xy}(f) = \frac{S_{xy}(f)}{\sqrt{S_x(f)S_y(f)}} [/math]

Znormalizowaną postać funkcji wzajemnej gęstości mocy widmowej nazywamy funkcją koherencji. 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).

Kilka słów o koherencji

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 C12 i C21.

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)