%% Time Frequency Transformation
% load bip_ref data
% apply bpf ([0.1 150])
% mirror data to avoid edge artifacts in TF transformation
% apply MV-Transformation
% save to disk

function TimeFrequencyTransformation_030422(subject,phase)
    tic
    fs = '\';
    dataDir = '..\Data';
    load('..\subject_list.mat');

    fileName = subNames{subject};
    disp([fileName ' - ' phase])
    cd([dataDir fs fileName])
    cd('bip_ref');

    %load data
    switch phase
        case 'enc'
            load data_enc_photo_bp.mat
        case 'ret'
            load data_ret_photo_bp.mat
    end

    %apply bandpassfilter
    cfg = [];
    cfg.bpfilter = 'yes'; 
    cfg.bpfreq = [0.1 150];
    data_prep = ft_preprocessing(cfg,data);  
    data = data_prep;

    %Reflect the data to avoid edge artifacts and get uniform
    %trial lengths
    Fs = data.fsample;    

    %flip data (mirror on both ends)
    disp('flipping')
    for trial = 1:length(data.trial)
        trialdata = data.trial{1,trial};
        trialtime = data.time{1,trial};

        %create reflected timeseries_data + corresponding time vectors
        %(each series is flipped and reflected on both ends
        reflected_trial = [fliplr(trialdata) trialdata fliplr(trialdata)];  % first segment goes to -1 seconds (idx = 1) % second segment flipped from stim onset (idx = 1001);

        %get new time array
        numTimes = length(trialtime);
        durTrial = (trialtime(end)-trialtime(1));
        refl_start = trialtime(1) - durTrial;
        reflected_time = [linspace(refl_start-1/Fs,trialtime(1)-1/Fs,numTimes) trialtime linspace(trialtime(end)+1/Fs,trialtime(end)+1/Fs+durTrial,numTimes)];

        %store reflected trials in new struct                                                                                                                                  
        flextrial{1,trial} = reflected_trial;
        flextime{1,trial} = reflected_time;
    end

    %%% check whether flipping was done properly
%     figure; subplot(1,2,1);
%     plot(data.time{1,1},data.trial{1,1}(1,:)); xlim([-6 8]);
%     subplot(1,2,2); plot(flextime{1,1},flextrial{1,1}(1,:));
    
    %store new data structure with reflected timeseries
    refl_data = {};
    refl_data = data;
    refl_data.trial = flextrial;
    refl_data.time = flextime;
    
    %set MV Transformation parameters
    disp('transforming')
    freqs = [1:150]; %frequencies to analyze [1:1:29 30:5:150] 
    cycles = [linspace(3,6,30) linspace(6,12,120)]; %linspaced cycles increasing with frequencies 
     % (from 3-6 for 1-30Hz range and from 6-12 for 31-150Hz range)

    cfg = [];
    cfg.method = 'wavelet';
    cfg.pad = 'nextpow2';
    cfg.output = 'pow';
    cfg.keeptrials = 'yes';
    cfg.trials = 'all';
    cfg.toi = -0.5:0.01:1.5;      %upper time-bound set 2s for all trials (some contain mirrored data)
    cfg.foi = freqs;              
    cfg.width = cycles;                                

    %compute morlet wavelet freqanalysis
    tic
    switch phase 
        case 'enc'
            disp(['TF - ENCODE - ' fileName])
            tfData150_enc = ft_freqanalysis(cfg,refl_data);
            data = tfData150_enc; clear tfData150_enc
        case 'ret'
            disp(['TF - RETRIEVE - ' fileName])
            tfData150_ret = ft_freqanalysis(cfg,refl_data);
            data = tfData150_ret; clear tfData150_ret
    end
    disp(['~ ' num2str(toc) ' seconds'])
    clear refl_data flextime flextrial
    
    %save
    cd([dataDir fs fileName])
    status = mkdir('tf_data'); cd('tf_data')
    disp('saving')

    if strcmp(phase,'enc')
        save('tfData150_enc.mat','data','-v7.3');
    else
        save('tfData150_ret.mat','data','-v7.3');
    end
end
