Abstract

MEG/EEG beamformer source imaging is a promising approach which can easily address spatiotemporal multi-dipole problems without a priori information on the number of sources and is robust to noise. Despite such promise, beamformer generally has weakness which is degrading localization performance for correlated sources and is requiring of dense scanning for covering all possible interesting (entire) source areas. Wide source space scanning yields all interesting area images, and it results in lengthy computation time. Therefore, an efficient source space scanning strategy would be beneficial in achieving accelerated beamformer source imaging. We propose a new strategy in computing beamformer to reduce scanning points and still maintain effective accuracy (good spatial resolution). This new strategy uses the distribution of correlation values between measurements and lead-field vectors. Scanning source points are chosen yielding higher RMS correlations than the predetermined correlation thresholds. We discuss how correlation thresholds depend on SNR and verify the feasibility and efficacy of our proposed strategy to improve the beamformer through numerical and empirical experiments. Our proposed strategy could in time accelerate the conventional beamformer up to over 40% without sacrificing spatial accuracy.

1. Introduction

Magnetoencephalography (MEG) and electroencephalography (EEG) are noninvasive imaging technologies which provide functional information about human brain dynamics by providing millisecond temporal images over the entire brain. These technologies have been widely used to diagnose epilepsy and forward neuroscience research. Particularly, MEG/EEG source localization estimates current sources from measured spatiotemporal data. Inherently, MEG/EEG source localization is mathematically ill-posed; that is, it has no unique solution and is very sensitive to noise. For a couple of decades, many researchers have tried to develop methods to deal with these difficulties in calculation, resulting in extensive research and commercialization of MEG/EEG source localization methods (see [1, 2] for a review).

In the early 1990s, beamformer techniques originated in the field of antenna signal processing [3]. Application in MEG/EEG source imaging soon followed. A beamformer represents a kind of linear spatial filter acting on spatial or spatiotemporal data within sensor space. It allows a neural signal produced only at designated source point to pass, filtering out signals originating from other source points. Therefore, even without a priori information on source quantity, beamformers can effectively image brain activities within a source space under the assumption that sources are uncorrelated [46]. Many varieties of beamformers have been investigated [710] and are roughly categorized into 2 classes: adaptive beamformers using measurement information, and nonadaptive beamformers independent of measurement information [11]. Among the variety of options, the minimum-variance (MV) beamformer has been most widely used and deeply investigated in regard to MEG/EEG source localization problems [5, 1217].

Recently, source imaging has been gaining more attention on continuous MEG/EEG (unaveraged) and single-trial MEG/EEG data in understanding rapidly changing brain dynamics [18, 19]. This understanding can better facilitate real-time brain activity monitoring, neurofeedback, brain computer interface (BCI) [2022], among others [23, 24]. Brain signals can usually be measured by means of MEG or EEG systems and their real-time interpretation can provide a variety of applications. Beamformer is a promising technique easily dealing with spatiotemporal multi-dipole source problems as well as being robust to noise. Regarding real-time source imaging, beamformer speed generally depends on the quantity of scanned source points, that is, the number of scanning points of interesting brain area. For this purpose, reducing beamformer scanning points (accelerating beamformer) without sacrifice of spatial resolution (accuracy) is greatly beneficial in real-time source imaging. In the present work, a procedure is proposed to accelerate beamformer as well as to maintain effective spatial resolution. Sensor measurements are composed of a (composite) linear combination of lead-field vectors (sensitivity of sensors to sources) at active points as well as noise [15]. For that reason, sensor measurements may produce relatively higher correlation values than correlation around inactive points. In the present paper, such reasoning is further investigated culminating in a proposal for reducing scanning points during beamformer source imaging.

This paper is an extended version of a short conference paper presented in BIOMAG 2010 [25] and is organized in sections. In Section 2, a conventional MV-beamformer is briefly explained. In Section 3 we follow with a theoretical proposal for reducing scanning regions and then discuss an effective strategy for computing the correlation distribution between measurement data and lead-field vectors. Next, both simulated and empirical experimental results are presented to verify the feasibility and efficacy of the proposed procedure. Lastly, we discuss other efforts in achieving accelerated beamformer.

2. Minimum-Variance (MV) Beamformer

