% this script does channel-level time-frequency analysis

definetrial       = true;
readdata          = true;
dofreq_notime     = true;
dofreq_mtmconvol  = true;
dofreq_wavelet    = true;
doplot            = true;
dofreq_singlechan = true;
writeflag         = true;

if ~exist('subj', 'var')
  error('specify a subject by calling subj = datainfo_subject(<number>)');
end

filename = fullfile(subj.outputpath, 'sensoranalysis', subj.name, sprintf('%s_data', subj.name));
load(filename, 'data')
%% create the trial definition for each run

if definetrial

  subj.trl = cell(6,1);
  for run_nr = 1:6

    cfg = [];
    cfg.dataset = subj.megfile{run_nr};
    cfg.readbids = 'no'; % it looks like BIDS, but is not complete

    % this is what we can see in the BIDS events.tsv file
    Famous      = [5 6 7];
    Unfamiliar  = [13 14 15];
    Scrambled   = [17 18 19];

    cfg.trialfun = 'ft_trialfun_general';
    cfg.trialdef.eventtype = 'STI101';
    cfg.trialdef.eventvalue = [Famous Unfamiliar Scrambled];
    cfg.trialdef.prestim = 0.5;
    cfg.trialdef.poststim = 1.2;

    cfg = ft_definetrial(cfg);

    % remember the trial definition for each run
    subj.trl{run_nr} = cfg.trl;

  end % for each run

end


%% read the data

if readdata

  rundata = cell(1,6);
  for run_nr = 1:6

    cfg         = [];
    cfg.dataset = subj.megfile{run_nr};
    cfg.trl     = subj.trl{run_nr};

    % MEG specific settings, minimal preprocessing %% NO PADDING, NO FILTER, EPOCH LENGTH IS NOT INCREASED!
    cfg.channel = 'MEG';
    cfg.demean  = 'yes';
    cfg.coilaccuracy = 0;
    data_meg    = ft_preprocessing(cfg);

    cfg            = [];
    cfg.resamplefs = 300;
    data_meg       = ft_resampledata(cfg, data_meg);

    % remember the data from this run
    rundata{run_nr} = data_meg;
  end % for each run

  data = ft_appenddata([], rundata{:});
  clear rundata;

  if writeflag
    filename = fullfile(subj.outputpath, 'sensoranalysis', subj.name, sprintf('%s_data', subj.name));
    save(filename, 'data');
  end

end

%% Frequency analysis not time-resolved

if dofreq_notime

  filename = fullfile(subj.outputpath, 'sensoranalysis', subj.name, sprintf('%s_data', subj.name));
  load(filename, 'data');

  cfg = [];
  cfg.method = 'mtmfft';
  cfg.taper = 'hanning';
  freq = ft_freqanalysis(cfg, data);

  sel = find(strcmp(freq.label, 'MEG0741')); % 84

  figure
  plot(freq.freq, freq.powspctrm(84,:));
  xlabel('frequency')
  ylabel('power')

  figure
  semilogy(freq.freq, freq.powspctrm(84,:))
  xlabel('frequency')
  ylabel('power')

  figure
  plot(freq.freq, 10*log10(freq.powspctrm(84,:)));
  xlabel('frequency')
  ylabel('power (dB)')
  grid on
  set(gca, 'XTick', 0:10:150)

  cfg = [];
  cfg.operation = '10*log10(x1)';
  cfg.parameter = 'powspctrm'; % this is represented by x, and x1 refers to the first (and only) input data structure
  freq_db = ft_math(cfg, freq);

  cfg = [];
  cfg.method = 'mtmfft';
  cfg.taper = 'hanning';
  cfg.keeptrials = 'yes';
  cfg.channel = 'MEG0741';
  freq_MEG0741 = ft_freqanalysis(cfg, data);

  imagesc(10*log10(squeeze(freq_MEG0741.powspctrm)))

end

%% Time-frequency analysis I - Hanning taper, fixed window length

