Elektroencefalografia/Metody analizy sygnałów EEG - analiza widmowa

Metody nieparametryczne

Transformacja Fouriera

Analiza wzrokowa sygnałów, jakkolwiek wciąż stosowana, szczególnie w praktyce klinicznej, nie wykorzystuje współczesnych możliwości analizy danych. Ponieważ duża część danych biomedycznych jest dzisiaj zbierana w sposób cyfrowy, możemy takie dane analizować przy użyciu komputerów. Obszerny dział wiedzy zwany analizą sygnałów dostarcza metod matematycznych służących do precyzyjnej i zaawansowanej oceny informacji, jaka zawarta jest w badanych zapisach. Warto podkreślić, że jakkolwiek wiele metod matematycznych analizy sygnałów było opracowanych już dość dawno, dopiero rozpowszechnienie się komputerów w laboratoriach umożliwiło ich praktyczne zastosowanie.

Analizując wzrokowo fragmenty zapisów EEG (i oczywiście także innego typu dane) możemy wyróżnić w nich struktury różnej postaci i kształtu. Szczególnie istotne są tzw. rytmy czyli składowe o postaci fal o pewnej określonej częstości. Ich kształty mogą być bardzo zróżnicowane w zależności od ilości takich składowych obserwowanych w sygnale, ich częstości i amplitud. Analiza wzrokowa wyróżnia w EEG kilka takich charakterystycznych rytmów; wiadomo też, że ich rola fizjologiczna jest inna. Dobrze więc byłoby móc umieć oddzielić poszczególne rytmy z sygnału i badać ich zachowanie dokładniej. Innymi słowy chcielibyśmy mieć funkcję zależną od częstości mówiącą o zawartości poszczególnych rytmów w sygnale. Funkcję taką nazywamy gęstością widmową mocy lub w skrócie widmem sygnału. Zagadnieniami wyznaczania takiej funkcji i opisem własności sygnału w zależności od częstości zajmuje się analiza widmowa.

Sygnały, które zbieramy w eksperymencie i których zapisy analizujemy wzrokowo są to wartości badanej przez nas wielkości w kolejnych chwilach czasu. Mówimy, że są to wartości rejestrowane (lub funkcje operujące) w dziedzinie czasu.

Do dalszych rozważań dobrze będzie zauważyć następujący fakt Wyobraźmy sobie, że mamy zmierzony nieskończony sygnał (czyli funkcję w dziedzinie czasu) — X(t), t = 0, ..., +∞. Jak sprawdzić czy nasz sygnał zawiera składową o częstości ω? Składowa taka ma postać funkcji sinusoidalnej, tak jak na przykład h(t) = sin(ωt).

Aby zbadać czy nasz sygnał zawiera poszukiwaną składową możemy policzyć iloczyn skalarny naszego sygnału z sygnałem sinusoidalnym. Iloczyn taki możemy uważać za ocenę korelacji tych dwóch sygnałów.

[math]\langle X(t), h(t) \rangle = \int_{0}^{\infty}X(t)h(t) dt= \int_{0}^{\infty}X(t)\sin(\omega t) dt[/math]

Możemy z takich iloczynów utworzyć następujące relacje

[math]\begin{array}{l}G(\mathrm{\omega})=\int_{0}^{\infty}X(t)\sin(\omega t) dt\\F(\omega)=\int_{0}^{\infty}X(t)\cos(\omega t) dt\end{array}[/math]

Relacje te nazywają się transformacjami sinus i kosinus Fouriera. Jak widać w rezultacie dostajemy dla naszego oryginalnego sygnału (funkcji X(t)) funkcje G(ω) lub F(ω) zwane transformatami sinus i cosinus Fouriera funkcji X(t). Transformaty są już funkcjami zależnymi od parametru ω czyli częstości. Mówimy, że są one wyrażone w dziedzinie częstości. Tak więc funkcji w dziedzinie czasu można przyporządkować jednoznacznie pewną inną funkcję operującą w dziedzinie częstości. Relacje te można również odwrócić i dla funkcji z dziedziny częstości uzyskać transformatę odwrotną produkującą funkcję w dziedzinie czasu (z dokładnością do pewnych współczynników po prawej stronie równań).

