(57)【特許請求の範囲】
【請求項１】
振幅・周波数が揺らぎ、ノイズを含む系において、高精度で位相及び振幅を推定する方法であって、
マルコフ確率場（MRF）のもとで確率伝播法による確率的推論手法を用い、瞬時の位相及び振幅を含むパラメータの周辺尤度を算出し、
その周辺尤度を最大化する第二種最尤推定法に基づいて、位相及び振幅のなめらかさを含むハイパーパラメータを算出し、
得られたハイパーパラメータのもとで、位相及び振幅を推定する
ことを特徴とする周期的情報抽出方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、振幅・周波数が揺らぎノイズを含む系において、高精度で周期的情報を抽出する方法に関する。
【背景技術】
【０００２】
神経活動など、振幅・周波数が揺らぎノイズを含む系において、高精度で位相及び振幅を推定したい場合がある。
例えば、α波は、8-13Hz の後頭の視覚野から発生する周期的活動であり、後頭由来の脳波の主要な成分である。近年、このα波の振幅が、working memory task 等のcognitive task のパフォーマンスに影響を与えること（非特許文献１〜２）や、刺激提示時のα波の位相が知覚に影響を与えること（非特許文献３〜５）が報告され、α波の位相及び振幅に注目が集まっている。
【０００３】
従来の脳波を用いた研究では、α波の瞬時における位相及び振幅を推定するために、主にヒルベルト変換（非特許文献６〜８）、ウェーブレット変換（非特許文献９〜１０） が用いられてきた。
【０００４】
しかし、脳波には多くのノイズが含まれていることから、単一試行からこれらの手法で推定したα波の瞬時における位相及び振幅の推定精度が悪いため、試行間平均をとることでノイズを軽減して議論してきた。
このように試行間平均をとることは、試行間で同期した活動しか抽出できないという問題点がある。脳の高次の活動は試行間で同期していない可能性が高く（非特許文献１１）、このような活動は単一試行で議論する必要がある。そして、そのためには瞬時における位相及び振幅の推定精度を向上させる必要がある。
【０００５】
従来の技術では、α波の位相や振幅を推定するには、その精度に限界があった。
【先行技術文献】
【非特許文献】
【０００６】
【非特許文献１】W. Klimesch: Int. J.Psychophysiol. 24 (1996) 61
【非特許文献２】H. Van Dijk, J.M. Schoffelen, R. Oostenveld, and O. Jensen:J. Neurosci. 28 (2008) 1816
【非特許文献３】N.A. Busch, J. Dubois, and R. VanRullen: J. Neurosci. 29 (2009) 7869
【非特許文献４】K.E. Mathewson, G. Gratton, M. Fabiani, D.M. Beck, and T. Ro: J.Neurosci. 29 (2009) 2725
【非特許文献５】F.J. Varela, A. Toro, E.R. John, and E.L. Schwartz: Neuropsychologia  19 (1981) 675
【非特許文献６】A. Matani, Y. Naruse, Y. Terazono, T. Iwasaki, N. Fujimaki,and T. Murata:to be published in IEEE Trans.Biomed. Eng
【非特許文献７】V.V. Nikulin, and T. Brismar: Neuroscience 137 (2006)647
【非特許文献８】V.V. Nikulin, K. Linkenkaer-Hansen, G. Nolte, S. Lemm,K.R. Muller,R.J.Ilmoniemi, and G. Curio: Eur. J. Neurosci.25 (2007) 3146
【非特許文献９】J.M. Palva, S. Palva, and K. Kaila: J. Neurosci. 25 (2005)3962
【非特許文献１０】Y. Naruse, A. Matani, Y. Miyawaki, and M. Okada: to be published in Hum.Brain. Mapp
【非特許文献１１】Arieli, A. Sterkin, A. Grinvald, and A. Aertsen: Science 273 (1996) 1868
【非特許文献１２】K. Linkenkaer-Hansen, V.V. Nikouline, J.M. Palva, and R.J. Ilmoniemi:J.Neurosci. 21 (2001) 1370
【非特許文献１３】J. Pearl: Probabilistic Reasoning in Intelligent Systems(Morgan Kaufmann,San Mateo, 1988)
【非特許文献１４】K. Tanaka, H. Shouno, M. Okada, and D.M. Titterington:J. Phys. A Math.Gen. 37 (2004) 8675
【非特許文献１５】K. Takiyama, K. Katahira, and M. Okada: J. Phys. Soc.Jpn. 78 (2009) 064003
【非特許文献１６】K. Watanabe, H. Tanaka, K. Miura, and M. Okada: IEICE Trans. Inf. Syst. E92d (2009) 1362
【非特許文献１７】C.M. Bishop: Pattern Recognition and Machine Learning (Spinger-Verlag,Berlin, 2006)
【非特許文献１８】F.H. Lopes da Silva, J.P. Pijn, D. Velis, and P.C.G. Nijssen:Int.J.Psychophysiol. 26 (1997) 237
【非特許文献１９】C.J. Stam, J.P.M. Pijn, P. Suffczynski, and F.H. Lopes da Silva: Clin.Neurophysiol. 110 (1999) 1801
【非特許文献２０】Y. Naruse, A. Matani, T. Hayakawa, and N. Fujimaki: Neuroimage 32 (2006) 1221
【非特許文献２１】S. Makeig, M. Westerfield, T.P. Jung, S. Enghoff, J.Townsend, E.Courchesne, and T.J. Sejnowski: Science 295 (2002) 690
【非特許文献２２】Mazaheri, and O. Jensen: Proc. Natl. Acad. Sci. U. S.A. 103 (2006) 2948
【非特許文献２３】M.L. Risner, C.J. Aura, J.E. Black, and T.J. Gawne: Neuroimage 45 (2009) 463
【非特許文献２４】W. Klimesch, P. Sauseng, and W. Gruber: Neuroimage 47 (2009) 5
【発明の概要】
【発明が解決しようとする課題】
【０００７】
そこで、本発明は、振幅・周波数が揺らぎノイズを含む系において、高精度で周期的情報を抽出する方法を提供することを課題とする。
【課題を解決するための手段】
【０００８】
上記課題を達成するため、本発明の周期的情報抽出方法は、振幅・周波数が揺らぎ、ノイズを含む系において、高精度で位相及び振幅を推定する方法であって、マルコフ確率場（MRF）のもとで確率伝播法による確率的推論手法を用い、瞬時の位相及び振幅を含むパラメータとその周辺尤度を算出し、その周辺尤度を最大化する第二種最尤推定法に基づいて、位相及び振幅のなめらかさを含むハイパーパラメータを算出し、得られたハイパーパラメータのもとで、瞬時の位相及び振幅を含むパラメータを推定することを特徴とする。
【０００９】
ここで、ハイパーパラメータに、中心周波数も含めてもよい。
【００１０】
ハイパーパラメータに、観測ノイズの強度も含も含めてもよい。
【００１１】
対象となる周期的情報としては、神経活動情報に好適に適用できる。
【００１２】
神経活動情報としては例えば、脳において観測されるα波が挙げられる。
【発明の効果】
【００１３】
本発明によると、脳活動のα波など、ゆらぎとノイズのある系において、位相及び振幅などの周期的情報を高精度で求められる。
【図面の簡単な説明】
【００１４】
【図１】脳波の実測データ(EEG) と観測データ(Observed data) のFourier spectrumの比較を示すグラフ
【図２】位相、振幅に対する確率の典型例を示すグラフ
【図３】観測値(Observed data)と本発明による推定値(MRF) と真値との間の絶対誤差の典型例を示すグラフ
【図４】観測値(Observed data)と本発明による推定値(MRF)と真値との間の平均絶対誤差の試行間平均値を示すグラフ
【図５】本発明(MRF)及びヒルベルト変換(HT)による位相推定値と真の位相値との間の絶対誤差の典型例を示すグラフ
【図６】本発明(MRF)及びヒルベルト変換(HT)による振幅推定値と真の振幅値との間の絶対誤差の典型例を示すグラフ
【図７】本発明(MRF)及びヒルベルト変換(HT) よる位相推定値と真の位相値との間の平均絶対誤差の試行間平均値を示すグラフ
【図８】本発明(MRF)及びヒルベルト変換(HT) よる振幅推定値と真の振幅値との間の平均絶対誤差の試行間平均値を示すグラフ
【発明を実施するための最良の形態】
【００１５】
以下に、本発明の実施形態を図面に基づいて説明する。実施形態は、前記特許文献など従来公知の技術を援用して適宜設計変更可能であり、また、計算の方法及び装置についても従来公知の技術を適宜適用可能である。
以下の実施例では、対象をα波とした例を挙げるが、本発明は、振幅・周波数が揺らぎノイズを含む系一般に適用可能である。
【００１６】
本発明では、マルコフ確率場(MRF) モデルを用い、従来の脳波解析で使われてきた手法よりも精度よく瞬時の位相及び振幅を推定する方法を提供する。観測されたα波は、特定の周波数の位相及び振幅があり、そこに観測ノイズが加わっているとモデル化することができる。
【００１７】
α波の位相は、各人により異なるα波中心周波数の周辺で時間とともになめらかに変化しながら進行し、その振幅も時間とともになめらかに変化している（非特許文献１２）。そこで、本発明では、特定のサンプリングポイントの瞬時の位相及び振幅は、その前後のサンプリングポイントからなめらかに変化しているという位相及び振幅に関するMRFモデルによりα波を記述する。
【００１８】
しかし、位相及び振幅の変化のなめらかさやα波の中心周波数は、被験者ごとに異なり、さらに、観測ノイズは実験ごとに異なっているため、定かではない。そこで、各サンプリングポイントにおける瞬時の位相及び振幅といったパラメータのみならず、位相及び振幅の変化のなめらかさ、α波の中心周波数、観測ノイズの強度といったハイパーパラメータも同時に推定する必要がある。
【００１９】
MRF の元では、確率的推論手法の一つである確率伝播法（非特許文献１３）を用いることで瞬時の位相及び振幅といったパラメータの推定及び周辺尤度を計算することが可能である（非特許文献１４〜１６）。そこで、本発明では、この周辺尤度を最大化する第二種最尤推定（非特許文献１７）に基づいて、位相及び振幅のなめらかさ、α波の中心周波数、観測ノイズの強度といったハイパーパラメータを推定する。そして、この推定したハイパーパラメータの元で瞬時の位相及び振幅といったパラメータを推定するアルゴリズムを構築する。
【００２０】
以下の実施例では、まず、脳波を模した人工データを用いて、本発明を用いて推定した時に観測ノイズを軽減できるかについて検証する。続いて、本発明を用いて推定した時の瞬時の位相及び振幅が、脳波データ解析で従来一般的に使われるヒルベルト変換よりも推定精度がよいことを示す。
【００２１】
脳波によって計測されるα波の時系列における各サンプリングポイントk の観測データyk は、サンプリングポイントk におけるα波の瞬時の位相及び振幅をそれぞれxk、ak とすると、そこに各時刻間で独立なノイズσ がのっているという次式で表すことができる。
【００２２】
【数１】
0005660436.tif 000002
【００２３】
ここで、ノイズをガウスノイズであると仮定すると、尤度関数は次式で表すことができる。
【００２４】
【数２】
0005660436.tif 000003
【００２５】
ここで、α はノイズ強度を表すハイパーパラメータである。
α波の周波数及び振幅は、時間とともになめらかに変化することが知られているので（非特許文献１１）、位相xk はα波中心周波数の周辺で時間とともになめらかに変化しながら進行し、また振幅ak も時間とともになめらかに変化する場合を考える。
このときの観測値y = {y1, y2, . . . , yN} から、振幅a = {a1, a2, . . . , aN}と、位相x = {x1, x2. 
.  . . , xN} を、確率伝播法により同時に推定する。
ここで、瞬時位相及び振幅の生成モデルは、次式とする。
【００２６】
【数３】
0005660436.tif 000004
【００２７】
【数４】
0005660436.tif 000005
【００２８】
ただし、I0(β) は、modified Bessel function of order zero であり、次式で表すことができる。
【００２９】
【数５】
0005660436.tif 000006
【００３０】
振幅ak は、0 以上であるため、規格化定数Zampは次式となる。
【００３１】
【数６】
0005660436.tif 000007
【００３２】
ここで、β、 γ、ω はそれぞれ、位相変化のなめらかさ、振幅変化のなめらかさ、α波中心周波数を表すハイパーパラメータである。
振幅と位相の事前分布を独立と仮定すると、次式となる。
【００３３】
【数７】
0005660436.tif 000008
【００３４】
マルコフ確率場を仮定していることから、観測データ数をN とすると、次式で表わされる。
【００３５】
【数８】
0005660436.tif 000009
【００３６】
また、周辺分布は次式で表わされる。
【００３７】
【数９】
0005660436.tif 000010
【００３８】
ここで、次式より、Z = p(y) とおくと、下式が得られる。
【００３９】
【数１０】
0005660436.tif 000011
【００４０】
【数１１】
0005660436.tif 000012
【００４１】
この式における重積分は一般に計算困難である。しかし、確率伝播法により以下で定義されるmessage ML、MR に関する漸化式を構築することで、数値積分をする上での振幅の分割数をna、位相の分割数をnp とすると2(N −1) 重積分の計算量をO(na(N-1)×np(N-1))からO((N−1)na2×(N−1)np2)へ削減することができる。
ここでML、MR を、次式とすると下式が得られる。
【００４２】
【数１２】
0005660436.tif 000013
【００４３】
【数１３】
0005660436.tif 000014
【００４４】
【数１４】
0005660436.tif 000015
【００４５】
0005660436.tif 000016
【００４６】
【数１５】
0005660436.tif 000017
【００４７】
次に、ML(ak, xk)、MR(ak, xk) の漸化式を構築する。
【００４８】
【数１６】
0005660436.tif 000018
【００４９】
【数１７】
0005660436.tif 000019
【００５０】
これらを用いることで、任意のk におけるML(ak, xk)、 MR(ak, xk)を計算できる。
また、ak、 xk を周辺化することで、次式として、規格化定数Z を求めることができる。
【００５１】
【数１８】
0005660436.tif 000020
【００５２】
0005660436.tif 000021
【００５３】
【数１９】
0005660436.tif 000022
【００５４】
0005660436.tif 000023
【００５５】
【数２０】
0005660436.tif 000024
【００５６】
【数２１】
0005660436.tif 000025
【００５７】
0005660436.tif 000026
【００５８】
0005660436.tif 000027
【００５９】
0005660436.tif 000028
【００６０】
図４は、観測値(Observed data) 及び、本発明による推定値(MRF)と真値との間の平均絶対誤差の試行間平均値を示すグラフである。
真値からの平均絶対誤差を観測データ及び本発明の推定値とで比較した結果、平均絶対誤差はすべての試行において観測データよりも推定値の方が小さかった(Wilcoxon signed rank test、 p < 0.0001)。このことから、本発明を用いて推定することで、観測ノイズを軽減できていることが示された。
【００６１】
次に、本発明による瞬時位相及び振幅の推定精度を調べた。
図５は、本発明(MRF)及びヒルベルト変換(HT)による位相推定値と真の位相値との間の絶対誤差の典型例を示すグラフであり、図６は、本発明(MRF)及びヒルベルト変換(HT)による振幅推定値と真の振幅値との間の絶対誤差の典型例を示すグラフである。
図７は、本発明(MRF)及びヒルベルト変換(HT) よる位相推定値と真の位相値との間の平均絶対誤差の試行間平均値を示すグラフであり、図８は、本発明(MRF)及びヒルベルト変換(HT) よる振幅推定値と真の振幅値との間の平均絶対誤差の試行間平均値を示すグラフである。
【００６２】
0005660436.tif 000029
【００６３】
【数２２】
0005660436.tif 000030
【００６４】
位相、振幅の平均絶対誤差はともにすべての試行において、本発明にて推定した方が小さかった(Wilcoxon signed rank test, p < 0.0001)。
【００６５】
以上のように、本発明では、α波を特定のサンプリングポイントの瞬時の位相及び振幅は、前後のサンプリングポイントからなめらかに変化するというMRF モデルにより記述し、これを確率伝播法を用いて瞬時位相及び振幅を推定し、周辺尤度を計算する手法を構築した。そして、周辺尤度を最大化することでハイパーパラメータを推定した。
尤度関数（数２）は、位相、振幅の非線形項を含む。しばしば非線形項を含む事後分布から隠れ変数を推定する際にはラプラス近似が用いられる（非特許文献１７）。しかし、ラプラス近似は単峰のガウス分布で事後分布を近似してしまうため、図２に示した多峰性の事後分布へ適用は、精度の低い推定値を与える可能性が高く、本発明にて構築したMRF モデルには適さない。
【００６６】
また、本発明では瞬時の位相及び振幅の推定にMPM 推定を用いたが、他にpoterior mean 推定(PM 推定) を用いることもできる。しかし、各瞬時の位相値及び振幅値に対する確率は多峰性を示すので、このような場合のPM 推定結果は、確率が高い部分の間の確率が低い部分を推定してしまう可能性があり、適切ではない。そのため、本発明ではMPM 推定を用いることが適切である。
従来、周波数が揺らいでいるときの脳波における瞬時の位相及び振幅の推定にヒルベルト変換が広く使われてきた。脳波には観測ノイズがのっているにもかかわらず、ヒルベルト変換ではそのノイズも信号と見なして瞬時の位相及び振幅を推定してしまっていたため、精度が悪かった。しかし、本発明では瞬時の位相及び振幅の推定と同時に観測ノイズも軽減していることから、ヒルベルト変換よりも精度よく瞬時位相及び振幅が推定できていると考えられる。
本発明を脳波の実データに適応することにより、より精度よく瞬時位相及び振幅の推定が可能となる。
【００６７】
近年、α波の振幅や位相がcognitive task のパフォーマンスや知覚に影響を与えることが報告されているが（非特許文献１〜５）、試行間平均をとった結果でしか議論できていないため、試行間で同期していない活動を検出できていない。さらに、α波は非線形振動子と見なすことが出来る可能性が示唆されているが（非特許文献１０、１８〜１９）、α波の位相が視覚刺激などの外部刺激により影響を受けるか否かについては結論が出ていない（非特許文献２０〜２４）。この原因も、単一試行における推定精度が悪いことに依存する。本発明によると、高い精度でα波の位相や振幅を推定することで、これらの問題を解決することができる可能性があり、脳波を用いた神経科学的研究の発展に寄与する。
【産業上の利用可能性】
【００６８】
本発明によると、脳活動のα波など、ゆらぎとノイズのある系において、位相及び振幅などの周期的情報を高精度でを求められるので、脳波を用いた神経科学的研究や工学的応用としてブレインマシンインターフェースなどに寄与するので、産業上利用価値が高い。
【図１】
0005660436.tif 000031
【図２】
0005660436.tif 000032
【図３】
0005660436.tif 000033
【図４】
0005660436.tif 000034
【図５】
0005660436.tif 000035
【図６】
0005660436.tif 000036
【図７】
0005660436.tif 000037
【図８】
0005660436.tif 000038