JP5852331B2 - Data processing apparatus and program - Google Patents

Description

  The present invention relates to a data processing apparatus and a program, and more particularly to noise removal.

  In devices that measure minute signals, such as a magnetoencephalograph (MEG), the signals are buried in various external noises, and thus noise removal is an issue. As a method for removing such noise, a method of averaging a plurality of measurement results has been conventionally known. Patent Document 1 discloses a technique for extracting a noise signal due to a magnetocardiogram by taking an addition average in synchronism with the generation time of a heartbeat with respect to a cerebral magnetic field signal.

JP 2009-195571 A

  However, in the noise removal method using the above-described average, the number of measurements must be increased as much as possible. However, depending on the object of measurement, it may be difficult to perform multiple measurements. In the case of a magnetoencephalograph or the like, the burden on the subject increases as the number of measurements increases.

  The present invention has been made in view of the above circumstances, and an object of the present invention is to provide a data processing apparatus that can estimate a signal even from a single measurement result.

  The present invention for solving the problems of the above-described conventional example is a data processing apparatus, including a receiving means for receiving measured data, and a noise calculation expression that includes at least one noise parameter and schematically represents noise. A means for determining the noise parameter based on a criterion determined using data obtained in the past, and a signal arithmetic expression schematically representing a signal included in the measured data, wherein at least one An estimation means for estimating the signal parameter based on a signal calculation formula including a signal parameter, measured data received by the receiving means, and the noise calculation formula, the estimated signal parameter, and the signal calculation And a means for generating and outputting an estimation result of a signal included in the measured data using the equation.

  According to the present invention, noise and a signal are expressed as a schematic arithmetic expression including a parameter, and the parameter is estimated, so that the signal can be estimated even from a single measurement result.

It is a block diagram showing an example of a data processing apparatus according to an embodiment of the present invention. It is a functional block diagram showing the example of the data processor which concerns on embodiment of this invention. It is a flowchart figure showing the example of a process by the data processor which concerns on embodiment of this invention. It is explanatory drawing showing the example of the waveform of the induced brain magnetic field obtained by the conventional method, and the isomagnetic-field diagram obtained from these waveforms. It is explanatory drawing showing the example of the time series data actually obtained including noise and the example of the isomagnetic field map obtained based on the said time series data. It is explanatory drawing showing the example of the waveform of the induced brain magnetic field which the data processor which concerns on the Example of this invention outputs, and the isomagnetic-field diagram obtained from these waveforms. It is explanatory drawing showing another example of the time series data actually obtained including a noise, and another example of the isomagnetic field diagram obtained based on the said time series data. It is explanatory drawing showing the example of the isomagnetic-field diagram obtained by estimating with the data processor which concerns on the Example of this invention. It is explanatory drawing showing the example of the coherence of the several isomagnetic field map obtained by estimating with the data processor which concerns on the Example of this invention, and the isomagnetic field map estimated by the conventional method. It is explanatory drawing showing the example of the isomagnetic-field diagram estimated by the conventional method.

  A data processing apparatus according to an embodiment of the present invention will be described with reference to the drawings. As illustrated in FIG. 1, the data processing apparatus 1 according to the present embodiment includes a control unit 11, a storage unit 12, an operation unit 13, an output unit 14, and an input interface unit 15. Hereinafter, an example in which data output from the magnetoencephalograph (MEG) (a series of time-varying data) is received via the input interface unit 15 will be described. However, the data processing apparatus 1 of the present embodiment is It is not limited to processing magnetoencephalographic data.

  The magnetoencephalographic data includes the data output by the magnetoencephalograph in the absence of the subject (noise-only data) and the data output by the magnetoencephalograph when the subject is in a resting state (in the absence of external stimulation) : Data called resting data) and data output by the magnetoencephalograph when an external stimulus is given to the subject (data in which noise, spontaneous brain magnetic field, and brain magnetic field induced by stimulation are accumulated: data during stimulation Called).

  The control unit 11 is a program control device such as a CPU (Central Processing Unit) and operates according to a program stored in the storage unit 12. The control unit 11 receives data measured by an external measuring device (here, for example, a magnetoencephalograph) input from the input interface unit 15. Moreover, this control part 11 determines a noise parameter based on the reference | standard defined using the data obtained in the past about the noise arithmetic expression which contains at least 1 noise parameter and represents noise typically.

  Further, the control unit 11 is a signal arithmetic expression that schematically represents a signal included in the measured data, and is based on the signal arithmetic expression including at least one signal parameter, the measured data, and the noise arithmetic expression. Then, the signal parameter is estimated, and the estimation result of the signal included in the measured data is generated and output using the estimated signal parameter and the signal arithmetic expression. Details of the processing of the control unit 11 will be described later.

  The storage unit 12 is a memory device or the like, and holds a program executed by the control unit 11. This program may be provided by being stored in a computer-readable recording medium such as a DVD-ROM (Digital Versatile Disc Read Only Memory) and copied to the storage unit 12. The storage unit 12 also operates as a work memory for the control unit 11.

  The operation unit 13 includes a keyboard, a mouse, and the like, accepts a user's instruction operation, and outputs the contents of the instruction operation to the control unit 11. The output unit 14 is a display, a printer, a network interface, or the like, and outputs information according to an instruction input from the control unit 11.

  The input interface unit 15 is a device that accepts data input, such as a USB (Universal Serial Bus) or a parallel port. In this example, the input interface unit 15 is connected to the magnetoencephalograph, accepts the data output by the magnetoencephalograph, and accepts the data. Data is output to the control unit 11.

