[FieldTrip] Advice on using LCMV beamformer

Jörn M. Horschig jorn at artinis.com
Fri Jul 17 09:42:08 CEST 2015

Hi Stephen,


Approach 1 and 3 would be fine, 2 would be not a wise thing to do, as you
get different noise estimates for the two periods (presumably). So beamform
all your data that will go into the final analysis. Instead of computing
some noise estimate you can also normalize the leadfields. According to
Robert that should also work fairly fine, though I haven’t tested that
myself. I would  go for beamforming the baseline and the activation period
(that’s because that’s what I already did earlier).


For getting the time courses, Eelke and me made an extended beamforming
tutorial based on the existing tutorials, which featured virtual channel
extraction. I just checked, and it can be found here:



So basically, yes, manually multiply the data with the filter. I guess that
should get you started. 


btw, not sure what you are doing wrong – maybe share the code you’re using?
Or check the extended beamforming tutorial, there ft_redefinetrial is used:


Hope this helps!








Jörn M. Horschig, PhD, Software Engineer

 <http://www.artinis.com/> Artinis Medical Systems  |  +31 481 350 980 


From: fieldtrip-bounces at science.ru.nl
[mailto:fieldtrip-bounces at science.ru.nl] On Behalf Of Stephen Whitmarsh
Sent: Friday, July 17, 2015 12:04 AM
To: fieldtrip, donders
Subject: [FieldTrip] Advice on using LCMV beamformer


Hi everyone,

I hope there is someone out there who has some time to help me with using
LCMV beamformer. I have no experience with it yet, although a little with
DICS, and understand they are at least conceptually very similar. However, I
am running into some questions. Before I go on, my data on sensor level
looks good and it works out with dipolefitting as well (attached).

I would like to beamform an early (unilaterally stimulated) sensory evoked
field. There is no conditional difference, I would just like to extract the
timecourse of some maximum voxels. In other words:

1) get a contralateral cluster or a number of voxels with maximum amplitude.

2) extract the timecourses of those clusters, then average them.

For step 1, it seem to have several options, either:

1a) calculate the filter on the whole timecourse (a common filter as would
be used in DICS)

1b) use that common filter to calculate the amplitude of the early component

1c) do the same for the baseline period

1d) subtract the baseline activation from the active period (or take a

Alternatively, I could:

2a) beamform the active period

2b) beamform the baseline period

2c) subtract the baseline activation from the active period (or take a

Or it seems, I might even:

3a) beamform the active period

3b) divide the active period by the projected noise

Do I get this correctly? Does anyone have an argument for one of these

Anyway, first practical problem is that when I plot either the whole, the
baseline or the activation period I get the same results (attached). A
difference therefore doesn't work, obviously. So this raises the question
how the selection of latency works in ft_sourceanalysis? I trust cfg.latency
works, but it doesn't seem to make a difference for me. At the same time,
using selectdata to select a time period doesn't make a difference either.
I'm obviously doing something wrong, but what?


Anyway, when I get this figured out, I still need to extract the
timecourses. How would I proceed doing this? I imagine the fastest way would
not be to loop ft_sourceanalysis for each timepoint (using a common filter I
would think). Should I 'manually' multiply the data with the filter?

It would be great to have some feedback on the general idea, and perhaps
someone can identify what I am doing wrong (script below). I would be happy
to finish the LCMV beamformer tutorial on the FT website once I get it
figured out, as I can imagine more people would have these questions.

All the best and thanks,


% timelock average + covariance

cfg = [];
cfg.covariance = 'yes';
cfg.covariancewindow = 'all';
cfg.trials = find(clean_data.trialinfo(:,2) == 1);
ERF_single{isubject} = ft_timelockanalysis(cfg,clean_data);
cfg.trials = find(clean_data.trialinfo(:,2) == 2);
ERF_double{isubject} = ft_timelockanalysis(cfg,clean_data);

% combine gradiometers (single stim)
cfg = [];
ERF_single_cmb{isubject} = ft_combineplanar(cfg,ERF_single{isubject});
ERF_double_cmb{isubject} = ft_combineplanar(cfg,ERF_double{isubject});