Beamformer techniques are categorized into two classes: one is adaptive and the other is nonadaptive. A nonadaptive spatial filter is independent of the measurement, but an adaptive spatial filter depends on the measurement. Among beamformers, the minimum-variance (MV) beamformer is superior in accuracy to others [11] and has been widely used in MEG/EEG source imaging [10]. MV beamformers can be scalar type and vector type. In the scalar type beamformers, the source orientation could be estimated to yield maximum source power (for details, see [10, 15]) or may be predetermined accordingly. On the contrary, the source orientation could be simultaneously estimated in the vector type beamformer (uses a set of three weights and each weight detects one component in Cartesian coordinate system) [4, 5, 26, 27]. Throughout this paper, MV vector type beamformer is used.

A vector type beamformer enables to estimate simultaneously source orientation and magnitude. A vector type spatial filter consists of a set of three weight vectors, 𝑤𝑥(𝐫), 𝑤𝑦(𝐫), and 𝑤𝑧(𝐫) depending on 𝑥,𝑦,and𝑧 components of the source vector, respectively. Denoting the weight matrix by 𝐖(𝐫)=[𝑤𝑥(𝐫);𝑤𝑦(𝐫);𝑤𝑧(𝐫)], vector-type beamformer is derived solving the following optimization: 𝐖(𝐫)=argmin𝐖𝑇𝐂𝐖,subjectto𝐖𝑇(𝐫)𝐋(𝐫)=𝐈,(2.1) where 𝐂=𝑚(𝑡)𝑚𝑇(𝑡)𝑡 is measurement covariance matrix estimated by time average, 𝐋(𝐫)=[𝑙𝑥(𝐫);𝑙𝑦(𝐫);𝑙𝑧(𝐫)] is the lead-field matrix representing the sensitivity of the whole sensor array to source activity at 𝐫, and 𝐈 represents an identity matrix. 𝑙𝜉(𝐫) is a lead-field vector of unit source activity oriented to 𝜉-axis at 𝐫. The weight matrix and output power of this vector type spatial filter are expressed as follows:𝐖(𝑟)=𝐂1𝐋𝐋(𝐫)𝑇(𝐫)𝐂1𝐋(𝐫)1,𝑄(𝐫,𝑡)2𝑡=𝐋𝑇(𝐫)𝐂𝟏𝐋(𝐫)1.(2.2)

3. Scanning Reduction Strategy

3.1. Definitions and Principles

In this section, a new strategy of reducing scanning points, thereby accelerating beamformer source imaging, is proposed without sacrificing accuracy. In beamformer source imaging, full source space scanning is necessary; sources are located on the brain’s cortical area and thus all scanning of the brain may be time intensive but achieves whole brain images. On the other hand, partial scanning accelerates beamformer source imaging. Evidently, reducing the number of scanning source points while keeping scanning resolution can accelerate the source imaging process without losing spatial resolution significantly. Our proposed idea is reducing the scanning region (eventually reduce the number of scanning points) by using the correlation between measurement and lead-field vectors, defined as𝑝(𝐫,𝑡)=𝑚(𝑡)𝑙(𝐫)(𝑚𝑡)𝑙(𝐫).(3.1) Here 𝑚(𝑡) and 𝑙(𝐫) represent the measurement vector at time 𝑡 and the lead-field vector at 𝐫, respectively. The inner product between 𝑚(𝑡) and 𝑙(𝐫) vector yields scalar value. At each source point 𝐫, the root mean square (RMS) correlation is computed as follows:𝑝rms(𝐫)=1𝐾𝐾𝑡=1𝑝(𝐫,𝑡)2,(3.2) where 𝐾 is the number of time points. In general, sensor measurements consist of a linear combination of (composite) lead-field vectors [15] at active points and noise:=𝑚(𝑡)=𝐋𝑠(𝑡)+𝑛(𝑡)𝑁𝑠𝑖=1𝑙𝐫𝑖𝑠𝑖(𝑡)+𝑛(𝑡)𝑛𝑎𝑖𝑎=1𝑙𝐫𝑖𝑎𝑠𝑖𝑎(𝑡)+𝑛(𝑡),𝑛𝑎𝑁𝑠.(3.3) Here 𝑠(𝑡)=[𝑠1(𝑡),𝑠2(𝑡),,𝑠𝑁𝑠(𝑡)]𝑇 is a source vector representing a signed source magnitude at each scanning point 𝐫𝑖 and time 𝑡. 𝐋=[𝑙(𝐫1),𝑙(𝐫2),,𝑙(𝐫𝑁𝑠)] is a matrix consisting of all lead-field vectors, and 𝑛(𝑡) is noise at time 𝑡, being commonly colored and independent of signal. Assuming that measurement is produced by a few significantly strong source activations (such case is common in BCI or spontaneous activity, evoked response activity, etc.), sources at most scanning points are inactive, yielding a negligibly small magnitude of 𝑠𝑖(𝑡), and only a few sources (𝑛𝑎 sources) are active, the corresponding 𝑠𝑖(𝑡) yields significant magnitudes. A correlation value between 𝑚(𝑡) and lead-field 𝑙(𝐫) at scanning point 𝐫 is expressed as𝑚(𝑡)𝑙(𝐫)=𝑛𝑎𝑖𝑎=1𝑠𝑖𝑎𝐫(𝑡)𝑙𝑖𝑎𝑙(𝐫)+𝑛(𝑡)𝑙(𝐫).(3.4)

