clear all;
close all;

% representative: subject 1, session 2, middle finger, second trial 

subject={'01','02','03','04','05','06','07','08','09','10','11','12','13','14','15','16','17','18','19','20'};

path_data=['D:\jiangxinyu\Open Access HDsEMG DataSet Analysis\hyser_dataset']; % change path_data to the location of the dataset
path_results='..\results\decomposition results\1DoF'; % change path_results to the location to save decomposition results

window_len=0.02; % 0.02 s sliding window to extract features
step_len=0.02; % 0.02 s step length of the sliding window to extract features
f_cutoff=10; % filter force data using a low pass filter with 10 Hz cut off frequency
fs_emg=2048;
fs_force=100;
R = 4; % Extension parameter
M = 300; % FastICA iteration number

threshold = 0.8; % threshold to select motor units

for i=1:20
    for session=1:2        
        
        force_finger_1dof=load_1dof(path_data,subject{1,i},session,'force');
        emg_finger_1dof=load_1dof(path_data,subject{1,i},session,'raw');
        
        mvc=get_mvc(path_data,subject{1,i},num2str(session));
        force_norm=normalize_force(force_finger_1dof,mvc);
        force_norm_preprocess=preprocess_force(force_norm,window_len,step_len,f_cutoff,fs_force,fs_emg);
        
        for u=1:5
            for v=1:3        
                
                SpikeTrainGoodSelect=[];
                SILSelect=[];
                for muscle_idx=1:2
                    i
                    session
                    u
                    v
                    muscle_idx
                    emg=emg_finger_1dof{u,v}(:,(muscle_idx-1)*128+1:muscle_idx*128);
                    
                    load([path_results,'/B_subject',subject{1,i},'_session',num2str(session),'_task',num2str(u),'_trial',num2str(v),'_R',num2str(R),'_M',num2str(M),'_muscle',num2str(muscle_idx),'.mat']);
                    load([path_results,'/SpikeTrain_subject',subject{1,i},'_session',num2str(session),'_task',num2str(u),'_trial',num2str(v),'_R',num2str(R),'_M',num2str(M),'_muscle',num2str(muscle_idx),'.mat']);
                    load([path_results,'/GoodIndex_subject',subject{1,i},'_session',num2str(session),'_task',num2str(u),'_trial',num2str(v),'_R',num2str(R),'_M',num2str(M),'_muscle',num2str(muscle_idx),'.mat']);
                    load([path_results,'/SpikeTrainGood_subject',subject{1,i},'_session',num2str(session),'_task',num2str(u),'_trial',num2str(v),'_R',num2str(R),'_M',num2str(M),'_muscle',num2str(muscle_idx),'.mat']);
                    load([path_results,'/sGood_subject',subject{1,i},'_session',num2str(session),'_task',num2str(u),'_trial',num2str(v),'_R',num2str(R),'_M',num2str(M),'_muscle',num2str(muscle_idx),'.mat']);
                    load([path_results,'/SIL_subject',subject{1,i},'_session',num2str(session),'_task',num2str(u),'_trial',num2str(v),'_R',num2str(R),'_M',num2str(M),'_muscle',num2str(muscle_idx),'.mat']);
                    
                    idx=find((SIL>threshold) & (SIL<0.99));
                    SpikeTrainGoodSelect=[SpikeTrainGoodSelect,SpikeTrainGood(:,idx)];
                    SILSelect=[SILSelect,SIL(idx)];
                    if(muscle_idx==1)
                        NumMU_muscle1(u,v,i,session)=size(SpikeTrainGood(:,idx),2);
                    else
                        NumMU_muscle2(u,v,i,session)=size(SpikeTrainGood(:,idx),2);
                    end
                end
                NumMU(u,v,i,session)=size(SpikeTrainGoodSelect,2);
                SILMean(u,v,i,session)=mean(SILSelect);
                if(size(SpikeTrainGoodSelect,2)~=0)
                    [FiringRate,TimeWin] = FiringRateEst(SpikeTrainGoodSelect,[1:NumMU(u,v,i,session)]',window_len,step_len,fs_emg,'On');
                    [b,a]= butter(8,10/(1/(TimeWin(2)-TimeWin(1))/2),'low'); 
                    FiringRate = filtfilt(b,a,double(FiringRate));
                    parameter = regress(force_norm_preprocess{u,v}(51:end-51,u),FiringRate(51:end-51,:));% remove the first and last 1 s signal
                    force_regress=FiringRate*parameter;
                    corr_mat=corrcoef([force_regress(51:end-51),force_norm_preprocess{u,v}(51:end-51,u)]);
                    correlation(u,v,i,session)=corr_mat(1,2);
                else
                    correlation(u,v,i,session)=NaN;
                end
            end
        end
    end
end
SIL_subject_mean=squeeze(mean(mean(mean(SILMean,1,'omitnan'),2,'omitnan'),4,'omitnan'));
correlation_subject_mean=squeeze(mean(mean(mean(correlation,1,'omitnan'),2,'omitnan'),4,'omitnan'));
NumMU_muscle1_subject_mean=squeeze(mean(mean(mean(NumMU_muscle1,1,'omitnan'),2,'omitnan'),4,'omitnan'));
NumMU_muscle2_subject_mean=squeeze(mean(mean(mean(NumMU_muscle2,1,'omitnan'),2,'omitnan'),4,'omitnan'));