sub = ['02';'04';'05';'06';'07';'09';'11';'13';'14';'15';'16';'17';'18';'19';'20';'21'] ; 
trialtype = ['3';'4';'5';'6'] ; 

Dir = 'C:\Users\EEG1\Documents\Data\FTtrials\';
mt = '.mat'; 
prefix = '1';

trialnum421 = ['02';'03';'04';'05';'06';'07';'08';'09';'10';'11';'12';'13';'14';'15';'16';'17';'18';'19';'20';'21';'22';'24'];
trialnum420 = ['01';'02';'03';'04';'05';'07';'08';'09';'10';'11';'12';'13';'14';'15';'16';'18';'19';'20';'21';'22';'23';'24';'25';'27';'28';'29';'30'];
trialnum419 = ['01';'02';'03';'04';'05';'06';'07';'08';'09';'10';'11';'12';'13';'14';'15';'16';'17';'18';'19';'20';'21';'22';'23';'24';'25';'26';'27';'28';'29';'30'];
trialnum418 = ['01';'02';'03';'04';'05';'06';'07';'08';'09';'10';'11';'12';'13';'14';'16';'17';'18';'19';'20';'21';'22';'23';'24';'25';'26';'27';'28';'29';'30'];
trialnum405 = ['01';'02';'03';'04';'05';'06';'07';'08';'09';'10';'12';'14';'15';'16';'22';'23';'24';'25';'26';'28';'29'];
trialnum404 = ['02';'03';'04';'05';'06';'07';'09';'10';'11';'12';'13';'14';'15';'16';'19';'21';'22';'23';'24';'25';'26';'27';'30'];
trialnum402 = ['02';'04';'05';'06';'07';'08';'09';'10';'14';'15';'16';'17';'18';'19';'20';'23';'25';'26';'28'];
trialnum406 = ['01';'02';'03';'04';'05';'06';'07';'08';'09';'10';'11';'12';'13';'14';'15';'16';'17';'18';'19';'20';'21';'22';'23';'24';'25';'26';'27';'28';'29';'30'];
trialnum407 = ['01';'02';'03';'04';'05';'06';'08';'09';'10';'11';'12';'14';'16';'17';'18';'19';'20';'21';'22';'23';'24';'25';'26';'27';'28';'30'];
trialnum409 = ['01';'02';'03';'04';'06';'09';'10';'11';'12';'13';'14';'15';'16';'18';'20';'21';'22';'27'];
trialnum411 = ['01';'02';'03';'04';'05';'06';'07';'08';'09';'10';'11';'12';'13';'14';'15';'16';'17';'18';'19';'20';'21';'22';'23';'24';'25';'26';'27';'28';'29';'30'];
trialnum413 = ['01';'02';'03';'04';'05';'06';'07';'08';'09';'11';'13';'14';'15';'16';'17';'20';'21';'22';'23';'24';'25';'26';'27';'29'];
trialnum414 = ['02';'03';'04';'05';'06';'07';'08';'10';'11';'12';'13';'14';'15';'16';'17';'18';'19';'21';'23';'24';'25';'28';'29';'30'];
trialnum415 = ['01';'02';'03';'04';'05';'06';'07';'08';'09';'10';'11';'13';'14';'15';'16';'17';'18';'19';'21';'22';'23';'24';'26';'27';'28';'29';'30'];
trialnum416 = ['01';'02';'03';'04';'05';'06';'07';'08';'09';'11';'12';'13';'14';'15';'17';'18';'19';'20';'21';'23';'24';'25';'27';'28';'29';'30'];
trialnum417 = ['01';'02';'03';'04';'05';'06';'07';'08';'10';'11';'12';'13';'14';'15';'16';'17';'18';'19';'21';'22';'23';'24';'25';'26';'27';'28';'29';'30'];

   j = length(trialnum421) ;
   l = length(trialnum420) ;
   m = length(trialnum419) ; 
   n = length(trialnum418) ; 
   o = length(trialnum417) ;
   p = length(trialnum416) ; 
   q = length(trialnum415) ; 
   r = length(trialnum414) ; 
   s = length(trialnum413) ; 
   t = length(trialnum411) ; 
   u = length(trialnum409) ; 
   v = length(trialnum407) ; 
   w = length(trialnum406) ; 
   x = length(trialnum405) ;
   y = length(trialnum404) ;