[math]\begin{array}{l}X(t)=\int_{0}^{\infty}G(\omega)\sin(\omega t) d\omega\\X(t)=\int_{0}^{\infty}F(\omega)\cos(\omega t) d\omega\end{array}[/math]

Jak można zauważyć, transformaty sinusowa czy kosinusowa nie są jeszcze dla nas poszukiwaną odpowiedzią na pytanie o zawartość składowych o określonej częstości w danym sygnale. Jeśli poszukiwana składowa nie pokrywa się akurat ani z funkcją sinus ani cosinus, ale występuje w oryginalnym sygnale z inną fazą, możemy źle ocenić jej wkład. Dlatego też bardziej ogólnym podejściem będzie policzenie tzw. całki Fouriera dającą w wyniku transformatę Fouriera:

[math]C(\omega)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{+\infty}X(t)\mathrm{e}^{-i\omega t} dt[/math]

Ponieważ e^−iωt = cos(ωt)−i sin(ωt), wyrażenie powyższe zawiera splot z kombinacją funkcji sinus i kosinus. W tym przypadku transformata Fouriera C(ω) jest funkcją o wartościach zespolonych (gdyż funkcja podcałkowa zawiera człon zespolony e^iωt), posiadających pewien moduł i fazę. Jest ona naszą poszukiwaną funkcją gęstości widmowej mocy. Moduł wartości tej funkcji mówi o „ilości” poszczególnych częstości w widmie sygnału czyli mocy względnej danej składowej. Faza funkcji C mówi o fazie względnej danej składowej.

Mamy też odwrotną transformatę Fouriera:

[math]X(t)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{+\infty}C(\omega)\mathrm{e}^{i\omega t} d\omega[/math]

Funkcja X i jej transformata opisują ten sam sygnał w różny sposób. Przykłady transformat Fouriera dla wybranych funkcji

X(t)	C(f )
1	δ(f )
δ(t)	1
eiat	δ(f–a/(2π))

Poszukiwane przez nas widmo sygnału, czyli zawartość w sygnale składowych o różnych częstościach najlepiej opisuje tzw. funkcja gęstości widmowej mocy określona poniższym wzorem:

<math>S(\omega)=C(\omega)\cdot \bar{C}(\omega)

Bliżej życia

Zauważmy, że całka Fouriera jest określona dla sygnałów nieskończonych: całkowanie przebiega w granicach od –∞ do +∞. Rejestrowane dane neurobiologiczne, które chcemy badać, są oczywiście skończone. Musimy więc rozszerzyć naszą teorię do badania sygnałów o skończonej długości. Pewną klasą sygnałów, które łatwo można w ten sposób analizować, są sygnały periodyczne. Ponieważ wartości takich sygnałów powtarzają się okresowo w czasie, zakładamy, że nasz sygnał poza czasem obserwacji możemy traktować jak nieskończony sygnał okresowy, przedłużając periodycznie posiadany sygnał poza okno obserwacji. Niestety, w przypadku danych neurobiologicznych dane praktycznie nigdy nie mają charakteru periodycznego, co więcej, są one przeważnie stochastyczne, czyli zawierają komponentę losową, sprawiającą, że ich wartości w czasie nie powtarzają się. Widmo takiego sygnału będzie więc zaburzone. Metoda transformacji Fouriera jest najczęściej wykorzystywaną metodą estymacji widma różnorodnych sygnałów. Jest tak być może dlatego, że istnieje FFT (Fast Fourier Transform, szybka transformata Fouriera) — efektywny algorytm komputerowy obliczania wartości tej transformaty (dla dyskretnych danych i w szczególnych przypadkach, co jednak w wielu sytuacjach okazuje się być wystarczające).