Independence from noise on a signal yields that the 2nd term on the right side of (3.4) is negligible and thus the correlation value is mainly affected by the 1st term. When a scanning point is around an active source point, that is, 𝐫 is very close or equal to 𝐫𝑖𝑎, then 𝑙(𝐫𝑖𝑎)𝑙(𝐫) gives a significant value. However, when the scanning point is away from all active points, 𝑙(𝐫𝑖𝑎) and 𝑙(𝐫) may be much less correlated, and thus 𝑙(𝐫𝑖𝑎)𝑙(𝐫) may yield a considerably smaller value in magnitude. Due to such arguments, the correlation value 𝑝(𝐫,𝑡) may be more significant at active points than others. Overall, RMS correlation values at active points are likely to give more significant magnitude than inactive points. Furthermore, it was observed in Section 4.1 that RMS correlations tended to have higher values around active than inactive points. From such reasoning and observations, the following strategy is proposed.(i) Correlation values 𝑝(𝐫,𝑡) for each scanning point 𝐫 in source space and time point 𝑡 are computed using (3.1).(ii)RMS correlation values 𝑝rms(𝐫) at each scanning source point 𝐫 are eventually computed using (3.2).(iii)Scanning regions are reduced by selecting source points with higher RMS correlation values than given correlation thresholds. Correlation thresholds are estimated empirically for various SNRs via the Monte Carlo simulation, while corresponding criteria are formulated in a least square sense (optimization) at a later time.(iv)Beamformer is performed at such reduced source region.

3.2. Correlation Computing Strategy

To reduce scanning points, RMS correlation should be computed at every point within the source space. Computing RMS correlations across the entire source space, including a plenty of corresponding time points, may be time intensive. For example, about 500 time samples (2 second-long duration for a 250 Hz acquisition system) are naturally needed to update brain dynamics or decode user intention in existing BCI systems. As shown in (3.1), correlation value computation is necessary for all time points, thereby prolonging overall computation. Downsampling is a possible strategy, but it risks losing important information. Instead, the use of eigen temporal window (called ETWB) is proposed in this work.

By projecting a measurement matrix onto eigen-temporal space while ignoring negligibly small eigenvalues and corresponding eigenvectors, the number of time points is dramatically reduced, thereby reducing correlation computation time. Such a procedure is formulated as follows.(i)Spatiotemporal measurement matrix 𝐌=[𝑚(𝑡1);𝑚(𝑡2);;𝑚(𝑡𝐾)] of size 𝑆×𝐾 can be decomposed by singular value decomposition (SVD) as follows: 𝑢𝐌=1𝑢𝐾𝜆1000000𝜆𝑆𝑣𝑇1𝑣𝑇𝑆=𝐔𝚺𝐕𝐓,𝜆1𝜆2𝜆𝑆.(3.5)(ii)Projection matrix 𝐕ETWB=[𝑣1;𝑣2;;𝑣𝐽] is constructed using 𝐽 significant temporal eigenvectors 𝑣𝑗 (determined by the significance of corresponding eigenvalues) of measurement matrix 𝐌. Threshold of eigenvalues or the number of eigenvalues 𝐽(𝐾) can be determined empirically.(iii)Projected measurement matrix 𝐌ETWB of size 𝑆×𝐽 is defined to compute as follows: 𝐌ETWB=𝐌𝐕ETWB.(3.6)

Reduction of the time dimensionality of the measurement matrix is thus ensured, thereby significantly reducing correlation computation time. Computing time of about 3 seconds is reduced to 0.7 seconds, which will be discussed later in this paper. The effectiveness of this strategy is demonstrated in the following section.