% load headmodel MEG
temp = load(['/home/stephen/analysis/metacognition/sourcemodel/pp_'
num2str(isubject) '_headmodel_singleshell.mat']);
headmodel_meg{isubject} = temp.headmodel_meg;

% headmodel_meg = ft_convert_units(headmodel_meg, 'cm');
gridtemp = load(['/home/stephen/analysis/metacognition/sourcemodel/pp_'
num2str(isubject) '_grid.mat']);

hdr             = ft_read_header(flist{1});
cfg             = [];
cfg.grad        = hdr.grad;
cfg.vol         = headmodel_meg{isubject};
cfg.grid        = gridtemp.grid;
cfg.channel     = {'MEGGRAD'};
cfg.normalize   = 'yes';
leadfield       = ft_prepare_leadfield(cfg);

% get the filter
cfg = [];
cfg.method              = 'lcmv';
cfg.grid                = leadfield;
cfg.vol                 = headmodel_meg{isubject};
cfg.keepfilter          = 'yes';
cfg.senstype            = 'meg';
cfg.channel             = 'meggrad';
cfg.lcmv.fixedori       = 'yes'; % use only right gradiometers
cfg.reducerank          = 2;
source_lcmv_single{isubject}   = ft_sourceanalysis(cfg,
source_lcmv_double{isubject}   = ft_sourceanalysis(cfg,

% project post stim latency with previous common filter

cfg = [];
cfg.latency = latency;
ERF_single_latency{isubject} = ft_selectdata(cfg,ERF_single{isubject});

cfg = [];
cfg.method              = 'lcmv';
cfg.grid                = leadfield;
cfg.vol                 = headmodel_meg{isubject};
cfg.grid.filter         = source_lcmv_single{isubject}.avg.filter;
cfg.channel             = 'meggrad';
cfg.senstype            = 'meg';
% cfg.latency             = latency;
cfg.lcmv.fixedori       = 'yes'; 
cfg.reducerank          = 2;
source_lcmv_single_latency{isubject} = ft_sourceanalysis(cfg,

% baseline

cfg = [];
cfg.latency = [-0.2 0];
ERF_single_baseline{isubject} = ft_selectdata(cfg,ERF_single{isubject});

cfg = [];
cfg.method              = 'lcmv';
cfg.grid                = leadfield;
cfg.vol                 = headmodel_meg{isubject};
cfg.grid.filter         = source_lcmv_single{isubject}.avg.filter;
cfg.channel             = 'meggrad';
cfg.senstype            = 'meg';
% cfg.latency             = latency;
cfg.lcmv.fixedori       = 'yes'; 
cfg.reducerank          = 2;
source_lcmv_single_baseline{isubject} = ft_sourceanalysis(cfg,

% get difference between baseline and component
source_lcmv_single_diff{isubject} = source_lcmv_single_latency{isubject};
source_lcmv_single_diff{isubject}.avg.pow =
source_lcmv_single_latency{isubject}.avg.pow -

template_mri =

cfg = [];
cfg.voxelcoord = 'no';
cfg.parameter = 'avg.pow';
cfg.interpmethod = 'nearest';
source_lcmv_single_int{isubject} = ft_sourceinterpolate(cfg,
source_lcmv_single{isubject}, template_mri);
source_lcmv_single_latency_int{isubject} = ft_sourceinterpolate(cfg,
source_lcmv_single_latency{isubject}, template_mri);
source_lcmv_single_diff_int{isubject} = ft_sourceinterpolate(cfg,
source_lcmv_single_diff{isubject}, template_mri);

cfg = [];
cfg.method = 'slice';
cfg.funparameter = 'avg.pow';
cfg.maskparameter = cfg.funparameter;
cfg.funcolorlim = [3e-22 5e-22];
% ft_sourceplot(cfg,source_lcmv_single_latency_int{isubject});
% ft_sourceplot(cfg,source_lcmv_single_int{isubject});

cfg = [];
cfg.method = 'slice';
cfg.funparameter = 'pow';
cfg.funcolorlim = [-1e-23 1e-23];



-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20150717/6a00e672/attachment-0002.html>

More information about the fieldtrip mailing list