Carica file su 'Script preprocessing EMG-EEG'

Script preprocessing EMG-EEG/M1-P15_peaks extraction.m

@@ -0,0 +1,174 @@
+%% Extract peaks
+clear all
+close all
+
+main_dir = fileparts(mfilename('fullpath')); % path of the script folder
+cd(main_dir)
+load('subjects.mat')
+Paths_file
+
+
+gen_dir = ''; %directory where TEP data is present
+
+
+for iSub=1:28
+    for block = 1:7
+        subjdir_filedef = [gen_dir subjects{iSub} '\' subjects{iSub} '_B_' sprintf('%01d', block) '\' '6_final_preprocessing_new_interp\'];
+        cd(subjdir_filedef)
+        filedefdir = dir(fullfile('*def.mat'));
+        filename = filedefdir.name;
+        load(filename)
+        if contains(filename, 'Mo_PA')
+            Mo_PA{iSub} = TEP_final2;
+        elseif contains(filename, 'Mo_AP')
+            Mo_AP{iSub} = TEP_final2;
+        elseif contains(filename, 'Mo_LM')
+            Mo_LM{iSub} = TEP_final2;
+        elseif contains(filename, 'Bi_PA_Contr')
+            Bi_PA_Contr{iSub} = TEP_final2;
+        elseif contains(filename, 'Bi_AP')
+            Bi_AP{iSub} = TEP_final2;
+        elseif contains(filename, 'Bi_LM')
+            Bi_LM{iSub} = TEP_final2;
+        else
+            Bi_PA{iSub} = TEP_final2;
+        end
+        clear TEP_final2
+    end
+end
+
+allcond = {'Bi_AP', 'Bi_LM', 'Bi_PA', 'Bi_PA_Contr', 'Mo_AP', 'Mo_LM', 'Mo_PA'};
+
+for c = 1:7
+    cond = allcond{c};
+    
+    if contains(cond, 'Mo_PA')
+        TEP_allsubj = Mo_PA;
+    elseif contains(cond, 'Mo_AP')
+        TEP_allsubj = Mo_AP;
+    elseif contains(cond, 'Mo_LM')
+        TEP_allsubj = Mo_LM;
+    elseif contains(cond, 'Bi_PA_Contr')
+        TEP_allsubj = Bi_PA_Contr;
+    elseif contains(cond, 'Bi_AP')
+        TEP_allsubj = Bi_AP;
+    elseif contains(cond, 'Bi_LM')
+        TEP_allsubj = Bi_LM;
+    else
+        TEP_allsubj = Bi_PA;
+    end
+    
+    cfg=[];
+    TEP_gavg = ft_timelockgrandaverage(cfg, TEP_allsubj{:})
+    cfg.keepindividual='yes';
+    TEP_gavg_ind = ft_timelockgrandaverage(cfg, TEP_allsubj{:})
+    
+    cfg = [];
+    cfg.channel={'F4' 'FC4'};
+    cfg.avgoverchan = 'yes';
+    P15 = ft_selectdata(cfg, TEP_gavg_ind);
+    
+    GAVG_all_peak=TEP_gavg_ind;
+    GAVG_all_peak.individual=cat(2,TEP_gavg_ind.individual, ...
+        zeros(length(TEP_gavg_ind.individual(:,1,1)), 1, length(TEP_gavg_ind.individual(1,1,:))), ...
+        zeros(length(TEP_gavg_ind.individual(:,1,1)), 1, length(TEP_gavg_ind.individual(1,1,:)))) %one extrachannel for the marker and one for the average
+    GAVG_all_peak.individual(:, (length(TEP_gavg_ind.label)+1), :) = P15.individual;
+    
+    GAVG_all_peak.label([1:length(TEP_gavg_ind.label)])=TEP_gavg_ind.label;
+    
+    GAVG_all_peak.label{length(TEP_gavg_ind.label)+1}='F4FC4';
+    GAVG_all_peak.label{length(TEP_gavg_ind.label)+2}='peakP15';
+    
+    % P15
+    channel_idx_p15 = [length(TEP_gavg_ind.label)+1];
+    
+    time_window = ones(length(subjects),2);
+    
+    for iI=1:length(TEP_gavg_ind.individual(:,1))
+        
+        % ADJUST WINDOWS FOR EACH SUBJ & CONDITIONS
+        original_window = [0.007 0.025];
+        late_window = [0.015 0.020];
+        early_window = [0.007 0.018];
+        verylate_window = [0.018 0.025];
+        
+        time_window(iI,:) = original_window;
+        
+        if contains(cond,'Bi_PA')
+            
+            if contains(cond, 'Bi_PA_Contr')
+                if iI == 4 || iI == 9
+                    time_window(iI,:) = late_window;
+                elseif iI == 7
+                    time_window(iI,:) = early_window;
+                end
+            else
+                
+                if iI == 4 || iI == 9
+                    time_window(iI,:) = late_window;
+                elseif iI == 23
+                    time_window(iI,:) = early_window;
+                end
+                
+            end      
+            
+        elseif cond == 'Bi_AP'
+            if iI == 23
+                time_window(iI,:) = verylate_window;
+            end
+            
+        elseif cond == 'Bi_LM'
+            if iI == 25
+                time_window(iI,:) = verylate_window;
+            end
+        end
+        
+        
+        samples_p15= find(GAVG_all_peak.time>=time_window(iI,1) & GAVG_all_peak.time<=time_window(iI,2));
+        
+        data_p15 = GAVG_all_peak.individual(iI, channel_idx_p15, samples_p15 ); % pull out the data for this channel and time window
+        
+        %amplitude
+        p_peak_a_p15(iI)=max(data_p15(:,:));
+        
+        %latency
+        p_peak_l_p15(iI) = GAVG_all_peak.time( samples_p15( find( data_p15(:,:)==max(data_p15(:,:)) )  ) );% get the peak latency in this window (max() can be replaced with min() for a negative peak)
+        
+        %calculate window around peak
+        p15_peak_l_min= find(GAVG_all_peak.time==p_peak_l_p15(iI))-25;
+        p15_peak_l_max= find(GAVG_all_peak.time==p_peak_l_p15(iI))+25;
+        p15_peak_l_window(iI,:)=GAVG_all_peak.time(p15_peak_l_min:p15_peak_l_max);
+        p15_peak_a_window(iI,:)=mean(GAVG_all_peak.individual(iI, channel_idx_p15, (p15_peak_l_min:p15_peak_l_max)),3);
+        
+        EEGm_max1(iI)=max(max(GAVG_all_peak.individual(iI,:,:)));
+        GAVG_all_peak.individual(iI, length(TEP_gavg_ind.label)+2, [samples_p15(1) samples_p15(length(samples_p15))])= max(EEGm_max1(iI));
+        GAVG_all_peak.individual(iI, length(TEP_gavg_ind.label)+2, samples_p15( find( data_p15(:,:)==max(data_p15(:,:))))-1:samples_p15( find( data_p15(:,:)==max(data_p15(:,:))))+1)= p_peak_a_p15(iI);
+        
+        
+        
+    end
+    
+    %plot (uncomment the following line after the time windows of all subjects have been adjusted
+    cfg = [];
+    cfg.latency = [-0.02 0.2];
+    tmp = ft_selectdata(cfg, GAVG_all_peak);
+    tmp.trial = {};
+    tmp.trial = GAVG_all_peak.individual;
+    tmp.dimord = 'trial_chan_time';
+    cfg.highlight = 'labels';
+    cfg.viewmode = 'butterfly';
+    cfg.layout = 'EEG1010_AZ.lay';
+    ft_databrowser(cfg, tmp);
+    title(cond)
+    
+    
+    
+    outdir = ''; %output directory for saving files
+    cd(outdir)
+    %%uncomment the following line after the time windows of all subjects have been adjusted
+    P15_amplat = cat(2, p_peak_l_p15', p15_peak_a_window);
+    save([outdir '\P15_peaks_allsubj_' cond '.txt'], 'P15_amplat', '-ascii')
+    
+    clearvars -except Bi_AP Bi_LM Bi_PA Bi_PA_Contr Mo_AP Mo_LM Mo_PA cond c allcond subjects
+    
+end

Script preprocessing EMG-EEG/MOBI_MEP_pipeline.m

@@ -0,0 +1,555 @@
+%%% Pipeline MEP Mobi
+
+clear all
+close all
+
+%main_dir = fileparts(mfilename('fullpath')); % path of the script folder
+cd(main_dir)
+%% paths for PLUG-INs
+
+addpath([main_dir '\PLUG-INs\SOUND']);%SOUND
+addpath([main_dir '\PLUG-INs\TESA1.1.1']); %TESA
+addpath([main_dir '\PLUG-INs\M_functions']); % specific functions for the script
+%path of the raw data
+%addpath('/raw_data');
+
+
+
+% initialize eeglab;
+addpath ('eeglab2022.1/') %% add path eeglab
+eeglab; close; %initialize eeg lab
+
+%initialize fieldtrip
+ft_defaults
+
+%% set folders for input and output: Y
+%%%% folder with data to be processed
+% datafolder_parent = %[main_dir '\mep_data'];  % parent data folder includes all the subdir of MEP dataset to be processed
+% mkdir(main_dir, '\analysis_MEP'); % build output folder
+% output = [main_dir, '\analysis_MEP'];
+
+%% load tables to transcribe them
+Load_tables_MEP
+cd(main_dir)
+
+%% Variables
+
+ref_emg=[2,4];
+epoching_long=[-0.2 0.5];% interval to epoch data, latency in sec relative to time-locking event
+epoching_long_VI=[-0.01 0.06];% interval to epoch data for visual inspection
+rmBase = [-100 -2]; % interval for baseline correction
+SR = 9600;
+%%step1.2
+interp = [0 3]; %interval to cut and interpolate TMS pulse for TESA method, time in ms related to the event
+interp_Win = [0.35,0.35];%times before and after artifact window for fitting cubic function(TESA) 
+HP_f = 1; % % HP filter frequency
+fftflag=1; % filterinf in the frequency domain, faster for high-order filters
+DS= 4800; %down sampling frequency
+tot_cond=7;
+
+%% 0. get data
+[subjects,n_subjects]= get_data_MOB_1(datafolder_parent,output,ref_emg,tot_cond);
+save('subjects.mat', 'subjects')
+save('n_subjects.mat', 'n_subjects')
+
+%% 0.1 HD OFF
+%uncomment these two lines if n_subjects already exist - and comment part
+%0. get data
+%load n_subjects.mat 
+%load subjects.mat 
+
+ %% 1
+for iSub = 1:n_subjects
+    load n_subjects.mat
+    load subjects.mat
+    subj =  subjects{iSub};
+    load([output '\' subj '\0_Data\Block_cond.mat'])
+    
+    for n_cond= 1:tot_cond
+               
+        name_condition = [subj '_' Block_cond{n_cond,2}];
+        name_folder= [subj '\' name_condition];
+        
+
+        mkdir([output, '\' name_folder '\1_Preprocessing\step1_1']);
+        current_output_folder = [output, '\' name_folder '\1_Preprocessing\step1_1'];
+    
+        if  isfile([ current_output_folder '\' Block_cond{n_cond,1} '_EMG_step1_1.set']) ...
+                        %&& global_overwrite==0 && prep1_1_overwrite==0
+                        warning(['skipping step1_1 for ' name_condition ' , file already present'])
+        else
+            load([output '\' subj '\0_Data\' subj '_' Block_cond{n_cond,2}]);
+            data = eval([Block_cond{n_cond,2}]);          
+
+            EMG = pop_importdata...
+            ('dataformat','matlab',...
+            'nbchan',height(data),...
+            'data', 'data',...
+            'setname', [subj '_rawmerged'],...
+            'srate',SR,...
+            'pnts',0,...
+            'xmin',0); 
+     
+            EMG = eeg_checkset( EMG );
+
+                % using the trigger identified from TMS artifact
+                % Read triggers from channel and then remove the channel and events
+                % different from trigger 
+            EMG = pop_chanevent(EMG, length(EMG.data(:,1)),'edge',Block_cond{n_cond,4},'edgelen',0,'delchan','on','typename', num2str(Block_cond{n_cond,5}));
+            EMG = pop_selectevent( EMG, 'type',str2num(Block_cond{n_cond,5}) ,'deleteevents','on');    
+            %per S06 blocco 4: EMG = pop_selectevent( EMG, 'type',3 ,'deleteevents','on');    
+            %S09 blocco 6: EMG = pop_selectevent( EMG, 'type',5 ,'deleteevents','on')
+            EMG = pop_saveset(EMG, 'filename',[Block_cond{n_cond,1} '_EMG_step1_1.set'], 'filepath',current_output_folder );
+         end
+
+        disp('Step 1.1 ended')
+  
+% load EMG data x condition
+% interpolation
+%filters
+% downsampling
+% epoching
+%% STEP 1.2 Remove TMS-pulse artifact and interpolate, HP filter, dwnsampling,epoching 
+
+        previous_folder= [output, '\' name_folder '\1_Preprocessing\step1_1'];
+        mkdir([output, '\' name_folder '\1_Preprocessing\step1_2']);
+        current_output_folder = [output, '\' name_folder '\1_Preprocessing\step1_2'];
+
+        if  isfile([ current_output_folder '\' Block_cond{n_cond,1} '_step1_2.set'])
+                warning(['skipping step1_2 for ' name_condition ' , file already present'])       
+        else
+            EMG = pop_loadset('filename',[Block_cond{n_cond,1} '_EMG_step1_1.set'],'filepath', previous_folder);
+
+            %output_folder = 'F:\Jack\MoBi\pipeline\analysis_MOBI_MEP';
+ 
+            % function to remove TMS pulse with TESA Method
+            % inputs: EEG, interpolation interval,times before and after artifact window for fitting cubic function 
+            EMG= remove_TMS_TESA_G(EMG,interp,interp_Win);
+            %plot with epoching and baseline correction only for visual purpose
+            myplot_F(EMG,[-100 300],[-300 300],['After interpolation TESA'],epoching_long,rmBase);
+            saveas(gcf, [current_output_folder '\after_interpolation_TESA'])
+            close
+
+            %%%%% High pass filter data 
+%             tic
+%             EMG =  high_pass_filter(EMG,HP_f,fftflag);
+%             toc
+            EMG = pop_eegfiltnew(EMG, 'locutoff',10, 'usefftfilt',1);
+            EMG = pop_eegfiltnew(EMG, 'hicutoff',2500, 'usefftfilt',1);
+            EMG = pop_eegfiltnew(EMG, 49.9,50.1,[],1,[],[],[],1);
+
+
+            %plot after filtering without baseline correction
+            myplot_F(EMG,[-100 300],[-500 500],['After filtering'],epoching_long);
+            saveas(gcf, [current_output_folder '\after_filtering'])
+            close
+
+            % DownSample 
+            EMG = pop_resample(EMG, DS);
+
+            myplot_F(EMG,[-100 300],[-500 500],['After resampling'],epoching_long);
+            saveas(gcf, [current_output_folder '\after_resampling'])
+            close
+
+            % Epoching data 
+            EMG = pop_epoch(EMG, {}, epoching_long);
+
+            %save figure
+            myplot_F(EMG,[-100 300],[-500 500],['After Filtering']);
+            saveas(gcf,[current_output_folder  '\after_filtering'])
+            close
+
+            savename = [Block_cond{n_cond,1} '_step1_2.set'];
+            EMG = pop_saveset( EMG, 'filename',savename,'filepath',current_output_folder);
+            
+        end
+        % save set file which is closer to MEPs time range for latency visual inspection
+        %if isfile([ current_output_folder '\' Block_cond{n_cond,1} '_step1_2_VisualInspection.set'])
+                %warning(['skipping step1_2 for ' name_condition ' , file already present']) 
+        %else  
+            EMG = pop_loadset('filename',[Block_cond{n_cond,1} '_step1_2.set'],'filepath',current_output_folder);
+            EMG2 = pop_epoch(EMG, {}, epoching_long_VI);
+            savename = [Block_cond{n_cond,1} '_step1_2_VisualInspection.set'];
+            EMG2 = pop_saveset( EMG2, 'filename',savename,'filepath',current_output_folder); 
+%         end
+        
+        disp('Step 1.2 ended')
+    end
+end
+
+
+%% STEP 2. ESTIMATE MEP AMPLITUDE x CONDITIONS
+%TW = cell2table(cell(n_subjects,7), 'VariableNames', {'Mo_PA','Mo_AP','Mo_LM','Bi_PA_Contr','Bi_AP','Bi_LM','Bi_PA'});
+
+
+for iSub = 1:n_subjects
+    load n_subjects.mat
+    load subjects.mat
+    subj =  subjects{iSub}
+    load([output '\' subj '\0_Data\Block_cond.mat'])
+    
+    for n_cond= 1:tot_cond
+               
+        name_condition = [subj '_' Block_cond{n_cond,2}];
+        name_folder= [subj '\' name_condition];
+        
+        previous_folder= [output, '\' name_folder '\1_Preprocessing\step1_2'];
+        mkdir([output, '\' name_folder '\2_MEP_AmplitudeP2P']);
+        current_output_folder = [output, '\' name_folder '\2_MEP_AmplitudeP2P'];
+    
+%         if  isfile([ previous_folder '\' Block_cond{n_cond,1} '_EMG_step1.set']) ...
+%                         %&& global_overwrite==0 && prep1_1_overwrite==0
+%                         warning(['skipping step2 for ' name_condition ' , file already present'])
+%         else
+          
+          EMG = pop_loadset('filename',[Block_cond{n_cond,1} '_step1_2.set'],'filepath',previous_folder);
+
+    condition = Block_cond{n_cond,1};
+    numtrials = length(EMG.data(1,1,:));
+    nomep_trials = 0;
+    handchan = 2; %RHAND
+       
+    p2p_mep = [];
+
+ for itrial = 1:numtrials
+%% 3.4) amplitude estimation
+   
+    a = EMG.data(handchan,:,itrial); %channel x time point x current trial
+
+    start_mep = 19;
+    end_mep = 77;
+    
+    inizio_mep = find(EMG.times>=start_mep,1,"first");
+    fine_mep = find(EMG.times>=end_mep,1,"first");
+    
+
+    z = a(inizio_mep:fine_mep);
+    
+    [minV,tpmin]= min(z);
+    [maxV,tpmax] = max(z); 
+
+    %time window
+    time_window = EMG.times(inizio_mep:fine_mep);
+    
+    %points in ms
+    tmin = time_window(tpmin);
+    tmax = time_window(tpmax);
+        
+    % peak2peak amplitude 
+    p2p = sprintf('%10.2f',(maxV - minV));
+    p2p(p2p==' ') =[];
+    p2p = str2num(p2p);
+    
+    %check trials without MEPs (exclude trial <50µV)
+    p2p50 = p2p;
+    if p2p50 <= 50
+        p2p50 = NaN;
+        nomep_trials = nomep_trials + 1;   
+    end
+    
+   
+    %amplitude x trial
+    p2p_mep(itrial) = p2p50;
+    p2p_all(itrial) = p2p;
+
+
+ end
+ 
+p2p_meptrials = p2p_mep;
+p2p_allamp = p2p_all;
+%save([current_output_folder '\' Block_cond{iSub,1} '_'  Block_cond{iSub,2}  '_p2p_meptrials.mat'],'p2p_meptrials'); 
+%save([current_output_folder '\' Block_cond{iSub,1} '_'  Block_cond{iSub,2} '_p2p_allamp.mat'],'p2p_allamp'); 
+save([current_output_folder '\p2p_meptrials.mat'],'p2p_meptrials'); 
+save([current_output_folder '\p2p_allamp.mat'],'p2p_allamp'); 
+
+
+
+ %waitfor(gcf)
+% output tables
+if contains(condition,"Mo_PA")
+         MEP_p2pAmp{iSub,1} = p2p_mep; 
+         nomep_sub(iSub, 1) = nomep_trials;
+         MEP_mean(iSub,1) = nanmean(p2p_mep);
+    elseif contains(condition,"Mo_AP")
+         MEP_p2pAmp{iSub,2} = p2p_mep;
+         nomep_sub(iSub, 2) = nomep_trials;
+         MEP_mean(iSub,2) = nanmean(p2p_mep);
+    elseif contains(condition,"Mo_LM")
+         MEP_p2pAmp{iSub,3} = p2p_mep;
+         nomep_sub(iSub, 3) = nomep_trials;
+         MEP_mean(iSub,3) = nanmean(p2p_mep);
+    elseif contains(condition,"Bi_PA_Contr")
+         MEP_p2pAmp{iSub,4} = p2p_mep;
+         nomep_sub(iSub, 4) = nomep_trials;
+         MEP_mean(iSub,4) = nanmean(p2p_mep);
+    elseif contains(condition,"Bi_AP")
+         MEP_p2pAmp{iSub,5} = p2p_mep;
+         nomep_sub(iSub, 5) = nomep_trials;
+         MEP_mean(iSub,5) = nanmean(p2p_mep);
+    elseif contains(condition,"Bi_LM")
+         MEP_p2pAmp{iSub,6} = p2p_mep;
+         nomep_sub(iSub, 6) = nomep_trials;
+         MEP_mean(iSub,6) = nanmean(p2p_mep);
+    elseif contains(condition,"Bi_PA")
+         MEP_p2pAmp{iSub,7} = p2p_mep;
+         nomep_sub(iSub, 7) = nomep_trials;
+         MEP_mean(iSub,7) = nanmean(p2p_mep);
+         
+         %aggiungo questo per S38  
+%     elseif contains(condition,"BI_PA")
+%         MEP_p2pAmp{iSub,7} = p2p_mep;
+%         nomep_sub(iSub, 7) = nomep_trials;
+%         MEP_mean(iSub,7) = nanmean(p2p_mep);
+end
+    
+%conditions in randomization order:
+MEP_p2pAmp_rc{iSub,n_cond} = p2p_mep;
+MEP_p2pAll_rc{iSub,n_cond} = p2p_all;%all amplitudes (also lower than 50uV)
+
+%     %end
+    end
+
+cd(output)
+TWMean(iSub,:) = array2table(MEP_mean(iSub,:),'VariableNames',{'Mo_PA','Mo_AP','Mo_LM','Bi_PA_Contr','Bi_AP','Bi_LM','Bi_PA'})
+TWP2P(iSub,:) = cell2table(MEP_p2pAmp(iSub,:),'VariableNames',{'Mo_PA','Mo_AP','Mo_LM','Bi_PA_Contr','Bi_AP','Bi_LM','Bi_PA'})
+TWNoM(iSub,:) = array2table(nomep_sub(iSub,:),'VariableNames',{'Mo_PA','Mo_AP','Mo_LM','Bi_PA_Contr','Bi_AP','Bi_LM','Bi_PA'})
+TQAMP(iSub,:) = cell2table(MEP_p2pAmp_rc(iSub,:),'VariableNames',{'B_1','B_2','B_3','B_4','B_5','B_6','B_7'})
+TQALL(iSub,:) = cell2table(MEP_p2pAll_rc(iSub,:),'VariableNames',{'B_1','B_2','B_3','B_4','B_5','B_6','B_7'})
+  
+end
+
+%aggiungo parentesi e aggiungo nel loop
+% TWMean = array2table(MEP_mean,'VariableNames',{'Mo_PA','Mo_AP','Mo_LM','Bi_PA_Contr','Bi_AP','Bi_LM','Bi_PA'})
+% TWP2P = cell2table(MEP_p2pAmp,'VariableNames',{'Mo_PA','Mo_AP','Mo_LM','Bi_PA_Contr','Bi_AP','Bi_LM','Bi_PA'})
+% TWNoM = array2table(nomep_sub,'VariableNames',{'Mo_PA','Mo_AP','Mo_LM','Bi_PA_Contr','Bi_AP','Bi_LM','Bi_PA'})
+% TQAMP = cell2table(MEP_p2pAmp_rc,'VariableNames',{'B_1','B_2','B_3','B_4','B_5','B_6','B_7'})
+% TQALL = cell2table(MEP_p2pAll_rc,'VariableNames',{'B_1','B_2','B_3','B_4','B_5','B_6','B_7'})
+
+
+save([output '\nomep_sub.mat'],'TWNoM'); %numero trial senza MEP
+save([output '\MEP_p2pAmp.mat'],'TWP2P');% ampiezza mep picco-picco %ordine: {'Mo_PA','Mo_AP','Mo_LM','Bi_PA_Contr','Bi_AP','Bi_LM','Bi_PA'}
+save([output '\MEP_mean.mat'],'TWMean'); %media dei trial x soggetto x condizione %ordine: {'Mo_PA','Mo_AP','Mo_LM','Bi_PA_Contr','Bi_AP','Bi_LM','Bi_PA'}
+save([output '\MEP_p2pAmp_randomorder.mat'],'TQAMP'); % {'B_1','B_2''B_3''B_4''B_5''B_6''B_7'}
+save([output '\MEP_p2pAll_randomorder.mat'],'TQALL'); 
+
+%% STEP 3. ESTIMATE MEP LATENCY x CONDITIONS
+
+% main_dir = 'C:\Users\delia\OneDrive\Documenti\MSc-cognitive neuroscience\Thesis\MEP';
+%datafolder_parent = 'C:\Users\delia\OneDrive\Documenti\MSc-cognitive neuroscience\Thesis\MEP'; 
+%'F:\Jack\MoBi\Data_MEP'; %[main_dir '\mep_data'];  % parent data folder includes all the subdir of MEP dataset to be processed
+
+
+for iSub = 1:n_subjects 
+    load n_subjects.mat
+    load subjects.mat
+    subj =  subjects{iSub};
+    load([output '\' subj '\0_Data\Block_cond.mat'])
+    
+    for n_cond = 1:tot_cond
+               
+        name_condition = [subj '_' Block_cond{n_cond,2}];
+        name_folder= [subj '\' name_condition];
+        
+        previous_folder= [output, '\' name_folder '\1_Preprocessing\step1_2'];
+        
+        %load amplitude mat file to get NaN amplitudes 
+        amplitude_folder= [output, '\' name_folder '\2_MEP_AmplitudeP2P'];
+        load([amplitude_folder,'\p2p_allamp.mat'])
+        
+        %create new folder for latency
+        mkdir([output, '\' name_folder '\3_MEP_Latency']);
+        current_output_folder = [output, '\' name_folder '\3_MEP_Latency'];
+        
+        %create new folder for plots
+        mkdir([output, '\' name_folder '\3_MEP_Latency\plots']);
+        plots_output_folder = [output, '\' name_folder '\3_MEP_Latency\plots'];
+    
+%         if  isfile([ previous_folder '\' Block_cond{n_cond,1} '_EMG_step1.set']) ...
+%                         %&& global_overwrite==0 && prep1_1_overwrite==0
+%                         warning(['skipping step2 for ' name_condition ' , file already present'])
+%         else
+
+        %load set file
+        EMG = pop_loadset('filename',[Block_cond{n_cond,1} '_step1_2.set'],'filepath',previous_folder);
+        %cut epochs from 10ms to 50ms
+        epoching_long=[0.01 0.05];% interval to epoch data, latency in sec relative to time-locking event
+        EMG = pop_epoch(EMG, {}, epoching_long);
+       
+        %variables
+        condition = Block_cond{n_cond,1};
+        numtrials = length(EMG.data(1,1,:));
+        handchan = 2;        
+        lat_mep = [];
+        lat_all = [];
+    
+    
+    for itrial = 1:numtrials
+
+data = EMG.data(handchan,:,itrial); %channel x time point x current trial
+data=data';
+time=EMG.times'; 
+
+% reverse data if the first peak is negative
+[MA, IA] = max(data);
+[MI, II] = min(data);
+if II<IA
+data=-data;
+%compute again IA after data has been reversed 
+[MA, IA] = max(data);
+end
+
+%find data from the beginning of the epoch to MEP's first peak 
+z = data(1:IA)
+
+% calculate derivative on a single trial
+ der1=diff([z(:); 0]); % to add one number at the end to have same samples as data
+
+%plot the derivative and MEP signal 
+%bo=cat(2, z, der1);
+%plot(bo)
+
+% to find longest vector of positive values in the derivative
+cons=der1>0; % to find  points of positive derivative (0 non positive; 1 positive)
+
+% Label each separate region/cluster of ones.
+[labeledX, numRegions] = bwlabel(cons);
+% Get lengths of each region
+props = regionprops(labeledX, 'Area', 'PixelList');
+regionLengths = [props.Area];
+k = find(regionLengths==max(regionLengths)); % to find which is the biggest cluster
+
+%if one max cluster is found, calculate MEP onset and plot it
+if length(k)==1
+    fprintf('Region #%d is %d elements long and starts at element #%d\n',...
+        k, regionLengths(k), props(k).PixelList(1,2));
+    start= props(k).PixelList(1,2); %index of start for longest cluster
+    time(start) % MEP onset
+    % to visualize 
+    MEPonset =zeros(size(z)); %to create a second channel
+    MEPonset(start)=100;   
+    % to plot figure
+    plot(cat(2, z, der1, MEPonset));    
+    %to save the plot
+    ax = gca;
+    numtrial=num2str(itrial);
+
+    %if the trial has no mep, name the image with "NO MEP"    
+    if p2p_allamp(itrial) <= 50
+        exportgraphics(ax,[plots_output_folder,'\latplot',numtrial,'_NO_MEP.jpg'])
+    elseif p2p_allamp(itrial) > 50
+        exportgraphics(ax,[plots_output_folder,'\latplot',numtrial,'.jpg'])
+    end
+    v=EMG.times(1:IA);
+    latency=v(start);
+    disp([ ' MEP latency in condition ', condition, ', trial ', num2str(itrial), ', is ' num2str(latency), ' ms' ])
+    latency50=latency;
+    %if it is not a MEP, save it as NaN
+    if p2p_allamp(itrial) <= 50
+        latency50 = NaN
+    end
+    
+        
+%if more than 1 max clusters are found, take into consideration the last one (in terms of time)      
+elseif length(k)>1
+    k=max(k);
+    fprintf('Region #%d is %d elements long and starts at element #%d\n',...
+        k, regionLengths(k), props(k).PixelList(1,2));
+    start= props(k).PixelList(1,2); %index of start for longest cluster
+    time(start) % MEP onset    
+    % to visualize
+    MEPonset =zeros(size(z)); %to create a second channel
+    MEPonset(start)=100;   
+    % to plot figure
+    plot(cat(2, z, der1, MEPonset));    
+    %to save the plot
+    ax = gca;
+    numtrial=num2str(itrial);
+    
+    %all images here are saved with "2P" (two peaks) to check them later on
+    if p2p_allamp(itrial) <= 50
+        exportgraphics(ax,[plots_output_folder,'\latplot',numtrial,'_NO_MEP_2P.jpg'])
+    elseif p2p_allamp(itrial) > 50
+        exportgraphics(ax,[plots_output_folder,'\latplot',numtrial,'_2P.jpg'])
+    end
+    v=EMG.times(1:IA);
+    latency=v(start);
+    disp([ ' MEP latency in condition ', condition, ', trial ', num2str(itrial), ', is ' num2str(latency), ' ms' ])
+    latency50=latency;    
+ 
+    if p2p_allamp(itrial) <= 50
+        latency50 = NaN
+    end    
+
+%if k is empty, the peak cannot be computed    
+elseif isempty(k)==1
+     warning(['skipping latency estimation for ' name_condition ' , missing minimum input arguments'])
+     latency=NaN; 
+     latency50=NaN;    
+end  
+
+    lat_mep(itrial)=latency50;
+    lat_all(itrial)=latency;
+    
+    end
+
+    lat_meptrials=lat_mep;
+    lat_alltrials=lat_all;
+    save([current_output_folder '\lat_meptrials.mat'],'lat_meptrials'); 
+    save([current_output_folder '\lat_alltrials.mat'],'lat_alltrials');
+    
+    if contains(condition,"Mo_PA")
+         MEP_lat{iSub,1} = lat_mep; 
+         %nomep_sub(iSub, 1) = nomep_trials;
+         MEAN_lat(iSub,1) = nanmean(lat_mep);
+    elseif contains(condition,"Mo_AP")
+         MEP_lat{iSub,2} = lat_mep;
+         %nomep_sub(iSub, 2) = nomep_trials;
+         MEAN_lat(iSub,2) = nanmean(lat_mep);
+    elseif contains(condition,"Mo_LM")
+         MEP_lat{iSub,3} = lat_mep;
+         %nomep_sub(iSub, 3) = nomep_trials;
+         MEAN_lat(iSub,3) = nanmean(lat_mep);
+    elseif contains(condition,"Bi_PA_Contr")
+         MEP_lat{iSub,4} = lat_mep;
+         %nomep_sub(iSub, 4) = nomep_trials;
+         MEAN_lat(iSub,4) = nanmean(lat_mep);
+    elseif contains(condition,"Bi_AP")
+         MEP_lat{iSub,5} = lat_mep;
+         %nomep_sub(iSub, 5) = nomep_trials;
+         MEAN_lat(iSub,5) = nanmean(lat_mep);
+    elseif contains(condition,"Bi_LM")
+         MEP_lat{iSub,6} = lat_mep;
+         %nomep_sub(iSub, 6) = nomep_trials;
+         MEAN_lat(iSub,6) = nanmean(lat_mep);
+    elseif contains(condition,"Bi_PA")
+         MEP_lat{iSub,7} = lat_mep;
+         %nomep_sub(iSub, 7) = nomep_trials;
+         MEAN_lat(iSub,7) = nanmean(lat_mep);
+         %aggiungo questo per S38
+    elseif contains(condition,"BI_PA")
+        MEP_lat{iSub,7} = lat_mep;
+        %nomep_sub(iSub, 7) = nomep_trials;
+        MEAN_lat(iSub,7) = nanmean(lat_mep);
+    end
+
+    %conditions in randomization order:
+MEP_lat_rc{iSub,n_cond} = lat_mep;
+MEP_latAll_rc{iSub,n_cond} = lat_all;%all amplitudes (also lower than 50uV)
+    
+    end
+
+cd(output)
+
+TOMeanLat(iSub,:) = array2table(MEAN_lat(iSub,:),'VariableNames',{'Mo_PA','Mo_AP','Mo_LM','Bi_PA_Contr','Bi_AP','Bi_LM','Bi_PA'})
+TOLat(iSub,:) = cell2table(MEP_lat(iSub,:),'VariableNames',{'Mo_PA','Mo_AP','Mo_LM','Bi_PA_Contr','Bi_AP','Bi_LM','Bi_PA'})
+TRLat(iSub,:) = cell2table(MEP_lat_rc(iSub,:),'VariableNames',{'B_1','B_2','B_3','B_4','B_5','B_6','B_7'})
+TRALLat(iSub,:) = cell2table(MEP_latAll_rc(iSub,:),'VariableNames',{'B_1','B_2','B_3','B_4','B_5','B_6','B_7'})
+     
+end
+
+
+save([output '\MEP_Lat.mat'],'TOLat');% latenza mep %ordine: {'Mo_PA','Mo_AP','Mo_LM','Bi_PA_Contr','Bi_AP','Bi_LM','Bi_PA'}
+save([output '\MEP_meanLat.mat'],'TOMeanLat'); %media dei trial x soggetto x condizione %ordine: {'Mo_PA','Mo_AP','Mo_LM','Bi_PA_Contr','Bi_AP','Bi_LM','Bi_PA'}
+save([output '\MEP_Lat_randomorder.mat'],'TRLat'); % {'B_1','B_2''B_3''B_4''B_5''B_6''B_7'}
+save([output '\MEP_AlLat_randomorder.mat'],'TRALLat');
+

Script preprocessing EMG-EEG/preprocessingTEP_pipeline.m

@@ -0,0 +1,426 @@
+ %%% Pipeline case gtec files
+
+clear all
+close all
+
+%main_dir = fileparts(mfilename('fullpath')); % path of the script folder
+
+%% paths for PLUG-INs functions
+
+addpath([main_dir '\PLUG-INs\SOUND']);%SOUND
+addpath([main_dir '\PLUG-INs\TESA1.1.1']); %TESA
+addpath([main_dir '\PLUG-INs\M_functions']); % specific functions for the script
+%path of the raw data
+%addpath('/raw_data');
+
+
+
+% initialize eeglab;
+%addpath ('eeglab2022.1/') %% add path eeglab
+%eeglab; close; %initialize eeg lab
+
+%initialize fieldtrip
+ft_defaults
+
+%% set folders for input and output: Y
+%%%% folder with data to be processed
+% datafolder_parent = [main_dir '\raw_data'];  % parent data folder includes all the subdir of all dataset to be processed
+% mkdir(main_dir, '\analysis_MOBI'); % build output folder
+% output = [main_dir, '\analysis_MOBI'];
+
+%% Variables
+
+%%control to overwrite steps
+global_overwrite=0 ;% 1 to overwrite all steps
+prep1_1_overwrite=0; % 1 to overwrite step 1.1
+prep1_2_overwrite=0; % 1 to overwrite step 1.2
+sound_overwrite=0; %1 to overwrite step 2 SOUND
+rej_overwrite=0; % 1 to overwrite step 3 artifact rejection
+ICA_overwrite=0; % 1 to overwrite step 4.1 ICA
+ICA_remove=0; %  1 to overwrite step 4.2 ICA selection components
+SSP_SIR_step1 =0; % 1 to overwrite step 5.1 SSP-SIR
+SSP_SIR_step2 =0; % 1 to overwrite step 5.2 SSP-SIR components selection
+
+%%step 1.1
+chan_loc_st = 'standard-10-5-cap385.elp';% standard 10-20 channel location
+load([main_dir '\PLUG-INs\M_functions\EEG_chanlocs_73' ]); chan_loc=EEG_chanlocs_73;%custom channel location
+ref_chan = 2;%row of the reference channel
+emg_refchan= [76, 78];
+emg_chan = [77, 79];
+eeg_channels=73;
+tot_cond=7;
+
+epoching_long=[-0.7 0.7];% interval to epoch data, latency in sec relative to time-locking event
+rmBase = [-200 -5]; % interval for baseline correction
+
+%%step1.2
+interp = [0 3]; %interval to cut and interpolate TMS pulse for TESA method, time in ms related to the event
+interp_Win = [0.35,0.35];%times before and after artifact window for fitting cubic function(TESA) 
+HP_f = 1; % % HP filter frequency
+fftflag=1; % filterinf in the frequency domain, faster for high-order filters
+DS= 4800; %down sampling frequency
+
+
+%%step2 SOUND
+ref_before_sound= 'FPz';
+lambda_value = 0.1; % SOUND parameter
+iter = 5;%SOUND parameter
+
+%%step 3 artifact rejection
+SD_first=5; %number of stnd dev to consider when computing joint probability and remove trial (first trial rejection)
+SD_last=5; %number of stnd dev to consider when computing joint probability and remove trial(last trial rejection)
+
+%%step 4.2 ICA
+threshold = [0 0;0 0; 0.5 1; 0 0; 0 0; 0 0; 0 0];% threshold to label eye components (50%)
+
+%%step 5 SSP-SIR
+SSP_SIR_window = [0,50];%the time window on which SSP-SIR is applied 
+limx=[-100 350];% xlim for figure with TEPs
+limy = [-30 30]; %ylim for figure with TEPs
+
+%%step 6 final preprocesing
+LP_f = 70;% frequency in Hz for lowpass filter
+%chan_loc2 = 'standard-10-5-cap385.elp';% standard 10-20 channel location
+old_ref = 'FPz';
+reref = {'TP9' 'TP10'};
+epoching_short = [-0.2 0.4];
+latency = [-0.1 0.4]; %epoch limit in final fieldtrip set
+ylimit=[-15 15];%ylimit for final figure
+
+
+%% 0. Get file ready for processing
+
+% % function to build the mat file row data
+[subjects,n_subjects]= get_data_MOBI(datafolder_parent,output,ref_chan,emg_refchan,tot_cond);
+save('subjects.mat', 'subjects')
+save('n_subjects.mat', 'n_subjects')
+
+%% 1. BASIC PREPROCESSING
+% Step 1.1
+% load set, add channel location, remove unwanted electrodes, define
+% reference, save the set
+
+%uncomment these two lines if n_sibjects.mat already exists (and skip/comment
+%section 0.):
+%load n_subjects.mat 
+%load subjects.mat 
+
+for iSub = 1:n_subjects
+    subj =  subjects{iSub}
+    load([output '\' subj '\0_Data\Block_cond.mat'])
+    
+    for n_cond= 1:tot_cond
+               
+        name_condition = [subj '_' Block_cond{n_cond,2}];
+        name_folder= [subj '\' name_condition];
+        
+        mkdir([output, '\' name_folder '\1_Preprocessing\step1_1']);
+        current_output_folder = [output, '\' name_folder '\1_Preprocessing\step1_1'];
+    
+        if  isfile([ current_output_folder '\' Block_cond{n_cond,1} '_step1_1.set']) ...
+                && isfile([current_output_folder '\' Block_cond{n_cond,1} '_EMG_step1_1.set'])...
+                && global_overwrite==0 && prep1_1_overwrite==0
+                warning(['skipping step1_1 for ' name_condition ' , file already present'])
+        else
+            %load row data
+            load([output '\' subj '\0_Data\' subj '_' Block_cond{n_cond,2}]);
+            data = eval([Block_cond{n_cond,2}]);          
+
+            % load raw data,discard EMG channels,read the trigger,load channel names
+            [EEG, EMG]= load_rawdata_MOBI(data,name_folder,Block_cond{n_cond,3},Block_cond{n_cond,5},Block_cond{n_cond,4},eeg_channels,chan_loc);
+
+            %plot before starting preprocessing signals. Epoching and baseline correction only
+            %for visual purpose
+            myplot_F(EEG,[-100 500],[-1000 45000],['Step 1.1 signals'],epoching_long);
+            saveas(gcf, [current_output_folder '\step1_1'] )
+            close
+            
+            EEG = pop_saveset(EEG, 'filename',[Block_cond{n_cond,1} '_step1_1'],'filepath', current_output_folder);
+            EMG = pop_saveset(EMG, 'filename',[Block_cond{n_cond,1} '_EMG_step1_1.set'], 'filepath',current_output_folder );
+
+        end
+
+        disp('Step 1.1 ended')
+   
+
+%% STEP 1.2 Remove TMS-pulse artifact and interpolate, HP filter, dwnsampling,epoching 
+        previous_folder= [output, '\' name_folder '\1_Preprocessing\step1_1'];
+        mkdir([output, '\' name_folder '\1_Preprocessing\step1_2']);
+        current_output_folder = [output, '\' name_folder '\1_Preprocessing\step1_2'];
+
+        if  isfile([ current_output_folder '\' Block_cond{n_cond,1} '_step1_2.set']) && global_overwrite==0 && prep1_2_overwrite ==0
+                warning(['skipping step1_2 for ' name_condition ' , file already present'])       
+        else
+            EEG = pop_loadset('filename',[Block_cond{n_cond,1} '_step1_1.set'],'filepath', previous_folder);
+
+            % function to remove TMS pulse with TESA Method
+            % inputs: EEG, interpolation interval,times before and after artifact window for fitting cubic function 
+           EEG= remove_TMS_TESA_G(EEG,interp,interp_Win);
+           %plot with epoching and baseline correction only for visual purpose
+           myplot_F(EEG,[-100 300],[-300 300],['After interpolation TESA'],epoching_long,rmBase);
+           saveas(gcf, [current_output_folder '\after_interpolation_TESA'])
+           close
+
+            %%%%% High pass filter data 
+            tic
+            EEG =  high_pass_filter(EEG,HP_f,fftflag);
+            toc
+
+            %plot after filtering without baseline correction
+            myplot_F(EEG,[-100 300],[-500 500],['After filtering'],epoching_long);
+            saveas(gcf, [current_output_folder '\after_filtering'])
+            close
+
+             % DownSample 
+            EEG = pop_resample(EEG, DS);
+
+            myplot_F(EEG,[-100 300],[-500 500],['After resampling'],epoching_long);
+            saveas(gcf, [current_output_folder '\after_resampling'])
+            close
+
+            % Epoching data 
+            EEG = pop_epoch(EEG, {}, epoching_long);
+
+            %save figure
+            myplot_F(EEG,[-100 300],[-500 500],['Step 1.2']);
+            saveas(gcf,[current_output_folder '\step1_2'])
+            close
+
+            savename = [Block_cond{n_cond,1} '_step1_2'];
+            EEG = pop_saveset( EEG, 'filename',savename,'filepath',current_output_folder);
+    
+        end
+        disp('Step 1.2 ended')
+
+
+%% STEP 2 SOUND
+% 2.1 - prepare for Sound
+        previous_folder= [output, '\' name_folder '\1_Preprocessing\step1_2'];
+        mkdir([output, '\' name_folder '\2_SOUND']);
+        current_output_folder = [output, '\' name_folder '\2_SOUND'];
+
+        if  isfile([ current_output_folder '\' Block_cond{n_cond,1} '_step2_2.set']) && global_overwrite==0 && sound_overwrite==0
+                warning(['skipping step 2 for ' name_condition ' , file already present'])
+        else
+            EEG = pop_loadset('filename',[Block_cond{n_cond,1} '_step1_2.set'],'filepath',previous_folder);
+
+            %function to construct leadfield and  to re-referencing 
+            [EEG,LFM_sphere]= before_SOUND_G(EEG,ref_before_sound,chan_loc_st);
+
+            % plot butterfly and save before SOUND     
+            myplot_F(EEG,[-100 300],[-50 50],['Before SOUND']);
+            saveas(gcf,[current_output_folder '\' name_condition '_before_SOUND'])
+            close 
+
+            %save dateset and leadfield matrix
+            %savename = [name_condition '_beforeSound'];
+            %EEG = pop_saveset(EEG, 'filename',savename,'filepath', current_output_folder);
+            % save leadfield
+            save([current_output_folder '\LFM_sphere.txt'], 'LFM_sphere', '-ascii'); 
+
+            disp('Step 2.1  pre-SOUND ended');
+
+            %%%% 2.2 SOUND
+
+            %function for SOUND algorithm
+            EEG_clean = SOUND_to_samples(EEG,LFM_sphere,lambda_value,iter);
+            delete(gcp)
+            close
+
+            EEG = EEG_clean;
+
+            % plot butterfly and save after SOUND
+            myplot_F(EEG,[-100 300],[-50 50],['After SOUND']);
+            saveas(gcf,[current_output_folder '\step2_2'])
+            close
+
+            %save set after SOUND
+            EEG.setname= [Block_cond{n_cond,1} '_afterSOUND'];
+            EEG = pop_saveset(EEG, 'filename',[Block_cond{n_cond,1} '_step2_2'],'filepath', current_output_folder);
+        end
+
+disp(' step 2.2 SOUND ended')
+
+%% 3 Automated Artifact rejection
+
+        previous_folder= [output, '\' name_folder '\2_SOUND'];
+        mkdir([output, '\' name_folder  '\3_Artifact_rejection']);
+        current_output_folder = [output, '\' name_folder  '\3_Artifact_rejection'];
+
+        if  isfile([ current_output_folder '\' Block_cond{n_cond,1}  '_step3.set']) && global_overwrite==0 && rej_overwrite==0
+                warning(['skipping step 3 for ' name_condition  ' , file already present'])
+
+        else
+            EEG = pop_loadset('filename',[Block_cond{n_cond,1}  '_step2_2.set'],'filepath',previous_folder);
+
+            % Automated Trial rejetion
+            [EEG,rejected_epochs] = reject_artifact(EEG,SD_first,SD_last);  
+
+            %save rejected epochs
+            save([current_output_folder '\rejected_epochs'], 'rejected_epochs');
+
+            %plot
+            myplot_F(EEG,[-100 300],[-50 50],['After trial rejction']);
+            saveas(gcf,[current_output_folder '\after_trialrejction'])
+            close
+
+            % save data after analysisG rejection:
+            EEG.setname= 'trialrejction';
+            EEG = pop_saveset(EEG, 'filename',[Block_cond{n_cond,1}  '_step3'],'filepath', current_output_folder);
+        end
+
+disp('step 3 Artifact rejection ended')
+
+%% 4 ICA 
+% 4.1 RUN ICA
+        previous_folder= [output, '\' name_folder  '\3_Artifact_rejection'];
+        mkdir([output,'\' name_folder '\4_ICA\4_1_RUNICA']);
+        current_output_folder = [output,'\' name_folder '\4_ICA\4_1_RUNICA'];
+
+        if  isfile([ current_output_folder '\' Block_cond{n_cond,1} '_step4_1.set']) && global_overwrite==0 && ICA_overwrite==0
+                warning(['skipping step 4.1 for ' name_condition ' , file already present'])
+
+        else
+                EEG = pop_loadset('filename',[Block_cond{n_cond,1} '_step3.set'],'filepath',  previous_folder);
+                
+                %EEG = pop_tesa_removedata( EEG,interp);
+
+                % Ica decomposition 
+                EEG = pop_runica(EEG, 'chanind', [1: size(EEG.data,1)], 'interupt','on');
+
+                %save set with ICA all components
+                EEG.setname= 'After ICA';
+                savename = [Block_cond{n_cond,1} '_step4_1'];
+                EEG = pop_saveset(EEG, 'filename',savename,'filepath', current_output_folder);    
+        end
+disp('Step 4.1 ICA ended')
+    end
+end
+%%  4.2 ICA - Remove ocular artifat with ICA
+for iSub = 1:n_subjects
+    subj =  subjects{iSub}
+    for n_cond =1:tot_cond
+        
+        load([output '\' subj '\0_Data\Block_cond.mat']);
+        name_condition = [subj '_' Block_cond{n_cond,2}];
+        name_folder= [subj '\' name_condition]
+
+            previous_folder= [output,'\' name_folder '\4_ICA\4_1_RUNICA'];
+            mkdir([output, '\' name_folder '\4_ICA\4_2_Remove_ocular_Artifact']);
+            current_output_folder = [output, '\' name_folder '\4_ICA\4_2_Remove_ocular_Artifact'];
+
+        if  isfile([ current_output_folder '\' Block_cond{n_cond,1} '_step4_2.set']) && global_overwrite==0 && ICA_remove==0
+                warning(['skipping step 4.2 for ' name_condition ' , file already present'])
+            else
+             EEG = pop_loadset('filename',[Block_cond{n_cond,1} '_step4_1.set'],'filepath',  previous_folder);
+             
+             Select_ICA_components       
+        end
+    disp('Step 4.2 ICA selection components ended')
+    end
+end
+%% 5. Use SSP-SIR to clean muscle artifacts
+%  SSP-SIR
+for iSub = 1:n_subjects
+    subj =  subjects{iSub}
+    for n_cond= 1:tot_cond
+        load([output '\' subj '\0_Data\Block_cond.mat']);
+        name_condition = [subj '_' Block_cond{n_cond,2}];
+        name_folder= [subj '\' name_condition]
+
+        previous_folder= [output, '\' name_folder '\4_ICA\4_2_Remove_ocular_Artifact'];
+        mkdir([output, '\' name_folder '\5_SSP-SIR\5_1_SSP-SIR calculations']);
+        current_output_folder = [output, '\' name_folder '\5_SSP-SIR\5_1_SSP-SIR calculations'];
+
+        if  size(dir([current_output_folder '/*.set' ]),1)== 6 && global_overwrite==0 && SSP_SIR_step1==0
+                warning(['skipping step 5.1 for ' name_condition ' , files already present'])
+        else
+            % re-load  LFM_sphere and EEG from previuos step
+            LFM_sphere= load([ output '\' name_folder '\2_SOUND\LFM_sphere.txt']);
+            EEG = pop_loadset('filename',[Block_cond{n_cond,1} '_step4_2.set'],'filepath',  previous_folder);    
+            
+            EEG = pop_tesa_removedata( EEG,interp);
+
+            % SSP-SIR with TEP with incremental number of removed components, up to 5
+            SSP_SIR_components(EEG,LFM_sphere,SSP_SIR_window, current_output_folder,name_condition,rmBase,limx,limy)
+        end
+
+disp('Step 5.1 SSP-SIR and TEPs with incremental number of removed components ended ')
+    end
+end
+  %%  5.2 SSP-SIR manual components selection
+for iSub = 1:n_subjects
+    subj =  subjects{iSub} 
+    for n_cond= 1:tot_cond
+        load([output '\' subj '\0_Data\Block_cond.mat']);
+        name_condition = [subj '_' Block_cond{n_cond,2}];
+        name_folder= [subj '\' name_condition]
+
+        previous_folder= [output, '\' name_folder '\5_SSP-SIR\5_1_SSP-SIR calculations'];
+        mkdir([output, '\' name_folder '\5_SSP-SIR\5_2_SSP-SIR_components']);
+        current_output_folder = [output,'\' name_folder '\5_SSP-SIR\5_2_SSP-SIR_components'];
+
+%         if  size(dir([current_output_folder '/*.mat' ]),1)==1 && global_overwrite==0 && SSP_SIR_step2==0
+%                      warning(['skipping step 5.2 for ' subj ' , files already present'])
+%         else
+            decide1_SSPSIR_comp_MOBI
+%         end
+disp('Step 5.2 SSP-SIR removed components ended ')
+    end
+end
+%% 6. Last step of processing after SSP-SIR
+for iSub = 1:n_subjects
+    subj =  subjects{iSub};
+    for n_cond= 1: tot_cond
+        load([output '\' subj '\0_Data\Block_cond.mat']);
+        name_condition = [subj '_' Block_cond{n_cond,2}];
+        name_folder= [subj '\' name_condition]
+        savename= [Block_cond{n_cond,1}];
+
+        previous_folder= [output,'\' name_folder '\5_SSP-SIR\5_2_SSP-SIR_components'];
+        mkdir([output, '\' name_folder '\6_final_preprocessing_new_interp']);
+        current_output_folder = [output,'\' name_folder '\6_final_preprocessing_new_interp'];
+        
+         if  isfile([current_output_folder '\' savename '_TEPs_def.set'])&& ...
+                isfile([current_output_folder '\' savename  '_TEPs_def.mat'])
+                     warning(['skipping step 6 for ' name_condition ' , files already present'])
+        else        
+            % to know the number of removed components in previuos step
+            load([previous_folder '\SSPSIR_delcomponents.mat'])
+            k = SSPSIR_delcomponents(1,1); 
+
+            % load the set after SSP-SIR with correct number of removed variables
+            EEG = pop_loadset('filename',[name_condition '_5SSP-SIR_' num2str(k) '.set'],'filepath', [output,'\' name_folder '\5_SSP-SIR\5_1_SSP-SIR calculations']);
+            EEG.comments=pop_comments(EEG.comments,'','SSP-SIR components evaluated by 2 scientists',1);
+
+            %EEG = pop_tesa_interpdata( EEG, 'cubic',interp_Win);
+
+            EEG= last_steps_MOBI_new(EEG,savename,LP_f,old_ref,chan_loc_st,rmBase,current_output_folder,epoching_short,reref);
+            %myplot_F(EEG,[-100 400],[-15 15],['Final']);
+
+            %artefact rejection 2 done in folder 6_final_preprocessing
+            load([output,'/' name_folder '/6_final_preprocessing/artrej2_rejepochs.mat'])
+            trialrej=zeros([1 size(EEG.data,3)]);
+            trialrej(artrej2_rejepochs)=1;
+            EEG = pop_rejepoch(EEG,trialrej,0);
+            
+            %remove data between 0 and 3
+            EEG = pop_tesa_removedata( EEG,interp);
+            myplot_F(EEG,[-100 400],[-20 40],['final remove']);
+            saveas(gcf,[current_output_folder '/final_remove'])
+            close
+
+            %save
+            EEG.setname= 'TEPs_def';
+            EEG = pop_saveset( EEG, 'filename',[savename '_TEPs_def.set'],'filepath',current_output_folder);
+
+            
+            % convert to fieldtrip structure and save figure TEPs def
+            TEP_final2=convert_to_fieldtrip(EEG,latency,current_output_folder,savename,ylimit);
+        
+            end
+disp('Step 6 ended ')
+    end
+    
+
+end