4. Results

4.1. Observation on Correlation Distribution

Correlation distribution between active and inactive points was investigated through numerical experiments in this section. A spherical homogeneous conductor head model with an 8.5 cm radius, centered at origin, was used for forward computing [28], while sensor geometry from a MEGVISION Yokogawa system with a 160-sensor gradiometer array was adopted. Two sources were assumed to be located within a box-shaped region defined by 8.5𝑥8.5, 8.5𝑦8.5, and 2.0𝑧8.5 (cm) with a total of 10,000 scanning points comprising equally spaced lattices. Such sources were located at (−3.7, −3.6, 5.0) and (3.3, 4.1, 4.3) cm, both directed to (1,0,0). Synthetic data was generated by adding white Gaussian noise to calculated sensor values computed through a forward model. A synthetic magnetic field with 500 time samples was generated with an SNR of 1. The SNR was defined by 𝐒F2/𝐍F2, where 𝐒 and 𝐍 are signal and noise matrices, respectively, and F is a Frobenius norm. Time courses for dipoles were provided as damping sinusoidals with different wavelengths; these sources were almost uncorrelated. A detailed configuration description is shown in Figure 1. Numerical experiments were conducted on a workstation (2x AMD Opteron CPU 2.3 GHz, 64 bit OS, and 64 GB RAM).

Correlation of full-source space with 10,000 scanning points was computed to illustrate the distribution using (3.2). In Figure 2, three different views of correlation distribution such as 𝑥𝑦-, 𝑥𝑧-, and 𝑦𝑧-plane projections are depicted. This 𝑥𝑦-plane projection map and others were generated by 𝑝rms(𝑖𝑥,𝑖𝑦)=max𝑖𝑧𝑝rms(𝑖𝑥,𝑖𝑦,𝑖𝑧) for each pixel point (𝑖𝑥,𝑖𝑦). Figure 2 shows that correlation values produced around active rather than inactive points were relatively higher. Therefore, beamformer scanning within regions with relatively high correlation values seemed to be adequate in reconstructing source information. It is obvious from this observation that choosing a scanning region in the proposed manner can be a good strategy in accelerating beamformer.

4.2. Correlation Threshold

In the previous section, correlation distribution was found to play a key role in getting a priori information on source locations. Intuitively, use of such correlation distribution can accelerate beamformer source imaging when beamformer scanning is confined to regions with relatively high correlation values. To apply such a concept, proper correlation thresholds should be predetermined in a reasonable way. In general, correlation distribution may influence various factors such as source locations, source magnitudes, SNR, and sensor geometry. For simplicity, correlation threshold criterion is assumed to depend only on SNR; estimation is thus possible with information easily obtainable or observable. Under such consideration, the following formulation on correlation threshold corrthresh(SNR) is proposed:corrthresh(SNR)=𝑚+𝜎a(SNR),𝑝𝑚=meanrms,()𝜎2𝑝=varrms.()(4.1)

Here 𝑚 and 𝜎 are the first- (mean) and second- (standard deviation) order statistics of RMS correlation values, which are easily estimated. Further, a(SNR) is a nonlinear function depending on SNR and should be determined with more caution. In order to estimate a(SNR) effectively, a Monte Carlo simulation study was executed. For a great number of simulated data with various SNRs (with real noise), proper a(SNR) values were computed and fitted to the given model in a least square sense—assuming that correlation values at active source points are higher than thresholds. From our experience, the following simple logarithmic model was ideally adopted in the present work:a(SNR)=max{𝛼log(𝛽SNR),0}.(4.2)

Correlation threshold criterion corrthresh(SNR) was estimated through the Monte Carlo simulation study. The Monte Carlo simulation was conducted with 200 randomly distributed dipole sources at random orientations across the entire brain. For each dipole (fully correlated source between lead-field vector 𝑙(𝐫) and measurement 𝑚(𝑡)), five different empirical noise realizations and 15 different SNRs (between 0.01 and 8) were generated by noise power control. Hence, a total of 15,000 single dipole problems were generated in estimating the correlation threshold. Estimated optimal parameters (𝛼 and 𝛽) for such ETWB strategy are as follows:a(SNR)max{0.214log(11.309SNR),0}.(4.3)

Empirical MEG noise used in the Monte Carlo simulation study was acquired under the following experimental paradigm.