if dofreq_mtmconvol

  filename = fullfile(subj.outputpath, 'sensoranalysis', subj.name, sprintf('%s_data', subj.name));
  load(filename, 'data');

  cfg        = [];
  cfg.method = 'mtmconvol';
  cfg.output = 'pow';
  cfg.foi    = 2.5:2.5:30;
  cfg.t_ftimwin = ones(1,numel(cfg.foi)).*0.4;
  cfg.taper  = 'hanning';
  cfg.toi    = (-0.8:0.05:1.3);
  cfg.pad    = 4;

  cfg.trials = ismember(data.trialinfo(:,1), Famous);
  freqlow_famous = ft_freqanalysis(cfg, data);

  cfg.trials = ismember(data.trialinfo(:,1), Unfamiliar);
  freqlow_unfamiliar = ft_freqanalysis(cfg, data);

  cfg.trials = ismember(data.trialinfo(:,1), Scrambled);
  freqlow_scrambled = ft_freqanalysis(cfg, data);

end

if doplot % Different figures from tutorial. Subject 15 used in tutorial?

    cfg = [];
    cfg.baseline     = [-0.6 -0.2];
    cfg.baselinetype = 'absolute';
    cfg.zlim         = [-5e-27 5e-27];
    cfg.showlabels   = 'yes';
    cfg.layout       = 'neuromag306mag_helmet.mat';
    figure; ft_multiplotTFR(cfg, freqlow_famous);

    cfg = [];
    cfg.baseline     = [-0.6 -0.2];
    cfg.baselinetype = 'absolute';
    cfg.maskstyle    = 'saturation';
    cfg.zlim         = [-5e-27 5e-27];
    cfg.channel      = 'MEG0731';
    cfg.layout       = 'neuromag306mag_helmet.mat';
    figure; ft_singleplotTFR(cfg, freqlow_famous);

    cfg = [];
    cfg.baseline     = [-0.6 -0.2];
    cfg.baselinetype = 'absolute';
    cfg.xlim         = [0.6 0.8];
    cfg.zlim         = [-5e-27 5e-27];
    cfg.ylim         = [4 6];
    cfg.marker       = 'on';
    cfg.layout       = 'neuromag306mag_helmet.mat';
    figure; ft_topoplotTFR(cfg, freqlow_famous);

end

% Ex1 - Relative or absolute baseline

% Ex2 - Plot the TFR of sensor MEG1921. How do you account for the increased power at ~100-200 ms? Hint: compare it to the ERFs.

% Ex3 - Visualize the TFR of the gradiometers. Use what you have learnt from the raw2erp tutorial to first combine the horizontal and planar gradient channels into a single estimate.


%% Time-frequency analysis II - Hanning taper, frequency-dependent window length

if dofreq_singlechan

  cfg              = [];
  cfg.output       = 'pow';
  cfg.channel      = 'MEG0741';
  cfg.method       = 'mtmconvol';
  cfg.taper        = 'hanning';
  cfg.foi          = 2:1:30;
  cfg.toi          = -0.8:0.05:1.5;
  cfg.trials = ismember(data.trialinfo(:,1), Famous); % DIFFERENT TRIAL SELECTION FROM TUTORIAL
  
  cfg.t_ftimwin    = 4./cfg.foi;   % 4 cycles per time window
  TFRhann4         = ft_freqanalysis(cfg, data);
  
  cfg.t_ftimwin    = 5./cfg.foi; % 5 cycles per time window
  TFRhann5         = ft_freqanalysis(cfg, data);
  
  cfg.t_ftimwin    = 7./cfg.foi;  % 7 cycles per time window
  TFRhann7         = ft_freqanalysis(cfg, data);

  cfg.t_ftimwin    = 10./cfg.foi; % 10 cycles per time window
  TFRhann10         = ft_freqanalysis(cfg, data);

end