ここにノイズは、要因ごとに複数の成分の和として表現されており、（１）式ではτ_tがベースラインドリフト成分、ｑ_tが交流電源成分、ｅ_tがベースラインドリフトや交流電源成分以外の定常ノイズ成分、ｎ_tは、このようにモデル化したことによるデータからの差を表すモデル化誤差成分であるものとしている。また本実施形態では、不定期に生じる測定値のとび（例えばグリッチノイズ）、眼球運動や心拍等の生体から生じるノイズはないものとする。 In the present embodiment, data input from the magnetoencephalograph changes every moment (time changes). The data m _t at time t is assumed to be represented by the sum of the target serving signal s _t and noise measurements. That is,

Here, the noise is expressed as the sum of a plurality of components for each factor. In equation (1), τ _t is the baseline drift component, q _t is the AC power supply component, and _et is other than the baseline drift or AC power supply component. The stationary noise component, n _t, is a modeling error component that represents a difference from the data resulting from modeling in this way. Further, in the present embodiment, it is assumed that there is no noise generated from a living body such as skips of measurement values that occur irregularly (for example, glitch noise), eye movements, and heartbeats.

と表す。 Each component of noise is specifically expressed by an arithmetic expression of the following form. Probabilistic trend model using smoothed prior distribution for time evolution of baseline drift component τ _t ,

It expresses.

Here, v _{τ, t} is a disturbance that causes a trend change, and is assumed to be a Gaussian distribution with an average zero variance σ _τ . In the expression (2) and the following description, time t + j means time after j × T using a predetermined unit time (eg, T) from t.

と表す。 Furthermore, the time evolution of q _t in equation (1) is calculated using an autoregressive model (AR model) determined from the sampling frequency and the AC frequency.

It expresses.

Here again, v _{q, t} is assumed to be a Gaussian distribution with mean zero variance σ _q . ΔT represents a discretized period determined from the sampling frequency and the frequency of the AC power supply.

とする。 Time evolution of stationary noise components e _t assumes AR model of order l. That is,

And

Again, it is assumed that v _{e, t} is a Gaussian distribution with mean zero variance σ _e . The modeling error component n _t is assumed to be a Gaussian distribution with a mean zero variance σ _n .

を用いる。このノイズ演算式に含まれるノイズパラメータを成分とするベクトルθ（モデルパラメータθベクトル）は、従って、

となる。 The control unit 11 of the present embodiment uses the sum of the expressions (1) to (4) as a noise calculation expression,

Is used. Therefore, the vector θ (model parameter θ vector) having the noise parameter as a component included in the noise calculation formula is

It becomes.

  The control unit 11 of the present embodiment functionally includes a data recording unit 21, a model learning unit 22, a noise parameter management unit 23, a signal estimation unit 24, and a signal output unit 25 as illustrated in FIG. Consists of including. Here, the data recording unit 21 accumulates and holds data received via the input interface unit 15 in the storage unit 12. In the present embodiment, the data recording unit 21, in accordance with an instruction input from the operation unit 13, converts received data into noise-only data, data including noise and a spontaneous brain magnetic field, noise, a spontaneous brain magnetic field, Then, the data is classified into any of the data including the induced cerebral magnetic field and stored in the storage unit 12.

  The model learning unit 22 refers to the data stored in the storage unit 12 by the data recording unit 21, and uses the data obtained in the past (data only for noise obtained in the absence of the subject) as a standard. Based on this, the component (each noise parameter) of the model parameter θ vector of equation (6) is determined.

Specifically, the model learning unit 22, the data m ₁ time series of only noise obtained in the past, m _2, ..., reads the m _t from the storage unit 12. Then, the vector θ that maximizes the probability p (m ₁ , m ₂ ,..., M _t | θ) that the data is observed in this order when the vector θ is given is obtained. That is, although the true probability model to generate a real observed variables m _t is unknown, the true probability model, and shall be given approximately by (5), approximate models the true probability model The parameter (vector θ component) of equation (5) that minimizes the KL (Calbach librarian) divergence (KL information amount), which is the pseudorange of

For example, the vector θ may be determined using the Akaike information amount criterion that minimizes the KL information amount. Therefore, the model learning unit 22 uses the logarithm log (p (m ₁ , m ₂ ,..., M _t | θ)) of the probability p (m ₁ , m ₂ ,..., M _t | θ) as a log likelihood. , Akaike Information Standard AIC,
AIC = −2 × p (m ₁ , m ₂ ,..., M _t | θ) + 2 × N
Calculate as Here, N represents the number of model parameters that can be freely set, and here is (l + 4), which is the dimension of θ.

The model learning unit 22 calculates the AIC value for each set of values that are a plurality of predetermined θ values. Specifically, a set of θ values predetermined here is
a ₁ , a ₂ ,..., a _l , σ _τ , σ _q , σ _e , σ _n
A plurality of grid points set on the vector space spanned by. Then, the θ that minimizes the calculated AIC value is selected from the plurality of values that are the predetermined θ. That is, the model learning unit 22 first searches the vector θ that minimizes the AIC value roughly globally by a grid search.

を求める。 The model learning unit 22 further optimizes the vector θ (which is assumed to be θ ₀ ) obtained by the grid search as an initial value by nonlinear optimization.

Ask for.

とする。ここで終了条件は、更新前のθと更新後のθとの距離が予め定めたしきい値を下回っているなどの広く使われる条件を適用すればよい。 Here, as an example of nonlinear optimization, there is a method of obtaining the log likelihood log (p (m ₁ , m ₂ ,..., M _t | θ)) by a Kalman filter. Using the vector θ ₀ obtained by the grid search as an initial value, the vector θ maximizing the logarithmic likelihood is repeatedly and sequentially updated until, for example, a predetermined termination condition is satisfied by the steepest descent method. Θ

