[FieldTrip] Beamformer and two different datasets
Tyler Grummett
tyler.grummett at flinders.edu.au
Fri Jul 11 04:17:27 CEST 2014
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-0001.html>
More information about the fieldtrip
mailing list