JA_TMS_2022/Pilot_analysis/JA_pilot_analysis_DEF_rep.m

%%
clear
general_directory = '/Users/b1037210/Desktop/Progetti/JA_Stage1/'
cd(general_directory)
addpath(general_directory)
%%
addpath(genpath('/Users/b1037210/Documents/MATLAB/fieldtrip-20161117'))

ft_defaults
%%

subjects = ...
    {'S01_pilot',...
    'S02_pilot',...
    'S03_pilot',...
    'S04_pilot',...
    'S05_pilot',...
    'S06_pilot',...
    'S07_pilot',...
    'S08_pilot',...
    'S09_pilot',...
    'S10_pilot',...
    'S11_pilot',...
    'S12_pilot'}


threshold = {[0.018 0.032],...
    [0.016 0.035],...
    [0.018 0.040],...
    [0.019 0.039],...
    [0.017 0.037],...
    [0.016 0.036],...
    [0.019 0.042],...
    [0.018 0.039],...
    [0.016 0.035],...
    [0.015 0.035],...
    [0.02 0.03],...
    [0.02 0.04]}



TMS_timings = {'-0.05', '0.0'};
final_table = table;
%%
for s = 1:length(subjects)
    
    %% cd directory
    cd([general_directory,subjects{s}])
    
    %% ...for each tms timing 
    for t = 1:length(TMS_timings)
        
        %% file loading
        file_tmp = (dir(strcat(subjects{s},'_', TMS_timings{t},'_TMS*')))
        
        
        if size(file_tmp,1)>1
        
        disp(strcat('concatenate ',file_tmp(1).name))
        [y SR MetaData] = GB_gtec_files_concatenation(file_tmp);
        
        else
        
        disp(strcat('loading ', file_tmp.name)) 
        load(file_tmp.name)

        end
        
        clear file_tmp
        %% name the trigger channel and time channel, then erase time from y
        
        % select trigger row and time row
        trigger = y(end,:); % trigger channel for glove touch detection
        time = y(1,:); % time-series in seconds
        
        % new data without time, so channels in the gtect settings and data
        % analysis have the same nnumber
        y = y(2:end,:);
        
        %% get triggers and sample points at which each event of interest happens
        
        % get were variations happen within trigger channel
        diff_trigger = diff(trigger);
        
        % get trigger codes wihtin trgger line
        trigger_numbers = unique(trigger);
        
        % find were the wash_out ends
        end_wash_out = find(diff_trigger == 1)+1;
        end_wash_out = end_wash_out(1);
        
        % for each code present, get the sample were the rising triggers are present
        for x = 1:length(trigger_numbers)
            
            % find rising triggers and add 1 sampling point (due to diff which removes one)
            tmp = find(diff_trigger == trigger_numbers(x))+1;
            
            % based on the code --> assign the tmp numbers to the correct variable
            switch trigger_numbers(x)
                
%                 case 14 % color press CF
%                     sp_slow_id = tmp(tmp>end_wash_out);
%                     
%                 case 15 % color press CF
%                     sp_fast_id = tmp(tmp>end_wash_out);
                
                case 16 % color press CF
                    sp_out_slow_id = tmp(tmp>end_wash_out);
                
                case 17 % color press CF
                    sp_out_fast_id = tmp(tmp>end_wash_out);
                    
