はじめまして
データサイエンス系のアプリケーションエンジニアとして働いています@mtaguchiと申します。
高校時代は生物が好きだったり、大学時代に脳波関連の研究をしていたりということもあり、そちら方面で何か世の中のお役に立てたらなと日々感じています。
一方で、数学はちょっと苦手な方なので、記事を書きながら勉強していきたいと思います。
万人ウケするネタは少なめかもしれませんが、主に備忘録的に書いていきますので応援のLGTMを頂けると嬉しいです。
今日は脳波の解析手法についてです
いきなりコアなネタ：Weighted Phase Lag Index (WPLI)の計算について載せます。先日ご縁があって調べたので。
平易な言葉で言うと、このWPLIという手法では、脳波電極間の結合の強さを計算します。脳の、どの部分とどの部分の結合の重みが大きいか、がわかるだけではなく、例えば重みの和を計算することで、何かのタスクをしている時に中心的な役割を果たしている脳の部位はどこか、という情報も得ることができます。この辺りはグラフ理論に繋がっていくわけですが...これについては後々触れられたらと思います。
さて、そもそも、もとの論文によると、Weighted Phase Lag IndexはPhase Lag Index (PLI) の改良版です[1]とのこと。PLIはいわゆる信号間の位相差を使ってカップリングを計算する手法[2]で0≦PLI≦1の範囲を取る（0:カップリング無し）けど、volume condunctionとノイズに対する感度、位相同期の変化を検出する能力が低いとのことで改良された手法がWPLIのようです。
式だけ引用させていただくと
まずPLIは参考文献[2]によると
PLI = |<sign[\Delta\phi(t_k)]>|
になっていて
WPLIは参考文献[3]によると
WPLI_{n_i,n_j,t} = |<\frac{|\sin(\Delta\phi_{n_i,n_j,\tau})|}{\sin(\Delta\phi_{n_i,n_j,\tau)}}>|
とのこと。
φは論文により添え字が違いますが、どちらも、各生信号をヒルベルト変換し、信号間の位相を計算させた値のようです。
上の2つの式を見る限り、解析を行う全時系列で、位相が0≦φ≦πの時はPLIもWPLIも同じ計算結果になる気がするのですがどうなんでしょう。。
とにもかくにも、WPLIの方がいろんな意味で外乱に強そうなメソッドなことはわかりました。以下は参考文献[3]の内容。
High-density electroencephalography (EEG) with active electrodes allows for monitoring of electrocortical dynamics during human walking but movement artifacts have the potential to dominate the signal. One potential method for recovering cognitive brain dynamics in the presence of gait-related artifact is the Weighted Phase Lag Index.
コーディング内容と結果の比較
恐らくこんな感じ。
論文からコードに落とすのって難しいですね。。。特にPLIは自信がないので、研究に使う際は本物の著者にご確認くださいね！
なお、hilbert関数は、デフォルトでは縦方向に時系列になってないとダメっぽいので、各列に各電極のデータが入っていて、行方向に、つまり、縦方向に時系列が伸びている状態にするようお気をつけください。
function y = pli(x)
[m,n] = size(x);
dif = zeros(m,1);
yh = hilbert(x);
pha = atan2(imag(yh./x),real(yh./x));
for ii = 1:n
for jj = 1:n
dif = pha(:,ii) - pha(:,jj);
y(ii,jj) = abs(mean(sign(dif)));
end
end
end
function y = wpli(x)
[m,n] = size(x);
dif = zeros(m,1);
yh = hilbert(x);
pha = atan2(imag(yh./x),real(yh./x));
for ii = 1:n
for jj = 1:n
dif = pha(:,ii) - pha(:,jj);
y(ii,jj) = abs(mean(abs(sin(dif))./sin(dif)));
end
end
end
ということで、上記の関数に、以下のような脳波を入れて計算させるとどうなるか計算させてみました。データはFieldTripのEEGサンプルデータSubjectEEG.zipを使用しており、今回は0.1秒分の脳波だけを使ってみました。（サンプリング周波数が500で、50点分のプロットとなります。）
結果はこんな感じです。
うーん、、やっぱりあまり変わらない、という感想ですが、近接行列（関数の出力y）の計算結果を見てみると値が少し違うので、π<φ<2πの成分が効いてきているのでしょうか。
ちなみに、WPLIを計算させる脳波の時間を1秒にのばしたところ、下記のような結果となりました。
カラーバーのスケールも変化しているので単純な比較はできませんが、より、結合が強い電極間と弱い電極間の差がくっきりと表れるようになったかと思います。ただし、PLIとWPLIの値は、こちらもあまり差がありませんね。
推奨されるデータ長など、メソッドがあるのかもしれません。少し疑問が残りますが、今日はここまでです。
あとがき
ヒルベルト変換
今回こちらの解析を行うにあたって、そもそもヒルベルト変換って何だったかなと思い、勉強に何日か要してしまいました。
色々なサイトを探しましたが、最終的に、自分のもやもやを1番多く解決してくださったのがこちらのスライドでした[4]。
長野日本無線株式会社様のどなたかと思いますが、お名前がわかりかねず、こちらにて御礼申し上げたいと思います。
上[4]にもある通り、
世の中に存在する信号は、すべて実信号である
実信号の観測において、信号の変化が、振幅の変化に起因するものか、位相の変化に起因するものかを判別できない
ということで、脳波で「位相」を見るときにはヒルベルト変換が必須なわけですね。昔、脳波の振幅にばかり目を向けていた自分にとっては、大きな一歩を踏み出せたような気持ちです。
さいごに
そうなると、コヒーレンスは？相関係数とは？など色々と気になってきたので、今後、また脳波を題材に数学的な紐解きも頑張ってみたいと思いますので、応援のLGTMを頂けるととても嬉しいです。
引き続き宜しくお願い致します。
参考文献
[1] Vinck, M., Oostenveld, R., Van Wingerden, M., Battaglia, F. and Pennartz, C.M., 2011. An improved index of phase-synchronization for electrophysiological data in the presence of volume-conduction, noise and sample-size bias. Neuroimage, 55(4), pp.1548-1565.
[2] Stam, C.J., Nolte, G. and Daffertshofer, A., 2007. Phase lag index: assessment of functional connectivity from multi channel EEG and MEG with diminished bias from common sources. Human brain mapping, 28(11), pp.1178-1193.
[3] Lau, T.M., Gwin, J.T., McDowell, K.G. and Ferris, D.P., 2012. Weighted phase lag index stability as an artifact resistant measure to detect cognitive EEG activity during locomotion. Journal of neuroengineering and rehabilitation, 9(1), pp.1-9
[4] Nagano Japan Radio Co.,Ltd., Understanding Analytic Signal & Digital Signal Processing for Software Defined Radio, Chronix RF/Microwave Seminar 2002 Tokyo & Kyoto, July 2002