Pracownia EEG/EEG wlasnosci EEG spoczynkowego: Różnice pomiędzy wersjami

Z Brain-wiki
(Nie pokazano 27 wersji utworzonych przez 3 użytkowników)
Linia 1: Linia 1:
 +
[[Pracownia EEG|Pracownia EEG]] / Własności EEG spoczynkowego: funkcja autokorelacji i widmo
  
 
=Estymacja funkcji autokowariancji, autokorelacji i koherencji sygnału.=
 
=Estymacja funkcji autokowariancji, autokorelacji i koherencji sygnału.=
Linia 11: Linia 12:
 
<math>
 
<math>
 
\gamma (\tau) = \mathrm{cov}(x(t),x(t-\tau ))=\mathrm{E}[(x(t)-\mu )(x(t-\tau )-\mu )]
 
\gamma (\tau) = \mathrm{cov}(x(t),x(t-\tau ))=\mathrm{E}[(x(t)-\mu )(x(t-\tau )-\mu )]
</math>
+
</math> &nbsp;&nbsp;&nbsp;&nbsp;(1)
 
</equation>
 
</equation>
  
Linia 19: Linia 20:
 
<math>
 
<math>
 
\mu = \mathrm{E}[x(t)]
 
\mu = \mathrm{E}[x(t)]
</math>
+
</math> &nbsp;&nbsp;&nbsp;&nbsp;(2)
 
</equation>
 
</equation>
  
W przypadku sygnałów ciągłych estymację tę można zapisać w poniższy sposób:
+
W przypadku sygnałów ciągłych estymatę tę można zapisać w poniższy sposób:
  
 
<equation id="uid81">
 
<equation id="uid81">
 
<math>
 
<math>
 
\gamma (\tau ) = \frac{1}{T}\int _0^{T}(x(t)-\mu )(x(t-\tau )-\mu )dt
 
\gamma (\tau ) = \frac{1}{T}\int _0^{T}(x(t)-\mu )(x(t-\tau )-\mu )dt
</math>
+
</math> &nbsp;&nbsp;&nbsp;&nbsp;(3)
 
</equation>
 
</equation>
  
Linia 35: Linia 36:
 
<math>
 
<math>
 
\gamma (k) = \frac{1}{N-1}\sum _{i=0}^{N-k}(x(i+k)-x_s)(x(i)-x_s)
 
\gamma (k) = \frac{1}{N-1}\sum _{i=0}^{N-k}(x(i+k)-x_s)(x(i)-x_s)
</math>
+
</math> &nbsp;&nbsp;&nbsp;&nbsp;(4)
 
</equation>
 
</equation>
  
Linia 43: Linia 44:
 
<math>
 
<math>
 
x_s = \frac{\sum _{i=0}^{N}x(i)}{N}
 
x_s = \frac{\sum _{i=0}^{N}x(i)}{N}
</math>
+
</math> &nbsp;&nbsp;&nbsp;&nbsp;(5)
 
</equation>
 
</equation>
  
Funkcja autokowariancji może osiągać dowolne wartości, dlatego aby można było porównać przebieg tej funkcji np. pomiędzy dwoma sygnałami, wprowadzono wersję znormalizowaną tej funkcji - ''funkcję autokorelacji''. Normalizacja ta wygląda następująco:
+
Funkcja autokowariancji może osiągać dowolne wartości, dlatego aby można było porównać przebieg tej funkcji np. pomiędzy dwoma sygnałami, wprowadzono wersję znormalizowaną tej funkcji &mdash; ''funkcję autokorelacji''. Normalizacja ta wygląda następująco:
  
 
<equation id="uid84">
 
<equation id="uid84">
 
<math>
 
<math>
 
\rho (k) = \frac{\gamma (\tau )}{\sigma^2}
 
\rho (k) = \frac{\gamma (\tau )}{\sigma^2}
</math>
+
</math> &nbsp;&nbsp;&nbsp;&nbsp;(6)
 
</equation>
 
</equation>
  
Linia 59: Linia 60:
 
<math>
 
