[FieldTrip] Help understanding the lcmv beamformer output

[Lorenzo Magazzini]

Hi Philipp,

Thanks so much for your answer, it does indeed solve the problem!

I thought the covariance matrix was no longer necessary after computation
of the (common) filter, but, as you showed by pointing to the code, this is
clearly not the case.

Now, I compared the results of
a) calculating power from the virtual sensors time-series (data multiplied
by the filter, estimated using the 'common' covariance), and
b) power taken from the output of ft_sourceanalysis, using two different
'condition-specific' covariance matrices.
Results are pretty identical in my case, apart from minor differences in
the power values.

I wonder if anyone could help me figure out if there is a mathematical
difference between these two approaches?

Anyway, many thanks again for your help!

Best wishes,

Lorenzo



> *From: *Philipp Ruhnau <mail at philipp-ruhnau.de>
> *Subject: **Re: [FieldTrip] Help understanding the lcmv beamformer output*
> *Date: *13 January 2016 at 16:46:09 GMT+1
> *To: *FieldTrip discussion list <fieldtrip at science.ru.nl>
>
> 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/7eba1bb3/attachment-0002.html>