dataFIC4 = [] ;
dataFIC4.sub = cell([1 16]);
sourceNAIInt4po = cell([1 16]);
sourceNAIInt4pe = cell([1 16]);
sourceNAIInt4pos = cell([1 16]);
sourceNAIInt4ze = cell([1 16]);
   
dataFIC4.cfg = [] ; 
dataFIC4.hdr = [] ; 
dataFIC4.trial = [] ; 
dataFIC4.time = [] ;
dataFIC4.grad = [] ;
dataFIC4.label = [] ; 
dataFIC4.hdr.Fs = 600 ; 
dataFIC4.hdr.nChans = 141 ; 
dataFIC4.hdr.nSamples = 5251 ; 
dataFIC4.hdr.nTrials = length(trialnum421) + length(trialnum420) + length(trialnum419) + length(trialnum418) + length(trialnum417) + length(trialnum416) + length(trialnum415) + length(trialnum414) + length(trialnum413) + length(trialnum411) + length(trialnum409) + length(trialnum407) + length(trialnum406) + length(trialnum405) + length(trialnum404) + length(trialnum402) ;  
dataFIC4.fsample = 600 ; 

for i=1:16
    dataFIC4.sub{i}.cfg = [] ; 
    dataFIC4.sub{i}.hdr = [] ; 
    dataFIC4.sub{i}.trial = [] ; 
    dataFIC4.sub{i}.time = [] ;
    dataFIC4.sub{i}.grad = [] ;
    dataFIC4.sub{i}.label = [] ; 
    dataFIC4.sub{i}.hdr.Fs = 600 ; 
    dataFIC4.sub{i}.hdr.nChans = 141 ; 
    dataFIC4.sub{i}.hdr.nSamples = 5251 ; 
    dataFIC4.sub{i}.fsample = 600 ; 
end

dataFIC4.sub{1}.hdr.nTrials = length(trialnum402);
dataFIC4.sub{2}.hdr.nTrials = length(trialnum404);
dataFIC4.sub{3}.hdr.nTrials = length(trialnum405);
dataFIC4.sub{4}.hdr.nTrials = length(trialnum406);
dataFIC4.sub{5}.hdr.nTrials = length(trialnum407);
dataFIC4.sub{6}.hdr.nTrials = length(trialnum409);
dataFIC4.sub{7}.hdr.nTrials = length(trialnum411);
dataFIC4.sub{8}.hdr.nTrials = length(trialnum413);
dataFIC4.sub{9}.hdr.nTrials = length(trialnum414);
dataFIC4.sub{10}.hdr.nTrials = length(trialnum415);
dataFIC4.sub{11}.hdr.nTrials = length(trialnum416);
dataFIC4.sub{12}.hdr.nTrials = length(trialnum417);
dataFIC4.sub{13}.hdr.nTrials = length(trialnum418);
dataFIC4.sub{14}.hdr.nTrials = length(trialnum419);
dataFIC4.sub{15}.hdr.nTrials = length(trialnum420);
dataFIC4.sub{16}.hdr.nTrials = length(trialnum421);