Splot

dla dwóch funkcji f(x) i g(x) możemy wprowadzić operację splotu tych funkcji

[math](f\ast g)(x)=\int_{0}^{x} f(\tau)g(x-\tau)d\tau[/math]

Podobną operację możemy wprowadzić nie tylko dla funkcji ciągłych, ale także dla ciągów liczb a_n i b_n

[math](a\ast b)_m=\sum_{n} a_nb_{m-n}[/math]

Operacja splotu jest niezwykle ważna w dziedzinie analizy sygnałów, gdyż zachodzi następujące twierdzenie o transformacie splotu

[math]\mathrm{F}(f\ast g)=\mathrm{F}(f)\cdot \mathrm{F}(g)[/math]

Oznacza to, że splot dwóch sygnałów w dziedzinie czasu transformuje się na iloczyn transformat tych sygnałów w dziedzinie częstości. Będziemy z tego faktu korzystać w dalszej części tego rozdziału.

Transformacja Z

Transformację Z możemy traktować jako analog transformacji Fouriera dla funkcji dyskretnych (znanych tylko w określonych, równo oddalonych chwilach czasu). Przekształca ona sygnał wejściowy X = {..., x₋₁, x₀, x₁, x₂, ...} na transformatę Z tego sygnału, operującą w tzw. dziedzinie Z, będącej dyskretnym analogiem dziedziny częstości.

Wyróżniamy dwustronną transformację Z:

[math] \mathrm{Z}(X)=\sum_{n=-\infty}^{\infty}x_n z^{-n} [/math]

oraz jednostronną transformację Z:

[math] \mathrm{Z}(X)=\sum_{n=0}^{\infty}x_n z^{-n} [/math]

Jeśli w tekście nie będziemy precyzować rodzaju stosowanej transformacji Z, będziemy mieć na myśli wersję jednostronną. W dodtaku, ponieważ operujemy ciągami współczynników o skończonej długości, górny zakres naszego sumowania będzie się kończył nie w nieskończoności, ale razem z ostatnim wyrazem ciągu.

Podstawiając z = e^iωΔt możemy wyrazić nasze transformaty od częstości ω, tak jak rozumieliśmy ją w przypadku transformacji Fouriera.

Filtry

Filtry cyfrowe są to funkcje, które aplikujemy do sygnału aby zmienić jego widmo. W dziedzinie częstości działanie filtru możemy przedstawić następującym równaniem

[math]Y(f)=H(f) \cdot X(f)[/math]

Działamy tutaj filtrem cyfrowym na sygnał X (tzw. wejście systemu), w wyniku tej operacji otrzymujemy (na tzw. wyjściu) przefiltrowany sygnał Y. Funkcję H nazywamy macierzą przejścia filtru.

Jeśli za funkcję X(f ) wybierzemy funkcję stałą o wartości 1, otrzymamy na wyjściu sygnał, którego względne moce poszczególnych składowych w częstości będą dokładnie takie jak wartości macierzy przejścia. Stosując transformację odwrotną (Fouriera lub Z) do wyrażenia H(f )·X(f ) otrzymamy w wyniku splot transformat odwrotnych h(t) i x(t). Dla funkcji X(f ) stałej, jej transformata odwrotna x(t) to funkcja delta Diraca δ(t). Transformatę odwrotną h(t) macierzy przejścia filtru nazywamy funkcją odpowiedzi impulsowej tego filtru. Opisuje ona zachowanie się filtru w przypadku pojawienia się na wejściu pojedynczego impulsu (delty Diraca).

Przekształcając równanie (Equation 12) możemy napisać

[math]H(f)=Y(f) / X(f)[/math]

Czyli funkcja H będzie mieć w ogólności postać ilorazu pewnych funkcji. Spośród wielu możliwych filtrów interesować nas będą filtry liniowe działające na sygnały próbkowane w czasie, czyli na dyskretne ciągi wartości. W ogólnym przypadku macierz przejścia takiego filtru ma postać (w dziedzinie Z — sygnały dyskretne):

