%% 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)