for i=1:length(trialnum421)
   %load the data for one individual, one trial type, one trial
   import = load([Dir prefix sub(16,:) trialtype(2) trialnum421(i,:) mt]);
   dataFIC4.trial{i} = import.avg(32:172,:) ; 
   dataFIC4.time{i} = import.time ; 
   dataFIC4.hdr.label = import.label(32:172)  ;
   dataFIC4.label = import.label(32:172)  ;
   dataFIC4.hdr.grad.chantype = import.grad.chantype(30:170) .' ; 
   dataFIC4.hdr.grad.chanpos = import.grad.chanpos(30:170,:) ;
   dataFIC4.hdr.grad.chanori = import.grad.chanori(30:170,:) ;
   dataFIC4.hdr.grad.tra = import.grad.tra(30:170,:) ;
   dataFIC4.hdr.grad.coilori = import.grad.coilori ; 
   dataFIC4.hdr.grad.coilpos = import.grad.coilpos ; 
   dataFIC4.hdr.grad.unit = import.grad.unit ; 
   dataFIC4.hdr.grad.label = import.grad.label(30:170) .' ; 
   dataFIC4.grad.chantype = import.grad.chantype(30:170) .' ; 
   dataFIC4.grad.chanpos = import.grad.chanpos(30:170,:) ;
   dataFIC4.grad.chanori = import.grad.chanori(30:170,:) ;
   dataFIC4.grad.tra = import.grad.tra(30:170,:) ;
   dataFIC4.grad.coilori = import.grad.coilori ; 
   dataFIC4.grad.coilpos = import.grad.coilpos ; 
   dataFIC4.grad.unit = import.grad.unit ; 
   dataFIC4.grad.label = import.grad.label(30:170) .' ; 
   
   dataFIC4.sub{16}.trial{i} = import.avg(32:172,:) ;
   dataFIC4.sub{16}.time{i} = import.time ;
   for i=1:16
    dataFIC4.sub{i}.hdr.label = import.label(32:172)  ;
    dataFIC4.sub{i}.label = import.label(32:172)  ;
    dataFIC4.sub{i}.hdr.grad.chantype = import.grad.chantype(30:170) .' ; 
    dataFIC4.sub{i}.hdr.grad.chanpos = import.grad.chanpos(30:170,:) ;
    dataFIC4.sub{i}.hdr.grad.chanori = import.grad.chanori(30:170,:) ;
    dataFIC4.sub{i}.hdr.grad.tra = import.grad.tra(30:170,:) ;
    dataFIC4.sub{i}.hdr.grad.coilori = import.grad.coilori ; 
    dataFIC4.sub{i}.hdr.grad.coilpos = import.grad.coilpos ; 
    dataFIC4.sub{i}.hdr.grad.unit = import.grad.unit ; 
    dataFIC4.sub{i}.hdr.grad.label = import.grad.label(30:170) .' ; 
    dataFIC4.sub{i}.grad.chantype = import.grad.chantype(30:170) .' ; 
    dataFIC4.sub{i}.grad.chanpos = import.grad.chanpos(30:170,:) ;
    dataFIC4.sub{i}.grad.chanori = import.grad.chanori(30:170,:) ;
    dataFIC4.sub{i}.grad.tra = import.grad.tra(30:170,:) ;
    dataFIC4.sub{i}.grad.coilori = import.grad.coilori ; 
    dataFIC4.sub{i}.grad.coilpos = import.grad.coilpos ; 
    dataFIC4.sub{i}.grad.unit = import.grad.unit ; 
    dataFIC4.sub{i}.grad.label = import.grad.label(30:170) .' ; 
   end
end

for i=1:length(trialnum420)
   %load the data for one individual, one trial type, one trial
   k = i + j ; 
   import = load([Dir prefix sub(15,:) trialtype(2) trialnum420(i,:) mt]);
   dataFIC4.trial{k} = import.avg(32:172,:) ; 
   dataFIC4.sub{15}.trial{i} = import.avg(32:172,:) ;
   dataFIC4.time{k} = import.time ;
   dataFIC4.sub{15}.time{i} = import.time ;
end

for i=1:length(trialnum419)
   %load the data for one individual, one trial type, one trial
   k = i + j + l; 
   import = load([Dir prefix sub(14,:) trialtype(2) trialnum419(i,:) mt]);
   dataFIC4.trial{k} = import.avg(32:172,:) ; 
   dataFIC4.sub{14}.trial{i} = import.avg(32:172,:) ;
   dataFIC4.time{k} = import.time ;
   dataFIC4.sub{14}.time{i} = import.time ;
end

