= Beamforming methods =
''Authors: Hui-Ling Chan, Francois Tadel, Sylvain Baillet''

The estimation of source distribution is an important step to understand the brain activity from EEG and MEG data. Dipole fitting, minimum norm estimation and beamformer are three commonly used methods. It has been proved that beamforming methods provide good spatial resolution. This tutorial will show how to apply beamforming methods to MEG data and obtain the statistic map of source activation.

We are going to use the protocol '''!TutorialRaw''' created in the [[http://neuroimage.usc.edu/brainstorm/Tutorials#Processing_continuous_recordings|introduction tutorials]]. If you have not followed these tutorials yet, please do it now.

<<TableOfContents(2,2)>>

== Introduction ==
Beamforming methods scan each targeted source position <<latex($\bf r$)>> and estimate the spatial filter <<latex(${\bf w}_{\bf r,q}$)>>. By multiplying with the MEG recordings <<latex(${\bf m}(t)$)>>, the spatial filter <<latex(${\bf w}_{\bf r,q}$)>> outputs the temporal waveform <<latex($y_{\bf r,q}(t)$)>>  of the dipole source at that position with the dipole orientation <<latex($\bf q$)>> as below:

<<latex($ y_{\bf r,\bf q}(t) = {{\bf w}_{\bf r,\bf q}}^{\rm T}{\bf m}(t)\enspace, $)>>

where 'T' indicates the transpose of a matrix or vector. The beamforming spatial filter can be vector-type or scalar-type.

==== Vector-type beamformer ====
For each position <<latex($\mathbf{r}$)>>, three orthogonal spatial filters <<latex($\mathbf{W}_{\mathbf{r}}=\left[\mathbf{w}_{\mathbf{r},\mathbf{q}_x},\mathbf{w}_{\mathbf{r},\mathbf{q}_y},\mathbf{w}_{\mathbf{r},\mathbf{q}_z}\right]$)>> are computed by applying the unit-gain constraint as well as the minimum norm and minimum variance criteria as below [Van  Veen et al., 1997]:

<<latex($ \hat{\mathbf{W}}_{\mathbf r}=(\mathbf{C}_{\rm s}+\alpha\mathbf{I})^{-1}\mathbf{L}_\mathbf{r}\left(\mathbf{L}_\mathbf{r}^{\rm T}(\mathbf{C}_{\rm s}+\alpha\mathbf{I})^{-1}\mathbf{L}_\mathbf{r}\right)^{-1} \enspace,$)>>

where <<latex(${\bf C}_{\rm s}$)>> is the covariance matrix of MEG recordings during window <<latex($T_{\rm_s}$)>>, <<latex($\mathbf{I}$)>> is the identity matrix, <<latex($\mathbf{L}_{\mathbf{r}}$)>> is the gain matrix for the dipole located at position <<latex($\bf{r}$)>>, and <<latex($\alpha$)>> is the regularization parameter which compromises the minimum norm and minimum variance criteria.

==== Scalar-type beamformer ====
For each position <<latex($\mathbf{r}$)>>, the source orientation <<latex($\mathbf{q}$)>> is first estimated to enable the spatial filter to output the source activity with maximum power or fitting some other criteria. The dipole orientation can be obtained by exhaustive search  [Vrba and Robinson, 2000] or analytical solution [Chen et al., 2006]. The estimated dipole orientation <<latex($\hat{\mathbf{q}}$)>> is then applied to calculate the spatial filter as follows [Vrba and Robinson, 2000]:

<<latex($ \hat{\mathbf{w}}_{\mathbf{r,q}}=(\mathbf{C}_{\rm  s}+\alpha\mathbf{I})^{-1}\mathbf{L}_\mathbf{r}\hat{\bf{q}}\left(\hat{\bf q}^{\rm  T}\mathbf{L}_\mathbf{r}^{\rm  T}(\mathbf{C}_{\rm  s}+\alpha\mathbf{I})^{-1}\mathbf{L}_\mathbf{r}\hat{\bf q}\right)^{-1} \enspace.$)>>

== LCMV beamformer ==
LCMV stands for '''linearly-constrained minimum variance'''. LCMV beamformer is vector-type beamformer. For each position <<latex($\mathbf{r}$)>>, this method calculates the neural activity index, which is  interpreted as the estimate of source to noise variance.

==== Neural activity index ====
For each position <<latex($\mathbf{r}$)>>, the variance of source activity during active state <<latex($T_{\rm a}$)>> is calculated as follows:

<<latex($ \rm{Var}_{\rm a}(\bf r) =trace \left \{\left[{\bf L}_{\bf r}^{\rm T}{\bf C}_{\rm a}{\bf L}_{\bf r}\right]^{-1}\right \}\enspace, $)>>

or

<<latex($ \rm{Var}_{\rm  a}(\bf r) =\lambda_1 \left \{\left[{\bf L}_{\bf r}^{\rm T}{\bf C}_{\rm  a}{\bf L}_{\bf r}\right]^{-1}\right \}\enspace, $)>>

where <<latex(${\bf C}_{\rm a}$)>> is the covariance matrix computed from MEG recordings during active state <<latex($T_{\rm a}$)>> and <<latex($\lambda_1\{\}$)>> indicates the maximum eigenvalue of the expression in braces.

When the location <<latex($\mathbf{r}$)>> is far from sensors, the elements of lead field matrix  <<latex($\mathbf{L}_{\bf r}$)>> are small. So the elements of <<latex($\left[{\bf L}_{\bf r}^{\rm T}{\bf C}_{\rm  a}{\bf L}_{\bf r}\right]^{-1}$)>> are generally large and the estimated variance for the deep source becomes large. When the location <<latex($\mathbf{r}$)>> is close to sensors, the elements of lead field matrix  <<latex($\mathbf{L}_{\bf r}$)>> are large. It results in small values of the elements of <<latex($\left[{\bf L}_{\bf r}^{\rm T}{\bf C}_{\rm  a}{\bf L}_{\bf r}\right]^{-1}$)>> and small estimated variance for the superficial source. To reduce the effect caused by the depth of source location, the estimated source variance <<latex($\rm{Var}_{\rm a}(\bf r)$)>> is normalized by the noise variance <<latex($\rm{Var}_{\rm c}(\bf r)$)>> as follows:

<<latex($ \rm{Var}_{\rm N}(\bf r) =\rm{trace} \left \{\left[{\bf L}_{\bf r}^{\rm T}{\bf C}_{\rm  a}{\bf L}_{\bf r}\right]^{-1}\right \} / \rm{trace} \left \{\left[{\bf L}_{\bf r}^{\rm T}{\bf C}_{\rm c}{\bf L}_{\bf r}\right]^{-1}\right \}\enspace, $)>>

or

<<latex($ \rm{Var}_{\rm  N}(\bf r) =\lambda_1 \left \{\left[{\bf L}_{\bf r}^{\rm T}{\bf  C}_{\rm  a}{\bf L}_{\bf r}\right]^{-1}\right \} / \lambda_1 \left  \{\left[{\bf L}_{\bf r}^{\rm T}{\bf C}_{\rm c}{\bf L}_{\bf  r}\right]^{-1}\right \}\enspace, $)>>

where <<latex(${\bf C}_{\rm c}$)>> is the covariance matrix computed from MEG recordings during control state <<latex($T_{\rm c}$)>>. The normalized variance <<latex($ \rm{Var}_{\rm N}(\bf r)$)>> is called neural activity index (NAI).

==== LCMV process ====
 1. Drag and drop all of the epochs in ''Subject01 / right / right | timeoffset (98 files)'' into ''Process1''.
  .

 1. Click on [RUN] and select the process "'''Sources > LCMV Beamformer'''".<<BR>><<BR>>''' {{attachment:LCMVmaxeig.png}} ''' <<BR>>     <<BR>>
 1. With this window you can select the method you want  to use to estimate the source variance and noise variance, the time window of active state, the regularization parameter, and the sensors you are going  to use for this estimation. You can edit  the following options:<<BR>>     <<BR>>
  * '''Comment''': This field contains what is going to be displayed in the database explorer.
  * '''Method''':  Please select ''Unconstrained (max eigenvalue)'' or ''Unconstrained (trace)''. The method  ''Unconstrained (max eigenvalue) ''calculates the variances of source activity and noise by using the maximum eigenvalue. The method ''Unconstrained (trace)'' calculates the variances of source activity and noise by using trace. The method ''Cortial Constrained''  use the orientation of cortical surface as dipole orientation. Note that the method'' Cortial Constrained ''is only available when the Head Model  provides the orientation for each grid location. The  unconstrained methods may give more smooth results than the cortical  constrained method.
  * '''Time window of active state''': Three parameters have to be set: ''Time range of interest'','' Active window size'', and ''Temporal resolution''. The neural activity index is computed using the variance of source activity during active state <<latex($T_a$)>>.  ''Active window size'' defines the duration of active state. Large active window size (at least larger than the number of sensors) is recommended. The time window of active state is going to be slided in ''Time range of interest'' with the interval set to be ''Temporal resolution''. If there are multiple active state windows, the time label of the neural activity index will be set as the middle time of each active state window, that is, <<latex($t_1,t_2,t_3,\ldots,t_{N-1},t_{N}$)>> in the following figure. <<BR>><<BR>>''' {{attachment:slidingwindow.png||height="219",width="437"}} '''<<BR>><<BR>>''' '''
  * '''Regularization parameter''':''' '''This value will be multiplied with the maximum eigenvalue of the covariance matrix computed from the recordings during ''Time range of interest''. The larger value of regularization parameter may give the more smooth results.
  * '''Sensors type''':  Modalities that are used for the reconstruction. Here we only have one  type of MEG sensors, so nothing to change.
  * Click on ''Run''.
 1. Two new files are available in the database explorer.<<BR>> <<BR>>''' '''
  * '''right: LCMV: spatial filter(Unconstr)|0.0_295.8ms''': This file saves the three orthogonal spatial filters for each dipole position. These spatial filters are normalized by the variance of noise. Note that there is one spatial filter for each dipole location if the method'' Cortically Constrained ''is used. <<BR>><<BR>>''' {{attachment:SpatialFilterright.png}} ''' <<BR>><<BR>>''' '''
  * '''right: LCMV: neural activity index(Unconstr)|0.0_295.8ms''': This file saves the neural activity index for every dipole position. <<BR>><<BR>>''' {{attachment:NAIright2.png}} '''<<BR>> <<BR>> Double click on this file. Right click on anywhere of the pop-out window. From the pop-out menu, select ''Colormap:Sources > Maximum: Custom...''. Set ''Minimum'' to be '''1.2''' and ''Maximum'' to be '''2.0'''. Click ''Ok''. <<BR>> <<BR>>''' {{attachment:Colormaplimit.png}} ''' <<BR>> <<BR>> Set the time to be anytime between 0 and 295.8 ms. Then click ''Surface'' tab on the right panel of database explorer. Set ''Data options > Amplitude'' to be '''50%'''.  <<BR>><<BR>>''' ''' {{attachment:NAIrightUnconstrained.png}} <<BR>><<BR>>The maps of neural activity index estimated using ''Unconstrained (trace)'' and ''Unconstrained (max eigenvalue)'' are shown in the left and right figures, respectively. Both maps display strong activations in the left sensory area.<<BR>><<BR>>''' '''

== Maximum constrast beamformer (MCB) ==
Maximum constrast beamformer (MCB) is scalar-type beamformer. For each position  <<latex(\bf r)>>, it provides an analytical solution of dipole orientation <<latex(\bf q)>>, which maximizes the constrast of beamformer outputs between active state and control state as bellow:

<<latex($\hat{\bf q}=\arg \max_{\bf q} \rm{Var}_{\rm  N}(\bf r) = \arg \max_{\bf q} \bf{w}_{\bf{r,q}}^{\rm T} {\bf C}_{\rm a} \bf{w}_{\bf{r,q}}  \left(\bf{w}_{\bf{r,q}}^{\rm T} {\bf C}_{\rm c} \bf{w}_{\bf{r,q}} \right)^{-1}\enspace,$)>>

where  <<latex(${\bf C}_{\rm a}$)>> is the covariance matrix computed from the measurements during active state and <<latex(${\bf C}_{\rm c}$)>> is the covariance matrix computed from the measurements during control state.

==== Solution of dipole orientation ====
The solution of spatial filter can be rewritten as

<<latex($ \mathbf{w}_{\mathbf{r,q}}=(\mathbf{C}_{\rm   s}+\alpha\mathbf{I})^{-1}\mathbf{L}_\mathbf{r}\bf{q}\left({\bf q}^{\rm  T}\mathbf{L}_\mathbf{r}^{\rm   T}(\mathbf{C}_{\rm   s}+\alpha\mathbf{I})^{-1}\mathbf{L}_\mathbf{r}{\bf q}\right)^{-1}  = {\bf A}_{\bf r}{\bf q}\left( {\bf q}^{\rm T}{\bf B_{r}q}\right)^{-1}\enspace,$)>>

where both <<latex(${\bf A}_{\bf r}=({\bf C}_{\rm s}+\alpha {\bf I})^{-1}{\bf L_r}$)>> and <<latex(${\bf B}_{\bf r}={\bf L_r}^{\rm T}{\bf A}_{\bf r}$)>> depend only on the dipole location <<latex($\bf r$)>>. Then the objective function for obtaining the dipole orientation <<latex(\bf q)>> can be rewritten as

<<latex($ \rm{Var}_{\rm  N}(\bf r) = {\bf q}^{\rm T}{\bf A}_{\bf  r}^{\rm T}{\bf C}_{\rm a}{\bf A}_{\bf r}{\bf q} \left({\bf q}^{\rm T}{\bf  A}_{\bf r}^{\rm T}{\bf C}_{\rm c}{\bf A}_{\bf r}{\bf q}   \right)^{-1}={\bf q}^{\rm T}{\bf P}_{\bf r}{\bf q} \left({\bf q}^{\rm T}{\bf Q}_{\bf r}{\bf q}   \right)^{-1}\enspace,$)>>

where both of the 3-by-3 matrices <<latex(${\bf P}_{\bf  r} = {\bf A}_{\bf  r}^{\rm T}{\bf C}_{\rm a}{\bf A}_{\bf r}$)>> and <<latex(${\bf Q}_{\bf  r}={\bf A}_{\bf  r}^{\rm T}{\bf C}_{\rm c}{\bf A}_{\bf r}$)>> does not depend on the dipole orientation <<latex($\bf q$)>>. The solution of dipole orientation <<latex($\hat{\bf q}$)>> is the eigenvector corresponding to the maximum eigenvalue of the matrix <<latex(${\bf Q}_{\bf  r}^{-1}{\bf P}_{\bf  r}$)>>.

==== Statistical mapping ====
After obtaining the dipole orientation <<latex($\hat{\bf{q}}$)>> and the spatial filter <<latex($\bf{\hat{w}_{r,q}}$)>> for dipole position <<latex(\bf r)>>, the F-statistic value at time <<latex($t$)>> can be obtained using the following formula:

<<latex(${\rm F}({\bf r},t)=\hat{{\bf w}}_{\bf{r,q}}^{\rm T} {\bf  C}_{T_t}\hat{\bf{w}}_{\bf{r,q}}  \left(\hat{\bf{w}}_{\bf{r,q}}^{\rm T} {\bf  C}_{\rm c} \hat{\bf{w}}_{\bf{r,q}} \right)^{-1} \enspace,$)>>

where <<latex(${\bf C}_{T_t}=\sum_{T_t}{\|{\bf m}(t)-(\sum_{T_t}{{\bf m}(t)/N_{\rm T}}) \|^{2}}/(N_{\rm T}-1)$)>> is the covariance matrix computed from the measurements during window <<latex($T_t$)>>, <<latex($N_{\rm T}$)>> is the size of <<latex($T_t$)>>, and the range of <<latex($T_t$)>> is  <<latex($t \pm 0.5N_{\rm T}$)>> . <<latex($T_t$)>> is a segment of active state <<latex($T_{\rm a}$)>>. If the number of <<latex($N_{\rm T}$)>> is small, the peak of temporal dynamics of F-statistic value may shift. In this case, it is more appropriate to replace the covariance matrix <<latex(${\bf  C}_{T_t}$)>> with

<<latex(${\bf P}_{T_t}=\sum_{T_t}{\|{\bf m}(t)\|^{2}}/(N_T-1)\enspace,$)>>

or

<<latex(${\bf Z}_{T_t}=\sum_{T_t}{\|{\bf m}(t)-(\sum_{T_{\rm B}}{{\bf m}(t)/N_{\rm B}}) \|^{2}}/(N_{\rm T}-1)\enspace,$)>>

where <<latex($T_{\rm B}$)>> represetns the window of baseline and its size is <<latex($N_{\rm B}$)>>. When <<latex(${\bf P}_{T_t}$)>> is used, the F-statistic value represents the estimate of source power to noise variance.

==== MCB process ====
 1. Drag and drop all of the epochs in ''Subject01 / right / right | timeoffset (98 files)'' into ''Process1''.

 1. Click on [RUN] and select the process "'''Sources > Maximum Contrast Beamformer'''".

 1. With this window you can select the dipole orientation you want to use, the time window of active state, the window size you are going to use for F-statistics computation,  the regularization parameter, and the sensors you are going  to use for  this estimation. You can edit  the following options:<<BR>>     <<BR>>
  * '''Comment''': This field contains what is going to be displayed in the database explorer.
  * '''Dipole orientation''':  Please select ''Unconstrained'', which calculates the dipole orientation using the maximum constrast criterion. The method ''Cortial Constrained''  use the orientation of cortical surface as dipole orientation. Note that the method'' Cortial Constrained ''is  only available when the Head Model  provides the orientation for each  grid location. The  unconstrained methods may give more smooth results  than the cortical  constrained method.
  * '''Time windows''': Three parameters have to be set: ''Time range of active window'','' F statistic window size'', and ''Temporal resolution''. The dipole orientation and spatial filter is estimated using the variance of source activity during active state <<latex($T_a$)>>. ''F statistic window size''  defines the length of the window for F-statistic computation, that is, <<latex($N_{\rm T}$)>>. The time window for computing F statistics is going to be slided in ''Time range of active state'' with the interval set to be ''Temporal resolution''.  If there are more than one F statistic windows, the time label of the F statistics will be set as the middle time of each window, that is, <<latex($t_1,t_2,t_3,\ldots,t_{N-1},t_{N}$)>> in the following figure. ''' '''<<BR>><<BR>>''' '''
  * '''Regularization parameter''':''' '''This value will be multiplied with the maximum eigenvalue of the covariance matrix computed from the recordings during ''Time range of interest''. The larger value of regularization parameter may give the more smooth results.
  * '''Sensors type''':  Modalities that are used for the reconstruction. Here we only have one  type of MEG sensors, so nothing to change.
  * '''Method for F statistic computation''':  Modalities that are used for the reconstruction. Here we only have one  type of MEG sensors, so nothing to change.
  * Click on ''Run''.

 1. Two new files are available in the database explorer.<<BR>> <<BR>>

== Beamformer-based correlation/coherence imaging ==
==== Dynamic imaging of coherent sources (DICS) ====
[Under construction]''' '''

==== Spatiotemporal imaging of linearly-related source components (SILSC) ====
[Under construction]''' '''

== References ==
 1. Van  Veen BD, Van Drongelen W, Yuchtman M, Suzuki A (1997)<<BR>>[[http://www.ncbi.nlm.nih.gov/pubmed/9282479|Localization of brain electrical activity  via  linearly constrained minimum variance spatial filtering]]<<BR>>IEEE   Transactions on Biomedical Engineering 44(9): 867-880.
 1. Vrba J, Robinson SE (2000)<<BR>>Differences between synthetic aperture  magnetometry (SAM) and linear beamformers<<BR>> Biomag 2000: 681-684.
 1. Chen YS, Cheng CY, Hsieh JC, Chen LF (2006)<<BR>>[[http://www.ncbi.nlm.nih.gov/pubmed/16941832|Maximum contrast beamformer for electromagnetic  mapping of brain activity]]<<BR>>IEEE Transactions on Biomedical Engineering  53(9): 1765-1774.

== Feedback ==
<<EmbedContent(http://neuroimage.usc.edu/bst/get_feedback.php?Tutorials/Beamformers)>>