Empirical MEG data was collected from a healthy male volunteer (24-year-old) who participated after appropriate informed consent was acquired. During spontaneous activity with eyes closed in a magnetically shielded room, an acquisition period of 120 seconds was launched. Data was digitized at 2 kHz with the online lowpass filter set at 500 Hz, postprocessing digital filter applied from 1 to 100 Hz, excluding 50 Hz due to electrical power conditions. Data was collected from a whole-head gradiometer system with 160 channels (MEGVISON Yokogawa system).

Our proposed threshold criterion was generated from a single dipole simulation study; thus other possible source configurations are not sure to be perfectly accounted for. Nevertheless, it was found that our proposal proved quite effective throughout the work.

4.3. Experimental Results: Simulated Data

In this section, investigation into the degree of acceleration and reduction of scanning points resulting from our proposed ETWB strategy was studied. For this purpose, beamformer and proposed scanning reduction strategy was integrated to apply to a simulated two dipole problem, which was generated in Section 4.1. For effective inversion of data covariance and higher output SNR [10], the inversion of data covariance matrix (𝐂+𝜀𝐈)1 in place of 𝐂1 was used, where 𝜀 and 𝐈 are a regularization factor and an identity matrix, respectively [29]. Regularization factor 𝜀 is estimated by 𝜀=max{𝜆1,𝜆2,,𝜆𝑆}105, where 𝜆𝑖,{𝑖=1,2,,𝑆} is the eigenvalue of the data covariance 𝐂. We tested beamformer using the ETWB strategy addressed in Section 3.1. Reconstructed source imaging by this method is illustrated in Figure 3. According to (4.1), the correlation threshold was properly chosen in reducing source scanning regions. In the ETWB strategy, significant eigenvectors were chosen explaining the amount of information, up to 70%, from measurement matrix 𝐌. Thus, approximately one fifth of 500 total time samples were selected to create an eigen-temporal projection matrix 𝐕ETWB.

Figure 3 shows which particular region has higher correlation values than threshold, and beamformer in the reduced region resulted in images being identical to the beamformer source imaging with the entire source space scanning. In this experiment, the number of scanning points was reduced to about 3,000, that is, about one third of the total original scanning points of 10,000. Total elapsed time of beamformer was about 0.53 seconds, even including correlation computation, while conventional beamformer time (full-brain scan) was about 1.03 seconds.

Source Imaging Effect on SNR
Conventional full-scan beamformer was compared with reduced scanning beamformer using the proposed strategy to determine how much scanning reduction was attained over different SNRs (0.1 or 4). The reduced scanning regions by ETWB strategy are illustrated in Figure 4. The same ETWB strategy as in previous experiments was used. For simulated data with SNR 0.1, Figure 4(a) shows the results from the reduced scanning region using the ETWB strategy and the reconstructed image from beamformer. Computation time for the conventional beamformer was about 1.03 seconds over the entire region of 10,000 points. Conversely, our proposed strategy took about 0.6 seconds, and scanning points were reduced up to about 4,400 points.
With regards to high SNR 4, the same experiment was conducted. Figure 4(b) shows the results from the reduced scanning region using the ETWB strategy for processing simulated data. In this case, computation time for beamformer using our proposed strategy was about 0.37 seconds, and scanning points were reduced to about 2,300 points. Such results show that our proposed strategy can accelerate beamformer while maintaining spatial resolution of source images with respect to SNRs. As previously mentioned in the initial simulation (SNR = 1), the number of significant eigenvectors needed for the ETWB strategy was about one fifth of the total time points.

Source Imaging Effect on Source Correlation
In this section, the performance of our proposed beamformer was investigated over different source correlations between two dipole sources. It is known that beamformer performance degrades for heavily correlated sources. To investigate how much our proposed method relatively effects source correlation, simulated data was generated under the same configurations as in Section 4.1, except for different source correlations such as 0.1, 0.6, and 0.9. To generate simulated data with different source correlations, the time course of the second dipole source was made using the formula: 𝑤2(𝑡)=(1𝜉)𝑤1(𝑡)+𝜉𝑤2(𝑡), where 𝑤1(𝑡) and 𝑤2(𝑡) represent the first and second time courses of the sources, respectively. The parameter 𝜉 controls the degree of the correlation between 𝑤1(𝑡) and 𝑤2(𝑡). Figure 5 shows the reduced scanning regions by our proposed strategy and reconstructed source imaging of beamformer. As correlation increases, sources get more broadly reconstructed; thereby spatial resolution is poorer. However, reduced scanning regions seem to be almost same; therefore, there is no significant effect on scanning reduction of correlation between sources.