for i=1:length(trialnum418)
   %load the data for one individual, one trial type, one trial
   k = i + j + l + m ; 
   import = load([Dir prefix sub(13,:) trialtype(2) trialnum418(i,:) mt]);
   dataFIC4.trial{k} = import.avg(32:172,:) ;
   dataFIC4.sub{13}.trial{i} = import.avg(32:172,:) ;
   dataFIC4.time{k} = import.time ;
   dataFIC4.sub{13}.time{i} = import.time ;
end

for i=1:length(trialnum417)
   %load the data for one individual, one trial type, one trial
   k = i + j + l + m + n ; 
   import = load([Dir prefix sub(12,:) trialtype(2) trialnum417(i,:) mt]);
   dataFIC4.trial{k} = import.avg(32:172,:) ; 
   dataFIC4.sub{12}.trial{i} = import.avg(32:172,:) ;
   dataFIC4.time{k} = import.time ;
   dataFIC4.sub{12}.time{i} = import.time ;
end

for i=1:length(trialnum416)
   %load the data for one individual, one trial type, one trial
   k = i + j + l + m + n + o ; 
   import = load([Dir prefix sub(11,:) trialtype(2) trialnum416(i,:) mt]);
   dataFIC4.trial{k} = import.avg(32:172,:) ; 
   dataFIC4.sub{11}.trial{i} = import.avg(32:172,:) ;
   dataFIC4.time{k} = import.time ;
   dataFIC4.sub{11}.time{i} = import.time ;
end

for i=1:length(trialnum415)
   %load the data for one individual, one trial type, one trial
   k = i + j + l + m + n + o + p ; 
   import = load([Dir prefix sub(10,:) trialtype(2) trialnum415(i,:) mt]);
   dataFIC4.trial{k} = import.avg(32:172,:) ; 
   dataFIC4.sub{10}.trial{i} = import.avg(32:172,:) ;
   dataFIC4.time{k} = import.time ;
   dataFIC4.sub{10}.time{i} = import.time ;
end

for i=1:length(trialnum414)
   %load the data for one individual, one trial type, one trial
   k = i + j + l + m + n + o + p + q; 
   import = load([Dir prefix sub(9,:) trialtype(2) trialnum414(i,:) mt]);
   dataFIC4.trial{k} = import.avg(32:172,:) ; 
   dataFIC4.sub{9}.trial{i} = import.avg(32:172,:) ;
   dataFIC4.time{k} = import.time ;
   dataFIC4.sub{9}.time{i} = import.time ;
end

for i=1:length(trialnum413)
   %load the data for one individual, one trial type, one trial 
   k = i + j + l + m + n + o + p + q + r; 
   import = load([Dir prefix sub(8,:) trialtype(2) trialnum413(i,:) mt]);
   dataFIC4.trial{k} = import.avg(32:172,:) ; 
   dataFIC4.sub{8}.trial{i} = import.avg(32:172,:) ;
   dataFIC4.time{k} = import.time ;
   dataFIC4.sub{8}.time{i} = import.time ;
end

for i=1:length(trialnum411)
   %load the data for one individual, one trial type, one trial
   k = i + j + l + m + n + o + p + q + r + s; 
   import = load([Dir prefix sub(7,:) trialtype(2) trialnum411(i,:) mt]);
   dataFIC4.trial{k} = import.avg(32:172,:) ; 
   dataFIC4.sub{7}.trial{i} = import.avg(32:172,:) ;
   dataFIC4.time{k} = import.time ;
   dataFIC4.sub{7}.time{i} = import.time ;
end

for i=1:length(trialnum409)
   %load the data for one individual, one trial type, one trial
   k = i + j + l + m + n + o + p + q + r + s + t; 
   import = load([Dir prefix sub(6,:) trialtype(2) trialnum409(i,:) mt]);
   dataFIC4.trial{k} = import.avg(32:172,:) ; 
   dataFIC4.sub{6}.trial{i} = import.avg(32:172,:) ;
   dataFIC4.time{k} = import.time ;
   dataFIC4.sub{6}.time{i} = import.time ;
end