And Here, as the termination condition, a widely used condition such as a distance between θ before update and θ after update being less than a predetermined threshold may be applied.

を、ノイズパラメータとして記憶部１２に格納する。また、このノイズパラメータ管理部２３は、信号推定部２４からノイズパラメータの要求を受けて、記憶部１２に格納したノイズパラメータを読み出し、信号推定部２４に対して出力する。 The noise parameter management unit 23 is obtained by calculation by the model learning unit 22

Is stored in the storage unit 12 as a noise parameter. The noise parameter management unit 23 receives a request for a noise parameter from the signal estimation unit 24, reads the noise parameter stored in the storage unit 12, and outputs the noise parameter to the signal estimation unit 24.

によって表される確率モデルで近似的に表されるものとする。ここにｒ_tは時刻ｔの自発脳磁場成分を意味し、δ波(約２Ｈｚ)，θ波(約４Ｈｚ)，α波(約１０Ｈｚ)，β波(約２０Ｈｚ)，low-γ波(約４０Ｈｚ)，high-γ波(約８０Ｈｚ)の成分から構成されるものとし、各成分を、

と表したものである。なお、ここでのθは、θ波に対応する自発脳磁場成分であることを意味するものとして便宜的に利用しているもので、パラメータのベクトルθとは異なる。 The signal estimation unit 24 according to the present embodiment includes a spontaneous brain magnetic field estimation unit 24a and an induced brain magnetic field estimation unit 24b. The spontaneous brain magnetic field estimation unit 24a reads out resting data (the brain magnetic field data obtained by measuring a subject in a resting state), which is stored in the storage unit 12 and includes noise and the spontaneous brain magnetic field. Here, the spontaneous brain magnetic field estimating unit 24a

Is approximately represented by a probability model represented by Here, r _t means a spontaneous brain magnetic field component at time t, and δ wave (about 2 Hz), θ wave (about 4 Hz), α wave (about 10 Hz), β wave (about 20 Hz), low-γ wave (about about 40 Hz), high-γ wave (about 80 Hz) components,

It is expressed. Here, θ is used for convenience to mean a spontaneous brain magnetic field component corresponding to the θ wave, and is different from the parameter vector θ.

  Each of these components is expressed by the following quasi-periodic vibration model.

ここでωは、δ，θ，α，β，γ_l，γ_hのそれぞれをとる。また、ｆ_ωは、成分ω（δ，θ，α，β，γ_l，γ_hのそれぞれをとる）の中心周波数であり、ｖ_t ^(ω)は、平均が０、分散がσ_ω ²であるガウス分布に属するとの条件、

を満足する確率的入力である。

Here, ω is δ, θ, α, β, γ _l , or γ _h . F _ω is the center frequency of the component ω (takes each of δ, θ, α, β, γ _l , and γ _h ), and v _t ^(ω) has an average of 0 and a variance of σ _ω ² . The condition of belonging to a Gaussian distribution,

Is a stochastic input that satisfies

となる。ここにおいて、左辺のθ_pのθはθベクトルを表すθであり、右辺中、σ_θ ²とあるθは、θ波に対応する自発脳磁場成分を表す便宜的なものである。 From the above, the θ vector that is a model parameter representing the component of the spontaneous brain magnetic field is

It becomes. Here, θ of θ _{p on} the left side is θ representing a θ vector, and θ, which is σ _θ ^{2 in the} right side, is a convenient one representing a spontaneous brain magnetic field component corresponding to the θ wave.

Spontaneous brain magnetic field estimation unit 24a, a storage unit resting time series containing the noise and spontaneous brain magnetic fields are stored in the 12 data m _1, m _2, ..., reads the m _t from the storage unit 12.

The spontaneous brain magnetic field estimation unit 24a requests the noise parameter management unit 23 for a noise parameter, and uses the noise parameter output by the noise parameter management unit 23 in response to this request and the noise (5). The time series estimation data ν ₁ , ν ₂ ,. Spontaneous brain magnetic field estimation unit 24a resting time series containing the noise and spontaneous brain magnetic field data m _1, m _2, ..., the corresponding time series estimation data of the noise from the data of the corresponding time m _t and subtracting the time of the data, 'obtaining the _i = m _{i -ν i,} data of the condition being assumed that there is no _{noise, m' m 1, m '} 2, ..., m' obtain _t.

Then, the spontaneous brain magnetic field estimating unit 24a gives the probability p (m ′ ₁ , m ′) that the resting data is actually observed in time series when the vector θ _p represented by the equation (10) is given. ₂ ,..., M ′ _t | θ _p ) is determined as the maximum vector θ _p . In this case as well, although the true probability model for generating the actual observed variable m ′ _t is unknown, the true probability model is approximately given by equation (7) and approximated to the true probability model. Since the parameter (7) (the component of the vector θ _{p in} (10)) that minimizes the KL (Cullback librarian) divergence (KL information amount), which is a pseudo-range with a typical model, is obtained. is there.

For example, the vector θ in equation (10) may be determined using the Akaike information amount criterion that minimizes the KL information amount. Therefore, the spontaneous brain magnetic field estimator 24a uses the logarithm log (p (m ′ ₁ , m ′ ₂ ,..., M ′ _t ) of the probability p (m ′ ₁ , m ′ ₂ ,..., M ′ _t | θ _p ). | Θ _p )) is the log likelihood, and the Akaike information criterion AIC is
AIC = −2 × p (m ′ ₁ , m ′ ₂ ,..., M ′ _t | θ _p ) + 2 × N
Calculate as Here, N represents the number of model parameters that can be freely set, and here is 6, which is the dimension of the vector θ _p in equation (10).