Source Imaging Effect on Source Depth and Extended Source Model
In general, as brain sources are activated far away from the sensor surface of the MEG/EEG system, it is dramatically less sensitive to those sources and localization accuracy gets worse. In this section, the performance of beamformer using our proposed strategy is investigated as sources gradually move away from the sensor surface. For this investigation, the same source configuration as in Figure 1 was used except for the depth of one source. Three source depths (𝑧-coordinate), 5.0, 2.9, and 0.8 cm, were considered. Figure 6 shows the reduced scanning regions by our proposed strategy over varying source depths. Obviously, as the source is located deeper, the source signal gets weaker; thus it is harder to reconstruct deep sources by our proposed beamformer. The reduced scanning region got slightly narrower since SNR gets smaller, but it is not significantly noticeable.
Further, the performance of our proposed beamformer was tested on extended source models. Extended sources were generated for three different source span diameters: 0, 2.0, and 5.0 cm. Here the source with a 0 cm span diameter means a dipole source. Extended sources were generated on the plane (𝑧=5.0 cm) parallel to the xy-plane. Each center of the extended sources is located on (−3.7, −3.6, 5.0) cm. All dipole sources within extended source regions were oriented to the same direction. Expectedly, there was no big variation on different span diameters, but extended sources with span diameter 5.0 cm seemed slightly broad spread. However, reduction of scanning regions was almost unchanged (not shown).

Comparison on Computing Time
So far, many kinds of simulated data over different SNRs, different source depths, different source correlations as well as different extended source spans have been tested to compare reduced scanning beamformer using our proposed strategy to conventional full-scan beamformer. Overall, our proposed beamformer shows its relative acceleration in time and yields comparable in accuracy to the conventional beamformer. Figure 7 depicts comparative computational time over numerous simulated problems, which were generated as follows.(i) For single dipole problem, 1,000 dipoles were randomly distributed within the spherical head model. For each dipole, 7 single dipole problems with different SNRs (between 0.01 and 8) were generated by adding white Gaussian noise, thus yielding 7,000 single dipole problems. Particularly, for different diameter of source simulation, 6 different diameters of source (between 0 and 50 mm) at each source were considered. Thus 6,000 single source problems were generated with keeping SNR of 1 in this case.(ii)For two-dipole problem, 1,000 pairs of two dipoles were randomly chosen among 1000 dipoles generated in single dipole problems. For each pair, two-dipole problems were generated with 7 different SNRs (between 0.01 and 8), 7 different source correlation coefficients (between 0 and 1), and 6 different interdistance between two sources (between 0 and 100 mm). Hence, total 294,000 two-dipole problems were generated to test.
For each problem, the number of scanning point and computational time were estimated. Evidently, computational time is linearly proportional to the reduced scanning points. In our proposed strategy, we found that averaged correlation computing time including SVD computation was about 0.15 seconds, which was added to computational time. The averaged computational time of the conventional beamformer was about 1.035 seconds. Each point in Figure 7 represents an averaged computing time with standard deviation over all corresponding simulated problems.
Figure 7(a) shows that our strategy speeds up as SNR gets higher. It means that reduction rate of scanning points is higher as SNR goes higher. Further, computational time has no noticeable difference as extent of source and source correlation vary (Figures 7(b) and 7(c)). However, computational time over inter-source distance seems slight difference (Figure 7(d)). As two sources get far, computational time increases slightly, but not substantial.
It is remarkable that correlation computation time did not change with respect to SNRs; however, we found that computation time of SVD varied slightly over SNRs. We guess that such small variation of SVD computing time may stem from different matrix characteristics. Total original scanning points amounted to 10,000 and the temporal size of the measurement matrix was 500. Figure 7(a) shows that beamformer utilizing our proposed ETWB strategy computed significantly quicker than the conventional beamformer, yielding about 30–60% improvement. These experiments demonstrate how the beamformer with our strategy was effective in accelerating beamformer source imaging.

4.4. Experimental Results: Empirical Data

In the present section, reduced scanning beamformer with our proposed strategy was applied to empirical MEG data. Four kinds of empirical MEG datasets were acquired on a whole-head gradiometer MEGVISION Yokogawa system with 160 channels—median nerve stimulations (right hand/left hand) and auditory stimulations (left ear/right ear). Experimental paradigm is detailed as follows.