[math]H(z)=B(z) / A(z)=\frac{b_0+b_1z^{-1}+\ldots+b_N z^{-N}}{1-(a_1z^{-1}+\ldots+a_M z^{-M})}[/math]

Licznik i mianownik wyrażenia na H(z) zawierają pewne wielomiany od zmiennej z^–1. Często określamy filtr przez podanie współczynników tych wielomianów: b₀, b₁,... b_N, a₀, a₁, ..., a_M.

Funkcja odpowiedzi impulsowej h(t) takiego filtru ma postać ciągu liczb. Jeśli A(z) ≡ 1, to funkcja h(t) przedstawia skończony ciąg współczynników, a filtr taki nazywany filtrem o skończonej odpowiedzi impulsowej (ang. FIR, finite impulse response). W pozostałych przypadkach h(t) jest nieskończonym ciągiem, a filtr nazywamy filtrem o nieskończonej odpowiedzi impulsowej (ang. IIR, infinite impulse response).

Często określa się tzw. rząd filtru. Jest to większa z liczb (M, N).

Filtry stosujemy głównie w celu usunięcia z sygnału pewnych zakresów częstości. Filtry dolnoprzepustowe mają za zadanie usuwać z sygnału wszystkie częstości powyżej pewnej częstości granicznej. Filtry górnoprzepustowe mają usuwać wszystkie częstości mniejsze niż częstość graniczna. Filtry pasmowe pozostawiają w sygnale (lub usuwają z sygnału) tylko określony zakres częstości, od dolnej częstości granicznej do górnej.

Istnieje bardzo rozbudowana teoria projektowania filtrów w zależności od konkretnych potrzeb. Efektywność tłumienia niepożądanych częstości, zachowanie się filtru w okolicy częstości granicznych i wpływ filtru na pozostawione w sygnale składowe są to parametry, które trzeba dobierać indywidualnie do rozwiązywanego problemu.

Metody parametryczne

Przedstawiona powyżej idea uzyskania widma sygnału metodą transformaty Fouriera nie jest jedyną możliwością w tym zakresie. W ogólności metody estymacji (oszacowania) widma danego sygnału możemy podzielić na dwie kategorie: nieparametryczne i parametryczne. Metody nieparametryczne bazują na wyznaczaniu widma bezpośrednio z wartości sygnału. Taką właśnie metodą jest transformacja Fouriera. Metody parametryczne polegają na założeniu pewnego modelu generacji posiadanych danych. Model powinien posiadać parametry, które dopasowujemy tak, aby badany sygnał jak najlepiej dawał się opisać wybranym modelem. Po dopasowaniu parametrów modelu dalsze wnioskowanie odbywa się już nie na badanym sygnale, ale na własnościach modelu.

W przypadku analizy EEG szeroko stosowany jest model autoregresyjny (AR). Zakłada on, że wartość sygnału w dowolnej chwili czasu t można wyznaczyć z pewnej liczby poprzednich wartości oraz z pewnej składowej czysto losowej E:

[math]X(t) = \sum_{j=1}^{p}A_jX(t-j)+E(t)[/math]

Model taki w większości wypadków bardzo dobrze opisuje sygnały EEG. Wynika to z jego własności i własności samego sygnału EEG. Otóż teoretyczne widmo takiego modelu ma postać pewnej liczby składowych o określonym zakresie częstości na tle szumowym co dobrze odpowiada „rytmom” zawartym w prawdziwym sygnale.

Procedura wyznaczania widma polega na dopasowaniu współczynników modelu A₁,... A_p tak, aby wariancja składowej szumowej była najmniejsza. Do tego celu służą specjalnie opracowane algorytmy, które można znaleźć w literaturze (na przykład metoda Yule’a-Walkera).