<math>
 
\sigma ^2 = \mathrm{E}[(x(t)-\mu )^2]
 
\sigma ^2 = \mathrm{E}[(x(t)-\mu )^2]
</math>
+
</math> &nbsp;&nbsp;&nbsp;&nbsp;(7)
 
</equation>
 
</equation>
  
Wariancję sygnału (<math>\gamma (0)=\sigma ^2</math>) można wyrazić przez funkcję autokowariancji dla przesunięcia <math>\tau =0</math>.  Wynika z tego, że funkcja korelacji przyjmuje wartości z zakresu <math>[-1, \, 1]</math>. Ostatecznie estymator funkcji autokorelacji można zapisać jak poniżej:
+
Wariancję sygnału (<math>\gamma (0)=\sigma ^2</math>) można wyrazić przez funkcję autokowariancji dla przesunięcia <math>\tau =0</math>.  Wynika z tego, że funkcja korelacji przyjmuje wartości z zakresu [&minus;1,&nbsp;1]. Ostatecznie estymator funkcji autokorelacji można zapisać jak poniżej:
  
 
<equation id="uid86">
 
<equation id="uid86">
 
<math>
 
<math>
 
\rho(k) = \frac{\gamma (k)}{\gamma (0)}
 
\rho(k) = \frac{\gamma (k)}{\gamma (0)}
</math>
+
</math> &nbsp;&nbsp;&nbsp;&nbsp;(8)
 
</equation>
 
</equation>
  
Linia 73: Linia 74:
  
  
Cechą charakterystyczną funkcji autokorelacji jest to, iż w przypadku sygnałów harmonicznych, przebieg funkcji ma charakter okresowy, z okresem takim samym o okres badanego sygnału. W przypadku szumu, funkcja autokorelacji ma kształt funkcji delta Diraca.
+
Cechą charakterystyczną funkcji autokorelacji jest to, iż w przypadku sygnałów harmonicznych, przebieg funkcji ma charakter okresowy, z okresem takim samym jak okres badanego sygnału. W przypadku szumu, funkcja autokorelacji ma kształt funkcji delta Diraca.
  
===Polecenie:===
+
==Polecenie:==
 
Zaimplementuj funkcję do obliczania funkcji korelacji zgodnie ze wzorem (4). Funkcja powinna przyjmować dwa wektory i maksymalne przesunięcie wzajemne tych wektorów, natomiast zwracać powinna wektor zawierający funkcję autokorelacji. Wywołanie przykładowe:
 
Zaimplementuj funkcję do obliczania funkcji korelacji zgodnie ze wzorem (4). Funkcja powinna przyjmować dwa wektory i maksymalne przesunięcie wzajemne tych wektorów, natomiast zwracać powinna wektor zawierający funkcję autokorelacji. Wywołanie przykładowe:
 
<source lang = python>
 
<source lang = python>
Linia 125: Linia 126:
 
-->
 
-->
  
 
+
==Zadanie 1:  Funkcje autokowariancji i autokorelacji==
 
+
W tym zadaniu posłużymy się sygnałami zarejestrowanymi w punkcie 3. poprzedniego ćwiczenia. Zaobserwuj, na którym kanale rytm alfa osiąga najwyższą wartość. Następnie zaimplementuj w Pythonie następujące kroki:
===Zadanie 1:  Funkcje autokowariancji i autokorelacji===
 
W tym zadaniu posłużymy się sygnałami zarejestrowanymi w punkcie 4. poprzedniego ćwiczenia. Zaobserwuj, na którym kanale rytm alfa osiąga najwyższą wartość. Następnie zaimplementuj w Pythonie następujące kroki:
 
  
 
# Wczytaj dane z wybranego kanału.
 
# Wczytaj dane z wybranego kanału.
 
# Oblicz funkcję autokorelacji dla sygnału zarejestrowanego w warunkach, gdy osoba badana siedziała z otwartymi oczami. Narysuj autokorelogram, to jest wykres wartości funkcji autokorelacji względem przesunięcia <math>\tau </math>. Oś <math>\tau </math> wyskaluj w sekundach.
 