The spontaneous brain magnetic field estimation unit 24a calculates the value of this AIC for each set of values that are a plurality of predetermined θ _p . Specifically, a set of θ _p values set in advance here is set on a vector space spanned by σ _ω ² (where ω is δ, θ, α, β, γ _l , γ _h ). A plurality of grid points can be used. Then, the set of θ values that minimizes the calculated AIC value is selected from the plurality of predetermined value sets as θ _p . That is, the spontaneous brain magnetic field estimation unit 24a first searches the vector θ _p that minimizes the AIC value roughly globally by a grid search.

を求める。 The spontaneous brain magnetic field estimation unit 24a further uses a vector θ _p (referred to as θ _p0 ) obtained by grid search as an initial value, and a parameter vector for the model of the spontaneous brain magnetic field optimized by nonlinear optimization.

Ask for.

とする。ここで終了条件は、更新前のθ_pと更新後のθ_pとの距離が予め定めたしきい値を下回っているなどの広く使われる条件を適用すればよい。 Here, as an example of nonlinear optimization, there is a method of obtaining the log likelihood log (p (m ′ ₁ , m ′ ₂ ,..., M ′ _t | θ _p )) using a Kalman filter. Using the vector θ _p0 obtained by the grid search as an initial value, the vector θ _p maximizing the log likelihood is repeatedly and sequentially updated until a predetermined termination condition is satisfied by the steepest descent method, for example. The obtained θ _p is

And Here termination conditions may be applied widely-used condition such as the distance between the theta _p and an updated previous theta _p is below a predetermined threshold value.

を、誘発脳磁場推定部２４ｂに出力する。 The spontaneous brain magnetic field estimation unit 24a thus obtained

Is output to the induced cerebral magnetic field estimation unit 24b.

のインパルス応答で近似できるものとする。 The induced cerebral magnetic field estimation unit 24b reads out the stimulation-time data (the cerebral magnetic field data obtained by measuring the subject in the stimulated state) including noise and the spontaneous brain magnetic field, which is stored in the storage unit 12. Here, the evoked cerebral magnetic field estimation unit 24b determines the evoked cerebral magnetic field as the following second-order linear model.

It can be approximated by the impulse response of.

The induced cerebral magnetic field estimation unit 24b of the present embodiment expresses the evoked cerebral magnetic field as a second-order stationary linear model having a stochastic input by the equation (11). In this equation (11), the model parameters are system input intensity and system parameters b ₁ and b ₂ that determine the response of the model to the system input.

Here, assuming that the evoked brain magnetic field can be approximated by an impulse response of a single linear model, the probability distribution of this stimulus (system input) is given as follows. That is, the evoked cerebral magnetic field estimation unit 24b includes the time series h ₁ , h ₂ ,... Of the impulse response h of the system having b ₁ and b ₂ as system parameters and the obtained stimulation time data m ₁ , m ₂ , Find the cross-correlation function with ... Then, the evoked brain magnetic field estimating unit 24b uses a binomial distribution B (2t, 0) centered on a temporal lag t (time t at which the evoked brain magnetic field appears) in which the probability distribution of the stimulus has the maximum absolute value of the cross-correlation. .5)
Suppose that Here, since time is discretized, a binomial distribution is adopted. The input intensity is a Gaussian distribution N (I, σ _I ² ) with the average value I being the value of the cross-correlation function at lag t.
Keep it as

The evoked brain magnetic field estimation unit 24b requests a noise parameter from the noise parameter management unit 23, and uses the noise parameter output by the noise parameter management unit 23 in response to the request and the equation (5) to calculate the noise. Time series estimation data ν ₁ , ν ₂ ,... Are obtained. The evoked cerebral magnetic field estimation unit 24b uses the parameters related to the spontaneous cerebral magnetic field input from the spontaneous cerebral magnetic field estimation unit 24a and the equation (7) to obtain time series estimation data g ₁ , g ₂ ,. obtain. The evoked field estimating unit 24b, stimulation time data m _1, m _2, ..., from the data of the corresponding time m _t, the time data and the corresponding time series estimation data of noise, when the spontaneous brain magnetic fields The data at the time corresponding to the sequence estimation data is subtracted to obtain m ″ _i = m _i −ν _i −g _i , and data in a state assumed to be free from noise and spontaneous brain magnetic field, m ″ ₁ , m ″ ₂ ,..., M ″ _t is obtained.

In the equation (11), data cannot be obtained directly from the model. Therefore, as described above, the stimulus probability distribution is expressed by an arithmetic expression using the Gaussian distribution (this includes the parameters b ₁ and b ₂ ). ), And the signal s _t at the time t of the evoked cerebral magnetic field estimated with respect to the stimulus is obtained by the stochastic difference equation using the expected value μ (t) at the time t and the variance σ ² in this arithmetic expression. and a signal s _t-1 at time t-1, previously obtained the time evolution of the computing equation of the state vector x includes in its components.

This arithmetic expression expresses the probability of time evolution (model of motion) p (s) from which the evoked brain magnetic field s _{t + 1} estimated at time _{t + 1} is obtained from the evoked brain magnetic field s _t estimated at times t and t−1. _{t + 1} | s _t , s _t-1 ).