for i=1:length(trialnum407)
   %load the data for one individual, one trial type, one trial
   k = i + j + l + m + n + o + p + q + r + s + t + u ; 
   import = load([Dir prefix sub(5,:) trialtype(2) trialnum407(i,:) mt]);
   dataFIC4.trial{k} = import.avg(32:172,:) ; 
   dataFIC4.sub{5}.trial{i} = import.avg(32:172,:) ;
   dataFIC4.time{k} = import.time ;
   dataFIC4.sub{5}.time{i} = import.time ;
end

for i=1:length(trialnum406)
   %load the data for one individual, one trial type, one trial
   k = i + j + l + m + n + o + p + q + r + s + t + u + v ; 
   import = load([Dir prefix sub(4,:) trialtype(2) trialnum406(i,:) mt]);
   dataFIC4.trial{k} = import.avg(32:172,:) ; 
   dataFIC4.sub{4}.trial{i} = import.avg(32:172,:) ;
   dataFIC4.time{k} = import.time ;
   dataFIC4.sub{4}.time{i} = import.time ;
end

for i=1:length(trialnum405)
   %load the data for one individual, one trial type, one trial
   k = i + j + l + m + n + o + p + q + r + s + t + u + v + w ; 
   import = load([Dir prefix sub(3,:) trialtype(2) trialnum405(i,:) mt]);
   dataFIC4.trial{k} = import.avg(32:172,:) ; 
   dataFIC4.sub{3}.trial{i} = import.avg(32:172,:) ;
   dataFIC4.time{k} = import.time ;
   dataFIC4.sub{3}.time{i} = import.time ;
end

for i=1:length(trialnum404)
   %load the data for one individual, one trial type, one trial
   k = i + j + l + m + n + o + p + q + r + s + t + u + v + w + x ; 
   import = load([Dir prefix sub(2,:) trialtype(2) trialnum404(i,:) mt]);
   dataFIC4.trial{k} = import.avg(32:172,:) ; 
   dataFIC4.sub{2}.trial{i} = import.avg(32:172,:) ;
   dataFIC4.time{k} = import.time ;
   dataFIC4.sub{2}.time{i} = import.time ;
end

for i=1:length(trialnum402)
   %load the data for one individual, one trial type, one trial
   k = i + j + l + m + n + o + p + q + r + s + t + u + v + w + x + y ; 
   import = load([Dir prefix sub(1,:) trialtype(2) trialnum402(i,:) mt]);
   dataFIC4.trial{k} = import.avg(32:172,:) ;
   dataFIC4.sub{1}.trial{i} = import.avg(32:172,:) ;
   dataFIC4.time{k} = import.time ;
   dataFIC4.sub{1}.time{i} = import.time ;
end