if doplot

  cfg              = [];
  cfg.baseline     = [-0.5 -0.1];
  cfg.baselinetype = 'relchange';
  cfg.maskstyle    = 'saturation';
  cfg.zlim         = [-1 1];
  cfg.channel      = 'MEG0741';
  cfg.interactive  = 'no';
  figure; ft_singleplotTFR(cfg, TFRhann4);
  figure; ft_singleplotTFR(cfg, TFRhann5);
  figure; ft_singleplotTFR(cfg, TFRhann7);
  figure; ft_singleplotTFR(cfg, TFRhann10);

end

%% Time-frequency analysis III - Multitapers
cfg        = [];
cfg.method = 'mtmconvol';
cfg.output = 'pow';
cfg.foi    = 30:5:80;
cfg.t_ftimwin = ones(1,numel(cfg.foi)).*0.2;
cfg.tapsmofrq = ones(1,numel(cfg.foi)).*10;
cfg.taper     = 'dpss';
cfg.toi       = (-0.8:0.05:1.3);
cfg.pad       = 4;

cfg.trials = ismember(data.trialinfo(:,1), Famous);
freqhigh_famous = ft_freqanalysis(cfg, data);

cfg.trials = ismember(data.trialinfo(:,1), Unfamiliar);
freqhigh_unfamiliar = ft_freqanalysis(cfg, data);

cfg.trials = ismember(data.trialinfo(:,1), Scrambled);
freqhigh_scrambled = ft_freqanalysis(cfg, data);

if doplot
    cfg = [];
    cfg.layout       = 'neuromag306mag_helmet.mat';
    cfg.baseline     = [-0.6 -0.2];
    cfg.baselinetype = 'relchange';
    cfg.zlim         = [-.2 .2];
    cfg.marker       = 'on';
    figure; ft_multiplotTFR(cfg, freqhigh_famous);
end

if writeflag
    filename = fullfile(subj.outputpath, 'sensoranalysis', subj.name, sprintf('%s_freq_mtmconvol', subj.name));
    save(filename, 'freqlow_famous', 'freqlow_unfamiliar', 'freqlow_scrambled', 'freqhigh_famous', 'freqhigh_unfamiliar', 'freqhigh_scrambled');
end

% Ex5 - Visualize the difference between 2 conditions.
cfg = [];
cfg.parameter = 'powspctrm';
cfg.operation = 'log10(x1)-log10(x2)';
freqhigh_contrast = ft_math(cfg, freqhigh_famous, freqhigh_scrambled);

%% Time-frequency analysis IV - Morlet wavelets

if dofreq_wavelet

  filename = fullfile(subj.outputpath, 'sensoranalysis', subj.name, sprintf('%s_data', subj.name));
  load(filename, 'data');

  cfg        = [];
  cfg.method = 'wavelet';
  cfg.output = 'pow';
  cfg.foi    = 1:1:60;
  cfg.width  = 7;
  cfg.toi    = (-0.8:0.05:1.3);
  cfg.pad    = 4;

  cfg.trials = ismember(data.trialinfo(:,1), Famous);
  freq_famous = ft_freqanalysis(cfg, data);

  cfg.trials = ismember(data.trialinfo(:,1), Unfamiliar);
  freq_unfamiliar = ft_freqanalysis(cfg, data);

  cfg.trials = ismember(data.trialinfo(:,1), Scrambled);
  freq_scrambled = ft_freqanalysis(cfg, data);

  if writeflag
    filename = fullfile(subj.outputpath, 'sensoranalysis', subj.name, sprintf('%s_freq_wavelet', subj.name));
    save(filename, 'freq_famous', 'freq_unfamiliar', 'freq_scrambled');
  end

end

if doplot

    cfg = [];
    cfg.baseline     = [-0.6 -0.1];
    cfg.baselinetype = 'absolute';
    cfg.zlim         = [-1e-25 1e-25];
    cfg.showlabels   = 'yes';
    cfg.layout       = 'neuromag306mag_helmet.mat';
    figure; ft_multiplotTFR(cfg, freq_famous)

end

% Ex6 - Adjust cfg.width and see how the TFRs change.