The evoked cerebral magnetic field estimation unit 24b generates an initial state p (x ₀ ) (this is induced) for a state vector x _t including signals s _t and s _t−1 represented by system parameters b ₁ and b ₂ as components. The probability that the state vector is x _t when the data observed at time t is m ″ _t using the above-mentioned probability of time evolution. (X _t | m ″ _t ) is calculated. The induced brain magnetic field estimation unit 24b of the present embodiment obtains this probability p (x _t | m ″ _t ) as follows.

が成立するので、これを変形して、

を得る。ここでｐ（ｍ″_t）は、

と表すことができ、またｐ（ｘ_t）は、時間発展の確率を用いて、

と表すことができるから、結局（１３）式は、

とすることができる。 That is, from Bayes' theorem,

Since this holds, transform this,

Get. Where p (m ″ _t ) is

And p (x _t ) can be expressed using the probability of time evolution,

Therefore, after all, the equation (13) becomes

It can be.

The equation (16), "When was _t, probability state vector is x _t density p (x _t | m" data observed at time t is m _t) is the model p observations (m " _t | x _t ), the motion model p (x _{t + 1} | x _t ), and the initial state p (x ₀ ).

The induced cerebral magnetic field estimation unit 24b obtains a vector θs that maximizes the probability p (x _t | m ″ _t ) when a vector θs including system parameters b ₁ and b ₂ is given. The estimation unit 24b determines the vector θs using the Akaike information criterion that minimizes the KL information amount, so that the evoked cerebral magnetic field estimation unit 24b determines the logarithm log () of the probability p (x _t | m ″ _t ). p (x _t | m ″ _t )) is a log likelihood, and the Akaike information criterion AIC is
AIC = −2 × p (x _t | m ″ _t ) + 2 × N
Calculate as Here, N represents the number of model parameters that can be freely set, and here is 2, which is the dimension of θs.

The evoked brain magnetic field estimation unit 24b calculates the value of this AIC for each set of values that are a plurality of predetermined θs. Specifically, a set of θs values predetermined here is
b ₁ , b ₂
A plurality of grid points set on the vector space spanned by. Then, θs that minimizes the calculated AIC value is selected from the plurality of predetermined value sets as θs. That is, the evoked cerebral magnetic field estimation unit 24b first searches the vector θs that minimizes the AIC value roughly globally by a grid search.

を求める。 The induced cerebral magnetic field estimation unit 24b further uses the vector θs (referred to as θs ₀ ) obtained by the grid search as an initial value, and is a vector optimized by nonlinear optimization.

Ask for.

とする。ここで終了条件は、更新前のθsと更新後のθsとの距離が予め定めたしきい値を下回っているなどの広く使われる条件を適用すればよい。 Here, as an example of nonlinear optimization, there is a method of obtaining the log likelihood log (p (x _t | m ″ _t )) by a Kalman filter. Using the vector θs ₀ obtained by the grid search as an initial value, log likelihood The vector θs that maximizes the degree is obtained by repeatedly and sequentially updating until a predetermined termination condition is satisfied by, for example, the steepest descent method, and the obtained θs is

And Here, as the termination condition, a widely used condition such as a distance between θs before the update and θs after the update is less than a predetermined threshold may be applied.

と、（１１）式で表されるモデルの式とを用いて、誘発脳磁場の時間変化波形を演算して出力する。 The signal output unit 25 is obtained by calculation by the induced brain magnetic field estimation unit 24b.

And the time-varying waveform of the evoked brain magnetic field is calculated and output using the model equation represented by equation (11).

のパラメータを得てもよい。このときには、複数のセンサの一つを注目センサとするとき、注目センサに係るパラメータを得るためには、ノイズのみのデータや安静時データ、刺激時データは、当該注目センサによって得られたものを用いてもよい。 Furthermore, the control unit 11 of the present embodiment performs the above calculation on each data obtained from a plurality of sensors (SQUIDs in the case of a magnetoencephalograph), and for each sensor.

May be obtained. In this case, when one of a plurality of sensors is used as a target sensor, in order to obtain parameters related to the target sensor, noise-only data, resting data, and stimulation data are obtained from the target sensor. It may be used.

  The control unit 11 uses the parameters obtained for each of the plurality of sensors and the equation (11) to calculate the time-varying waveform of the evoked brain magnetic field at each sensor position and create an isomagnetic field diagram. May be output. Since this isomagnetic field diagram creation process is widely known, a detailed description thereof will be omitted here.

[Calculation by Monte Carlo method]
By the way, the calculation of the expression (16) in the control unit 11 is difficult because the integral of the right-hand side denominator (integral with the product of the prior probability and the posterior probability in Bayesian estimation) cannot be obtained as an analytic solution. Become. Therefore, the control unit 11 according to the present embodiment may include a plurality of arithmetic elements and perform calculations as follows.

Here, the plurality of arithmetic elements may use arithmetic elements such as GPU (Graphics Processing Unit) instead of the CPU. When a plurality of arithmetic elements can be used in this way, the control unit 11 causes each arithmetic element to randomly select values of x _t and x _t−1 from within the integration range, and the integration variable selected at random. The integrand values obtained by substituting the values of (if there are a plurality of values, randomly selected from the integration range) into the integrand are calculated in parallel. Then, the control unit 11 accumulates a plurality of integrand values obtained by the respective arithmetic elements and performs integration by the Monte Carlo method.

[Timing for obtaining noise-only data]
In the above description, the noise-only data used for the determination of the noise parameter may be obtained at any time as long as it is in the past from the time of calculation. For example, it may be immediately before obtaining stimulus data (data obtained prior to measurement immediately before a subject who actually measures stimulus data enters), or obtained after the subject leaves after measuring stimulus data. It may be the data obtained.

