Laboratorium EEG/AR 1

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

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)