Gdy przyjrzymy się równaniu opisującemu model AR, po przekształceniu do poniższej postaci:

[math]E(t)=\sum_{j=0}^{p}A_jX(t-j)[/math]

zauważamy, że jest to tak naprawdę splot współczynników modelu i wartości sygnału.

Stosując transformację Z obu stron równania (Equation 16) przechodzimy z równaniem do dziedziny Z:

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

Dostrzegamy w tym równaniu działanie filtru o macierzy przejścia H = A⁻¹ na sygnał szumowy E. Tak więc model AR możemy przedstawić jako filtr liniowy o macierzy przejścia

[math]\mathbf{H}(z)=\mathbf{A}^{-1}(z)=\frac{1}{\mathbf{A}(z)}=\frac{1}{\mathbf{I}-A_1 z^{-1}-A_2 z^{-2}-\ldots -A_p z^{-p}}[/math]

Widmo mocy otrzymujemy więc zgodnie ze znaną już relacją (wzór Equation 6):

[math]\mathbf{S}=\mathbf{X}\cdot \bar{\mathbf{X}}^{\mathrm{T}} =\mathbf{HE}\cdot \bar{\mathbf{E}}^{\mathrm{T}}\bar{\mathbf{H}}^{\mathrm{T}}=\mathbf{HV\bar{ H}^{\mathrm{T}}}[/math]

Analiza danych wielokanałowych

Mówiliśmy do tej pory o sygnałach w postaci pewnego zbioru wartości X(t). Jeśli dokonujemy pomiaru jednej wartości w poszczególnych chwilach czasu, nasze wartości X(t) są to po prostu liczby. Możemy jednak obserwować i rejestrować wartości jednocześnie z wielu źródeł sygnału. Dzieje się tak w przypadku zapisu EEG z wielu elektrod. W tej sytuacji wartości X(t) są (dla każdej chwili czasu t) wektorami o długości odpowiadającej liczbie obserwowanych źródeł, np. k: X(t) = (X₁(t), X₂(t), ..., X_k(t)), t = (t₁, t₂, ..., t_n).

Oczywiście możemy prowadzić analizę każdego z rejestrowanych sygnałów osobno. Musimy jednak zdawać sobie sprawę z faktu, że wielokanałowy zestaw danych zawiera jeszcze dodatkową informację o współzależności sygnałów między sobą. Analizując sygnały osobno, informację tę całkowicie pomijamy. Aby ją wydobyć, musimy użyć specjalnych funkcji, zaprojektowanych do badania związków między sygnałami.

Większość wzorów z poprzednich rozdziałów ny. analizy widmowej może być w łatwy sposób rozszerzona na przypadek danych wielokanałowych poprzez traktowanie wartości badanego sygnału X(t) jako wielkości wektorowej (wzory z rozdziału o analizie parametrycznej od razu są napisane przy tym założeniu). Widmo mocy wielokanałowego zestawu danych jest w takim wypadku macierzą o wymiarze k×k. Na przekątnej tej macierzy znajdują się widma własne (autowidma) każdego z sygnałów osobno, a poza przekątną mamy widma wzajemne (krosswidma), mówiące o stopniu zależności sygnałów między sobą w dziedzinie częstości.

[math]\mathbf{S}(f)=\left( \begin{array}{cccc} S_{11}(f)&S_{12}(f)&\ldots&S_{1k}(f)\\ S_{21}(f)&S_{22}(f)&\ldots&S_{2k}(f)\\ \vdots& \vdots& \ddots&\vdots \\ S_{k1}(f)&S_{k2}(f)&\ldots&S_{kk}(f)\\ \end{array} \right)[/math]

Koherencje

Widma wzajemne mówią nam o tym, na ile składowe o określonych częstościach przebiegają „wspólnie” dla dwóch sygnałów. Moduł widma wzajemnego mówi na o wspólnej amplitudzie danej składowej, a faza o ich spójnym wzajemnym przesunięciu w fazie w danej częstości. Wadą tej wielkości jest jej zależność od mocy całkowitej sygnału. Wprowadza się więc znormalizowaną miarę współzależności sygnałów w dziedzinie częstości nazywaną koherencją (zwyczajną):