[Example of obtaining spontaneous brain magnetic field by Bayesian estimation]
In the description so far, the spontaneous brain magnetic field estimation unit 24a has the probability p (m ′) that the resting data is observed in the order of the measured time series when the vector θ represented by the equation (10) is given. ₁ , m ′ ₂ ,..., M ′ _t | θ) is obtained as a vector θ. However, the present embodiment is not limited to this.

That is, the spontaneous brain magnetic field estimation unit 24a uses, as the state vector x _t at time t, data g _t , g _t−1 ,... Representing the spontaneous brain magnetic field before time t and data representing each noise before time t. τ _t , q _t , q _t-1 ,…, q _{t-T + 2} , _et ,…, et _{t-l + 1} , n _t …
x _t = [τ _t , q _t , q _t−1 ,…, q _{t−T + 2} , e _t ,…, et _{l−l + 1} , n _t , g _t , g _t−1 ,…] ^T
Is generated. The shoulder T represents transposition.

Here, how much noise data and spontaneous brain magnetic field data are included from the data before the time t depends on the noise, and the number is sufficient to predict the corresponding noise at the time t + 1. It also depends on the model adopted, for example, in terms of the spontaneous brain magnetic field, there may be a case of g _t only one.

The spontaneous brain magnetic field estimator 24a has a probability p (x _t | m ′ _t when the vector θ of equation (10) is given by the same operation as the equation (16) performed by the induced brain magnetic field estimator 24b. ) Is maximized. Specifically, the spontaneous brain magnetic field estimation unit 24a determines the vector θ using the Akaike information amount criterion that minimizes the KL information amount. Thus spontaneous brain magnetic field estimation unit 24a, the probability _{p (x t | m 't} ) of the logarithm _{log (p (x t | m} ' t)) of the log likelihood, the Akaike information criterion AIC,
AIC = −2 × p (x _t | m ′ _t ) + 2 × N
Calculate as Here, N represents the number of model parameters that can be freely set, and here is 6 which is the dimension of θ.

  The spontaneous brain magnetic field estimation unit 24a may find θ that minimizes the AIC by a method that combines the already described grid search and nonlinear optimization.

  In this case, the calculation of the equation (16) can also be performed by the Monte Carlo method.

[Example of including spontaneous brain magnetic field in Bayesian estimation]
Further, in the description so far, the evoked brain magnetic field estimation unit 24b uses the parameters related to the spontaneous brain magnetic field input from the spontaneous brain magnetic field estimation unit 24a and the equation (7), and the time series estimation data g of the spontaneous brain magnetic field g _1, g _2, to give ... a stimulation time of data m _1, m _2, ..., from each corresponding time data of m _t, the time-series estimation data of noise corresponding time data and, of spontaneous brain magnetic fields By subtracting the corresponding time data of the time series estimation data, m ″ _i = m _i −ν _i −g _i is obtained, and m ″ ₁ , data in a state assumed to have no noise and spontaneous brain magnetic field, m " ₂ , ..., m" _t .

However, the present embodiment is not limited to this. That evoked field estimating unit 24b, stimulation time data m _1, m _2, ..., from the data of the corresponding time m _t, and subtracting the data of the corresponding time of the time-series estimation data of noise, m ″ _I = m _i −ν _i is obtained, data “m ″ ₁ , m ″ ₂ ,..., M ″ _t in a state assumed to be free of noise is obtained, and processing is performed as follows. Good.

In this case, the induced cerebral magnetic field estimation unit 24b includes data g _t , g _t−1 ,... Representing spontaneous brain magnetic fields before time t in the state vector x _t at time t.
x _t = [s _t , s _t-1 , g _t , g _t-1 , ...] ^T
And The shoulder T represents transposition.

Here, how much the spontaneous brain magnetic field is included from the data before time t is a number sufficient to predict g _{t + 1} . It also depends on the model adopted, there may be a case, for example, g _t only one.

を求めることができる。 The calculation of the expression (16) performed by the evoked brain magnetic field estimation unit 24b can be performed in the same manner as described above even if the state vector components increase in this way. Just as I said

Can be requested.

[Example of including noise in Bayesian estimation]
Furthermore, in the description up to this point, evoked field estimating unit 24b, the time-series estimation data g _1, g ₂ of spontaneous brain magnetic field, ... was obtained, and stimulation when data m _1, m _2, ..., each of m _t By subtracting the corresponding time data of the time series estimation data of the noise and the corresponding time data of the time series estimation data of the spontaneous brain magnetic field from the corresponding time data, m ″ _i = m _i −ν _i -G _i was obtained, and data m ″ ₁ , m ″ ₂ ,..., M ″ _t were obtained in a state where there was no noise and no spontaneous brain magnetic field.

However, the present embodiment is not limited to this. That evoked field estimating unit 24b, stimulation time data m _1, m _2, ..., from the data of the corresponding time m _t, the time-series estimation data g _1, g ₂ of spontaneous brain magnetic field, ... corresponding time To obtain m ″ _i = m _i −g _i, and obtain data m ″ ₁ , m ″ ₂ ,..., M ″ _t assuming that there is no spontaneous brain magnetic field. The processing may be performed as follows.

In this case, the evoked cerebral magnetic field estimation unit 24b adds the data τ _t , q _t , q _t−1 ,..., Q _{t−T + 2} , representing each noise before the time t to the state vector x _t at the time t. including e _t ,…, e _{t-l + 1} , n _t …
x _t = [s _t , s _t-1 , τ _t , q _t , q _t-1 ,…, q _{t-T + 2} , e _t ,…, e _{t-l + 1} , n _t ] ^T
And The shoulder T represents transposition.

  Here, how much data of each noise is included from the data before time t depends on the noise, and is a number sufficient to predict the corresponding noise at time t + 1. This depends on the model used.

