[FieldTrip] Beamformer and two different datasets

[Tyler Grummett]

Hello fieldtrip,


As mentioned in my previous email, I had success at calculating beamformer with one dataset but not with another.


The dropbox link to dataset1 is:

https://www.dropbox.com/s/2nyps8pph7xszf0/Dataset1.mat


The dropbox link to dataset2 is:

https://www.dropbox.com/s/pkmkdv871y4w67z/Dataset2.mat



In the datasets are structured in the following way:

datasetx.data

datasetx.timelock

datasetx.vol

datasetx.sourcemodel

datasetx.grid

datasetx.virtualchans

datasetx.sourcemodel2


source wasnt included as it will make the file too big.


The following code was used:


-------------------------------------------------------------

%% timelock data
cfg = [];
cfg.channel = 'EEG';
cfg.vartrllength = 2;
cfg.covariance = 'yes';
cfg.covariancewindow = 'all';
cfg.keeptrials = 'yes';
timelock = ft_timelockanalysis(cfg, data);


-------------------------------------------------------------

%% create headmodel
% segment MRI (return probabilistic tissue maps of gray/white/csf
% compartments
cfg = [];
cfg.write = 'no';
cfg.coordsys = 'spm';
cfg.output = { 'scalp', 'skull', 'brain'};
segmentedmri = ft_volumesegment(cfg, mri);

cfg = [];
cfg.method = 'iso2mesh';
cfg.numvertices = 10000;
bnd = ft_prepare_mesh( cfg, segmentedmri);

% fix mesh
[ bnd( 1).pnt, bnd( 1).tri] = meshresample( bnd( 1).pnt, bnd( 1).tri, 1000/size( bnd( 1).pnt, 1));
[ bnd( 2).pnt, bnd( 2).tri] = meshresample( bnd( 2).pnt, bnd( 2).tri, 2000/size( bnd( 2).pnt, 1));
[ bnd( 3).pnt, bnd( 3).tri] = meshresample( bnd( 3).pnt, bnd( 3).tri, 3000/size( bnd( 3).pnt, 1));
for ii = 1:size( bnd),
    [ bnd( ii).pnt, bnd( ii).tri] = meshcheckrepair( bnd( ii).pnt, bnd( ii).tri, 'dup');
    [ bnd( ii).pnt, bnd( ii).tri] = meshcheckrepair( bnd( ii).pnt, bnd( ii).tri, 'isolated');
    [ bnd( ii).pnt, bnd( ii).tri] = meshcheckrepair( bnd( ii).pnt, bnd( ii).tri, 'deep');
    [ bnd( ii).pnt, bnd( ii).tri] = meshcheckrepair( bnd( ii).pnt, bnd( ii).tri, 'meshfix');
end

% calculate headmodel % reordered to brain skull scalp
cfg = [];
cfg.method = 'bemcp'; %openmeeg doesnt work with multiple output from ft_volumesegment
vol = ft_prepare_headmodel(cfg, bnd);
clear bnd

-------------------------------------------------------------
%% calculate sourcemodel
cfg = [];
cfg.mri = mri;
cfg.vol = vol;
cfg.grid.warpmni = 'yes';
cfg.grid.template = template.sourcemodel;
cfg.grid.nonlinear = 'yes';
cfg.moveinward = 1; % actually uses vol mesh
cfg.inwardshift = 0; % needs to be expressed to work with moveinward
​cfg.elec = timelock.elec;
sourcemodel = ft_prepare_sourcemodel( cfg);

-------------------------------------------------------------
%% beamformer calculation

% compute lead field
cfg = [];
cfg.elec = timelock.elec;
cfg.vol = vol;
cfg.grid = sourcemodel;
cfg.reducerank = 3; % 3 for EEG, 2 for MEG
cfg.backproject = 'yes';
cfg.normalize = 'yes'; % if you are not contrasting the activity of interest again another condition or baseline time-window
grid = ft_prepare_leadfield( cfg, timelock);

% Source Analysis: without contrasting condition
cfg = [];
cfg.channel = 'EEG';
cfg.method = 'lcmv';
cfg.grid = grid;
cfg.vol = vol;
cfg.keepfilter = 'yes';
cfg.lcmv.fixedori = 'yes'; % project on axis of most variance using SVD
source = ft_sourceanalysis( cfg, timelock);

-------------------------------------------------------------
%% map beamformer source locations onto an anatomical label in an atlas
cfg = [];
cfg.interpmethod = 'nearest';
cfg.parameter = 'tissue';
sourcemodel2 = ft_sourceinterpolate( cfg, atlas, sourcemodel);

-------------------------------------------------------------
%% compute virtual channels
% start building vchan
vchan = [];
label = lower( atlas.tissuelabel);
label = label( 1:90);
vchan.time = data.time;
vchan.fsample = data.fsample;
Ntr = numel( data.trial);
vchan.trial = cell( 1, Ntr);

% find sensor names and indices
chans = ft_channelselection( 'EEG', data.label);
chans = match_str( data.label, chans);

count = 1;
tic
for i = 1:numel( label),
    atlas_sources = find( sourcemodel2.tissue == i);
    ai = ismember( atlas_sources, find( sourcemodel.inside));
    bregion_sources = atlas_sources( ai); clear atlas_sources
    if isempty( bregion_sources),
        continue;
    end
    for f = 1:numel( bregion_sources),
        source_inx = bregion_sources( f);
        dipole_data = cell( 1, Ntr);

        % multiply spatial filter (3,Nchan) by the original data
        if isempty( source.avg.filter{ source_inx}), continue; end
        for tr = 1:Ntr,
            dipole_data{ tr} = source.avg.filter{ source_inx} * data.trial{ tr}(chans,:);
        end

        % concatenate data, run svd on data, multiple data by the
        % orientation of the dipole in which it is strongest
        time_series = cat( 2, dipole_data{ :});
        [ U1, ~, ~] = svd( time_series, 'econ'); % u is the spatial decomposition, v the temporal and s the eigenvalues         along diagonal
        for tr = 1:Ntr,
            %             tt.trial{ tr}( f, :) = U1( :, 1)' * dipole_data{ tr};
            tt.trial{ tr}( f, :) = dipole_data{ tr};
        end
        clear source_inx dipole_data U1 timeseries
    end

    % mean channels with brain region
    for tr = 1:Ntr,
        vchan.trial{ tr}( i, :) = mean( tt.trial{ tr});
    end

    % include position and power for each source
    vchan.label( count) = label( i);
    fprintf( 'created virtual channel %d\n', count);
    count = count + 1;
    clear tt U S sv si temp_data bregion_sources bregion_source
end

cfg = [];
vchan = ft_preprocessing( cfg, vchan);

-------------------------------------------------------------​


I will greatly appreciate the help once again. As beamformer is the basically the key element of my Phd I really want it to get it working.


Kind regards,


Tyler


*************************

Tyler Grummett ( BBSc, BSc(Hons I))
PhD Candidate
Brain Signals Laboratory
Flinders University
Rm 5A301
Ext 66124
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20140711/9b913899/attachment.html>