[FieldTrip] Advice on using LCMV beamformer
Stephen Whitmarsh
stephen.whitmarsh at gmail.com
Fri Jul 17 11:06:53 CEST 2015
Hi everyone. Thanks so much for your responses so far. They are very
clarifying.
Anyway, after browsing the mailinglist archive I found two related Q&A:
This clears up the lack of differences, because the covariance is
calculated during ft_timelockanalysis, and changing the latency/timewindow
later has no effect:
http://mailman.science.ru.nl/pipermail/fieldtrip/2009-February/001965.html
When it comes to the common filter etc. this conversation was somewhat
helpful, at least how to compute it:
http://mailman.science.ru.nl/pipermail/fieldtrip/2013-January/006120.html
Although these problems are solved now and I have "results", they are still
confusing. They do not look totally out-there, but they are certainly not
unequivocal. After trying out different options I see that the covariance
timewindow of the common filter makes a big difference. It seems I get the
best (or at least better) results when this is restricted to the
post-stimulus window in the time-period of component I am interested in.
See attachment - blue brains.
However, before I go any further I am somewhat on my toes because
projecting the power of the active or baseline period seems strangely
shifted to the anterior. See attachment - red brains.
Your comments have been very helpful thinking about the common filter. I
actually now wonder whether it is a good idea at all in my case. One other
thing could try is to somehow calculate a common filter by concatenating
the baseline period and the active component period before doing
calculating the covariance.
Or, alternatively, I could calculate the common filter over the whole
post-stimulus interval (of say 300ms instead of 30ms of my component of
interest), similarly as Johanna has done here:
http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2899153/#R28
Best and thanks,
Stephen
p.s. I attached my script.
On 17 July 2015 at 09:46, Jörn M. Horschig <jorn at artinis.com> wrote:
> whoopsy, saw your code just now - so I take that back ;)
>
>
>
> *--*
>
>
>
> *Jörn M. Horschig, PhD*, Software Engineer
>
> Artinis Medical Systems <http://www.artinis.com/> | +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
> ratio)
>
> 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
> ratio)
>
> 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
> option?
>
> 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,
>
> Stephen
>
>
> % 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,
> ERF_single{isubject});
> source_lcmv_double{isubject} = ft_sourceanalysis(cfg,
> ERF_double{isubject});
>
> % 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,
> ERF_single_latency{isubject});
>
> % 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,
> ERF_single_baseline{isubject});
>
> % 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 -
> source_lcmv_single_baseline{isubject}.avg.pow;
>
> template_mri =
> ft_read_mri(['/opt/fieldtrip/external/spm8/templates/T1.nii']);
>
> 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});
> ft_sourceplot(cfg,source_lcmv_single_latency_int{isubject});
>
> cfg = [];
> cfg.method = 'slice';
> cfg.funparameter = 'pow';
> cfg.funcolorlim = [-1e-23 1e-23];
> ft_sourceplot(cfg,source_lcmv_single_diff_int{isubject});
>
>
>
>
>
>
> _______________________________________________
> fieldtrip mailing list
> fieldtrip at donders.ru.nl
> http://mailman.science.ru.nl/mailman/listinfo/fieldtrip
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20150717/7b7da298/attachment-0002.html>
-------------- next part --------------
function metacognition_lcmv(isubject,covlatency,latency,baseline)
addpath /home/stephen/analysis/metacognition/scripts/
addpath /home/stephen/analysis/metacognition/trialfuns/
% addpath /opt/fieldtrip/
% ft_defaults
metacognition_subjectinfo;
% get subject filename details
flist = metacognition_filelist(isubject);
% load raw cleaned data (fixing some changes of filenames)
temp = load(['/home/stephen/analysis/metacognition/preprocessed/pp_' num2str(isubject) '_data_stim.mat']);
if isfield(temp,'data')
clean_data = temp.data;
end
if isfield(temp,'clean_data')
clean_data = temp.clean_data;
end
clear temp
% get header
hdr = ft_read_header(flist{1});
clean_data.grad = hdr.grad;
clean_data.elec = hdr.elec;
% baselinecorrect
cfg = [];
cfg.demean = 'yes';
cfg.baselinewindow = [-0.3 0];
clean_data = ft_preprocessing(cfg,clean_data);
% 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 using COV of whole timeseries
cfg = [];
cfg.covariance = 'yes';
cfg.covariancewindow = covlatency;
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);
cfg = [];
cfg.method = 'lcmv';
cfg.grid = leadfield;
cfg.vol = headmodel_meg{isubject};
cfg.lcmv.keepfilter = 'yes';
cfg.lcmv.lambda = '5%';
cfg.senstype = 'meg';
cfg.channel = 'meggrad';
cfg.lcmv.fixedori = 'yes'; % use only right gradiometers
cfg.reducerank = 2;
source_lcmv_single{isubject} = ft_sourceanalysis(cfg, ERF_single{isubject});
% source_lcmv_double{isubject} = ft_sourceanalysis(cfg, ERF_double{isubject});
% project post stim latency with previous common filter
cfg = [];
cfg.covariance = 'yes';
cfg.covariancewindow = latency;
cfg.trials = find(clean_data.trialinfo(:,2) == 1);
ERF_single_active{isubject} = ft_timelockanalysis(cfg,clean_data);
% cfg.trials = find(clean_data.trialinfo(:,2) == 2);
% ERF_double_active{isubject} = ft_timelockanalysis(cfg,clean_data);
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_active{isubject} = ft_sourceanalysis(cfg, ERF_single_active{isubject});
% project baseline latency with previous common filter
cfg = [];
cfg.covariance = 'yes';
cfg.covariancewindow = baseline;
cfg.trials = find(clean_data.trialinfo(:,2) == 1);
ERF_single_baseline{isubject} = ft_timelockanalysis(cfg,clean_data);
% cfg.trials = find(clean_data.trialinfo(:,2) == 2);
% ERF_double_baseline{isubject} = ft_timelockanalysis(cfg,clean_data);
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, ERF_single_baseline{isubject});
% get difference between baseline and component
source_lcmv_single_diff{isubject} = source_lcmv_single_active{isubject};
source_lcmv_single_diff{isubject}.avg.pow = source_lcmv_single_active{isubject}.avg.pow - source_lcmv_single_baseline{isubject}.avg.pow;
template_mri = ft_read_mri(['/opt/fieldtrip/external/spm8/templates/T1.nii']);
cfg = [];
cfg.voxelcoord = 'no';
cfg.parameter = 'pow';
cfg.interpmethod = 'nearest';
% source_lcmv_single_int{isubject} = ft_sourceinterpolate(cfg, source_lcmv_single{isubject}, template_mri);
source_lcmv_single_active_int{isubject} = ft_sourceinterpolate(cfg, source_lcmv_single_active{isubject}, template_mri);
source_lcmv_single_baseline_int{isubject} = ft_sourceinterpolate(cfg, source_lcmv_single_baseline{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 = '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});
ft_sourceplot(cfg,source_lcmv_single_active_int{isubject});
ft_sourceplot(cfg,source_lcmv_single_baseline_int{isubject});
cfg = [];
cfg.method = 'slice';
cfg.funparameter = 'pow';
% cfg.maskparameter = cfg.funparameter;
% cfg.funcolorlim = [-1e-23 1e-23];
ft_sourceplot(cfg,source_lcmv_single_diff_int{isubject});
cfg = [];
cfg.method = 'surface';
cfg.funparameter = 'pow';
cfg.maskparameter = cfg.funparameter;
cfg.projmethod = 'nearest';
cfg.surffile = 'surface_white_both.mat';
cfg.surfdownsample = 10;
% cfg.funcolorlim = [-1e-23 1e-23];
ft_sourceplot(cfg,source_lcmv_single_diff_int{isubject});
view([90 15])
1 == 2
-------------- next part --------------
A non-text attachment was scrubbed...
Name: SEF_lcmv_diff_active.jpg
Type: image/jpeg
Size: 28082 bytes
Desc: not available
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20150717/7b7da298/attachment-0006.jpg>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: SEF_lcmv_diff_baseline.jpg
Type: image/jpeg
Size: 26299 bytes
Desc: not available
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20150717/7b7da298/attachment-0007.jpg>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: SEF_lcmv_diff_from_baseline.jpg
Type: image/jpeg
Size: 26823 bytes
Desc: not available
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20150717/7b7da298/attachment-0008.jpg>
More information about the fieldtrip
mailing list