# Oblicz funkcję autokorelacji dla sygnału zarejestrowanego w warunkach, gdy osoba badana siedziała z otwartymi oczami. Narysuj autokorelogram, to jest wykres wartości funkcji autokorelacji względem przesunięcia <math>\tau </math>. Oś <math>\tau </math> wyskaluj w sekundach.
# Powtórz krok 2, tym razem za sygnału zebranego w warunkach czuwania z zamkniętymi oczami.
+
# Powtórz krok 2, tym razem dla sygnału zebranego w warunkach czuwania z zamkniętymi oczami.
 
# Porównaj autokorelogramy.
 
# Porównaj autokorelogramy.
  
==Związek autokorelacji z widmem sygnału==
+
=Związek autokorelacji z widmem sygnału=
===Wstęp===
+
==Wstęp==
 
Zgodnie z twierdzeniem Chinczyna, z którym zapoznaliśmy się na wykładzie z [[Twierdzenie_Wienera-Chinczyna|Analizy Sygnałów]], widmową gęstość mocy sygnału można policzyć jako transformatę Fouriera funkcji autokowariancji:
 
Zgodnie z twierdzeniem Chinczyna, z którym zapoznaliśmy się na wykładzie z [[Twierdzenie_Wienera-Chinczyna|Analizy Sygnałów]], widmową gęstość mocy sygnału można policzyć jako transformatę Fouriera funkcji autokowariancji:
  
 
<equation id="uid93">
 
<equation id="uid93">
 
<math>
 
<math>
S(f) = \int _{-\infty }^{\infty }\gamma (\tau )e^{-2\pi i f \tau}d\tau </math>
+
S(f) = \int _{-\infty }^{\infty }\gamma (\tau )e^{-2\pi i f \tau}d\tau </math> &nbsp;&nbsp;&nbsp;&nbsp;(9)
 
</equation>
 
</equation>
  
Linia 256: Linia 255:
 
-->
 
-->
  
 
+
==Polecenie ==
=====Polecenie 1=====
+
Zaimplementuj funkcję obliczającą transformację Fouriera dyskretyzując wzór (9) dla zadanego wektora częstości <tt>f</tt> i zadanej częstości próbkowania sygnału (tutaj: 10).
Zaimplementuj funkcję obliczającą transformatę Fouriera dyskretyzując wzór (9) dla zadanego wektora częstości <tt>f</tt> i zadanej częstości próbkowania sygnału (tutaj 10.0):
+
<!--
 
 
 
Wywołanie przykładowe:
 
Wywołanie przykładowe:
 
<source lang = python>
 
<source lang = python>
Linia 269: Linia 267:
 
</source>
 
</source>
  
Powinno dać:
+
 
[  3.15975012e-16 +5.19678720e-16j  1.05325004e-16 +3.51083347e-16j
 
  -4.56408351e-16 -2.10650008e-16j  4.91516686e-16 +1.58113883e+00j
 
  -1.40433339e-16 -7.02166694e-17j  0.00000000e+00 +0.00000000e+00j
 
  -1.40433339e-16 +7.02166694e-17j  4.91516686e-16 -1.58113883e+00j
 
  -4.56408351e-16 +2.10650008e-16j  1.05325004e-16 -3.51083347e-16j]
 
  
 
Natomiast wywołanie:
 
Natomiast wywołanie:
Linia 289: Linia 282:
  
 
[[Plik:Fourier_test.png]]  
 
[[Plik:Fourier_test.png]]  
 +
-->
 +
 +
 
<!-- <tt>*</tt> -->
 
<!-- <tt>*</tt> -->
 
<!--
 
<!--
 
import numpy as np
 
import numpy as np
 
import pylab as py
 
import pylab as py
 +
 +
 +
 +
