[FieldTrip] Help understanding the lcmv beamformer output

[Philipp Ruhnau]

Dear Lorenzo and others,
I’ve just looked into this, can reproduce the issue (same result with two different conditions/activity-baseline) 

The issue arises because power is calculated based on the filter and the covariance of the input data (following VanVeen or a slight change for precomputed filters), yet in the example below no covariance is precomputed for activity/baseline timelock structures (only for the common dataset) and in that case fieldtrip creates an identity matrix for both activity and baseline (around line 713 in ft_sourceanalysis). 
there is actually a warning about this, yet just from this warning (No covariance matrix found - will assume identity covariance matrix (mininum-norm solution)) it might to hard to get that this explains the results:

in essence from what I understand, the input data is ignored in the formula (for those interested this is line 278 in beamformer_lcmv) to compute the time series power, because the covariance is an identity matrix, and only the filter determines what ends up in the .pow field. and as the filter is the same for activity and baseline, the difference between the .pow output is zero.

now, computing covariances for the to-be-beamed segments (add a cfg.covariance = ‘yes’ to the timelock step for baseline and activity in the example) produces different power values (just tested with a simple auditory ERF dataset). 

I suppose this is desired behavior (the identity matrix covariance), otherwise I would suggest to produce an error. 

how well the ERF/timeseries is reconstructed of course depends on the data (time window) you use to create the filters. still it seems to me using the .mom field or creating virtual sensors yields better (in a sense of more focused) results. I still wonder why virtual sensors/.mom averages and .pow results are quite different, maybe one of the more knowledgeable ft developers could comment on at some point, but this might exceed the limit of this list…

cheers
philipp 

> On 12 Jan 2016, at 15:46, Lorenzo Magazzini <lorenzomagazzini at gmail.com> wrote:
> 
> Hi,
> 
> I hope someone can help me to understand the use of the lcmv beamformer, because either I am using it incorrectly, or I misunderstood the meaning of its output.
> 
> I am trying to get source-level power estimates of 'baseline' and 'active' time-windows. I first compute a common filter, and then apply it separately to the two conditions. I thought the power estimates would be stored, voxel-wise, in the '.avg.pow' field of the ft_sourceanalysis output. However, the values I get in the '.pow' field for the analysis of 'baseline' and 'active' windows are identical. Furthermore, the '.pow' values are identical to the '.noise' values when using cfg.projectnoise = 'yes';.
> 
> I have checked if the sensor-level data that passed on to ft_sourceanalysis is different between baseline and active - it is.
> 
> Also, if I do the following...
> 1. reconstruct virtual sensors by applying the common filter separately to baseline and active, 
> 2. average the single-trial (virtual sensor) power estimates, 
> 3. replace the '.pow' value for each voxel, 
> ...then the contrast between baseline and active is correct.
> 
> Can anyone please help me figure out what I'm doing wrong here? Below I copied part of the code.
> 
> Many thanks in advance!
> 
> Lorenzo
> 
> 
> %compute the leadfield
> cfg                 = [];
> cfg.channel         = {'MEG'};
> cfg.grid            = sourcemodel;
> cfg.vol             = hdm;
> cfg.grad            = grad;
> cfg.normalize       = 'yes';
> leadfield = ft_prepare_leadfield(cfg);
> 
> %compute covariance matrix
> cfg = [];
> cfg.channel = {'MEG'};
> cfg.removemean = 'no';
> cfg.covariance = 'yes';
> cfg.covariancewindow = [-1.5 1.5];
> tlck = ft_timelockanalysis(cfg, data);
> 
> %get the common spatial filter
> cfg = [];
> cfg.method          = 'lcmv';
> cfg.grid            = leadfield;
> cfg.vol             = hdm;
> cfg.grad            = grad;
> cfg.lcmv.fixedori   = 'yes';
> cfg.lcmv.keepfilter = 'yes';
> cfg.lcmv.projectnoise = 'yes';
> cfg.lcmv.lambda     = '5%';
> src = ft_sourceanalysis(cfg, tlck);
> 
> %store the filter
> wts = src.avg.filter;
> 
> %separate baseline and active windows of bandpass-filtered data
> cfg = [];
> cfg.toilim = [-1.5 -0.3];
> data_bsln = ft_redefinetrial(cfg, data);
> cfg = [];
> cfg.toilim = [0.3 1.5];
> data_actv = ft_redefinetrial(cfg, data);
> 
> %timelock analysis of baseline and active
> cfg             = [];
> cfg.channel     = {'MEG'};
> cfg.keeptrials  = 'yes';
> tlck_bsln = ft_timelockanalysis(cfg, data_bsln);
> tlck_actv = ft_timelockanalysis(cfg, data_actv);
> 
> %get source-power estimates of bsln and actv separately using the pre-computed common filter (wts)
> cfg                 = [];
> cfg.method          = 'lcmv';
> cfg.vol             = hdm;
> cfg.grid            = leadfield; 
> cfg.grid.filter     = wts;
> % cfg.keepfilter      = 'no';
> cfg.lcmv.fixedori   = 'yes';
> cfg.projectnoise    = 'yes';
> cfg.lambda          = '5%';
> cfg.keeptrials      = 'yes';
> src_bsln = ft_sourceanalysis(cfg, tlck_bsln);
> src_actv = ft_sourceanalysis(cfg, tlck_actv);
> 
> 
> _______________________________________________
> 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/20160113/74387d30/attachment.html>