データ処理装置、及びプログラム
【課題】単一の測定結果からでも信号の推定を可能とするデータ処理装置を提供する
【解決手段】測定されたデータを受け入れ、少なくとも一つのノイズパラメータを含み、ノイズを模式的に表すノイズ演算式について、過去に得られたデータを用いて定められる基準に基づいて、ノイズパラメータを決定する。また測定されたデータに含まれる信号を模式的に表す信号演算式であって、少なくとも一つの信号パラメータを含む信号演算式と、受入手段が受け入れる測定されたデータと、ノイズ演算式とに基づいて、信号パラメータを推定し、測定されたデータに含まれる信号の推定結果を生成して出力する。
【解決手段】測定されたデータを受け入れ、少なくとも一つのノイズパラメータを含み、ノイズを模式的に表すノイズ演算式について、過去に得られたデータを用いて定められる基準に基づいて、ノイズパラメータを決定する。また測定されたデータに含まれる信号を模式的に表す信号演算式であって、少なくとも一つの信号パラメータを含む信号演算式と、受入手段が受け入れる測定されたデータと、ノイズ演算式とに基づいて、信号パラメータを推定し、測定されたデータに含まれる信号の推定結果を生成して出力する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、データ処理装置、及びプログラムに係り、特にノイズの影響除去に関する。
【背景技術】
【０００２】
脳磁計（ＭＥＧ）等、微小信号を計測する装置では、外来の種々のノイズに信号が埋もれてしまうため、ノイズの除去が課題となっている。かかるノイズ除去の一手法として、従来より、複数回の測定結果を加算平均する方法が知られている。また特許文献１には、脳磁場信号に対して、心拍の発生時刻に同期させて加算平均を取ることにより、心磁によるノイズ信号を抽出する技術が開示されている。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特開２００９−１９５５７１号公報
【発明の概要】
【発明が解決しようとする課題】
【０００４】
しかしながら、上記加算平均を用いるノイズの除去方法では、測定回数をなるべく多くしなければならない。ところが測定の対象によっては多数回の測定が困難である場合もある。また脳磁計等の場合、測定回数を増大させるほど被験者の負担が大きくなってしまう。
【０００５】
本発明は上記実情に鑑みて為されたもので、単一の測定結果からでも信号の推定を可能とするデータ処理装置を提供することを、その目的の一つとする。
【課題を解決するための手段】
【０００６】
上記従来例の問題点を解決するための本発明は、データ処理装置であって、測定されたデータを受け入れる受入手段と、少なくとも一つのノイズパラメータを含み、ノイズを模式的に表すノイズ演算式について、過去に得られたデータを用いて定められる基準に基づいて、前記ノイズパラメータを決定する手段と、前記測定されたデータに含まれる信号を模式的に表す信号演算式であって、少なくとも一つの信号パラメータを含む信号演算式と、前記受入手段が受け入れる測定されたデータと、前記ノイズ演算式とに基づいて、前記信号パラメータを推定する推定手段と、前記推定された信号パラメータと、前記信号演算式とを用いて、測定されたデータに含まれる信号の推定結果を生成し、出力する手段と、を含むこととしたものである。
【発明の効果】
【０００７】
本発明によれば、ノイズ及び信号を、パラメータを含んだ模式的な演算式として表現し、そのパラメータを推定するようにしたことで、単一の測定結果からでも信号の推定を可能とした。
【図面の簡単な説明】
【０００８】
【図１】本発明の実施の形態に係るデータ処理装置の例を表す構成ブロック図である。
【図２】本発明の実施の形態に係るデータ処理装置の例を表す機能ブロック図である。
【図３】本発明の実施の形態に係るデータ処理装置による処理例を表すフローチャート図である。
【図４】従来の方法で得られた誘発脳磁場の波形と、これらの波形から得られた等磁界線図の例を表す説明図である。
【図５】ノイズを含めて実際に得られる時系列データの例と、当該時系列データを基に得られた等磁界線図の例を表す説明図である。
【図６】本発明の実施例に係るデータ処理装置が出力する誘発脳磁場の波形と、これらの波形から得られた等磁界線図の例を表す説明図である。
【図７】ノイズを含めて実際に得られる時系列データの別の例と、当該時系列データを基に得られた等磁界線図の別の例を表す説明図である。
【図８】本発明の実施例に係るデータ処理装置で推定して得た等磁界線図の例を表す説明図である。
【図９】本発明の実施例に係るデータ処理装置で推定して得た複数の等磁界線図と、従来の方法で推定された等磁界線図とのコヒーレンスの例を表す説明図である。
【図１０】従来の方法によって推定された等磁界線図の例を表す説明図である。
【発明を実施するための形態】
【０００９】
本発明の実施の形態に係るデータ処理装置を図面を参照しつつ説明する。本実施の形態のデータ処理装置１は、図１に例示するように、制御部１１、記憶部１２、操作部１３、出力部１４、及び入力インタフェース部１５を含んで構成される。以下では、具体的に、脳磁計（ＭＥＧ）の出力するデータ（時間変化する一連のデータ）を入力インタフェース部１５を介して受け入れる例について説明するが、本実施の形態のデータ処理装置１は、脳磁計のデータを処理するものに限られない。
【００１０】
脳磁計のデータには、被験者のいない状態で脳磁計が出力するデータ（ノイズのみのデータ）と、被験者が安静な状態にあるときに脳磁計が出力するデータ（外的刺激を受けていない状態：安静時データと呼ぶ）と、被験者に対して外的刺激を与えたときに脳磁計が出力するデータ（ノイズ、自発脳磁場及び、刺激による誘発脳磁場が累算されたデータ：刺激時データと呼ぶ）とがある。
【００１１】
制御部１１は、ＣＰＵ（Central Processing Unit）等のプログラム制御デバイスであり、記憶部１２に格納されたプログラムに従って動作する。この制御部１１は、入力インタフェース部１５から入力される、外部の測定装置（ここでは例えば脳磁計とする）により測定されたデータを受け入れる。また、この制御部１１は、少なくとも一つのノイズパラメータを含み、ノイズを模式的に表すノイズ演算式について、過去に得られたデータを用いて定められる基準に基づいて、ノイズパラメータを決定する。
【００１２】
さらに制御部１１は、測定されたデータに含まれる信号を模式的に表す信号演算式であって、少なくとも一つの信号パラメータを含む信号演算式と、測定されたデータと、ノイズ演算式とに基づいて、信号パラメータを推定し、当該推定した信号パラメータと、信号演算式とを用いて、測定されたデータに含まれる信号の推定結果を生成し、出力する。この制御部１１の詳しい処理の内容は後に説明する。
【００１３】
記憶部１２は、メモリデバイス等であり、制御部１１によって実行されるプログラムを保持している。このプログラムは、例えばＤＶＤ−ＲＯＭ（Digital Versatile Disc Read Only Memory）等のコンピュータ可読な記録媒体に格納されて提供され、この記憶部１２に複写されたものであってもよい。また、この記憶部１２は制御部１１のワークメモリとしても動作する。
【００１４】
操作部１３は、キーボードやマウスなどを含み、利用者の指示操作を受け入れて、制御部１１に当該指示操作の内容を出力する。出力部１４は、ディスプレイや、プリンタ、ネットワークインタフェース等であり、制御部１１から入力される指示に従って、情報を出力する。
【００１５】
入力インタフェース部１５は、ＵＳＢ（Universal Serial Bus）やパラレルポート等のデータの入力を受け入れるデバイスであり、ここでの例では脳磁計に接続されて、脳磁計が出力するデータを受け入れ、当該受け入れたデータを制御部１１に出力する。
【００１６】
本実施の形態では、脳磁計から入力されるデータは時々刻々と変化する（時間変化する）。そして時刻ｔにおけるデータｍtが、測定の対象たる信号ｓtとノイズとの和により表現されているものとする。すなわち、
【数１】
ここにノイズは、要因ごとに複数の成分の和として表現されており、（１）式ではτtがベースラインドリフト成分、ｑtが交流電源成分、ｅtがベースラインドリフトや交流電源成分以外の定常ノイズ成分、ｎtは、このようにモデル化したことによるデータからの差を表すモデル化誤差成分であるものとしている。また本実施形態では、不定期に生じる測定値のとび（例えばグリッチノイズ）、眼球運動や心拍等の生体から生じるノイズはないものとする。
【００１７】
ノイズの各成分は、具体的に次のような形の演算式で表されるものとする。ベースラインドリフト成分τtの時間発展を、平滑化事前分布を用いた確率的トレンドモデルを採用し、
【数２】
と表す。
【００１８】
ここでvτ ,tはトレンドの変化を起こす外乱であり、平均ゼロ分散στのガウス分布と仮定しておく。なお、この（２）式、及び以下の説明において時刻ｔ＋ｊは、ｔから予め定めた単位時間（例えばＴとする）を用いてｊ×Ｔだけ後の時刻を意味するものである。
【００１９】
さらに（１）式におけるｑtの時間発展を、サンプリング周波数と交流周波数とから決まる自己回帰モデル（ARモデル）を用いて、
【数３】
と表す。
【００２０】
ここでも同様に、vq ,tは平均ゼロ分散σqのガウス分布と仮定する。またΔTはサンプリング周波数と交流電源の周波数から決定される、離散化された周期を表す。
【００２１】
定常ノイズ成分ｅtの時間発展は次数ｌのＡＲモデルを仮定する。すなわち、
【数４】
とする。
【００２２】
これにおいても、ｖe ,tは平均ゼロ分散σeのガウス分布と仮定する。モデル化誤差成分ｎtは、平均ゼロ分散σnのガウス分布と仮定しておく。
【００２３】
本実施の形態の制御部１１は、ノイズ演算式として、（１）から（４）式の和、
【数５】
を用いる。このノイズ演算式に含まれるノイズパラメータを成分とするベクトルθ（モデルパラメータθベクトル）は、従って、
【数６】
となる。
【００２４】
本実施の形態の制御部１１は、機能的には、図２に例示するように、データ記録部２１、モデル学習部２２、ノイズパラメータ管理部２３、信号推定部２４、及び信号出力部２５を含んで構成される。ここにデータ記録部２１は、入力インタフェース部１５を介して受け入れられるデータを記憶部１２に蓄積して保持する。本実施の形態では、このデータ記録部２１は、操作部１３から入力される指示に従い、受け入れたデータを、ノイズのみのデータ、ノイズと自発脳磁場とを含んだデータ、ノイズ、自発脳磁場、及び誘発脳磁場を含んだデータのいずれかに分類して記憶部１２に格納する。
【００２５】
モデル学習部２２は、データ記録部２１によって記憶部１２に蓄積されたデータを参照し、過去に得られたデータ（被験者のいない状態で得られたノイズのみのデータ）を用いて定められる基準に基づいて、（６）式のモデルパラメータθベクトルの成分（各ノイズパラメータ）を決定する。
【００２６】
具体的に、このモデル学習部２２は、過去に得られたノイズのみの時系列のデータｍ1，ｍ2，…，ｍtを記憶部１２から読み出す。そしてベクトルθが与えられたときにデータがこの順に観測される確率ｐ（ｍ1，ｍ2，…，ｍt｜θ）を最大とするベクトルθを求める。すなわち、現実の観測変数ｍtを生成する真の確率モデルは未知であるものの、当該真の確率モデルを、近似的に（５）式で与えられるものとし、真の確率モデルと近似的なモデルとの擬距離であるＫＬ（カルバック・ライブラー）ダイバージェンス（ＫＬ情報量）を最小とする（５）式のパラメータ（ベクトルθの成分）を求めることとするのである。
【００２７】
例えば上記ＫＬ情報量を最小とする赤池情報量基準を用いてベクトルθを決定すればよい。このためモデル学習部２２は、上記確率ｐ（ｍ1，ｍ2，…，ｍt｜θ）の対数log（ｐ（ｍ1，ｍ2，…，ｍt｜θ））を対数尤度とし、赤池情報量基準ＡＩＣを、
ＡＩＣ＝−２×ｐ（ｍ1，ｍ2，…，ｍt｜θ）＋２×Ｎ
として演算する。ここにＮは自由に設定できるモデルパラメータの数を表し、ここではθの次元である、（ｌ＋４）となる。
【００２８】
モデル学習部２２は、このＡＩＣの値を、複数の予め定めたθである値の組ごとに演算する。具体的にここで予め定めたθの値の組は、
ａ1，ａ2，…，ａl，στ，σq，σe，σn
で張られるベクトル空間上に設定された複数の格子点とすることができる。そして演算したＡＩＣの値が最小となるθを、上記複数の予め定めたθである値の組のうちから選択する。つまりモデル学習部２２は、ＡＩＣの値を最小とするベクトルθを、まずはグリッドサーチにより大域的に粗く探索する。
【００２９】
モデル学習部２２は、さらにグリッドサーチにより求めたベクトルθ（θ0とする）を初期値として、非線形最適化により最適化されたベクトル
【数７】
を求める。
【００３０】
ここで非線形最適化の一例としては、上記対数尤度log（ｐ（ｍ1，ｍ2，…，ｍt｜θ））をカルマンフィルタにより求める方法がある。上記グリッドサーチにより求めたベクトルθ0を初期値として、対数尤度を最大とするベクトルθを、例えば最急降下法により予め定めた終了条件を満足するまで繰り返し逐次的に更新して求め、当該求めたθを、
【数７】
とする。ここで終了条件は、更新前のθと更新後のθとの距離が予め定めたしきい値を下回っているなどの広く使われる条件を適用すればよい。
【００３１】
ノイズパラメータ管理部２３は、モデル学習部２２が演算により求めた
【数７】
を、ノイズパラメータとして記憶部１２に格納する。また、このノイズパラメータ管理部２３は、信号推定部２４からノイズパラメータの要求を受けて、記憶部１２に格納したノイズパラメータを読み出し、信号推定部２４に対して出力する。
【００３２】
本実施の形態の信号推定部２４は、自発脳磁場推定部２４ａと、誘発脳磁場推定部２４ｂとを含んで構成される。自発脳磁場推定部２４ａは、記憶部１２に格納されている、ノイズと自発脳磁場とを含んだ安静時データ（安静状態にある被験者を測定した脳磁場のデータ）を読み出す。ここで自発脳磁場推定部２４ａは、自発脳磁場が
【数８】
によって表される確率モデルで近似的に表されるものとする。ここにｒtは時刻ｔの自発脳磁場成分を意味し、δ波(約２Ｈｚ)，θ波(約４Ｈｚ)，α波(約１０Ｈｚ)，β波(約２０Ｈｚ)，low-γ波(約４０Ｈｚ)，high-γ波(約８０Ｈｚ)の成分から構成されるものとし、各成分を、
【数９】
と表したものである。なお、ここでのθは、θ波に対応する自発脳磁場成分であることを意味するものとして便宜的に利用しているもので、パラメータのベクトルθとは異なる。
【００３３】
また、これらの成分の各々はいずれも、次の準周期振動モデルで表現されるものとする。
【００３４】
【数１０】
ここでωは、δ，θ，α，β，γl，γhのそれぞれをとる。また、ｆωは、成分ω（δ，θ，α，β，γl，γhのそれぞれをとる）の中心周波数であり、ｖt(ω)は、平均が０、分散がσω2であるガウス分布に属するとの条件、
【数１１】
を満足する確率的入力である。
【００３５】
以上より、自発脳磁場の成分を表すモデルパラメータであるθベクトルは、
【数１２】
となる。ここにおいて、左辺のθpのθはθベクトルを表すθであり、右辺中、σθ2とあるθは、θ波に対応する自発脳磁場成分を表す便宜的なものである。
【００３６】
自発脳磁場推定部２４ａは、記憶部１２に格納されているノイズと自発脳磁場とを含んだ時系列の安静時データｍ1，ｍ2，…，ｍtを記憶部１２から読み出す。
【００３７】
また自発脳磁場推定部２４ａは、ノイズパラメータ管理部２３に対してノイズパラメータを要求し、この要求に応答してノイズパラメータ管理部２３が出力するノイズパラメータと（５）式とを用いて、ノイズの時系列推定データν1，ν2，…を得る。自発脳磁場推定部２４ａは、ノイズと自発脳磁場とを含んだ時系列の安静時データｍ1，ｍ2，…，ｍtのそれぞれ対応する時刻のデータからノイズの時系列推定データの対応する時刻のデータを差引きして、ｍ′i＝ｍi−νiを求め、ノイズがないと仮定される状態のデータ、ｍ′1，ｍ′2，…，ｍ′tを得る。
【００３８】
そして自発脳磁場推定部２４ａは、（１０）式で表されるベクトルθpが与えられたときに安静時データが実際に測定された時系列順に観測される確率ｐ（ｍ′1，ｍ′2，…，ｍ′t｜θp）を最大とするベクトルθpを求める。この場合も、現実の観測変数ｍ′tを生成する真の確率モデルは未知であるものの、当該真の確率モデルを、近似的に（７）式で与えられるものとし、真の確率モデルと近似的なモデルとの擬距離であるＫＬ（カルバック・ライブラー）ダイバージェンス（ＫＬ情報量）を最小とする（７）式のパラメータ（（１０）式のベクトルθpの成分）を求めることとするのである。
【００３９】
例えば上記ＫＬ情報量を最小とする赤池情報量基準を用いて（１０）式のベクトルθを決定すればよい。このため自発脳磁場推定部２４ａは、上記確率ｐ（ｍ′1，ｍ′2，…，ｍ′t｜θp）の対数log（ｐ（ｍ′1，ｍ′2，…，ｍ′t｜θp））を対数尤度とし、赤池情報量基準ＡＩＣを、
ＡＩＣ＝−２×ｐ（ｍ′1，ｍ′2，…，ｍ′t｜θp）＋２×Ｎ
として演算する。ここにＮは自由に設定できるモデルパラメータの数を表し、ここでは（１０）式のベクトルθpの次元である、６となる。
【００４０】
自発脳磁場推定部２４ａは、このＡＩＣの値を、複数の予め定めたθpである値の組ごとに演算する。具体的にここで予め定めたθpの値の組は、σω2（ただしωはδ，θ，α，β，γl，γhのそれぞれをとる）で張られるベクトル空間上に設定された複数の格子点とすることができる。そして演算したＡＩＣの値が最小となるθの値の組を、上記複数の予め定めたθpとしての値の組のうちから選択する。つまり自発脳磁場推定部２４ａは、ＡＩＣの値を最小とするベクトルθpを、まずはグリッドサーチにより大域的に粗く探索する。
【００４１】
自発脳磁場推定部２４ａは、さらにグリッドサーチにより求めたベクトルθp（θp0とする）を初期値として、非線形最適化により最適化された自発脳磁場のモデルのためのパラメータのベクトル
【数１３】
を求める。
【００４２】
ここで非線形最適化の一例としては、上記対数尤度log（ｐ（ｍ′1，ｍ′2，…，ｍ′t｜θp））をカルマンフィルタにより求める方法がある。上記グリッドサーチにより求めたベクトルθp0を初期値として、対数尤度を最大とするベクトルθpを、例えば最急降下法により予め定めた終了条件を満足するまで繰り返し逐次的に更新して求め、当該求めたθpを、
【数１３】
とする。ここで終了条件は、更新前のθpと更新後のθpとの距離が予め定めたしきい値を下回っているなどの広く使われる条件を適用すればよい。
【００４３】
自発脳磁場推定部２４ａは、こうして求めた
【数１３】
を、誘発脳磁場推定部２４ｂに出力する。
【００４４】
誘発脳磁場推定部２４ｂは、記憶部１２に格納されている、ノイズと自発脳磁場とを含んだ刺激時データ（刺激を受けた状態にある被験者を測定した脳磁場のデータ）を読み出す。ここで誘発脳磁場推定部２４ｂは、誘発脳磁場を次の２次の線型モデル
【数１４】
のインパルス応答で近似できるものとする。
【００４５】
本実施の形態の誘発脳磁場推定部２４ｂは、誘発脳磁場を、確率的入力を持つ２次の定常線形モデルとして（１１）式で表現しておく。この（１１）式においてモデルパラメータはシステム入力強度と、システム入力に対するモデルの応答を決めるシステムパラメータｂ1，ｂ2である。
【００４６】
ここでは単一の線型モデルのインパルス応答で誘発脳磁場を近似できるものとして、この刺激（システム入力）の確率分布を次のようにして与える。すなわち、誘発脳磁場推定部２４ｂは、システムパラメータとしてｂ1，ｂ2をもつシステムのインパルス応答ｈの時系列ｈ1，ｈ2，…と、得られている刺激時データｍ1，ｍ2，…との相互相関関数を求める。そして誘発脳磁場推定部２４ｂは、刺激の確率分布が、上記相互相関の絶対値が最大となる時間的ラグｔ（誘発脳磁場が現れる時刻ｔ）を中心とする二項分布
Ｂ（２ｔ，０．５）
であるとする。ここでは時間が離散化されているので二項分布とした。また、入力強度はラグｔでの相互相関関数の値を平均値Iとするガウス分布
Ｎ（Ｉ，σI2）
としておく。
【００４７】
誘発脳磁場推定部２４ｂは、ノイズパラメータ管理部２３に対してノイズパラメータを要求し、この要求に応答してノイズパラメータ管理部２３が出力するノイズパラメータと（５）式とを用いて、ノイズの時系列推定データν1，ν2，…を得る。誘発脳磁場推定部２４ｂは、また自発脳磁場推定部２４ａから入力される自発脳磁場に関するパラメータと（７）式とを用いて、自発脳磁場の時系列推定データｇ1，ｇ2，…を得る。そして誘発脳磁場推定部２４ｂは、刺激時データｍ1，ｍ2，…，ｍtのそれぞれ対応する時刻のデータから、ノイズの時系列推定データの対応する時刻のデータ及び、自発脳磁場の時系列推定データの対応する時刻のデータを差引きして、ｍ″i＝ｍi−νi−ｇiを求め、ノイズ及び自発脳磁場がないと仮定される状態のデータ、ｍ″1，ｍ″2，…，ｍ″tを得る。
【００４８】
（１１）式では、モデルから直接にデータが得られるものではないので、上述のように刺激の確率分布を、ガウス分布を用いた演算式（これには上記パラメータｂ1，ｂ2が含まれる）で表現しておき、この演算式における時刻ｔでの期待値μ（ｔ）と分散σ2とを用いた確率差分方程式により、刺激に対する推定される誘発脳磁場の時刻ｔでの信号ｓt及び時刻ｔ−１での信号ｓt-1を、その成分に含む状態ベクトルｘの時間発展の演算式を得ておく。
【００４９】
この演算式は、時刻ｔ，ｔ−１において推定される誘発脳磁場ｓtから、時刻ｔ＋１において推定される誘発脳磁場ｓt+1が得られる時間発展の確率（運動のモデル）ｐ（ｓt+1｜ｓt，ｓt-1）を表す。
【００５０】
誘発脳磁場推定部２４ｂは、システムパラメータｂ1，ｂ2を用いて表される信号ｓt，ｓt-1を成分として含む状態ベクトルｘtについて、初期状態ｐ（ｘ0）（これは誘発脳磁場が発生していないものとして既知である）と、上記時間発展の確率とを用いて、時刻ｔで観測されるデータがｍ″tであったとき、状態ベクトルがｘtである確率ｐ（ｘt｜ｍ″t）を演算する。本実施の形態の誘発脳磁場推定部２４ｂは、この確率ｐ（ｘt｜ｍ″t）を、次のように求める。
【００５１】
すなわち、ベイズの定理より、
【数１５】
が成立するので、これを変形して、
【数１６】
を得る。ここでｐ（ｍ″t）は、
【数１７】
と表すことができ、またｐ（ｘt）は、時間発展の確率を用いて、
【数１８】
と表すことができるから、結局（１３）式は、
【数１９】
とすることができる。
【００５２】
この（１６）式は、時刻ｔで観測されるデータがｍ″tであったとき、状態ベクトルがｘtである確率密度ｐ（ｘt｜ｍ″t）が、観測のモデルｐ（ｍ″t｜ｘt）と、運動のモデルｐ（ｘt+1｜ｘt）と、初期状態ｐ（ｘ0）とから推定できることを表す。
【００５３】
誘発脳磁場推定部２４ｂは、システムパラメータｂ1，ｂ2を含むベクトルθsが与えられたときに確率ｐ（ｘt｜ｍ″t）を最大にするベクトルθsを求める。具体的に誘発脳磁場推定部２４ｂは、ＫＬ情報量を最小とする赤池情報量基準を用いてベクトルθsを決定する。このため誘発脳磁場推定部２４ｂは、上記確率ｐ（ｘt｜ｍ″t）の対数log（ｐ（ｘt｜ｍ″t））を対数尤度とし、赤池情報量基準ＡＩＣを、
ＡＩＣ＝−２×ｐ（ｘt｜ｍ″t）＋２×Ｎ
として演算する。ここにＮは自由に設定できるモデルパラメータの数を表し、ここではθsの次元である、２となる。
【００５４】
誘発脳磁場推定部２４ｂは、このＡＩＣの値を、複数の予め定めたθsである値の組ごとに演算する。具体的にここで予め定めたθsの値の組は、
ｂ1，ｂ2
で張られるベクトル空間上に設定された複数の格子点とすることができる。そして演算したＡＩＣの値が最小となるθsを、上記複数の予め定めたθsとしての値の組のうちから選択する。つまり誘発脳磁場推定部２４ｂは、ＡＩＣの値を最小とするベクトルθsを、まずはグリッドサーチにより大域的に粗く探索する。
【００５５】
誘発脳磁場推定部２４ｂは、さらにグリッドサーチにより求めたベクトルθs（θs0とする）を初期値として、非線形最適化により最適化されたベクトル
【数２０】
を求める。
【００５６】
ここで非線形最適化の一例としては、上記対数尤度log（ｐ（ｘt｜ｍ″t））をカルマンフィルタにより求める方法がある。上記グリッドサーチにより求めたベクトルθs0を初期値として、対数尤度を最大とするベクトルθsを、例えば最急降下法により予め定めた終了条件を満足するまで繰り返し逐次的に更新して求め、当該求めたθsを、
【数２０】
とする。ここで終了条件は、更新前のθsと更新後のθsとの距離が予め定めたしきい値を下回っているなどの広く使われる条件を適用すればよい。
【００５７】
信号出力部２５は、誘発脳磁場推定部２４ｂが演算により得た
【数２０】
と、（１１）式で表されるモデルの式とを用いて、誘発脳磁場の時間変化波形を演算して出力する。
【００５８】
さらに本実施の形態の制御部１１は、以上の演算を、複数のセンサ（脳磁計の場合はSQUIDs）から得られた各データについて行い、センサごとに
【数２０】
のパラメータを得てもよい。このときには、複数のセンサの一つを注目センサとするとき、注目センサに係るパラメータを得るためには、ノイズのみのデータや安静時データ、刺激時データは、当該注目センサによって得られたものを用いてもよい。
【００５９】
そして制御部１１は、複数のセンサの各々に対応して得られたパラメータと（１１）式とを用い、各センサの位置における誘発脳磁場の時間変化波形を演算し、等磁界線図を作成して出力してもよい。この等磁界線図の作成処理は広く知られているので、ここでの詳細な説明を省略する。
【００６０】
［モンテカルロ法による演算］
ところで制御部１１における（１６）式の演算では、右辺分母の積分（ベイズ推定における、事前確率と事後確率との積を被積分関数とする積分）が解析的な解が得られないので困難になる。そこで本実施の形態の制御部１１は、複数の演算素子を含み、次のように演算を行うこととしてもよい。
【００６１】
ここで複数の演算素子は、ＣＰＵでなく、ＧＰＵ（Graphics Processing Unit）等の演算素子を用いてもよい。このように複数の演算素子を利用できる場合、制御部１１は、各演算素子に対し、積分範囲内からランダムにｘtとｘt-1との値を選択させ、当該ランダムに選択した積分変数の値（複数ある場合はそれぞれ積分範囲からランダムに選択する）を被積分関数に代入して得た被積分関数値を、並列的に演算させる。そして制御部１１は、各演算素子で得られた複数の被積分関数値を累算してモンテカルロ法により積分を遂行する。
【００６２】
［ノイズのみのデータを得るタイミング］
また、以上の説明において、ノイズパラメータの決定に用いた、ノイズのみのデータは、演算を行う時点より過去であれば、いつ得られたものであってもよかった。例えば刺激データを得る直前（刺激データを実際に測定する被験者が入る直前に、測定に先立って得られたデータ）であってもよいし、刺激データを測定した後、被験者が退席してから得られたデータであってもよい。
【００６３】
［自発脳磁場をベイズ推定で求める例］
ここまでの説明において、自発脳磁場推定部２４ａは、（１０）式で表されるベクトルθが与えられたときに、安静時データが、測定された時系列順に観測される確率ｐ（ｍ′1，ｍ′2，…，ｍ′t｜θ）を最大とするベクトルθを求めていた。しかしながら、本実施の形態はこれに限られない。
【００６４】
すなわち、自発脳磁場推定部２４ａは、時刻ｔでの状態ベクトルｘtとして、時刻ｔ以前の自発脳磁場を表すデータｇt，ｇt-1，…と、時刻ｔ以前の各ノイズを表すデータτt，ｑt,ｑt-1, …, ｑt-T+2, ｅt, …, ｅt-l+1，ｎt…とを含めた
xt=[τt，ｑt,ｑt-1, …, ｑt-T+2, ｅt, …, ｅt-l+1，ｎt，ｇt，ｇt-1，…]T
を生成する。なお、肩のＴは、転置を表す。
【００６５】
ここで、各ノイズのデータや自発脳磁場のデータを、時刻ｔよりどれだけ前のデータから含めるかは、ノイズによって異なり、それぞれ時刻ｔ＋１における対応するノイズを予測するに十分な数とする。これは採用するモデルによっても異なり、例えば自発脳磁場でいえば、ｇt一つだけの場合もあり得る。
【００６６】
自発脳磁場推定部２４ａは、誘発脳磁場推定部２４ｂが行う（１６）式の演算と同様の演算により、（１０）式のベクトルθが与えられたときに確率ｐ（ｘt｜ｍ′t）を最大にするベクトルθを求める。具体的に自発脳磁場推定部２４ａは、ＫＬ情報量を最小とする赤池情報量基準を用いてベクトルθを決定する。このため自発脳磁場推定部２４ａは、上記確率ｐ（ｘt｜ｍ′t）の対数log（ｐ（ｘt｜ｍ′t））を対数尤度とし、赤池情報量基準ＡＩＣを、
ＡＩＣ＝−２×ｐ（ｘt｜ｍ′t）＋２×Ｎ
として演算する。ここにＮは自由に設定できるモデルパラメータの数を表し、ここではθの次元である、６となる。
【００６７】
そして自発脳磁場推定部２４ａは、このＡＩＣを最小とするθを、既に述べたグリッドサーチや非線形最適化を組み合わせた方法で見出すこととしてもよい。
【００６８】
この際の（１６）式の演算においても、モンテカルロ法による演算を行うことができる。
【００６９】
［自発脳磁場もベイズ推定に含める例］
さらに、ここまでの説明では、誘発脳磁場推定部２４ｂは、自発脳磁場推定部２４ａから入力される自発脳磁場に関するパラメータと（７）式とを用いて、自発脳磁場の時系列推定データｇ1，ｇ2，…を得て、刺激時データｍ1，ｍ2，…，ｍtのそれぞれ対応する時刻のデータから、ノイズの時系列推定データの対応する時刻のデータ及び、自発脳磁場の時系列推定データの対応する時刻のデータを差引きして、ｍ″i＝ｍi−νi−ｇiを求め、ノイズ及び自発脳磁場がないと仮定される状態のデータ、ｍ″1，ｍ″2，…，ｍ″tを得ていた。
【００７０】
しかしながら本実施の形態はこれに限られない。すなわち誘発脳磁場推定部２４ｂは、刺激時データｍ1，ｍ2，…，ｍtのそれぞれ対応する時刻のデータから、ノイズの時系列推定データの対応する時刻のデータを差引きして、ｍ″i＝ｍi−νiを求め、ノイズがないと仮定される状態のデータ、ｍ″1，ｍ″2，…，ｍ″tを得ておき、次のように処理を遂行してもよい。
【００７１】
この場合、誘発脳磁場推定部２４ｂは、時刻ｔでの状態ベクトルｘtに、時刻ｔ以前の自発脳磁場を表すデータｇt，ｇt-1，…を含めて、
xt=[st,st-1，ｇt，ｇt-1，…]T
とする。なお、肩のＴは、転置を表す。
【００７２】
ここで、自発脳磁場を時刻ｔよりどれだけ前のデータから含めるかは、ｇt+1を予測するに十分な数とする。これは採用するモデルによっても異なり、例えばｇt一つだけの場合もあり得る。
【００７３】
誘発脳磁場推定部２４ｂが行う（１６）式の演算は、このように状態ベクトルの成分が増えても、既に説明した方法と同様に行うことができ、従って誘発脳磁場推定部２４ｂは、既に述べたと同様にして
【数２０】
を求めることができる。
【００７４】
［ノイズもベイズ推定に含める例］
さらに、ここまでの説明では、誘発脳磁場推定部２４ｂは、自発脳磁場の時系列推定データｇ1，ｇ2，…を得て、刺激時データｍ1，ｍ2，…，ｍtのそれぞれ対応する時刻のデータから、ノイズの時系列推定データの対応する時刻のデータ及び、自発脳磁場の時系列推定データの対応する時刻のデータを差引きして、ｍ″i＝ｍi−νi−ｇiを求め、ノイズ及び自発脳磁場がないと仮定される状態のデータ、ｍ″1，ｍ″2，…，ｍ″tを得ていた。
【００７５】
しかしながら本実施の形態はこれに限られない。すなわち誘発脳磁場推定部２４ｂは、刺激時データｍ1，ｍ2，…，ｍtのそれぞれ対応する時刻のデータから、自発脳磁場の時系列推定データｇ1，ｇ2，…の対応する時刻のデータを差引きして、ｍ″i＝ｍi−ｇiを求め、自発脳磁場がないと仮定される状態のデータ、ｍ″1，ｍ″2，…，ｍ″tを得ておき、次のように処理を遂行してもよい。
【００７６】
この場合、誘発脳磁場推定部２４ｂは、時刻ｔでの状態ベクトルｘtに、時刻ｔ以前の各ノイズを表すデータτt，ｑt,ｑt-1, …, ｑt-T+2, ｅt, …, ｅt-l+1，ｎt…を含め、
xt=[st, st-1，τt，ｑt, ｑt-1, …, ｑt-T+2,ｅt, …, ｅt-l+1，ｎt]T
とする。なお、肩のＴは、転置を表す。
【００７７】
ここで、各ノイズのデータを、時刻ｔよりどれだけ前のデータから含めるかは、ノイズによって異なり、それぞれ時刻ｔ＋１における対応するノイズを予測するに十分な数とする。これは採用するモデルによっても異なる。
【００７８】
誘発脳磁場推定部２４ｂが行う（１６）式の演算は、このように状態ベクトルの成分が増えても、既に説明した方法と同様に行うことができ、従って誘発脳磁場推定部２４ｂは、既に述べたと同様にして
【数２０】
を求めることができる。
【００７９】
［ノイズも自発脳磁場もベイズ推定に含める例］
さらに誘発脳磁場推定部２４ｂは、刺激時データｍ1，ｍ2，…，ｍtのそれぞれを、そのまま上記ｍ″1，ｍ″2，…，ｍ″tとして、次のように処理を遂行してもよい。
【００８０】
すなわち誘発脳磁場推定部２４ｂは、時刻ｔでの状態ベクトルｘtに、時刻ｔ以前の各ノイズを表すデータτt，qt,ｑt-1, …, ｑt-T+2, ｅt, …, ｅt-l+1，ｎt…を含め、
xt=[st,st-1，τt，ｑt, ｑt-1, …, ｑt-T+2,ｅt, …, ｅt-l+1，ｎt，ｇt，ｇt-1，…]T
とする。なお、肩のＴは、転置を表す。
【００８１】
ここで、各ノイズのデータを、時刻ｔよりどれだけ前のデータから含めるかは、ノイズによって異なり、それぞれ時刻ｔ＋１における対応するノイズを予測するに十分な数とする。これは採用するモデルによっても異なる。
【００８２】
また誘発脳磁場推定部２４ｂは、時刻ｔ以前の自発脳磁場を表すデータｇt，ｇt-1，…も時刻ｔでの状態ベクトルｘtに含める。ここで、自発脳磁場を時刻ｔよりどれだけ前のデータから含めるかは、ｇt+1を予測するに十分な数とする。これは採用するモデルによっても異なり、例えばｇt一つだけの場合もあり得る。
【００８３】
誘発脳磁場推定部２４ｂが行う（１６）式の演算は、このように状態ベクトルの成分が増えても、既に説明した方法と同様に行うことができ、従って誘発脳磁場推定部２４ｂは、既に述べたと同様にして
【数２０】
を求めることができる。
【００８４】
［グリッドサーチにおけるＡＲモデルの変数変換］
さらにここまでの説明では定常ノイズ成分ｅtの時間発展を表すＡＲモデルにおいては、安定性の条件が見にくいので、次の方法により繰り返し変数変換を行って、１次のＡＲモデルに変換してもよい。
【００８５】
すなわち、Ｌ次のＡＲモデルをＬ−１次のＡＲモデルに変換する方法として、
【数２１】
を用いる。ただし、ｉ＝１，２，…，ｌ−１である。ここにａの肩にある添字のｌは、ｌ次の変数であることを意味する。この（１７）式による変換を、１次となるまで繰り返して行うと、各変換の段階で、ａjjが得られる（ここにｊ＝１，２，…ｌ）。すると、安定なＡＲモデルとなるための必要十分条件が、
【数２２】
と得られる。そこで、グリッドサーチにおいては、この（１８）式で表される範囲を定義域として、この定義域を均等に分割した格子点を用いればよい。
【００８６】
［他のモデル］
なお、一般に活動部位や刺激の種類に応じて、誘発脳磁場モデルの入力の確率分布は異なると考えられるうえ、線形モデルに限らず非線形なモデルが適当な場合が想定される。また本実施の形態は脳磁計のデータに限らず、一般的な測定データの処理に用いることができるものである。そこで利用者がそれぞれの仮説に従い、測定対象となるデータのモデル（ここでの例では誘発脳磁場のモデル）の確率分布を作成してもよい。本実施の形態の処理は、種々の確率分布を仮定した場合でも実行可能なものである。
【００８７】
また、自発脳磁場のモデルについても、ここでは（８）式によって表されるものとしたが、振幅の時間変化が大きい場合には、確率的入力項のパワーが時間的に変化するものと仮定し、当該パワーをパラメータを含む何らかの関数形で記述し、当該パラメータをモデルパラメータに含めて処理を行ってもよい。さらに、被験者の各自発脳磁場帯域のピーク周波数が典型的なものと大きく異なる場合は、中心周波数もモデルパラメータに含めてもよい。これらの場合も、上述の処理と同様にモデルパラメータの一つとして推定を行うことができる。
【００８８】
［動作］
本実施の形態は、以上の構成を有し、次のように動作する。すなわち、本実施の形態のデータ処理装置１は、ノイズのみのデータ、安静時データ、刺激時データを受け入れて記憶部１２に保持しており、図３に例示するように、まずノイズのみのデータを用いて予め定めたノイズのモデルに含まれるノイズパラメータを推定する（Ｓ１）。
【００８９】
またデータ処理装置１は、安静時データに含まれるノイズを、上記推定したノイズパラメータと予め定めたノイズのモデルとを用いて除去し、ノイズを除去した安静時データを参照して、予め定めた安静時データのモデルに含まれるパラメータを推定する（Ｓ２）。
【００９０】
そしてデータ処理装置１は、刺激時データから、処理Ｓ１で推定したノイズパラメータと予め定めたノイズのモデルとを用いて除去し、刺激時データに含まれる自発脳磁場を、処理Ｓ２で推定したパラメータと、予め定めた安静時データのモデルとを用いて除去する。そしてデータ処理装置１は、こうしてノイズ及び自発脳磁場を除去したデータを参照して、予め定めた刺激時データのモデル（誘発脳磁場のモデル）に含まれるパラメータを推定する（Ｓ３）。
【００９１】
データ処理装置１は、処理Ｓ３で推定したパラメータと誘発脳磁場のモデルとを用いて、誘発脳磁場のデータを推定して出力する（Ｓ４）。
【実施例】
【００９２】
本発明の一実施例について説明する。ここでは被験者に、「ボタンを押してもらう」という刺激を与えたときの誘発脳磁場を得たときの例について述べる。
【００９３】
図４は、ＬＰＦ（Low Pass Filter）を経た１６９回の測定データを加算平均した、従来の方法で得られた誘発脳磁場の波形（複数のセンサの信号を重ね合わせて表している、以下も同様）と、これらの波形から得られた等磁界線図である。図５は、ノイズを含めて実際に得られる時系列データの例と、当該時系列データを基に得られた等磁界線図である。この図５のようなデータを１６９回得て加算平均したものが図４に示した従来の結果である。
【００９４】
図６は、この図５に示したデータ（単一測定のデータ）に基づき、上記本実施の形態で説明した処理を経て本発明の実施例に係るデータ処理装置１が出力する誘発脳磁場の波形（複数のセンサの信号を重ね合わせて表している）と、これらの波形から得られた等磁界線図である。
【００９５】
図４に示した従来の解析図と、図６に示した本実施例の解析図とを比較すると、略一致していることが理解される。
【００９６】
図７は、自発脳磁場の強度が大きい場合のノイズを含めて実際に得られる時系列データの例と、当該時系列データを基に得られた等磁界線図である。この図７に例示したようなデータを複数回得て加算平均した、従来の方法によって推定された等磁界線図と、この図７に例示したような複数のデータのそれぞれから、本実施例のデータ処理装置１で推定して得た複数の等磁界線図（図８に例示する）とのコヒーレンスを演算して表したものを図９に示す。また、従来の方法で推定された等磁界線図を図１０に示す。図９，図１０から見られるように、図１０で得られた等磁界線図のうち、信号の強度が大きい領域でのコヒーレンスは、図９に示されるように０．８を超えている。つまり、本実施例のデータ処理装置１での推定は、ロバストネスの面からも十分であることが理解される。
【符号の説明】
【００９７】
１ データ処理装置、１１ 制御部、１２ 記憶部、１３ 操作部、１４ 出力部、１５ 入力インタフェース、２１ データ記録部、２２ モデル学習部、２３ ノイズパラメータ管理部、２４ 信号推定部、２５ 信号出力部。
【技術分野】
【０００１】
本発明は、データ処理装置、及びプログラムに係り、特にノイズの影響除去に関する。
【背景技術】
【０００２】
脳磁計（ＭＥＧ）等、微小信号を計測する装置では、外来の種々のノイズに信号が埋もれてしまうため、ノイズの除去が課題となっている。かかるノイズ除去の一手法として、従来より、複数回の測定結果を加算平均する方法が知られている。また特許文献１には、脳磁場信号に対して、心拍の発生時刻に同期させて加算平均を取ることにより、心磁によるノイズ信号を抽出する技術が開示されている。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特開２００９−１９５５７１号公報
【発明の概要】
【発明が解決しようとする課題】
【０００４】
しかしながら、上記加算平均を用いるノイズの除去方法では、測定回数をなるべく多くしなければならない。ところが測定の対象によっては多数回の測定が困難である場合もある。また脳磁計等の場合、測定回数を増大させるほど被験者の負担が大きくなってしまう。
【０００５】
本発明は上記実情に鑑みて為されたもので、単一の測定結果からでも信号の推定を可能とするデータ処理装置を提供することを、その目的の一つとする。
【課題を解決するための手段】
【０００６】
上記従来例の問題点を解決するための本発明は、データ処理装置であって、測定されたデータを受け入れる受入手段と、少なくとも一つのノイズパラメータを含み、ノイズを模式的に表すノイズ演算式について、過去に得られたデータを用いて定められる基準に基づいて、前記ノイズパラメータを決定する手段と、前記測定されたデータに含まれる信号を模式的に表す信号演算式であって、少なくとも一つの信号パラメータを含む信号演算式と、前記受入手段が受け入れる測定されたデータと、前記ノイズ演算式とに基づいて、前記信号パラメータを推定する推定手段と、前記推定された信号パラメータと、前記信号演算式とを用いて、測定されたデータに含まれる信号の推定結果を生成し、出力する手段と、を含むこととしたものである。
【発明の効果】
【０００７】
本発明によれば、ノイズ及び信号を、パラメータを含んだ模式的な演算式として表現し、そのパラメータを推定するようにしたことで、単一の測定結果からでも信号の推定を可能とした。
【図面の簡単な説明】
【０００８】
【図１】本発明の実施の形態に係るデータ処理装置の例を表す構成ブロック図である。
【図２】本発明の実施の形態に係るデータ処理装置の例を表す機能ブロック図である。
【図３】本発明の実施の形態に係るデータ処理装置による処理例を表すフローチャート図である。
【図４】従来の方法で得られた誘発脳磁場の波形と、これらの波形から得られた等磁界線図の例を表す説明図である。
【図５】ノイズを含めて実際に得られる時系列データの例と、当該時系列データを基に得られた等磁界線図の例を表す説明図である。
【図６】本発明の実施例に係るデータ処理装置が出力する誘発脳磁場の波形と、これらの波形から得られた等磁界線図の例を表す説明図である。
【図７】ノイズを含めて実際に得られる時系列データの別の例と、当該時系列データを基に得られた等磁界線図の別の例を表す説明図である。
【図８】本発明の実施例に係るデータ処理装置で推定して得た等磁界線図の例を表す説明図である。
【図９】本発明の実施例に係るデータ処理装置で推定して得た複数の等磁界線図と、従来の方法で推定された等磁界線図とのコヒーレンスの例を表す説明図である。
【図１０】従来の方法によって推定された等磁界線図の例を表す説明図である。
【発明を実施するための形態】
【０００９】
本発明の実施の形態に係るデータ処理装置を図面を参照しつつ説明する。本実施の形態のデータ処理装置１は、図１に例示するように、制御部１１、記憶部１２、操作部１３、出力部１４、及び入力インタフェース部１５を含んで構成される。以下では、具体的に、脳磁計（ＭＥＧ）の出力するデータ（時間変化する一連のデータ）を入力インタフェース部１５を介して受け入れる例について説明するが、本実施の形態のデータ処理装置１は、脳磁計のデータを処理するものに限られない。
【００１０】
脳磁計のデータには、被験者のいない状態で脳磁計が出力するデータ（ノイズのみのデータ）と、被験者が安静な状態にあるときに脳磁計が出力するデータ（外的刺激を受けていない状態：安静時データと呼ぶ）と、被験者に対して外的刺激を与えたときに脳磁計が出力するデータ（ノイズ、自発脳磁場及び、刺激による誘発脳磁場が累算されたデータ：刺激時データと呼ぶ）とがある。
【００１１】
制御部１１は、ＣＰＵ（Central Processing Unit）等のプログラム制御デバイスであり、記憶部１２に格納されたプログラムに従って動作する。この制御部１１は、入力インタフェース部１５から入力される、外部の測定装置（ここでは例えば脳磁計とする）により測定されたデータを受け入れる。また、この制御部１１は、少なくとも一つのノイズパラメータを含み、ノイズを模式的に表すノイズ演算式について、過去に得られたデータを用いて定められる基準に基づいて、ノイズパラメータを決定する。
【００１２】
さらに制御部１１は、測定されたデータに含まれる信号を模式的に表す信号演算式であって、少なくとも一つの信号パラメータを含む信号演算式と、測定されたデータと、ノイズ演算式とに基づいて、信号パラメータを推定し、当該推定した信号パラメータと、信号演算式とを用いて、測定されたデータに含まれる信号の推定結果を生成し、出力する。この制御部１１の詳しい処理の内容は後に説明する。
【００１３】
記憶部１２は、メモリデバイス等であり、制御部１１によって実行されるプログラムを保持している。このプログラムは、例えばＤＶＤ−ＲＯＭ（Digital Versatile Disc Read Only Memory）等のコンピュータ可読な記録媒体に格納されて提供され、この記憶部１２に複写されたものであってもよい。また、この記憶部１２は制御部１１のワークメモリとしても動作する。
【００１４】
操作部１３は、キーボードやマウスなどを含み、利用者の指示操作を受け入れて、制御部１１に当該指示操作の内容を出力する。出力部１４は、ディスプレイや、プリンタ、ネットワークインタフェース等であり、制御部１１から入力される指示に従って、情報を出力する。
【００１５】
入力インタフェース部１５は、ＵＳＢ（Universal Serial Bus）やパラレルポート等のデータの入力を受け入れるデバイスであり、ここでの例では脳磁計に接続されて、脳磁計が出力するデータを受け入れ、当該受け入れたデータを制御部１１に出力する。
【００１６】
本実施の形態では、脳磁計から入力されるデータは時々刻々と変化する（時間変化する）。そして時刻ｔにおけるデータｍtが、測定の対象たる信号ｓtとノイズとの和により表現されているものとする。すなわち、
【数１】
ここにノイズは、要因ごとに複数の成分の和として表現されており、（１）式ではτtがベースラインドリフト成分、ｑtが交流電源成分、ｅtがベースラインドリフトや交流電源成分以外の定常ノイズ成分、ｎtは、このようにモデル化したことによるデータからの差を表すモデル化誤差成分であるものとしている。また本実施形態では、不定期に生じる測定値のとび（例えばグリッチノイズ）、眼球運動や心拍等の生体から生じるノイズはないものとする。
【００１７】
ノイズの各成分は、具体的に次のような形の演算式で表されるものとする。ベースラインドリフト成分τtの時間発展を、平滑化事前分布を用いた確率的トレンドモデルを採用し、
【数２】
と表す。
【００１８】
ここでvτ ,tはトレンドの変化を起こす外乱であり、平均ゼロ分散στのガウス分布と仮定しておく。なお、この（２）式、及び以下の説明において時刻ｔ＋ｊは、ｔから予め定めた単位時間（例えばＴとする）を用いてｊ×Ｔだけ後の時刻を意味するものである。
【００１９】
さらに（１）式におけるｑtの時間発展を、サンプリング周波数と交流周波数とから決まる自己回帰モデル（ARモデル）を用いて、
【数３】
と表す。
【００２０】
ここでも同様に、vq ,tは平均ゼロ分散σqのガウス分布と仮定する。またΔTはサンプリング周波数と交流電源の周波数から決定される、離散化された周期を表す。
【００２１】
定常ノイズ成分ｅtの時間発展は次数ｌのＡＲモデルを仮定する。すなわち、
【数４】
とする。
【００２２】
これにおいても、ｖe ,tは平均ゼロ分散σeのガウス分布と仮定する。モデル化誤差成分ｎtは、平均ゼロ分散σnのガウス分布と仮定しておく。
【００２３】
本実施の形態の制御部１１は、ノイズ演算式として、（１）から（４）式の和、
【数５】
を用いる。このノイズ演算式に含まれるノイズパラメータを成分とするベクトルθ（モデルパラメータθベクトル）は、従って、
【数６】
となる。
【００２４】
本実施の形態の制御部１１は、機能的には、図２に例示するように、データ記録部２１、モデル学習部２２、ノイズパラメータ管理部２３、信号推定部２４、及び信号出力部２５を含んで構成される。ここにデータ記録部２１は、入力インタフェース部１５を介して受け入れられるデータを記憶部１２に蓄積して保持する。本実施の形態では、このデータ記録部２１は、操作部１３から入力される指示に従い、受け入れたデータを、ノイズのみのデータ、ノイズと自発脳磁場とを含んだデータ、ノイズ、自発脳磁場、及び誘発脳磁場を含んだデータのいずれかに分類して記憶部１２に格納する。
【００２５】
モデル学習部２２は、データ記録部２１によって記憶部１２に蓄積されたデータを参照し、過去に得られたデータ（被験者のいない状態で得られたノイズのみのデータ）を用いて定められる基準に基づいて、（６）式のモデルパラメータθベクトルの成分（各ノイズパラメータ）を決定する。
【００２６】
具体的に、このモデル学習部２２は、過去に得られたノイズのみの時系列のデータｍ1，ｍ2，…，ｍtを記憶部１２から読み出す。そしてベクトルθが与えられたときにデータがこの順に観測される確率ｐ（ｍ1，ｍ2，…，ｍt｜θ）を最大とするベクトルθを求める。すなわち、現実の観測変数ｍtを生成する真の確率モデルは未知であるものの、当該真の確率モデルを、近似的に（５）式で与えられるものとし、真の確率モデルと近似的なモデルとの擬距離であるＫＬ（カルバック・ライブラー）ダイバージェンス（ＫＬ情報量）を最小とする（５）式のパラメータ（ベクトルθの成分）を求めることとするのである。
【００２７】
例えば上記ＫＬ情報量を最小とする赤池情報量基準を用いてベクトルθを決定すればよい。このためモデル学習部２２は、上記確率ｐ（ｍ1，ｍ2，…，ｍt｜θ）の対数log（ｐ（ｍ1，ｍ2，…，ｍt｜θ））を対数尤度とし、赤池情報量基準ＡＩＣを、
ＡＩＣ＝−２×ｐ（ｍ1，ｍ2，…，ｍt｜θ）＋２×Ｎ
として演算する。ここにＮは自由に設定できるモデルパラメータの数を表し、ここではθの次元である、（ｌ＋４）となる。
【００２８】
モデル学習部２２は、このＡＩＣの値を、複数の予め定めたθである値の組ごとに演算する。具体的にここで予め定めたθの値の組は、
ａ1，ａ2，…，ａl，στ，σq，σe，σn
で張られるベクトル空間上に設定された複数の格子点とすることができる。そして演算したＡＩＣの値が最小となるθを、上記複数の予め定めたθである値の組のうちから選択する。つまりモデル学習部２２は、ＡＩＣの値を最小とするベクトルθを、まずはグリッドサーチにより大域的に粗く探索する。
【００２９】
モデル学習部２２は、さらにグリッドサーチにより求めたベクトルθ（θ0とする）を初期値として、非線形最適化により最適化されたベクトル
【数７】
を求める。
【００３０】
ここで非線形最適化の一例としては、上記対数尤度log（ｐ（ｍ1，ｍ2，…，ｍt｜θ））をカルマンフィルタにより求める方法がある。上記グリッドサーチにより求めたベクトルθ0を初期値として、対数尤度を最大とするベクトルθを、例えば最急降下法により予め定めた終了条件を満足するまで繰り返し逐次的に更新して求め、当該求めたθを、
【数７】
とする。ここで終了条件は、更新前のθと更新後のθとの距離が予め定めたしきい値を下回っているなどの広く使われる条件を適用すればよい。
【００３１】
ノイズパラメータ管理部２３は、モデル学習部２２が演算により求めた
【数７】
を、ノイズパラメータとして記憶部１２に格納する。また、このノイズパラメータ管理部２３は、信号推定部２４からノイズパラメータの要求を受けて、記憶部１２に格納したノイズパラメータを読み出し、信号推定部２４に対して出力する。
【００３２】
本実施の形態の信号推定部２４は、自発脳磁場推定部２４ａと、誘発脳磁場推定部２４ｂとを含んで構成される。自発脳磁場推定部２４ａは、記憶部１２に格納されている、ノイズと自発脳磁場とを含んだ安静時データ（安静状態にある被験者を測定した脳磁場のデータ）を読み出す。ここで自発脳磁場推定部２４ａは、自発脳磁場が
【数８】
によって表される確率モデルで近似的に表されるものとする。ここにｒtは時刻ｔの自発脳磁場成分を意味し、δ波(約２Ｈｚ)，θ波(約４Ｈｚ)，α波(約１０Ｈｚ)，β波(約２０Ｈｚ)，low-γ波(約４０Ｈｚ)，high-γ波(約８０Ｈｚ)の成分から構成されるものとし、各成分を、
【数９】
と表したものである。なお、ここでのθは、θ波に対応する自発脳磁場成分であることを意味するものとして便宜的に利用しているもので、パラメータのベクトルθとは異なる。
【００３３】
また、これらの成分の各々はいずれも、次の準周期振動モデルで表現されるものとする。
【００３４】
【数１０】
ここでωは、δ，θ，α，β，γl，γhのそれぞれをとる。また、ｆωは、成分ω（δ，θ，α，β，γl，γhのそれぞれをとる）の中心周波数であり、ｖt(ω)は、平均が０、分散がσω2であるガウス分布に属するとの条件、
【数１１】
を満足する確率的入力である。
【００３５】
以上より、自発脳磁場の成分を表すモデルパラメータであるθベクトルは、
【数１２】
となる。ここにおいて、左辺のθpのθはθベクトルを表すθであり、右辺中、σθ2とあるθは、θ波に対応する自発脳磁場成分を表す便宜的なものである。
【００３６】
自発脳磁場推定部２４ａは、記憶部１２に格納されているノイズと自発脳磁場とを含んだ時系列の安静時データｍ1，ｍ2，…，ｍtを記憶部１２から読み出す。
【００３７】
また自発脳磁場推定部２４ａは、ノイズパラメータ管理部２３に対してノイズパラメータを要求し、この要求に応答してノイズパラメータ管理部２３が出力するノイズパラメータと（５）式とを用いて、ノイズの時系列推定データν1，ν2，…を得る。自発脳磁場推定部２４ａは、ノイズと自発脳磁場とを含んだ時系列の安静時データｍ1，ｍ2，…，ｍtのそれぞれ対応する時刻のデータからノイズの時系列推定データの対応する時刻のデータを差引きして、ｍ′i＝ｍi−νiを求め、ノイズがないと仮定される状態のデータ、ｍ′1，ｍ′2，…，ｍ′tを得る。
【００３８】
そして自発脳磁場推定部２４ａは、（１０）式で表されるベクトルθpが与えられたときに安静時データが実際に測定された時系列順に観測される確率ｐ（ｍ′1，ｍ′2，…，ｍ′t｜θp）を最大とするベクトルθpを求める。この場合も、現実の観測変数ｍ′tを生成する真の確率モデルは未知であるものの、当該真の確率モデルを、近似的に（７）式で与えられるものとし、真の確率モデルと近似的なモデルとの擬距離であるＫＬ（カルバック・ライブラー）ダイバージェンス（ＫＬ情報量）を最小とする（７）式のパラメータ（（１０）式のベクトルθpの成分）を求めることとするのである。
【００３９】
例えば上記ＫＬ情報量を最小とする赤池情報量基準を用いて（１０）式のベクトルθを決定すればよい。このため自発脳磁場推定部２４ａは、上記確率ｐ（ｍ′1，ｍ′2，…，ｍ′t｜θp）の対数log（ｐ（ｍ′1，ｍ′2，…，ｍ′t｜θp））を対数尤度とし、赤池情報量基準ＡＩＣを、
ＡＩＣ＝−２×ｐ（ｍ′1，ｍ′2，…，ｍ′t｜θp）＋２×Ｎ
として演算する。ここにＮは自由に設定できるモデルパラメータの数を表し、ここでは（１０）式のベクトルθpの次元である、６となる。
【００４０】
自発脳磁場推定部２４ａは、このＡＩＣの値を、複数の予め定めたθpである値の組ごとに演算する。具体的にここで予め定めたθpの値の組は、σω2（ただしωはδ，θ，α，β，γl，γhのそれぞれをとる）で張られるベクトル空間上に設定された複数の格子点とすることができる。そして演算したＡＩＣの値が最小となるθの値の組を、上記複数の予め定めたθpとしての値の組のうちから選択する。つまり自発脳磁場推定部２４ａは、ＡＩＣの値を最小とするベクトルθpを、まずはグリッドサーチにより大域的に粗く探索する。
【００４１】
自発脳磁場推定部２４ａは、さらにグリッドサーチにより求めたベクトルθp（θp0とする）を初期値として、非線形最適化により最適化された自発脳磁場のモデルのためのパラメータのベクトル
【数１３】
を求める。
【００４２】
ここで非線形最適化の一例としては、上記対数尤度log（ｐ（ｍ′1，ｍ′2，…，ｍ′t｜θp））をカルマンフィルタにより求める方法がある。上記グリッドサーチにより求めたベクトルθp0を初期値として、対数尤度を最大とするベクトルθpを、例えば最急降下法により予め定めた終了条件を満足するまで繰り返し逐次的に更新して求め、当該求めたθpを、
【数１３】
とする。ここで終了条件は、更新前のθpと更新後のθpとの距離が予め定めたしきい値を下回っているなどの広く使われる条件を適用すればよい。
【００４３】
自発脳磁場推定部２４ａは、こうして求めた
【数１３】
を、誘発脳磁場推定部２４ｂに出力する。
【００４４】
誘発脳磁場推定部２４ｂは、記憶部１２に格納されている、ノイズと自発脳磁場とを含んだ刺激時データ（刺激を受けた状態にある被験者を測定した脳磁場のデータ）を読み出す。ここで誘発脳磁場推定部２４ｂは、誘発脳磁場を次の２次の線型モデル
【数１４】
のインパルス応答で近似できるものとする。
【００４５】
本実施の形態の誘発脳磁場推定部２４ｂは、誘発脳磁場を、確率的入力を持つ２次の定常線形モデルとして（１１）式で表現しておく。この（１１）式においてモデルパラメータはシステム入力強度と、システム入力に対するモデルの応答を決めるシステムパラメータｂ1，ｂ2である。
【００４６】
ここでは単一の線型モデルのインパルス応答で誘発脳磁場を近似できるものとして、この刺激（システム入力）の確率分布を次のようにして与える。すなわち、誘発脳磁場推定部２４ｂは、システムパラメータとしてｂ1，ｂ2をもつシステムのインパルス応答ｈの時系列ｈ1，ｈ2，…と、得られている刺激時データｍ1，ｍ2，…との相互相関関数を求める。そして誘発脳磁場推定部２４ｂは、刺激の確率分布が、上記相互相関の絶対値が最大となる時間的ラグｔ（誘発脳磁場が現れる時刻ｔ）を中心とする二項分布
Ｂ（２ｔ，０．５）
であるとする。ここでは時間が離散化されているので二項分布とした。また、入力強度はラグｔでの相互相関関数の値を平均値Iとするガウス分布
Ｎ（Ｉ，σI2）
としておく。
【００４７】
誘発脳磁場推定部２４ｂは、ノイズパラメータ管理部２３に対してノイズパラメータを要求し、この要求に応答してノイズパラメータ管理部２３が出力するノイズパラメータと（５）式とを用いて、ノイズの時系列推定データν1，ν2，…を得る。誘発脳磁場推定部２４ｂは、また自発脳磁場推定部２４ａから入力される自発脳磁場に関するパラメータと（７）式とを用いて、自発脳磁場の時系列推定データｇ1，ｇ2，…を得る。そして誘発脳磁場推定部２４ｂは、刺激時データｍ1，ｍ2，…，ｍtのそれぞれ対応する時刻のデータから、ノイズの時系列推定データの対応する時刻のデータ及び、自発脳磁場の時系列推定データの対応する時刻のデータを差引きして、ｍ″i＝ｍi−νi−ｇiを求め、ノイズ及び自発脳磁場がないと仮定される状態のデータ、ｍ″1，ｍ″2，…，ｍ″tを得る。
【００４８】
（１１）式では、モデルから直接にデータが得られるものではないので、上述のように刺激の確率分布を、ガウス分布を用いた演算式（これには上記パラメータｂ1，ｂ2が含まれる）で表現しておき、この演算式における時刻ｔでの期待値μ（ｔ）と分散σ2とを用いた確率差分方程式により、刺激に対する推定される誘発脳磁場の時刻ｔでの信号ｓt及び時刻ｔ−１での信号ｓt-1を、その成分に含む状態ベクトルｘの時間発展の演算式を得ておく。
【００４９】
この演算式は、時刻ｔ，ｔ−１において推定される誘発脳磁場ｓtから、時刻ｔ＋１において推定される誘発脳磁場ｓt+1が得られる時間発展の確率（運動のモデル）ｐ（ｓt+1｜ｓt，ｓt-1）を表す。
【００５０】
誘発脳磁場推定部２４ｂは、システムパラメータｂ1，ｂ2を用いて表される信号ｓt，ｓt-1を成分として含む状態ベクトルｘtについて、初期状態ｐ（ｘ0）（これは誘発脳磁場が発生していないものとして既知である）と、上記時間発展の確率とを用いて、時刻ｔで観測されるデータがｍ″tであったとき、状態ベクトルがｘtである確率ｐ（ｘt｜ｍ″t）を演算する。本実施の形態の誘発脳磁場推定部２４ｂは、この確率ｐ（ｘt｜ｍ″t）を、次のように求める。
【００５１】
すなわち、ベイズの定理より、
【数１５】
が成立するので、これを変形して、
【数１６】
を得る。ここでｐ（ｍ″t）は、
【数１７】
と表すことができ、またｐ（ｘt）は、時間発展の確率を用いて、
【数１８】
と表すことができるから、結局（１３）式は、
【数１９】
とすることができる。
【００５２】
この（１６）式は、時刻ｔで観測されるデータがｍ″tであったとき、状態ベクトルがｘtである確率密度ｐ（ｘt｜ｍ″t）が、観測のモデルｐ（ｍ″t｜ｘt）と、運動のモデルｐ（ｘt+1｜ｘt）と、初期状態ｐ（ｘ0）とから推定できることを表す。
【００５３】
誘発脳磁場推定部２４ｂは、システムパラメータｂ1，ｂ2を含むベクトルθsが与えられたときに確率ｐ（ｘt｜ｍ″t）を最大にするベクトルθsを求める。具体的に誘発脳磁場推定部２４ｂは、ＫＬ情報量を最小とする赤池情報量基準を用いてベクトルθsを決定する。このため誘発脳磁場推定部２４ｂは、上記確率ｐ（ｘt｜ｍ″t）の対数log（ｐ（ｘt｜ｍ″t））を対数尤度とし、赤池情報量基準ＡＩＣを、
ＡＩＣ＝−２×ｐ（ｘt｜ｍ″t）＋２×Ｎ
として演算する。ここにＮは自由に設定できるモデルパラメータの数を表し、ここではθsの次元である、２となる。
【００５４】
誘発脳磁場推定部２４ｂは、このＡＩＣの値を、複数の予め定めたθsである値の組ごとに演算する。具体的にここで予め定めたθsの値の組は、
ｂ1，ｂ2
で張られるベクトル空間上に設定された複数の格子点とすることができる。そして演算したＡＩＣの値が最小となるθsを、上記複数の予め定めたθsとしての値の組のうちから選択する。つまり誘発脳磁場推定部２４ｂは、ＡＩＣの値を最小とするベクトルθsを、まずはグリッドサーチにより大域的に粗く探索する。
【００５５】
誘発脳磁場推定部２４ｂは、さらにグリッドサーチにより求めたベクトルθs（θs0とする）を初期値として、非線形最適化により最適化されたベクトル
【数２０】
を求める。
【００５６】
ここで非線形最適化の一例としては、上記対数尤度log（ｐ（ｘt｜ｍ″t））をカルマンフィルタにより求める方法がある。上記グリッドサーチにより求めたベクトルθs0を初期値として、対数尤度を最大とするベクトルθsを、例えば最急降下法により予め定めた終了条件を満足するまで繰り返し逐次的に更新して求め、当該求めたθsを、
【数２０】
とする。ここで終了条件は、更新前のθsと更新後のθsとの距離が予め定めたしきい値を下回っているなどの広く使われる条件を適用すればよい。
【００５７】
信号出力部２５は、誘発脳磁場推定部２４ｂが演算により得た
【数２０】
と、（１１）式で表されるモデルの式とを用いて、誘発脳磁場の時間変化波形を演算して出力する。
【００５８】
さらに本実施の形態の制御部１１は、以上の演算を、複数のセンサ（脳磁計の場合はSQUIDs）から得られた各データについて行い、センサごとに
【数２０】
のパラメータを得てもよい。このときには、複数のセンサの一つを注目センサとするとき、注目センサに係るパラメータを得るためには、ノイズのみのデータや安静時データ、刺激時データは、当該注目センサによって得られたものを用いてもよい。
【００５９】
そして制御部１１は、複数のセンサの各々に対応して得られたパラメータと（１１）式とを用い、各センサの位置における誘発脳磁場の時間変化波形を演算し、等磁界線図を作成して出力してもよい。この等磁界線図の作成処理は広く知られているので、ここでの詳細な説明を省略する。
【００６０】
［モンテカルロ法による演算］
ところで制御部１１における（１６）式の演算では、右辺分母の積分（ベイズ推定における、事前確率と事後確率との積を被積分関数とする積分）が解析的な解が得られないので困難になる。そこで本実施の形態の制御部１１は、複数の演算素子を含み、次のように演算を行うこととしてもよい。
【００６１】
ここで複数の演算素子は、ＣＰＵでなく、ＧＰＵ（Graphics Processing Unit）等の演算素子を用いてもよい。このように複数の演算素子を利用できる場合、制御部１１は、各演算素子に対し、積分範囲内からランダムにｘtとｘt-1との値を選択させ、当該ランダムに選択した積分変数の値（複数ある場合はそれぞれ積分範囲からランダムに選択する）を被積分関数に代入して得た被積分関数値を、並列的に演算させる。そして制御部１１は、各演算素子で得られた複数の被積分関数値を累算してモンテカルロ法により積分を遂行する。
【００６２】
［ノイズのみのデータを得るタイミング］
また、以上の説明において、ノイズパラメータの決定に用いた、ノイズのみのデータは、演算を行う時点より過去であれば、いつ得られたものであってもよかった。例えば刺激データを得る直前（刺激データを実際に測定する被験者が入る直前に、測定に先立って得られたデータ）であってもよいし、刺激データを測定した後、被験者が退席してから得られたデータであってもよい。
【００６３】
［自発脳磁場をベイズ推定で求める例］
ここまでの説明において、自発脳磁場推定部２４ａは、（１０）式で表されるベクトルθが与えられたときに、安静時データが、測定された時系列順に観測される確率ｐ（ｍ′1，ｍ′2，…，ｍ′t｜θ）を最大とするベクトルθを求めていた。しかしながら、本実施の形態はこれに限られない。
【００６４】
すなわち、自発脳磁場推定部２４ａは、時刻ｔでの状態ベクトルｘtとして、時刻ｔ以前の自発脳磁場を表すデータｇt，ｇt-1，…と、時刻ｔ以前の各ノイズを表すデータτt，ｑt,ｑt-1, …, ｑt-T+2, ｅt, …, ｅt-l+1，ｎt…とを含めた
xt=[τt，ｑt,ｑt-1, …, ｑt-T+2, ｅt, …, ｅt-l+1，ｎt，ｇt，ｇt-1，…]T
を生成する。なお、肩のＴは、転置を表す。
【００６５】
ここで、各ノイズのデータや自発脳磁場のデータを、時刻ｔよりどれだけ前のデータから含めるかは、ノイズによって異なり、それぞれ時刻ｔ＋１における対応するノイズを予測するに十分な数とする。これは採用するモデルによっても異なり、例えば自発脳磁場でいえば、ｇt一つだけの場合もあり得る。
【００６６】
自発脳磁場推定部２４ａは、誘発脳磁場推定部２４ｂが行う（１６）式の演算と同様の演算により、（１０）式のベクトルθが与えられたときに確率ｐ（ｘt｜ｍ′t）を最大にするベクトルθを求める。具体的に自発脳磁場推定部２４ａは、ＫＬ情報量を最小とする赤池情報量基準を用いてベクトルθを決定する。このため自発脳磁場推定部２４ａは、上記確率ｐ（ｘt｜ｍ′t）の対数log（ｐ（ｘt｜ｍ′t））を対数尤度とし、赤池情報量基準ＡＩＣを、
ＡＩＣ＝−２×ｐ（ｘt｜ｍ′t）＋２×Ｎ
として演算する。ここにＮは自由に設定できるモデルパラメータの数を表し、ここではθの次元である、６となる。
【００６７】
そして自発脳磁場推定部２４ａは、このＡＩＣを最小とするθを、既に述べたグリッドサーチや非線形最適化を組み合わせた方法で見出すこととしてもよい。
【００６８】
この際の（１６）式の演算においても、モンテカルロ法による演算を行うことができる。
【００６９】
［自発脳磁場もベイズ推定に含める例］
さらに、ここまでの説明では、誘発脳磁場推定部２４ｂは、自発脳磁場推定部２４ａから入力される自発脳磁場に関するパラメータと（７）式とを用いて、自発脳磁場の時系列推定データｇ1，ｇ2，…を得て、刺激時データｍ1，ｍ2，…，ｍtのそれぞれ対応する時刻のデータから、ノイズの時系列推定データの対応する時刻のデータ及び、自発脳磁場の時系列推定データの対応する時刻のデータを差引きして、ｍ″i＝ｍi−νi−ｇiを求め、ノイズ及び自発脳磁場がないと仮定される状態のデータ、ｍ″1，ｍ″2，…，ｍ″tを得ていた。
【００７０】
しかしながら本実施の形態はこれに限られない。すなわち誘発脳磁場推定部２４ｂは、刺激時データｍ1，ｍ2，…，ｍtのそれぞれ対応する時刻のデータから、ノイズの時系列推定データの対応する時刻のデータを差引きして、ｍ″i＝ｍi−νiを求め、ノイズがないと仮定される状態のデータ、ｍ″1，ｍ″2，…，ｍ″tを得ておき、次のように処理を遂行してもよい。
【００７１】
この場合、誘発脳磁場推定部２４ｂは、時刻ｔでの状態ベクトルｘtに、時刻ｔ以前の自発脳磁場を表すデータｇt，ｇt-1，…を含めて、
xt=[st,st-1，ｇt，ｇt-1，…]T
とする。なお、肩のＴは、転置を表す。
【００７２】
ここで、自発脳磁場を時刻ｔよりどれだけ前のデータから含めるかは、ｇt+1を予測するに十分な数とする。これは採用するモデルによっても異なり、例えばｇt一つだけの場合もあり得る。
【００７３】
誘発脳磁場推定部２４ｂが行う（１６）式の演算は、このように状態ベクトルの成分が増えても、既に説明した方法と同様に行うことができ、従って誘発脳磁場推定部２４ｂは、既に述べたと同様にして
【数２０】
を求めることができる。
【００７４】
［ノイズもベイズ推定に含める例］
さらに、ここまでの説明では、誘発脳磁場推定部２４ｂは、自発脳磁場の時系列推定データｇ1，ｇ2，…を得て、刺激時データｍ1，ｍ2，…，ｍtのそれぞれ対応する時刻のデータから、ノイズの時系列推定データの対応する時刻のデータ及び、自発脳磁場の時系列推定データの対応する時刻のデータを差引きして、ｍ″i＝ｍi−νi−ｇiを求め、ノイズ及び自発脳磁場がないと仮定される状態のデータ、ｍ″1，ｍ″2，…，ｍ″tを得ていた。
【００７５】
しかしながら本実施の形態はこれに限られない。すなわち誘発脳磁場推定部２４ｂは、刺激時データｍ1，ｍ2，…，ｍtのそれぞれ対応する時刻のデータから、自発脳磁場の時系列推定データｇ1，ｇ2，…の対応する時刻のデータを差引きして、ｍ″i＝ｍi−ｇiを求め、自発脳磁場がないと仮定される状態のデータ、ｍ″1，ｍ″2，…，ｍ″tを得ておき、次のように処理を遂行してもよい。
【００７６】
この場合、誘発脳磁場推定部２４ｂは、時刻ｔでの状態ベクトルｘtに、時刻ｔ以前の各ノイズを表すデータτt，ｑt,ｑt-1, …, ｑt-T+2, ｅt, …, ｅt-l+1，ｎt…を含め、
xt=[st, st-1，τt，ｑt, ｑt-1, …, ｑt-T+2,ｅt, …, ｅt-l+1，ｎt]T
とする。なお、肩のＴは、転置を表す。
【００７７】
ここで、各ノイズのデータを、時刻ｔよりどれだけ前のデータから含めるかは、ノイズによって異なり、それぞれ時刻ｔ＋１における対応するノイズを予測するに十分な数とする。これは採用するモデルによっても異なる。
【００７８】
誘発脳磁場推定部２４ｂが行う（１６）式の演算は、このように状態ベクトルの成分が増えても、既に説明した方法と同様に行うことができ、従って誘発脳磁場推定部２４ｂは、既に述べたと同様にして
【数２０】
を求めることができる。
【００７９】
［ノイズも自発脳磁場もベイズ推定に含める例］
さらに誘発脳磁場推定部２４ｂは、刺激時データｍ1，ｍ2，…，ｍtのそれぞれを、そのまま上記ｍ″1，ｍ″2，…，ｍ″tとして、次のように処理を遂行してもよい。
【００８０】
すなわち誘発脳磁場推定部２４ｂは、時刻ｔでの状態ベクトルｘtに、時刻ｔ以前の各ノイズを表すデータτt，qt,ｑt-1, …, ｑt-T+2, ｅt, …, ｅt-l+1，ｎt…を含め、
xt=[st,st-1，τt，ｑt, ｑt-1, …, ｑt-T+2,ｅt, …, ｅt-l+1，ｎt，ｇt，ｇt-1，…]T
とする。なお、肩のＴは、転置を表す。
【００８１】
ここで、各ノイズのデータを、時刻ｔよりどれだけ前のデータから含めるかは、ノイズによって異なり、それぞれ時刻ｔ＋１における対応するノイズを予測するに十分な数とする。これは採用するモデルによっても異なる。
【００８２】
また誘発脳磁場推定部２４ｂは、時刻ｔ以前の自発脳磁場を表すデータｇt，ｇt-1，…も時刻ｔでの状態ベクトルｘtに含める。ここで、自発脳磁場を時刻ｔよりどれだけ前のデータから含めるかは、ｇt+1を予測するに十分な数とする。これは採用するモデルによっても異なり、例えばｇt一つだけの場合もあり得る。
【００８３】
誘発脳磁場推定部２４ｂが行う（１６）式の演算は、このように状態ベクトルの成分が増えても、既に説明した方法と同様に行うことができ、従って誘発脳磁場推定部２４ｂは、既に述べたと同様にして
【数２０】
を求めることができる。
【００８４】
［グリッドサーチにおけるＡＲモデルの変数変換］
さらにここまでの説明では定常ノイズ成分ｅtの時間発展を表すＡＲモデルにおいては、安定性の条件が見にくいので、次の方法により繰り返し変数変換を行って、１次のＡＲモデルに変換してもよい。
【００８５】
すなわち、Ｌ次のＡＲモデルをＬ−１次のＡＲモデルに変換する方法として、
【数２１】
を用いる。ただし、ｉ＝１，２，…，ｌ−１である。ここにａの肩にある添字のｌは、ｌ次の変数であることを意味する。この（１７）式による変換を、１次となるまで繰り返して行うと、各変換の段階で、ａjjが得られる（ここにｊ＝１，２，…ｌ）。すると、安定なＡＲモデルとなるための必要十分条件が、
【数２２】
と得られる。そこで、グリッドサーチにおいては、この（１８）式で表される範囲を定義域として、この定義域を均等に分割した格子点を用いればよい。
【００８６】
［他のモデル］
なお、一般に活動部位や刺激の種類に応じて、誘発脳磁場モデルの入力の確率分布は異なると考えられるうえ、線形モデルに限らず非線形なモデルが適当な場合が想定される。また本実施の形態は脳磁計のデータに限らず、一般的な測定データの処理に用いることができるものである。そこで利用者がそれぞれの仮説に従い、測定対象となるデータのモデル（ここでの例では誘発脳磁場のモデル）の確率分布を作成してもよい。本実施の形態の処理は、種々の確率分布を仮定した場合でも実行可能なものである。
【００８７】
また、自発脳磁場のモデルについても、ここでは（８）式によって表されるものとしたが、振幅の時間変化が大きい場合には、確率的入力項のパワーが時間的に変化するものと仮定し、当該パワーをパラメータを含む何らかの関数形で記述し、当該パラメータをモデルパラメータに含めて処理を行ってもよい。さらに、被験者の各自発脳磁場帯域のピーク周波数が典型的なものと大きく異なる場合は、中心周波数もモデルパラメータに含めてもよい。これらの場合も、上述の処理と同様にモデルパラメータの一つとして推定を行うことができる。
【００８８】
［動作］
本実施の形態は、以上の構成を有し、次のように動作する。すなわち、本実施の形態のデータ処理装置１は、ノイズのみのデータ、安静時データ、刺激時データを受け入れて記憶部１２に保持しており、図３に例示するように、まずノイズのみのデータを用いて予め定めたノイズのモデルに含まれるノイズパラメータを推定する（Ｓ１）。
【００８９】
またデータ処理装置１は、安静時データに含まれるノイズを、上記推定したノイズパラメータと予め定めたノイズのモデルとを用いて除去し、ノイズを除去した安静時データを参照して、予め定めた安静時データのモデルに含まれるパラメータを推定する（Ｓ２）。
【００９０】
そしてデータ処理装置１は、刺激時データから、処理Ｓ１で推定したノイズパラメータと予め定めたノイズのモデルとを用いて除去し、刺激時データに含まれる自発脳磁場を、処理Ｓ２で推定したパラメータと、予め定めた安静時データのモデルとを用いて除去する。そしてデータ処理装置１は、こうしてノイズ及び自発脳磁場を除去したデータを参照して、予め定めた刺激時データのモデル（誘発脳磁場のモデル）に含まれるパラメータを推定する（Ｓ３）。
【００９１】
データ処理装置１は、処理Ｓ３で推定したパラメータと誘発脳磁場のモデルとを用いて、誘発脳磁場のデータを推定して出力する（Ｓ４）。
【実施例】
【００９２】
本発明の一実施例について説明する。ここでは被験者に、「ボタンを押してもらう」という刺激を与えたときの誘発脳磁場を得たときの例について述べる。
【００９３】
図４は、ＬＰＦ（Low Pass Filter）を経た１６９回の測定データを加算平均した、従来の方法で得られた誘発脳磁場の波形（複数のセンサの信号を重ね合わせて表している、以下も同様）と、これらの波形から得られた等磁界線図である。図５は、ノイズを含めて実際に得られる時系列データの例と、当該時系列データを基に得られた等磁界線図である。この図５のようなデータを１６９回得て加算平均したものが図４に示した従来の結果である。
【００９４】
図６は、この図５に示したデータ（単一測定のデータ）に基づき、上記本実施の形態で説明した処理を経て本発明の実施例に係るデータ処理装置１が出力する誘発脳磁場の波形（複数のセンサの信号を重ね合わせて表している）と、これらの波形から得られた等磁界線図である。
【００９５】
図４に示した従来の解析図と、図６に示した本実施例の解析図とを比較すると、略一致していることが理解される。
【００９６】
図７は、自発脳磁場の強度が大きい場合のノイズを含めて実際に得られる時系列データの例と、当該時系列データを基に得られた等磁界線図である。この図７に例示したようなデータを複数回得て加算平均した、従来の方法によって推定された等磁界線図と、この図７に例示したような複数のデータのそれぞれから、本実施例のデータ処理装置１で推定して得た複数の等磁界線図（図８に例示する）とのコヒーレンスを演算して表したものを図９に示す。また、従来の方法で推定された等磁界線図を図１０に示す。図９，図１０から見られるように、図１０で得られた等磁界線図のうち、信号の強度が大きい領域でのコヒーレンスは、図９に示されるように０．８を超えている。つまり、本実施例のデータ処理装置１での推定は、ロバストネスの面からも十分であることが理解される。
【符号の説明】
【００９７】
１ データ処理装置、１１ 制御部、１２ 記憶部、１３ 操作部、１４ 出力部、１５ 入力インタフェース、２１ データ記録部、２２ モデル学習部、２３ ノイズパラメータ管理部、２４ 信号推定部、２５ 信号出力部。
【特許請求の範囲】
【請求項１】
測定されたデータを受け入れる受入手段と、
少なくとも一つのノイズパラメータを含み、ノイズを模式的に表すノイズ演算式について、過去に得られたデータを用いて定められる基準に基づいて、前記ノイズパラメータを決定する手段と、
前記測定されたデータに含まれる信号を模式的に表す信号演算式であって、少なくとも一つの信号パラメータを含む信号演算式と、前記受入手段が受け入れる測定されたデータと、前記ノイズ演算式とに基づいて、前記信号パラメータを推定する推定手段と、
前記推定された信号パラメータと、前記信号演算式とを用いて、測定されたデータに含まれる信号の推定結果を生成し、出力する手段と、
を含むデータ処理装置。
【請求項２】
請求項１記載のデータ処理装置であって、
前記受入手段は、測定されたデータとして時間変化するデータを受け入れ、
前記推定手段は、時刻ｔにおけるデータの値ｍtから、時刻ｔの信号の値を含む情報の確率密度を与える演算式として、当該情報がｘtのときに測定されるデータの値ｍtの確率密度である事後確率と、前記信号の値を含む情報の初期値ｘ0、及び当該情報の確率密度の時間発展を表す模式的な演算式の積により表現される事前確率とを用いたベイズ推定により前記信号パラメータを推定する推定手段であるデータ処理装置。
【請求項３】
請求項２記載のデータ処理装置であって、
複数の演算装置要素を含み、
前記推定手段は、前記ベイズ推定における、事前確率と事後確率との積を被積分関数とする積分を行うにあたり、積分範囲内からランダムに選択される積分変数の値を代入したときの前記被積分関数の値の演算を、複数の演算装置要素に分散的に行わせ、モンテカルロ法により前記積分を遂行するデータ処理装置。
【請求項４】
請求項１から３のいずれか一項に記載のデータ処理装置であって、
前記ノイズパラメータを決定する手段は、前記ノイズ演算式に含まれる少なくとも一つのノイズパラメータについては、データの実際の測定に先立って、前記受入手段が受け入れるデータを用いて定められる基準に基づいて決定するデータ処理装置。
【請求項５】
コンピュータを、
測定されたデータを受け入れる受入手段と、
少なくとも一つのノイズパラメータを含み、ノイズを模式的に表すノイズ演算式について、過去に得られたデータを用いて定められる基準に基づいて、前記ノイズパラメータを決定する手段と、
前記測定されたデータに含まれる信号を模式的に表す信号演算式であって、少なくとも一つの信号パラメータを含む信号演算式と、前記受入手段が受け入れる測定されたデータと、前記ノイズ演算式とに基づいて、前記信号パラメータを推定する推定手段と、
前記推定された信号パラメータと、前記信号演算式とを用いて、測定されたデータに含まれる信号の推定結果を生成し、出力する手段と、
として機能させるプログラム。
【請求項１】
測定されたデータを受け入れる受入手段と、
少なくとも一つのノイズパラメータを含み、ノイズを模式的に表すノイズ演算式について、過去に得られたデータを用いて定められる基準に基づいて、前記ノイズパラメータを決定する手段と、
前記測定されたデータに含まれる信号を模式的に表す信号演算式であって、少なくとも一つの信号パラメータを含む信号演算式と、前記受入手段が受け入れる測定されたデータと、前記ノイズ演算式とに基づいて、前記信号パラメータを推定する推定手段と、
前記推定された信号パラメータと、前記信号演算式とを用いて、測定されたデータに含まれる信号の推定結果を生成し、出力する手段と、
を含むデータ処理装置。
【請求項２】
請求項１記載のデータ処理装置であって、
前記受入手段は、測定されたデータとして時間変化するデータを受け入れ、
前記推定手段は、時刻ｔにおけるデータの値ｍtから、時刻ｔの信号の値を含む情報の確率密度を与える演算式として、当該情報がｘtのときに測定されるデータの値ｍtの確率密度である事後確率と、前記信号の値を含む情報の初期値ｘ0、及び当該情報の確率密度の時間発展を表す模式的な演算式の積により表現される事前確率とを用いたベイズ推定により前記信号パラメータを推定する推定手段であるデータ処理装置。
【請求項３】
請求項２記載のデータ処理装置であって、
複数の演算装置要素を含み、
前記推定手段は、前記ベイズ推定における、事前確率と事後確率との積を被積分関数とする積分を行うにあたり、積分範囲内からランダムに選択される積分変数の値を代入したときの前記被積分関数の値の演算を、複数の演算装置要素に分散的に行わせ、モンテカルロ法により前記積分を遂行するデータ処理装置。
【請求項４】
請求項１から３のいずれか一項に記載のデータ処理装置であって、
前記ノイズパラメータを決定する手段は、前記ノイズ演算式に含まれる少なくとも一つのノイズパラメータについては、データの実際の測定に先立って、前記受入手段が受け入れるデータを用いて定められる基準に基づいて決定するデータ処理装置。
【請求項５】
コンピュータを、
測定されたデータを受け入れる受入手段と、
少なくとも一つのノイズパラメータを含み、ノイズを模式的に表すノイズ演算式について、過去に得られたデータを用いて定められる基準に基づいて、前記ノイズパラメータを決定する手段と、
前記測定されたデータに含まれる信号を模式的に表す信号演算式であって、少なくとも一つの信号パラメータを含む信号演算式と、前記受入手段が受け入れる測定されたデータと、前記ノイズ演算式とに基づいて、前記信号パラメータを推定する推定手段と、
前記推定された信号パラメータと、前記信号演算式とを用いて、測定されたデータに含まれる信号の推定結果を生成し、出力する手段と、
として機能させるプログラム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【公開番号】特開２０１２−２４９８１７（Ｐ２０１２−２４９８１７Ａ）
【公開日】平成２４年１２月２０日（２０１２．１２．２０）
【国際特許分類】
【出願番号】特願２０１１−１２４６８８（Ｐ２０１１−１２４６８８）
【出願日】平成２３年６月２日（２０１１．６．２）
【出願人】（５０４１３７９１２）国立大学法人 東京大学 (1,942)
【Ｆターム（参考）】
【公開日】平成２４年１２月２０日（２０１２．１２．２０）
【国際特許分類】
【出願日】平成２３年６月２日（２０１１．６．２）
【出願人】（５０４１３７９１２）国立大学法人 東京大学 (1,942)
【Ｆターム（参考）】
[ Back to top ]