def koreluj_pol(x,y,max_tau):
 +
    x = x - np.mean(x)
 +
    y = y - np.mean(y)
 +
    N= len(x)
 +
    cor = np.zeros(max_tau+1)
 +
    cor[0] = np.sum(x[:]*y[:])
 +
    for i in range(1,max_tau+1):
 +
        cor[i] = np.sum(x[i:]*y[:-i])
 +
       
 +
    cor = cor /(N)
 +
    return cor
 +
 +
def koreluj(x,y,max_tau):
 +
    cor = np.zeros(2*max_tau+1)   
 +
    cor[max_tau:] = koreluj_pol(x,y,max_tau)
 +
    tmp = koreluj_pol(y,x,max_tau)
 +
    cor[:max_tau+1] = tmp[-1::-1]
 +
    return cor 
 +
   
 +
 +
def fourier_chin(x, FF, FS):
 +
    ak = koreluj(x,x,len(x)-1)
 +
    M = np.zeros(len(FF),dtype='complex')
 +
    for i,f in enumerate( FF):
 +
        for tau in range(len(ak)):
 +
    M[i] += ak[tau]*np.exp(-2*np.pi*1j*f*( tau- len(x))/FS)
 +
    #M/=np.sqrt(len(ak))
 +
    return M
 +
   
 
def fourier(ak, FF, FS):
 
def fourier(ak, FF, FS):
 +
 
 
     M = np.zeros(len(FF),dtype='complex')  
 
     M = np.zeros(len(FF),dtype='complex')  
 
     for i,f in enumerate( FF):
 
     for i,f in enumerate( FF):
Linia 300: Linia 329:
 
     M/=np.sqrt(len(ak))
 
     M/=np.sqrt(len(ak))
 
     return M
 
     return M
t= np.arange(0,1,0.1)
+
t= np.arange(0,1,0.05)
 
x = np.sin(2*np.pi*2*t)
 
x = np.sin(2*np.pi*2*t)
 
f = np.arange(-5,5,0.01)
 
f = np.arange(-5,5,0.01)
X = fourier(x,f,10.0)
+
 
 +
X = fourier_chin(x,f,10.0)
 +
Xf = fourier(x,f,10.0)
 
print X
 
print X
py.plot(f,np.abs(X))
+
py.plot(f,np.abs(X), f,np.abs(Xf)**2)
 
py.show()
 
py.show()
  
 
-->
 