for i=1:16
   %time windows of interest
    cfg = [];      %empty configuration                                     
    cfg.toilim = [-2.25 -1.25];          %3.125             
    dataPre4 = ft_redefinetrial(cfg, dataFIC4.sub{i});
    cfg.toilim = [0.5 6.25];                       
    dataPost4 = ft_redefinetrial(cfg, dataFIC4.sub{i});
    
    %calculating the cross spectral density matrix
    cfg = [];
    cfg.method    = 'mtmfft';
    cfg.output    = 'powandcsd';
    cfg.tapsmofrq = 4;
    cfg.foilim    = [11.5 12.5]; %change freq here
    freqPre4 = ft_freqanalysis(cfg, dataPre4);
    cfg = [];
    cfg.method    = 'mtmfft';
    cfg.output    = 'powandcsd';
    cfg.tapsmofrq = 4;
    cfg.foilim    = [11.5 12.5]; %change freq here
    freqPost4 = ft_freqanalysis(cfg, dataPost4);
    
    % forward model mri
    mri = ft_read_mri('C:\Users\EEG1\Downloads\mni_colin27_1998_nifti\colin27_t1_tal_lin.nii');
    mri.coordsys = 'ctf' ; 
    cfg = [];
    cfg.write      = 'no';
    [segmentedmri] = ft_volumesegment(cfg, mri);
    %load('C:\Users\EEG1\Documents\Data\segmentedmri.mat');
    cfg = [];
    cfg.method = 'singleshell';
    headmodel = ft_prepare_headmodel(cfg, segmentedmri);
    
    % compute lead field 
    cfg                 = [];
    elec = ft_read_sens('C:\Users\EEG1\Documents\Data\head coordinates\head_coordinates_avg\head_vef_20090408_01-export_cc_HC.sfp');
    cfg.elec = elec; 
    cfg.grad            = freqPost4.grad;
    cfg.headmodel       = headmodel;
    cfg.reducerank      = 2;
    cfg.channel         = 'MEG';
    cfg.grid.resolution = 1;   % use a 3-D grid with a 1 cm resolution
    cfg.grid.unit       = 'cm';
    cfg.senstype = 'MEG';
    %cfg.normalize = 'yes';
    [grid] = ft_prepare_leadfield(cfg);
    
   
    % source analysis without contrast
    cfg              = []; 
    cfg.method       = 'dics';
    cfg.frequency    = 12;  %change freq here
    cfg.grid         = grid; 
    cfg.headmodel    = headmodel;
    cfg.dics.projectnoise = 'yes';
    cfg.dics.lambda       = 0;
    sourcePost_nocon = ft_sourceanalysis(cfg, freqPost4); 
    sourcePre_nocon = ft_sourceanalysis(cfg, freqPre4);
    save sourcePost_nocon sourcePost_nocon
    save sourcePre_nocon sourcePre_nocon
    load sourcePost_nocon
    load sourcePre_nocon
    mri = ft_read_mri('C:\Users\EEG1\Downloads\mni_colin27_1998_nifti\colin27_t1_tal_lin.nii') ;
    mri = ft_volumereslice([], mri);

    cfg            = [];
    cfg.downsample = 2;
    cfg.parameter = 'avg.pow';
    sourcePostInt_nocon  = ft_sourceinterpolate(cfg, sourcePost_nocon , mri);
    sourcePreInt_nocon  = ft_sourceinterpolate(cfg, sourcePre_nocon , mri);
    cfg              = [];
    cfg.method       = 'slice';
    cfg.funparameter = 'avg.pow';
    %%%ft_sourceplot(cfg, sourcePostInt_nocon);
    %%%ft_sourceplot(cfg, sourcePreInt_nocon);
    cfg = [];
    cfg.method        = 'ortho';
    cfg.funparameter  = 'avg.pow';
    cfg.maskparameter = cfg.funparameter;
    %cfg.funcolorlim   = [0.0 1.7];
    %cfg.opacitylim    = [0.0 1.7]; 
    cfg.opacitymap    = 'rampup';  
    %%%ft_sourceplot(cfg, sourcePostInt_nocon);
    %%%ft_sourceplot(cfg, sourcePreInt_nocon);
    
    sourceNAIpo4 = sourcePost_nocon;
    sourceNAIpo4.avg.pow = sourcePost_nocon.avg.pow ./ sourcePost_nocon.avg.noise;
    sourceNAIpe4 = sourcePre_nocon;
    sourceNAIpe4.avg.pow = sourcePre_nocon.avg.pow ./ sourcePre_nocon.avg.noise;
 
    cfg = [];
    cfg.downsample = 2;
    cfg.parameter = 'avg.pow';
    sourceNAIInt4po{i} = ft_sourceinterpolate(cfg, sourceNAIpo4 , mri); % add {i}
    sourceNAIInt4pe{i} = ft_sourceinterpolate(cfg, sourceNAIpe4 , mri); % add {i}
    
    cfg = [];
    cfg.method        = 'slice';
    cfg.funparameter  = 'avg.pow';
    cfg.maskparameter = cfg.funparameter;
    cfg.funcolorlim   = [3.0 12.0];
    cfg.opacitylim    = [3.0 12.0]; 
    cfg.opacitymap    = 'rampup';  
    %%%ft_sourceplot(cfg, sourceNAIInt4po);
    
    % source analysis with contrast
    dataAll4 = ft_appenddata([], dataPre4, dataPost4);
    cfg = [];
    cfg.method    = 'mtmfft';
    cfg.output    = 'powandcsd'; 
    cfg.tapsmofrq = 4;
    cfg.foilim    = [11.5 12.5]; % change freq here
    freqAll3 = ft_freqanalysis(cfg, dataAll4);
    cfg              = [];
    cfg.method       = 'dics';
    cfg.frequency    = 12; %change freq here
    cfg.grid         = grid;
    cfg.headmodel    = headmodel;
    cfg.dics.projectnoise = 'yes';
    cfg.dics.lambda       = '5%';
    cfg.dics.keepfilter   = 'yes';
    cfg.dics.realfilter   = 'yes';
    sourceAll = ft_sourceanalysis(cfg, freqAll3);
    cfg.grid.filter = sourceAll.avg.filter;
    sourcePre_con  = ft_sourceanalysis(cfg, freqPre4 );
    sourcePost_con = ft_sourceanalysis(cfg, freqPost4);
    save sourcePre_con sourcePre_con 
    save sourcePost_con sourcePost_con
    sourceDiff = sourcePost_con;
    sourceDiff.avg.pow = (sourcePost_con.avg.pow - sourcePre_con.avg.pow) ./ sourcePre_con.avg.pow;
    mri = ft_read_mri('C:\Users\EEG1\Downloads\mni_colin27_1998_nifti\colin27_t1_tal_lin.nii') ;
    mri = ft_volumereslice([], mri) ;
    cfg            = [];
    cfg.downsample = 2;
    cfg.parameter  = 'avg.pow';
    sourceDiffInt  = ft_sourceinterpolate(cfg, sourceDiff , mri); 
    %plot the power ratios
    cfg = [];
    cfg.method        = 'slice';
    cfg.funparameter  = 'avg.pow';
    cfg.maskparameter = cfg.funparameter;
    cfg.funcolorlim   = [-1.0 -0.1];
    cfg.opacitylim    = [-1.0 -0.1]; 
    cfg.opacitymap    = 'rampup';  
    %%%ft_sourceplot(cfg, sourceDiffInt);
    cfg = [];
    cfg.method        = 'ortho';
    cfg.funparameter  = 'avg.pow';
    cfg.maskparameter = cfg.funparameter;
    cfg.funcolorlim   = [-1.0 -0.1];
    cfg.opacitylim    = [-1.0 -0.1]; 
    cfg.opacitymap    = 'rampup';  
    %%%ft_sourceplot(cfg, sourceDiffInt);
    