A healthy male volunteer (24-year-old) participated in the MEG measurement after appropriate informed consent. First, the somatosensory electrostimulation median nerve of his right/left hand was stimulated and all measurements were done in a magnetically shielded room. During measurement, the subject was instructed to keep his eyes closed. A 2.9 (right)/3.5 (left) mA current and a 2 Hz sampling were applied with 0.3-millisecond current pulses. Interstimulus-interval (ISI) was randomized with 0.5 s duration. Second, the auditory stimulation of his right/left ear was stimulated with eyes closed. A 80 dB sound pressure and a 2 Hz sampling were applied with 40 ms plateau with 10 ms rise and falls. Also, ISI was randomized with a 2 s duration (random 50%, 1–3 s). Data were digitized at 2 kHz, lowpass filtered at 500 Hz as well as postprocessed via a digital filter at 1–100 Hz (median nerve stimulation)/1–50 Hz (auditory stimulation), excluding 50 Hz due to electrical power conditions. In the case of median nerve stimulation, a total 399 (right-hand)/418 (left-hand) single trials were acquired. Total 73 (right-ear)/70 (left-ear) single trials were obtained in the case of auditory simulations.

Each collected single trial data (unaveraged) with 0–250 millisecond (time-locked, 500 time samples) time window after stimulation onset were analyzed. Such analysis can better facilitate real-time brain activity monitoring for neurofeedback and brain computer interface (BCI) [23], among others. For single trial analysis, the reduced scanning beamformer with our proposed strategy and the conventional full-scan beamformer were used. We used entire the source region consisting of a total of 8,640 scanning points and a temporal window with 500 samples corresponding to a duration of 250 milliseconds. The correlation threshold was accordingly determined through formula defined in (4.1). Eigen-temporal projection matrix 𝐕ETWB was estimated as having explained 70% of the measurement matrix information. To compare source reconstruction performance, all of each single trial was reconstructed and their source information (number of scanning points × number of time samples × number of trials) was averaged over the trials. Also, total computing time elapsed for single trial analysis was measured over different datasets and two kinds of beamformers. Figure 8 shows the averaged-over-trial-reconstructed source power maps of both conventional full-scan beamformer and reduced scanning beamformer with our proposed strategy for median nerve stimulation data (right and left). To verify the capability of our proposed reduced scanning beamformer, we attempted single-trial analysis for empirical data. Both beamformers were applied to compare and averaged source power maps over trial were generated by considering zero power out of the reduced scanning region. Since each trial has different reduced scanning region based on correlation threshold, averaged-over-trial source power map was seen in the whole region.

Beamformer source imaging for averaged median nerve stimulation data yielded focused source activity in the upper contralateral posterior central sulcus (not shown here), which is relevant to existing literature (see [30] and therein). As illustrated in Figure 8, strong source activity appeared in the upper contralateral area and most mildly active sources seemed scattered around strong source activity, which is relevant to beamformer results (not shown) for averaged median nerve stimulation data. Evidently, both beamformers showed almost identical results, thereby reduced scanning beamformer being comparable in accuracy to the conventional one. However, significant improvement in computing time was achieved in the reduced scanning beamformer with our proposed strategy, and comparative total elapsed computing time for each experiment was tabulated in Table 1.

These results tell that overall improvement in computing time of reduced scanning beamformer with our proposed strategy is over 40%. In conclusion, our proposed reduced scanning beamformer is feasible and very applicable in accelerating the conventional beamformer.

5. Discussion

5.1. Other Possible Accelerating Strategies

In this work, one beamformer using ETWB strategy was proposed to reduce the beamformer scanning regions by computing correlation distribution, thereby accelerating beamformer. Similar to the ETWB strategy, another strategy is possibly proposed by doing eigen-temporal projection and in addition doing eigen-spatial projection. We call this strategy eigen-spatiotemporal window-based method (we call it ESTWB). The discarding of relatively insignificant sensor information via the ESTWB strategy seems more efficient in relation to computation time while the ETWB strategy will inherently provide a more accurate correlation distribution as well as yielding more reliable beamformer source imaging. In addition to ETWB and ESTWB strategies, temporal downsampling is another idea; however, it possibly loses useful information.

Computing correlation in (3.1) and (3.2) amounts to a spatial matched filter when normalization term of measurement vector 𝑚(𝑡) in (3.1) is dropped. So, it is not surprising that correlation behaves like another spatial filter, thereby yielding rough less accurate but fast spatial imaging. This information is used for prescreening to determine final scanning region of more elegant spatial filter, that is, MV beamformer. Such argument gives more general perspective in developing speeding-up strategies of scanning methods. Any hybrid between a fast scanning method and an accurate scanning method can be possible when they can yield synergy effect due to reasonable combination strategies.