-->
  
 
+
==Zadanie 2: Związek autokorelacji z widmem sygnału==
 
 
===Zadanie 2: Związek autokorelacji z widmem sygnału===
 
 
Oblicz gęstość widmową mocy sygnału zarejestrowanego w trakcie czuwania z zamkniętymi oczami, korzystając z twierdzenia Chinczyna oraz [[Nieparametryczne_widmo_mocy#Metoda_Welcha | metodą Welcha]].
 
Oblicz gęstość widmową mocy sygnału zarejestrowanego w trakcie czuwania z zamkniętymi oczami, korzystając z twierdzenia Chinczyna oraz [[Nieparametryczne_widmo_mocy#Metoda_Welcha | metodą Welcha]].
Znajdź częstość rytmu <math>\alpha</math> dla osoby, która była badana.
+
Znajdź częstość rytmu &alpha; dla osoby, która była badana.
  
==Funkcja kowariancji i korelacji==
+
==Funkcja kowariancji (wzajemnej)==
=====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:
 
W celu scharakteryzowania zależności wzajemnej dwóch sygnałów losowych, stosuje się funkcję kowariancji, zdefiniowaną w następujący sposób:
Linia 337: Linia 365:
 
</equation>
 
</equation>
  
W przypadku sygnałów ciągłych estymację tę można zapisać w poniższy sposób:
+
W przypadku sygnałów ciągłych estymatę tę można zapisać w poniższy sposób:
  
 
<equation id="uid100">
 
<equation id="uid100">
Linia 365: Linia 393:
 
<equation id="uid103">
 
<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}
+
\rho (\tau) = \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>
 
</math>
 
</equation>
 
</equation>
Linia 371: Linia 399:
 
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.
 
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.
  
 
+
=====Polecenie =====
 
 
====Zadanie 3: 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).  
 
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:
 
Przykładowe wywołanie:
Linia 409: Linia 435:
 
</source>
 
</source>
 
{{hidden end}} -->
 
{{hidden end}} -->
 +
===Zadanie 3===
 
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.
 
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.
  
Linia 429: Linia 456:
 
<li> Obliczenie funkcji  korelacji wzajemnej dla sygnałów policzonych w punkcie A.
 
<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> 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> 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> 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.
 
<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>
 
</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.
+
Procedura opisana powyżej ma jednak '''zasadniczą 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ń''.
 
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.  
+
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 ''p'' 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):
+
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:
  
 
<ol type="A">
 
<ol type="A">
 
<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> 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> 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>    Zapamiętanie maksymalnej wartości bezwzględnej funkcji korelacji z punktu B (maksimum bierzemy po wszystkich przesunięciach i po wszystkich parach kanałów; dla funkcji autokorelacji, ze względu na jej normalizację do 1 dla zerowego przesunięcia, tam maksymalnych wartości poszukujemy dla przesunięć innych niż 0).
 
<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> 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>    Wyznaczenie 95 centyla rozkładu wartości maksymalnych.
Linia 453: Linia 480:
 
[[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.]]
 
[[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.]]
  
==Wzajemna gęstość widmowa sygnałów i koherencja==
+
==Wzajemna gęstość widmowa sygnałów==
--> to chyba przeniesiemy do Labu
+
 
====Zadanie 4: 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:
 
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:
  
Linia 484: Linia 509:
 
</equation>
 
</equation>
  
 +
<!--
 
Podobnie jak w przypadku funkcji autokorelacji i korelacji wzajemnej, funkcję wzajemnej gęstości mocy widmowej można znormalizować:
 
Podobnie jak w przypadku funkcji autokorelacji i korelacji wzajemnej, funkcję wzajemnej gęstości mocy widmowej można znormalizować:
  
Linia 494: Linia 520:
 
Znormalizowaną postać funkcji wzajemnej gęstości mocy widmowej nazywamy funkcją ''koherencji''.  
 
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).
 
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).
 +
-->
 +
 +
===Zadanie 4===
 +
Zaimplementuj funkcję obliczającą wzajemną gęstość widmową dla pary kanałów.
 +
<!--Niech argumentami tej funkcji będą dwa wektory zawierające sygnały, zakres częstości, częstość próbkowania. -->
 +
Oblicz i narysuj macierz gęstości widmowych (własnych i wzajemnych) dla kolejnych par kanałów (tych samych co w zadaniu 3). Wyniki zaprezentuj w postaci kwadratowej macierzy rysunków. Ponieważ są to funkcje zespolone, dobrze jest zaprezentować osobno ich wartość i fazę. Uzyskane wartości bezwzględne narysuj nad przekątną tej macierzy, a fazę pod przekątną.
  
=====Kilka słów o koherencji=====
+
===Zadanie 5===
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>
+
Przygotuj sygnał dwukanałowy, w którym jako pierwszy sygnał wybierz fragment sygnału EEG (z danych zebranych wcześniej) o długości 2000 próbek, a jako drugiego sygnału użyj tego samego fragmentu EEG, ale opóźnionego o wybraną liczbę (1 - 5) próbek. Oblicz widma wzajemne tych sygnałów i zaprezentuj ich fazy na rysunku. Na podstawie tych widm znajdź wartość przesunięcia czasowego tych sygnałów.
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=====
+
Obliczenia powtórz w przypadku, gdy do drugiego sygnału dodany będzie szum o wariancji równej 0,25 wariancji oryginalnego sygnału.
Zaimplementuj funkcję obliczającą koherencję dla pary kanałów.
 
<!--Niech argumentami tej funkcji będą dwa wektory zawierające sygnały, zakres częstości, częstość próbkowania. -->
 
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)
 

Wersja z 07:59, 22 paź 2019

Pracownia EEG / Własności EEG spoczynkowego: funkcja autokorelacji i widmo

Estymacja funkcji autokowariancji, autokorelacji i koherencji sygnału.

Wstęp

Z funkcjami tymi spotkaliśmy się już na zajęciach z analizy sygnałów.

Funkcja autokowariancji sygnału charakteryzuje liniową zależność wartości tego sygnału w danej określonej chwili czasu od wartości (tego samego sygnału) w innej chwili. W przypadku stacjonarnych procesów stochastycznych, przebieg tej funkcji nie zależy od czasu. Oznacza to, że obliczając funkcję autokorelacji sygnału pomiędzy chwilą czasu [math]x(t)[/math] i [math]x(t+\tau )[/math] otrzymamy tę samą wartość, jak dla przypadku obliczania funkcji autokorelacji pomiędzy momentami [math]x(t + T)[/math] i [math]x(t + T+\tau )[/math], gdzie [math]T[/math] to dowolny przedział czasu. Innymi słowy, funkcja autokorelacji procesu stacjonarnego zależy tylko od odstępu czasu pomiędzy próbkami [math]\tau[/math], dla którego jest wyznaczana, a nie od konkretnej chwili czasu. Odrębną klasę sygnałów stanowią procesy niestacjonarne, w przypadku których funkcja autokorelacji będzie zależeć od czasu [math]t[/math] w którym jest obliczana. Estymator funkcji autokowariancji uzyskuje się poprzez obliczanie iloczynów wartości sygnału [math]x[/math] w chwilach czasu [math]t[/math] czyli [math]x(t)[/math] i wartości sygnału [math]x[/math] w chwili czasu t+τ czyli [math]x(t+\tau)[/math] i uśredniając wartości iloczynów po czasie [math]T[/math]:

[math] \gamma (\tau) = \mathrm{cov}(x(t),x(t-\tau ))=\mathrm{E}[(x(t)-\mu )(x(t-\tau )-\mu )] [/math]     (1)

gdzie:

[math] \mu = \mathrm{E}[x(t)] [/math]     (2)

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

[math] \gamma (\tau ) = \frac{1}{T}\int _0^{T}(x(t)-\mu )(x(t-\tau )-\mu )dt [/math]     (3)

natomiast dla sygnałów dyskretnych jako:

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

gdzie:

[math] x_s = \frac{\sum _{i=0}^{N}x(i)}{N} [/math]     (5)

Funkcja autokowariancji może osiągać dowolne wartości, dlatego aby można było porównać przebieg tej funkcji np. pomiędzy dwoma sygnałami, wprowadzono wersję znormalizowaną tej funkcji — funkcję autokorelacji. Normalizacja ta wygląda następująco:

[math] \rho (k) = \frac{\gamma (\tau )}{\sigma^2} [/math]     (6)

gdzie:

[math] \sigma ^2 = \mathrm{E}[(x(t)-\mu )^2] [/math]     (7)

Wariancję sygnału ([math]\gamma (0)=\sigma ^2[/math]) można wyrazić przez funkcję autokowariancji dla przesunięcia [math]\tau =0[/math]. Wynika z tego, że funkcja korelacji przyjmuje wartości z zakresu [−1, 1]. Ostatecznie estymator funkcji autokorelacji można zapisać jak poniżej:

[math] \rho(k) = \frac{\gamma (k)}{\gamma (0)} [/math]     (8)

Funkcję autokorelacji estymuje się w celu określenia, w jakim stopniu wartości sygnału w danej chwili czasu wpływają na wartości sygnału w kolejnych chwilach czasu. Ma to kluczowe znaczenie przy rozpoznawaniu rodzaju procesów fizycznych odpowiedzialnego za generowanie sygnału. Funkcja ta zawsze mam maksimum dla przesunięcia [math]\tau =0[/math].


Cechą charakterystyczną funkcji autokorelacji jest to, iż w przypadku sygnałów harmonicznych, przebieg funkcji ma charakter okresowy, z okresem takim samym jak okres badanego sygnału. W przypadku szumu, funkcja autokorelacji ma kształt funkcji delta Diraca.

Polecenie:

Zaimplementuj funkcję do obliczania funkcji korelacji zgodnie ze wzorem (4). Funkcja powinna przyjmować dwa wektory i maksymalne przesunięcie wzajemne tych wektorów, natomiast zwracać powinna wektor zawierający funkcję autokorelacji. Wywołanie przykładowe:

a = np.array([1,2,3])
print koreluj(a,a,2)

powinno dać wynik:

[-0.5  0.   1.   0.  -0.5]


Zadanie 1: Funkcje autokowariancji i autokorelacji

W tym zadaniu posłużymy się sygnałami zarejestrowanymi w punkcie 3. poprzedniego ćwiczenia. Zaobserwuj, na którym kanale rytm alfa osiąga najwyższą wartość. Następnie zaimplementuj w Pythonie następujące kroki:

  1. Wczytaj dane z wybranego kanału.
  2. Oblicz funkcję autokorelacji dla sygnału zarejestrowanego w warunkach, gdy osoba badana siedziała z otwartymi oczami. Narysuj autokorelogram, to jest wykres wartości funkcji autokorelacji względem przesunięcia [math]\tau [/math]. Oś [math]\tau [/math] wyskaluj w sekundach.
  3. Powtórz krok 2, tym razem dla sygnału zebranego w warunkach czuwania z zamkniętymi oczami.
  4. Porównaj autokorelogramy.

Związek autokorelacji z widmem sygnału

Wstęp

Zgodnie z twierdzeniem Chinczyna, z którym zapoznaliśmy się na wykładzie z Analizy Sygnałów, widmową gęstość mocy sygnału można policzyć jako transformatę Fouriera funkcji autokowariancji:

[math] S(f) = \int _{-\infty }^{\infty }\gamma (\tau )e^{-2\pi i f \tau}d\tau [/math]     (9)

gdzie:

  • [math]f[/math] — częstość
  • [math]S(f)[/math] — gęstość widmowa mocy


Polecenie

Zaimplementuj funkcję obliczającą transformację Fouriera dyskretyzując wzór (9) dla zadanego wektora częstości f i zadanej częstości próbkowania sygnału (tutaj: 10).


Zadanie 2: Związek autokorelacji z widmem sygnału

Oblicz gęstość widmową mocy sygnału zarejestrowanego w trakcie czuwania z zamkniętymi oczami, korzystając z twierdzenia Chinczyna oraz metodą Welcha. Znajdź częstość rytmu α dla osoby, która była badana.

Funkcja kowariancji (wzajemnej)

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 estymatę 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 (\tau) = \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.

Polecenie

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]

Zadanie 3

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 zasadniczą 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 p 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:

    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; dla funkcji autokorelacji, ze względu na jej normalizację do 1 dla zerowego przesunięcia, tam maksymalnych wartości poszukujemy dla przesunięć innych niż 0).
    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

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]


Zadanie 4

Zaimplementuj funkcję obliczającą wzajemną gęstość widmową dla pary kanałów. Oblicz i narysuj macierz gęstości widmowych (własnych i wzajemnych) dla kolejnych par kanałów (tych samych co w zadaniu 3). Wyniki zaprezentuj w postaci kwadratowej macierzy rysunków. Ponieważ są to funkcje zespolone, dobrze jest zaprezentować osobno ich wartość i fazę. Uzyskane wartości bezwzględne narysuj nad przekątną tej macierzy, a fazę pod przekątną.

Zadanie 5

Przygotuj sygnał dwukanałowy, w którym jako pierwszy sygnał wybierz fragment sygnału EEG (z danych zebranych wcześniej) o długości 2000 próbek, a jako drugiego sygnału użyj tego samego fragmentu EEG, ale opóźnionego o wybraną liczbę (1 - 5) próbek. Oblicz widma wzajemne tych sygnałów i zaprezentuj ich fazy na rysunku. Na podstawie tych widm znajdź wartość przesunięcia czasowego tych sygnałów.

Obliczenia powtórz w przypadku, gdy do drugiego sygnału dodany będzie szum o wariancji równej 0,25 wariancji oryginalnego sygnału.