end
    
    
    % run statistics over subjects %
cfg=[];
cfg.dim         = sourceNAIInt4po{1}.dim;
cfg.method      = 'analytic';
cfg.statistic   = 'ft_statfun_depsamplesT';
cfg.parameter   = 'pow';
%cfg.correctm    = 'cluster';
%cfg.numrandomization = 1000;
cfg.alpha       = 0.05; % note that this only implies single-sided testing
cfg.tail        = 0;

%sourceNAIInt4ze = sourceNAIInt4pe ; 
%sourceNAIInt4ze.pow = zeros(2097152,1);
%for i=1:16
%    sourceNAIInt4ze{1,i} = sourceNAIInt4pe ;
%    sourceNAIInt4ze{1,i}.pow = zeros(2097152,1);
%    %sourceNAIInt4ze{1,i}.inside = zeros(128,128,128);
%    sourceNAIInt4pos{1,i} = sourceNAIInt4po;
%    sourceNAIInt4pes{1,i} = sourceNAIInt4pe;
%end
nsubj=length(sub);
cfg.design(1,:) = [1:nsubj 1:nsubj];
cfg.design(2,:) = [ones(1,nsubj) ones(1,nsubj)*2];
cfg.uvar        = 1; % row of design matrix that contains unit variable (in this case: subjects)
cfg.ivar        = 2; % row of design matrix that contains independent variable (the conditions)

stat = ft_sourcestatistics(cfg, sourceNAIInt4pe{:}, sourceNAIInt4po{:});