Furthermore, there exist other strategies to speed up scanning methods. One can consider downsampling in spatial scanning to reduce computation effort. It is commonly applicable in scanning methods like beamformer. When our strategy would be applied together, scanning methods should be boosted easily. Another easy strategy is to increase computational resources such as parallel computing to use many personal computers at the same time. This provides powerful computing ability; thus tens to hundreds of times improvement would be possible. Certainly, such a strategy demands a relatively high cost; thus, it may be not easily affordable to most investigators. It is noted that our proposed strategy of simple consideration of correlation distribution can be achieved on a single personal computer.

5.2. Applicability of Our Proposed Strategy

As discussed in Section 4.1, simple description and observation of correlation distribution around active and inactive sources motivated us to conduct this kind of investigation. For a few strong brain sources being well separated, our strategy is favorably relevant. However, for many sources, the correlation threshold may become lower due to more mixed source intercorrelations; thus, scanning points will be slightly reduced and our proposed correlation threshold criterion may not be effective in such cases. Nevertheless, real-time brain activity monitoring systems (or BCI systems) are usually designed to find relatively rapidly changing source information (or discriminative information among different conditions); thus, a few strong sources may be accountable and helpful for providing such information.

In real-time monitoring or BCI systems, seamless monitoring or control decision is of great importance. In reality, most systems have some delay coming from computation implicitly required to do any defined action; however, such delay is reasonably small or one may design the system to have delay constant or within an affordable bound; thus it looks real-time system to user. Particularly, source imaging analysis introduced in these systems may not avoid significant delay incurring from intense computation. However, well-desiged system to hide such significant delay may be possible. For example, assuming that every 1-second-long window data is collected to analyze and source imaging requires 1 second, such system may be designed that data is analyzed to give final result to the system during system being collected next 1 second-long data; then system continues to do the same procedure. This system can achieve seamless procedure (it looks real-time system from user) with constant delay of 1 second. From this perspective on real-time monitoring or BCI systems, the faster source imaging can achieve the higher monitoring scan or decision rate. Therefore, our proposed strategy yielding about 40% speed-up of source imaging is well applicable in this reasoning.

On the other hand, in BCI systems, one could define a region-of-interest (ROI) manually [23, 24] or by elegant procedures [31, 32]. Source imaging on such ROI could increase the spatial specificity or enhance SNR in source space [31]. Furthermore, it is possible to achieve far more accelerated source imaging than whole brain source imaging. To achieve real-time source activity monitoring or source information-based BCI [20, 2224], naturally higher spatial resoluted imaging within ROI is beneficial in such development. However, it requires more computation time. Our proposed strategy can be still applicable on ROI in the same manner. Furthermore, single-trial analysis using fast beamformer has been beneficial in developing many applications [18, 19]. Accordingly, our proposed method was applied to single-trial data, such as empirical median nerve stimulation and auditory stimulation, as discussed in Section 4.4. It revealed several dynamic activation regions, usually not visible during averaged data analysis (not shown here). Such single-trial analysis and online BCI application will be investigated in more detail and reported in a subsequent paper.

5.3. Other Issues

Through the Monte Carlo simulation, our correlation threshold criterion, depending on SNR, was formulated in a least square sense. This criterion was dependent on sensor configuration (geometry); so it should be reestimated under other sensor configurations. Factors other than SNR may have some effects on correlation thresholds. Accordingly, more thorough research is currently under investigation.

Even though the conventional vector-type MV-beamformer among the beamformer variants was adopted in this work, scalar-type or other nonadaptive beamformers are similarly applicable. Furthermore, this strategy is very straightforward for EEG or simultaneous MEG/EEG data application, and subsequent work on simultaneous MEG/EEG beamformers [33] is currently under investigation.

Acknowledgments

This work was supported by a KRF grant (KRF-331-2008-1-D00768), NRF grant (NRF-2010-32A-B00283), the PLSI supercomputing resources of KISTI, the BioImaging Research Center at GIST, and the National IT Industry Promotion Agency (NIPA-2011-C1090-1131-0006). The authors appreciate Dr. Haruta at Yokogawa Electric Corp for his assistance with MEG data acquisition.