datadir = '../data';  % CHANGE THIS FOR THE CORRECT LOCATION OF THE DATA
subj    = 15;         % CHANGE THIS NUMBER FOR EACH SUBJECT

%% read the data from all separate runs

% this will contain the runs for a single subject
rundata = {};

for run=1:6
  trialdef = fullfile(datadir, sprintf('Sub%02d', subj), 'MEEG', 'Trials', sprintf('run_%02d_trldef.txt', run));
  dataset  = fullfile(datadir, sprintf('Sub%02d', subj), 'MEEG', sprintf('run_%02d_sss.fif', run));
  
  [begsample, endsample, offset, trialtype] = textread(trialdef, '%d%d%d%s');
  
  trialcode = nan(size(trialtype));
  trialcode(strcmp(trialtype, 'Famous'))      = 1;
  trialcode(strcmp(trialtype, 'Unfamiliar'))  = 2;
  trialcode(strcmp(trialtype, 'Scrambled'))   = 3;
  
  % construct the trial definition matrix, usually done with FT_DEFINETRIAL
  trl = [begsample(:) endsample(:) offset(:) trialcode(:)];
  
  cfg         = [];
  cfg.dataset = dataset;
  cfg.trl     = trl;
  
  % MEG specific settings
  cfg.channel = 'MEG';
  cfg.demean  = 'yes';
  data_meg = ft_preprocessing(cfg);
  
  % EEG specific settings
  cfg.channel    = 'EEG';
  cfg.demean     = 'yes';
  cfg.reref      = 'yes';
  cfg.refchannel = 'all'; % average reference
  data_eeg = ft_preprocessing(cfg);
  
  % settings for all other channels
  cfg.channel = {'all', '-MEG', '-EEG'};
  cfg.demean  = 'no';
  cfg.reref   = 'no';
  data_other = ft_preprocessing(cfg);
  
  cfg = [];
  cfg.resamplefs = 300;
  data_meg   = ft_resampledata(cfg, data_meg);
  data_eeg   = ft_resampledata(cfg, data_eeg);
  data_other = ft_resampledata(cfg, data_other);
  
  %% append the different channel sets into a single structure
  
  rundata{run} = ft_appenddata(cfg, data_meg, data_eeg, data_other);
  clear data_meg data_eeg data_other
  
end % for each run


%% append the 6 runs into a single structure
data = ft_appenddata(cfg, rundata{:});


%% compute the overall average and condition-specific averages

cfg = [];
cfg.trials = find(data.trialinfo==1);
avg_Famous = ft_timelockanalysis(cfg, data);
cfg.trials = find(data.trialinfo==2);
avg_Unfamiliar = ft_timelockanalysis(cfg, data);

cfg.trials = find(data.trialinfo==3);
avg_Scrambled = ft_timelockanalysis(cfg, data);

cfg.trials = find(data.trialinfo==1 | data.trialinfo==2);
avg_Faces = ft_timelockanalysis(cfg, data);

cfg = [];
% cfg.layout = 'neuromag306all.lay';
cfg.layout = 'neuromag306mag.lay';
ft_multiplotER(cfg, avg_Faces, avg_Scrambled);


%% compute the difference between faces and

cfg = [];
cfg.parameter = 'avg';
cfg.operation = 'x1-x2';
avg_Faces_vs_Scrambled = ft_math(cfg, avg_Faces, avg_Scrambled);

cfg        = [];
% cfg.layout = 'neuromag306all.lay';
cfg.layout = 'neuromag306mag.lay';
ft_multiplotER(cfg, avg_Faces_vs_Scrambled);

cfg           = [];
cfg.layout    = 'neuromag306mag.lay';
cfg.colorbar  = 'yes';
figure; ft_movieplotER(cfg, avg_Faces_vs_Scrambled);

% for saving to disk
prefix = sprintf('Sub%02d', subj);

%% save the raw data to disk
% save([prefix '_raw'], 'data');

%% save the averages to disk
prefix = sprintf('Sub%02d', subj);
% save([prefix '_avg_Famous'],     'avg_Famous');
% save([prefix '_avg_Unfamiliar'], 'avg_Unfamiliar');
% save([prefix '_avg_Scrambled'],  'avg_Scrambled');
% save([prefix '_avg_Faces'],      'avg_Faces');
% save([prefix '_avg_Faces_vs_Scrambled'], 'avg_Faces_vs_Scrambled');

%% look at the analysis history
cfg           = [];
cfg.filename  = [prefix '_avg_Faces_vs_Scrambled.html'];
ft_analysispipeline(cfg, avg_Faces_vs_Scrambled);