[math]K_{ij}(f)=\frac{S_{ij}(f)}{\sqrt{S_{ii}(f)S_{jj}(f)}}[/math]

Miara ta zbudowana jest z odpowiednich elementów macierzy widmowej. Element macierzy K jest odpowiadającym mu elementem macierzy S podzielonym przez odpowiednie widma własne. W taki sposób moduł koherencji przyjmuje wartości z przedziału [0, 1]. Wartość 0 oznacza brak związku, a wartość 1 identyczność badanych sygnałów (dla danej częstości).

Koherencja jest miarą szeroko stosowaną w analizie najróżnorodniejszych sygnałów. Jest łatwa w użyciu i daje natychmiast wgląd w sytuację: czy dane dwa sygnały mają ze sobą coś wspólnego? Warto zauważyć, że aby koherencja miała wysoką wartość nie wystarczy, aby dwa sygnały zawierały składową o takiej samej częstości. Składowe te w obu sygnałach muszą być ze sobą związane, na przykład posiadać niezmienne (słabo zmienne) przesunięcie fazowe. Dlatego też czasami mówi się, że koherencja opisuje liniowy związek faz sygnałów.

Związki przyczynowe

Istnienie związku między sygnałami może oznaczać, że jeden z nich jest źródłem informacji dla drugiego. Czy możemy określić kierunek wzajemnego wpływu sygnałów? Możemy na przykład zbadać fazę koherencji — kierunek przesunięcia fazowego powinien wskazywać, który sygnał był wcześniejszy niż drugi. W praktyce jednak metoda taka działa tylko dla bardzo prostych sygnałów. Dla danych biologicznych, o charakterze stochastycznym, jest ona praktycznie bezużyteczna. Z drugiej strony, zagadnienie to jest trudne, dlatego w literaturze istnieje wiele propozycji funkcji opisujących związki przyczynowe między sygnałami, działających lepiej lub gorzej w różnych sytuacjach.

Jednym z możliwych podejść jest oparcie definicji funkcji opisującej związki przyczynowe o macierz przejścia modelu autoregresyjnego. Jeśli przyjrzymy się równaniu (Equation 17), zauważymy, że nasz sygnał X w dziedzinie częstości otrzymujemy z transformaty sygnału szumowego E (którego widmo nie zależy od częstości), na który działa macierz przejścia H. Tak więc całość zależności między sygnałami z całego zestawu, również te kierunkowe, zawarte są w tej macierzy. Na tej idei bazuje defincja skierowanej funkcji przejścia (ang. directed transfer function, DTF). Jej nienormalizowana wersja (NDTF) jest zdefiniowana na podstawie odpowiednich elenetów macierzy przejścia H; transmisja z kanału j do i opisana jest jako:

[math]\mathrm{NDTF}(j \rightarrow i)=\theta_{ij}^2(f)=|H_{ij}(f)|^2[/math]

Wartość 0 oznacza brak transmisji z kanału j do i w częstości f. Funkcja powyższa może przyjmować dowolnie wartości. Możemy wprowadzić normalizację tej wartości w taki sposób, że będziemy opisywać stosunek transmisji z kanału j do kanału i do sumy transmisji do kanału i ze wszystkich kanałów. Taka postać funkcji DTF była zaproponowana w pracy []:

[math]\mathrm{DTF}(j \rightarrow i)=\gamma_{ij}^2(f)=\frac{|H_{ij}(f)|^2}{\sum_{m=1}^{k}{|H_{im}(f)|^2}}[/math]

W tym przypadku wartość 0 funkcji oznacza brak transmisji. Maksymalną wartością jest zaś 1 oznaczające, że do kanału i informacja dopływa wyłącznie z kanału j.