を求めることができる。 The calculation of the expression (16) performed by the evoked brain magnetic field estimation unit 24b can be performed in the same manner as described above even if the state vector components increase in this way. Just as I said

Can be requested.

[Example of including noise and spontaneous brain magnetic field in Bayesian estimation]
Further evoked field estimating unit 24b, stimulation time data m _1, m _2, ..., each of m _t, as the _{m "1, m" 2,} ..., as m _"t, perform processing as follows May be.

That is, the evoked cerebral magnetic field estimation unit 24b adds the data τ _t , q _t , q _t−1 ,..., Q _{t−T + 2} , et representing the noise before the time t to the state vector x _t at the time _t. ,…, Including e _{t-l + 1} , n _t …
_{x t = [s t, s} t-1, τ t, q t, q t-1, ..., q t-T + 2, e t, ..., e t-l + 1, n t, g t, g _t-1 , ...] ^T
And The shoulder T represents transposition.

  Here, how much data of each noise is included from the data before time t depends on the noise, and is a number sufficient to predict the corresponding noise at time t + 1. This depends on the model used.

The induced brain magnetic field estimation unit 24b also includes data g _t , g _t−1 ,... Representing the spontaneous brain magnetic field before time t in the state vector x _t at time t. Here, how much the spontaneous brain magnetic field is included from the data before time t is a number sufficient to predict g _{t + 1} . It also depends on the model adopted, there may be a case, for example, g _t only one.

を求めることができる。 The calculation of the expression (16) performed by the evoked brain magnetic field estimation unit 24b can be performed in the same manner as described above even if the state vector components increase in this way. Just as I said

Can be requested.

[AR model variable conversion in grid search]
In yet AR model representing the time evolution of the stationary noise component e _t in the previous sections, the stability condition is hard to see, performs iteration variable converted by the following method, be converted into first-order AR model Good.

を用いる。ただし、ｉ＝１，２，…，ｌ−１である。ここにａの肩にある添字のｌは、ｌ次の変数であることを意味する。この（１７）式による変換を、１次となるまで繰り返して行うと、各変換の段階で、ａ_j ^jが得られる（ここにｊ＝１，２，…ｌ）。すると、安定なＡＲモデルとなるための必要十分条件が、

と得られる。そこで、グリッドサーチにおいては、この（１８）式で表される範囲を定義域として、この定義域を均等に分割した格子点を用いればよい。 That is, as a method of converting an L-order AR model into an L-1 order AR model,

Is used. However, i = 1, 2,..., L-1. Here, the subscript l on the shoulder of a means an l-order variable. When the conversion according to the equation (17) is repeated until it becomes the first order, a _j ^j is obtained at each conversion stage (where j = 1, 2,... L). Then, the necessary and sufficient conditions to become a stable AR model are

And obtained. Therefore, in the grid search, the range represented by the equation (18) may be used as a defined area, and lattice points obtained by equally dividing the defined area may be used.

[Other models]
In general, the probability distribution of the input of the evoked brain magnetic field model is considered to differ depending on the active site and the type of stimulus, and it is assumed that not only a linear model but also a non-linear model is appropriate. The present embodiment is not limited to magnetoencephalographic data, but can be used for processing general measurement data. Therefore, the user may create a probability distribution of a model of data to be measured (in this example, an induced brain magnetic field model) according to each hypothesis. The processing of the present embodiment can be executed even when various probability distributions are assumed.

  In addition, the model of the spontaneous brain magnetic field is also expressed here by equation (8). However, when the time change of the amplitude is large, it is assumed that the power of the stochastic input term changes with time. Then, the power may be described in some form of function including a parameter, and the process may be performed by including the parameter in the model parameter. Further, if the peak frequency of each subject's spontaneous brain magnetic field band is significantly different from the typical one, the center frequency may also be included in the model parameter. In these cases, estimation can be performed as one of the model parameters as in the above-described processing.

