Hi Fieldtrippers,<br><br>I have a slight issue that I was hoping someone can help with. I have some CTF data, that I was trying to perform DICS beamforming and then source interpolation on. I am trying to beamform to the source of gamma frequency in the right hemisphere, comparing pre and post stimuli. Heres the basic script:<br>
<br><b>People can skip the first part of code, I included it as possibly useful information for this plead <br></b><br>%%%%% Data has already been preprocessed and timelocked%%%%%%<br><br> load('processed_data_500Hz.mat')<br>
load('timelock_data_500Hz.mat')<br><br>%%%% following the beamforming tutorial %%%%%%%<br><br>%%%%%%%%%%%%%split and combine%%%%%%%%%%%%% <br> <br> cfg = []; <br>
cfg.toilim = baseline; <br> cfg.channel=channel<br> total_bsl = ft_redefinetrial(cfg, data);<br> <br> <br> cfg = []; <br> cfg.toilim = interest; <br>
cfg.channel=channel<br> total_interest = ft_redefinetrial(cfg, data);<br> <br> cfg = [];<br> total_cmb = ft_appenddata(cfg, total_bsl, total_interest);<br> <br> total_cmb.trialinfo = [zeros(length(total_bsl.trial), 1); ones(length(total_interest.trial), 1)];<br>
<br><br> %%%%%%%%%%%%%cross spectrial density for dics and tf for common filter%%%%%%%%%%%%%%<br> <br>cfg = [];<br>cfg.method = 'mtmfft';<br>cfg.output = 'powandcsd';<br>cfg.tapsmofrq = 25;<br>
cfg.foilim = 55;<br>freqPre_total = ft_freqanalysis(cfg, total_bsl);<br><br>cfg = [];<br>cfg.method = 'mtmfft';<br>cfg.output = 'powandcsd';<br>cfg.tapsmofrq = 25;<br>cfg.foilim = 55;<br>freqPost_total = ft_freqanalysis(cfg, total_interest);<br>
<br>cfg = [];<br>cfg.method = 'mtmfft';<br>cfg.output = 'powandcsd';<br>cfg.tapsmofrq = 25;<br>cfg.foilim = 55;<br>freqALL_total = ft_freqanalysis(cfg, total_cmb);<br><br><br>%%%%%%%%%%%%%%%%%%%%%%%% calc leadfield %%%%%%%%%%%%%%%%%%%%%%%%%%%<br>
<br>% <br>T=ft_read_mri('R128_V2.mri')<br><br> cfg.fiducial.nas = T.hdr.fiducial.mri.nas<br> cfg.fiducial.lpa = T.hdr.fiducial.mri.lpa<br> cfg.fiducial.rpa = T.hdr.fiducial.mri.rpa<br> cfg.method = 'fiducial'<br>
cfg.coordsys = 'ctf' <br> <br> [mri] = ft_volumerealign(cfg, T)<br> <br>cfg = [];<br>cfg.coordsys = 'ctf'; % the MRI is expressed in the CTF coordinate system, see below<br>segmentedmri = ft_volumesegment(cfg, mri);<br>
<br>cfg = [];<br>cfg.method = 'singleshell';<br>hdm = ft_prepare_headmodel(cfg, segmentedmri);<br><br> vol = ft_convert_units(hdm, 'cm');<br><br>%%%%%%%%%%%%%%%<br><br><br>cfg = [];<br>
cfg.grad = freqPost_total.grad;<br>cfg.vol = vol;<br>cfg.reducerank = 2;<br>cfg.channel = 'MR'<br>[grid] = ft_prepare_leadfield(cfg);<br><br><br><br><b>Upto this point everything works. My volume conduction model, grads, and grid looks good when verified visually...<br>
<br>I then compute the common filter and apply it to the beamforming of each condition</b><br><br><br>%%%%%%%%%%%%%%%%%%%%% computer common filter%%%%%%%%%%%<br><br>cfg = [];<br>cfg.method = 'dics';<br>
cfg.grad = freqALL_total.grad;<br>cfg.keeptrials = 'yes';<br>cfg.frequency = freqALL_total.freq;<br>cfg.grid = grid;<br>cfg.vol = vol;<br>cfg.dics.projectnoise = 'yes';<br>
cfg.dics.lambda = '5%';<br>cfg.dics.keepfilter = 'yes';<br>cfg.dics.realfilter = 'yes';<br>cfg.channel = 'MR'<br>sourceAll = ft_sourceanalysis(cfg, freqALL_total);<br><br>
%%%%%%%%%%%%%%%%%% source for bsl + interest %%%%%%%%%%%%%<br><br>cfg.grid.filter = sourceAll.avg.filter;<br>cfg.grid = grid;<br>sourcePre_con = ft_sourceanalysis(cfg, freqPre_total );<br>sourcePost_con = ft_sourceanalysis(cfg, freqPost_total);<br>
<br><br>%%%%%%%%%% Calc difference %%%%%%%%%%%%%%%%%<br><br>sourceDiff = sourcePost_con;<br>sourceDiff.avg.pow = (sourcePost_con.avg.pow - sourcePre_con.avg.pow) ./ sourcePre_con.avg.pow;<br><br><br><br><b>Heres where it falls apart, because I try to source interpolate</b><br>
<br><br>%%%%%%%%%%%%%%%%%%%%% graph time %%%%%%%%%%%%%%%%%%<br><br>% slicemri = ft_volumereslice([], mri);<br>% % mri = ft_volumereslice([],segmentedmri) %%<b> didnt appear to help with the current problem<br></b><br>cfg = [];<br>
cfg.downsample = 2;<br>cfg.parameter = 'avg.pow';<br>cfg.coordsys = 'ctf';<br>cfg.channel = 'MR'<br>source_diffInt = ft_sourceinterpolate(cfg, sourceDiff , mri); &&& I also tried segmentedmri, but no difference<br>
<br>It works until it reachs line 237; tmp = interpmat*tmp; where there is the matrixs are mismatch (interpmat is a 16777216(256x256x256) x642 matrix, and tmp is a 1 * 642 matrix. From my basic skills I gather that tmp is sourceDiff.avg.pow (the power difference for each leadfield point). interpmat is the interpolation matrix return by<br>
<br>interpmat = interp_ungridded(functional.pos <b>(position of the leadfield locations)</b>, pos (mri voxels), 'projmethod', cfg.interpmethod (<b>nearest</b>), 'sphereradius', cfg.sphereradius) (0.<b>05</b>);<br>
<br>what I get out of this is a 16777216(<b>256x256x256</b> voxels) x642 <b>(leadfield position)</b> matrix<b>. </b>it then undergoes interpmat(~anatomical.inside(:), :) = 0. My <b>anatomical.inside is all voxels (all of the 256 x 256 x 256 matrix= 1), as a result of the defaults in ft_checkdata</b>.<br>
<br>When I try to find the minimum value in this matrix I get (via min(interpmat) = All zero sparse: 1-by-642<br><br>I think what is supposed to happen is that the interpmat is a matrix that relates meg to mri space, and as such then mutlipling with a vector that represents the power in meg space, I can see how that would get points of power in mri space. Though this must be wrong since the math is not allowed.... I cannot see where I went off the beaten track. It seems as if other people have not had this problem. <br>
<br>Sorry for the bother,<br><br>Russ Port<br><br>Graduate Student @ Upenn<br><br><br>P.S. thanks to anyone who tries to help in advance...<br><br><br><br><br><br><br><br><br>