function out = UpsampleTimecourse(timefolder,GLMinput,regions,Data)
% Load the time stamps we've made in Matlab for FSL
% Load the time courses
% Upsample the time courses and cut them according to time stamps for all
% phases of the experiment that we need (e.g. first decision)

%% Presets for upsampling
TR=1.5;   %from scan
upsample=10; %upsampling factor
norm_timec=1; %normalise time course
%nrphases = 4;

%% Load the timing flags
for isub=1:length(Data.SubjectIDs)
    onset.subj{isub}.phases{1}=GLMinput.SelfDec.Onsets{isub};
    onset.subj{isub}.phases{2}=GLMinput.OtherDec.Onsets{isub};
    onset.subj{isub}.phases{3}=GLMinput.Sampling.Onsets{isub};
    onset.subj{isub}.phases{4}=GLMinput.Outcome.Onsets{isub};
end

%% Load, cut (according to timings as for FSL) and upsample the time course file

for iregion = 1:length(regions)
    for isub= 1:length(Data.SubjectIDs)
        
        % Load and upsample timecourse
        %timetxt = ['S' Data.SubjectIDs{isub} '_' regions{iregion} '.txt']; %file name for the timecourse
        
        timetxt = [Data.SubjectIDs{isub} '_' regions{iregion} '.txt']; %file name for the timecourse

        
        loadpathtimc = fullfile(timefolder,timetxt);
        timecourseT=load(loadpathtimc); %load the timecourse (BOLD responses at the measured times, each measurement is TR apart from one before)
        timecourse = timecourseT - mean(timecourseT); % remove the mean
        if norm_timec %normalise the timecourse
            timecourse=zscore(timecourse);
        end
        
        % Upsample the time course
        t=0:length(timecourse)-1;
        t_ups=0:(1/upsample):length(timecourse)-1;
        ts_ups=spline(t,timecourse,t_ups); %upsample the timecourse (i.e. the bold values)
        ts_upsDerivative=diff([0 ts_ups]);
        % Cut the time course into the phases I need
        %%%%%%%
        % For DECISION PHASE
        window = 12; % epoch window duration in seconds
        TCShift=0; % in case you don't want to start event log at zero but later.
        pseudosamples=round((window./TR)*upsample); % number of pseudosamples in epoch
        
        % Get the onsets
        onsets=round(((onset.subj{isub}.phases{1}(:,1))./TR-1)*upsample)+1;% for ech event in phase, get the time value in the new framework
        %onsets are in the new timeframe, the beginning of ts_ups is at time =0 (as the timing file i have read in was already like that)
        % the '..TR-xx' means that the epoch starts xx seconds before timing flag
        tmpmat=repmat(onsets,1,pseudosamples)+repmat(0:pseudosamples-1,size(onsets,1),1);
        
        ps= ts_ups(tmpmat); %ts_ups is the new timecourse, ps is same time course, but cut into epochs
        %psDerivative=ts_upsDerivative(tmpmat);
        
        % Store all the upsampled and epoched time courses together
        
        %%%% Self Dec
        out.SelfDec.Region(iregion).v{isub} = ps;
        out.SelfDec.Region(iregion).name    = regions{iregion};
        out.SelfDec.info                    = 'Time course time locked to the SelfDec phase';
        out.SelfDec.pseudosamples           = pseudosamples;
        out.SelfDec.window                  = window;
        out.SelfDec.TCShift                 = TCShift;
        out.SelfDec.TR                      = TR;
        out.SelfDec.upsample                = upsample;
               
       
        %%%%% Other Dec
        window=6;TCShift=3;
        pseudosamples=round((window./TR)*upsample);
        onsets=round(((onset.subj{isub}.phases{2}(:,1)+TCShift)./TR-1)*upsample)+1;
        tmpmat=repmat(onsets,1,pseudosamples)+repmat(0:pseudosamples-1,size(onsets,1),1);
        ps= ts_ups(tmpmat);
        out.OtherDec.Region(iregion).v{isub} = ps;
        out.OtherDec.Region(iregion).name    = regions{iregion};
        out.OtherDec.info                    = 'Time course time locked to the onset of the OtherDec phase';
        out.OtherDec.pseudosamples           = pseudosamples;
        out.OtherDec.window                  = window;
        out.OtherDec.TCShift                 = TCShift;
        out.OtherDec.TR                      = TR;
        out.OtherDec.upsample                = upsample;

        
         %%%%% Sampling         
        window=12;TCShift=0;
        pseudosamples=round((window./TR)*upsample);
        onsets=round(((onset.subj{isub}.phases{3}(:,1))./TR-1)*upsample)+1;
        tmpmat=repmat(onsets,1,pseudosamples)+repmat(0:pseudosamples-1,size(onsets,1),1);
        ps= ts_ups(tmpmat);
        out.Sampling.Region(iregion).v{isub} = ps;
        out.Sampling.Region(iregion).name    = regions{iregion};
        out.Sampling.info                    = 'Time course time locked to the onset of the Sampling phase';
        out.Sampling.pseudosamples           = pseudosamples;
        out.Sampling.window                  = window;
        out.Sampling.TCShift                 = TCShift;        
        out.Sampling.TR                      = TR;
        out.Sampling.upsample                = upsample;

          %%%%% Outcome         
        window=8;TCShift=0;
        pseudosamples=round((window./TR)*upsample);
        onsets=round(((onset.subj{isub}.phases{4}(:,1))./TR-1)*upsample)+1;
        tmpmat=repmat(onsets,1,pseudosamples)+repmat(0:pseudosamples-1,size(onsets,1),1);
        ps= ts_ups(tmpmat);
        out.Outcome.Region(iregion).v{isub} = ps;
        out.Outcome.Region(iregion).name    = regions{iregion};
        out.Outcome.info                    = 'Time course time locked to the onset of the Outcome phase';
        out.Outcome.pseudosamples           = pseudosamples;
        out.Outcome.window                  = window;
        out.Outcome.TCShift                 = TCShift;
        out.Outcome.TR                      = TR;
        out.Outcome.upsample                = upsample;

    end
end