%                 case 18 % color press CF
%                     sp_out_ok_catch_id = tmp(tmp>end_wash_out);
%                      
%                 case 19 % color press CF
%                     sp_out_failed_catch_id = tmp(tmp>end_wash_out);
%                 
                    
                    
                    
                    
                case 20 % color press CF
                    sp_color_CF_press = tmp(tmp>end_wash_out);
                    
                case 21 % color lift CF
                    sp_color_CF_lift = tmp(tmp>end_wash_out);
                    
                case 22 % color catch CF
                    sp_color_CF_catch = tmp(tmp>end_wash_out);
                    
                case 30 % color press SBJ
                    sp_color_SBJ_press = tmp(tmp>end_wash_out);
                    
                case 31 % color_lift SBJ
                    sp_color_SBJ_lift = tmp(tmp>end_wash_out);
                    
                case 40 % CF hitting
                    sp_CF_hitting = tmp(tmp>end_wash_out);
                    
                case 41 % SBJ hitting
                    sp_SBJ_hitting = tmp(tmp>end_wash_out);
                    
                case 128 % TMS pulse
                    sp_TMS = tmp(tmp>end_wash_out);
            end
            
        end
        
        %%
        
        % check dimensions
        size(sp_color_CF_press)
        size(sp_color_CF_lift)
        size(sp_color_CF_catch)
        size(sp_TMS)
        
        %
        [CF_colors_idx] = sortrows(vertcat(sp_color_CF_press',sp_color_CF_lift',sp_color_CF_catch'));
        CF_colors_order = CF_colors_idx;
        
        % assign codes to sample points corresponding to press lift and catch
        CF_colors_order(ismember(CF_colors_idx, sp_color_CF_press)) = 1;
        CF_colors_order(ismember(CF_colors_idx, sp_color_CF_lift)) = -1;
        CF_colors_order(ismember(CF_colors_idx, sp_color_CF_catch)) = 0;
        
        %%
        

        size(sp_out_slow_id)
        size(sp_out_fast_id)

         
  [trials_idx] = sortrows(vertcat(sp_out_slow_id',sp_out_fast_id'))

 
         
       %%
       for ex = 1:length(trials_idx)
          
           [v_abs(ex) p_abs(ex)] = min(abs(CF_colors_idx - trials_idx(ex)))
          ddd = CF_colors_idx - trials_idx(ex);
           exclude(ex) = p_abs(ex);
           is_exclude_positive(ex) = ddd(exclude(ex))
       end
       
       clear v_abs
       clear ddd
       clear is_exclude_positive
       
       
        
        %% EMG preprocessing
        
        ECD_right = y(2,:);
        FDS_right = y(4,:);
        ECD_response = y(6,:);
        
        
        
        Fsample = SR;
        hp_filter = 10
        lp_filter = 2500
        
   
        ECD_right_filt = ft_preproc_highpassfilter(ECD_right, Fsample, hp_filter,1);
        FDS_right_filt = ft_preproc_highpassfilter(FDS_right, Fsample, hp_filter,1);
        
        ECD_right_filt = ft_preproc_lowpassfilter(ECD_right_filt, Fsample, lp_filter,1);
        FDS_right_filt = ft_preproc_lowpassfilter(FDS_right_filt, Fsample, lp_filter,1);
        
        
        ECD_right_filt = ft_preproc_dftfilter(ECD_right_filt, Fsample);
        FDS_right_filt = ft_preproc_dftfilter(FDS_right_filt, Fsample);
        
        
 %% take away the decay from MEPs
        
cfg = [];
cfg.start_time = 0.005;
cfg.end_time = 0.15;
cfg.SR = SR;
cfg.sp_trigger = sp_TMS;

[ECD_sweep_cuts_tmp, ECD_sweep_time, ECD_sample_points] = gb_cutting_sweeps(cfg, ECD_right_filt);
[FDS_sweep_cuts_tmp, FDS_sweep_time, FDS_sample_points] = gb_cutting_sweeps(cfg, FDS_right_filt);




        %% get MEPs p2p amplitudes
        cfg = [];
        cfg.start_time = threshold{s}(1);
        cfg.end_time = threshold{s}(2);
        cfg.SR = SR;
        cfg.sp_trigger = sp_TMS;
        
        [ECD_MEP_amplitudes] = gb_MEP_p2p(cfg, ECD_right_filt);
        [FDS_MEP_amplitudes] = gb_MEP_p2p(cfg, FDS_right_filt);
        
        
%% fit the decay and remove it

cfg = [];
cfg.start_time = -0.15;
cfg.end_time = -0.002;
cfg.SR = SR;
cfg.sp_trigger = sp_TMS;

[ECD_sweep_cuts_bkgrnd_tmp, ECD_sweep_time_bkgrnd, ECD_sample_points_bkgrnd] = gb_cutting_sweeps(cfg, ECD_right_filt)

[FDS_sweep_cuts_bkgrnd_tmp, FDS_sweep_time_bkgrnd, FDS_sample_points_bkgrnd] = gb_cutting_sweeps(cfg, FDS_right_filt)

%%
for x = 1:length(ECD_sweep_cuts_bkgrnd_tmp)
[~, ~, resECD{x}] = fit(ECD_sweep_time_bkgrnd',ECD_sweep_cuts_bkgrnd_tmp{x}','poly3'); % time in ms
ECD_sweep_cuts_bkgrnd{x} =resECD{x}.residuals';

[~, ~, resFDS{x}] = fit(FDS_sweep_time_bkgrnd',FDS_sweep_cuts_bkgrnd_tmp{x}','poly3'); % time in ms
FDS_sweep_cuts_bkgrnd{x} =resFDS{x}.residuals';

end
clear resECD resFDS

%%
for x = 1:length(ECD_sample_points_bkgrnd)
    
ECD_right_filt(ECD_sample_points_bkgrnd{x}) = ECD_sweep_cuts_bkgrnd{x};
FDS_right_filt(FDS_sample_points_bkgrnd{x}) = FDS_sweep_cuts_bkgrnd{x};

end


        %% any contraction befor TMS?
        
        cfg = [];
        cfg.start_time = -0.05;
        cfg.end_time = -0.002;
        cfg.SR = SR;
        cfg.sp_trigger = sp_TMS;
        cfg.analysis = 'p2p';
        
        [ECD_MEP_bkgrnd] = gb_background_activity(cfg, ECD_right_filt)
        [FDS_MEP_bkgrnd] = gb_background_activity(cfg, FDS_right_filt)
        
        
        %% muscle contraction before TMS?
        
        MEP_excl_ECD_idx = find(ECD_MEP_bkgrnd>200) % microvolt
        MEP_excl_FDS_idx = find(FDS_MEP_bkgrnd>200) % microvolt
        
        
        %% EXCLUDE MEP<50 p2p microVolts
        
        MEP_excl2_ECD_idx = find(ECD_MEP_amplitudes<50) % microvolt
       
        
         %% EXCLUDE out slow and out fast trials
        
      ECD_MEP_amplitudes(exclude) = nan;
     
        
        %% exclude trials with background contraction
        
        ECD_MEP_amplitudes(union(MEP_excl_ECD_idx, MEP_excl_FDS_idx)) = nan;
    
        
        ECD_MEP_amplitudes(MEP_excl2_ECD_idx) = nan;
    
 
        %% median amplitude based on color trigger         
        summary.ECD_press{s,t} = nanmedian((ECD_MEP_amplitudes(CF_colors_order == 1)));
        summary.ECD_lift{s,t} = nanmedian((ECD_MEP_amplitudes(CF_colors_order == -1)));
        %%
              
cfg = [];
cfg.start_time = -0.015;
cfg.end_time = 0.15;
cfg.SR = SR;
cfg.sp_trigger = sp_TMS;

[bbb.ECD_sweep_PLOT_CHECKS, bbb.ECD_sweep_time_PLOT_CHECKS, bbb.ECD_sample_points_PLOT_CHECKS] = gb_cutting_sweeps(cfg, ECD_right_filt);
    
        %%
        
        trial = (1:length(CF_colors_order))';
        subject_rep = repmat(subjects(s), [length(CF_colors_order),1]);
        TMS_timings_rep = repmat(TMS_timings(t), [length(CF_colors_order),1]);
        
        subject_table{t} = table(subject_rep,trial, ECD_MEP_amplitudes,...
            ECD_MEP_bkgrnd,...
            FDS_MEP_bkgrnd,...
            TMS_timings_rep,...
            CF_colors_order)
        
        
        
        
        
        
        clearvars -except s t subjects TMS_timings threshold subject_table general_directory final_table summary
        
        %%
       
    end
    
    final_subject_table = vertcat(subject_table{:})
    
    
     output_filename = strcat('final_subject_table_',subjects{s},'_' , '.mat');
     save(output_filename,'final_subject_table')
     
     
     final_table = vertcat(final_table, final_subject_table)
     
     clearvars -except s t subjects TMS_timings threshold general_directory final_table summary
     
     
end


%%
cd(general_directory)
filename = 'final_table_for_R_10-2500_def_o1_rep.xlsx';
writetable(final_table,filename,'Sheet',1)