load template_hdm.mat
load template_grad.mat
load hdm.mat

% 1. create grid for voxel of interest
	cfg                 = [];
	cfg.hdmfile         = hdm;
	cfg.inwardshift     = 0.55;
	cfg.grid.xgrid      = -0.6;
	cfg.grid.ygrid      = -9.6;
	cfg.grid.zgrid      = 0;

	virt_grid = prepare_dipole_grid(cfg, template_hdm, grad); 
	virt_grid.dim       = [1 1 1];

% 2. load data

load dataALL

% 3. preprocess data set

    cfg                     = [];
    cfg.channel             = {'MEG'};
    cfg.dftfilter           = 'yes';
    cfg.dftfreq             = [50 100 150];
    cfg.blc                 = 'no';                             
    preprocessed            = preprocessing(cfg, dataALL);

% 4. redefine trials to pre and post stimulus time ranges

    cfg                     = [];
    cfg.toilim              = [-2 0];
    cfg.minlength           = 0.3;
    
    dataPre = redefinetrial(cfg, preprocessed);
    
    cfg                     = [];
    cfg.toilim              = [0 2];
    cfg.minlength           = 0.3;
    
    dataPost = redefinetrial(cfg, preprocessed);
  
% 5. do timelock analysis on certain time windows
   
	cfg                     = [];
	cfg.trials              = trialvis_pre;
	cfg.keeptrials          = 'yes';
	cfg.blcovariance        = 'no';
	cfg.vartrllength        = 2;
	cfg.covariance          = 'yes';
	cfg.covariancewindow    = [-2 0];

	tlckpre = timelockanalysis(cfg, dataPre);

	cfg.trials              = trialvis_post;
	cfg.covariancewindow    = [0 2.0];

	tlckpost = timelockanalysis(cfg, dataPost);

% 6. do source analysis on results of timelock analysis

	cfg                     = [];
	cfg.method              = 'lcmv';
	cfg.grid                = virt_grid;
	cfg.vol                 = template_hdm;
	cfg.grad                = grad;
	cfg.lambda              = '5%';
	cfg.projectnoise        = 'yes';
	cfg.keeptrials          = 'yes';
	cfg.rawtrial            = 'yes';

	sourcepre  = sourceanalysis(cfg, tlckpre);
	sourcepost  = sourceanalysis(cfg, tlckpost);

% 7. make a structure resembling a preprocessed file and remove NaNs in trials 

	datapost                    = [];
	datapost.fsample            = 999.999;
	datapost.cfg                = sourcepost.cfg;
	datapost.label              = {'dip1', 'dip2', 'dip3'};
		
	for i = 1:numel(sourcepost.trial)
    		datapost.time{i}        = sourcepost.time;
    		datapost.trial{i}       = sourcepost.trial(i).mom{1,1}; % something with source.trial(i).mom{1} and mom{2};
	end

	for h = 1:numel(datapost.trial)
		[ind1 ind2]                 = find(isnan(datapost.trial{h}));
		datapost.trial{h}(:,ind2)   = [];
		datapost.time{h}(:,ind2)    = [];
	end

% same done for pre data

% 8. calculate TFR over results of source analysis

	cfg = [];
	cfg.trials              = 'all';
	cfg.output              = 'pow';
	cfg.keeptrials          = 'no';

	cfg.method              = 'mtmconvol';
	cfg.toi                 = -2:.05:0; 
	cfg.foi                 = 1:2:150; 
	cfg.t_ftimwin           = 0.2 * ones(numel(cfg.foi));
	cfg.tapsmofrq           = 10 * ones(numel(cfg.foi));

	TFRpre                  = freqanalysis(cfg,datapre); 

	cfg.toi                 = 0:.05:2; 
	TFRpost                 = freqanalysis(cfg,datapost); 

% 9. calculate difference between pre and post stimulation 
	
	TFRDiff                 = TFRpost;

	TFRDiff.powspctrm       = (TFRpost.powspctrm - TFRpre.powspctrm)...
                             ./ TFRpre.powspctrm;
                   
	TFRDiff2                = TFRpost;
	TFRDiff2.powspctrm      = TFRpost.powspctrm ./ TFRpre.powspctrm;                      

% 10. plot results

	cfg = [];
	cfg.channel = 'dip1';
	cfg.baseline = 'no';
	cfg.interactive = 'yes';
	figure;
	singleplotTFR(cfg,TFRDiff)
	figure;
	singleplotTFR(cfg,TFRDiff2)