[Operation]
The present embodiment has the above configuration and operates as follows. That is, the data processing apparatus 1 of the present embodiment accepts noise-only data, resting data, and stimulation data and holds them in the storage unit 12, and as illustrated in FIG. Is used to estimate a noise parameter included in a predetermined noise model (S1).

  Further, the data processing apparatus 1 removes noise included in the resting data using the estimated noise parameter and a predetermined noise model, and refers to the resting data from which the noise has been removed to determine the predetermined data. The parameters included in the resting data model are estimated (S2).

  Then, the data processing device 1 removes from the stimulation data using the noise parameter estimated in the process S1 and a predetermined noise model, and the spontaneous brain magnetic field included in the stimulation data is the parameter estimated in the process S2. And a predetermined resting data model. Then, the data processing apparatus 1 refers to the data from which the noise and the spontaneous brain magnetic field have been removed in this way, and estimates parameters included in a predetermined model at the time of stimulation (induced brain magnetic field model) (S3).

  The data processing device 1 estimates and outputs evoked cerebral magnetic field data using the parameters estimated in process S3 and the evoked cerebral magnetic field model (S4).

  An embodiment of the present invention will be described. Here, we describe an example when the subject is given an evoked brain magnetic field when a stimulus “getting a button pressed” is given.

  FIG. 4 shows the waveform of the induced cerebral magnetic field obtained by the conventional method obtained by averaging the measurement data of 169 times passed through the LPF (Low Pass Filter) (the signals of a plurality of sensors are overlapped). The same), and isomagnetic field diagrams obtained from these waveforms. FIG. 5 is an example of time-series data actually obtained including noise, and an isomagnetic field diagram obtained based on the time-series data. The conventional result shown in FIG. 4 is obtained by averaging 169 times of data as shown in FIG.

  FIG. 6 shows the induced brain magnetic field output from the data processing apparatus 1 according to the example of the present invention through the processing described in the present embodiment based on the data (single measurement data) shown in FIG. It is a magnetic field diagram obtained from waveforms (represented by superimposing signals from a plurality of sensors) and these waveforms.

  When the conventional analysis diagram shown in FIG. 4 is compared with the analysis diagram of the present embodiment shown in FIG.

  FIG. 7 is an example of time series data actually obtained including noise when the intensity of the spontaneous brain magnetic field is large, and an isomagnetic field diagram obtained based on the time series data. From the isomagnetic field diagram estimated by the conventional method obtained by averaging the data as illustrated in FIG. 7 a plurality of times and averaging, and the plurality of data as illustrated in FIG. FIG. 9 shows a calculation result of coherence with a plurality of isomagnetic field diagrams (illustrated in FIG. 8) obtained by estimation by the data processing apparatus 1. FIG. 10 shows an isomagnetic field diagram estimated by the conventional method. As can be seen from FIGS. 9 and 10, the coherence in the region where the signal intensity is large in the isomagnetic field diagram obtained in FIG. 10 exceeds 0.8 as shown in FIG. 9. That is, it is understood that the estimation by the data processing device 1 of the present embodiment is sufficient from the viewpoint of robustness.

  DESCRIPTION OF SYMBOLS 1 Data processing device, 11 Control part, 12 Storage part, 13 Operation part, 14 Output part, 15 Input interface, 21 Data recording part, 22 Model learning part, 23 Noise parameter management part, 24 Signal estimation part, 25 Signal output part .

Claims (5)

Receiving means for receiving measured time-series resting data including noise and spontaneous brain magnetic field, and time-series stimulation data including noise, spontaneous brain magnetic field and evoked brain magnetic field ;
The noise parameter included in the noise calculation formula that schematically represents the noise using at least one noise parameter is determined based on a criterion determined using only noise data obtained in the past, and a time series of noise Means for obtaining estimated data ;
Means for subtracting the time series estimation data of the noise from the data at the corresponding time of the resting data to obtain data in a state assumed to be free of noise;
When the signal included in the spontaneous brain magnetic field is schematically expressed as a sum of the components of the spontaneous brain magnetic field, the model parameter representing the component of the spontaneous brain magnetic field is the data assumed to be free of noise. Based on a model parameter representing the component of the spontaneous brain magnetic field and a model schematically representing the signal contained in the spontaneous brain magnetic field as a sum of the components of the spontaneous brain magnetic field. An estimation means for obtaining time series estimation data of
The time series estimation data of the noise and the time series estimation data of the spontaneous brain magnetic field are subtracted from the data of the corresponding time from the stimulation time data, in a state where there is no noise and no spontaneous brain magnetic field. Means for obtaining and outputting data;
Including a data processing apparatus. The data processing apparatus according to claim 1, wherein
When the resting data at time t is mt, the estimating means has a probability distribution p where the state vector including data representing the spontaneous brain magnetic field before time t and data representing noise before time t is xt. When the model parameter that maximizes (xt | mt) is a probability p (xt) that the state vector is xt at time t, and the state vector is xt at time t, the spontaneous brain magnetic field The probability distribution p (mt | xt) of resting data at time t represented by a model that approximately represents the signal included in the model as a sum of spontaneous magnetic field components is used as the posterior probability, and the prior probability and posterior A data processing apparatus which is estimation means for estimating by Bayesian estimation using an integration result having a product of probabilities as an integrand . A data processing apparatus according to claim 2, wherein
Including a plurality of arithmetic unit elements;
It said estimating means, in the Bayesian estimation when performing integration to the integrand a product of the prior probability and the a posteriori probability, wherein when randomly assigns the value of the integration variables selected from the integration range A data processing apparatus for performing the integration by a Monte Carlo method by performing a calculation of a value of an integrand on a plurality of arithmetic element elements in a distributed manner. A data processing device according to any one of claims 1 to 3,
The means for obtaining the time series estimation data of the noise uses data received by the receiving means as the data of only the noise prior to actual measurement of the resting data or stimulation data. Computer
Receiving means for receiving measured time-series resting data including noise and spontaneous brain magnetic field, and time-series stimulation data including noise, spontaneous brain magnetic field and evoked brain magnetic field ;
The noise parameter included in the noise calculation formula that schematically represents the noise using at least one noise parameter is determined based on a criterion determined using only noise data obtained in the past, and a time series of noise Means for obtaining estimated data ;
Means for subtracting the time series estimation data of the noise from the data at the corresponding time of the resting data to obtain data in a state assumed to be free of noise;
When the signal included in the spontaneous brain magnetic field is schematically expressed as a sum of the components of the spontaneous brain magnetic field, the model parameter representing the component of the spontaneous brain magnetic field is the data assumed to be free of noise. Based on a model parameter representing the component of the spontaneous brain magnetic field and a model schematically representing the signal contained in the spontaneous brain magnetic field as a sum of the components of the spontaneous brain magnetic field. An estimation means for obtaining time series estimation data of
The time series estimation data of the noise and the time series estimation data of the spontaneous brain magnetic field are subtracted from the data of the corresponding time from the stimulation time data, in a state where there is no noise and no spontaneous brain magnetic field. Means for obtaining and outputting data;
Program to function as.

Images