%Code for analysis of Study 1 dataset related to 
% Worley et al. The PiNe box: development and validation of 
% an electronic device to time-lock multimodal responses to sensory stimuli in hospitalised infants 
%
% The associated data set can also be downloaded as a supplementary file
% Please cite this paper if you use the code.

%% Identify the data - this is found as S1 Dataset
fpath='~/pine_test_study1/';
files=dir(strcat(fpath,'*_waves.csv'));

ip_epoch_store=[]; file_no_store=[]; %assign variables

%loop through files and store the IP signal
for i=1:length(files)
    
fid = fopen(strcat(fpath,files(i).name),'r');
H=textscan(fid,'%s%s%s%s%s%s%s%s%s%s',1,'Delimiter',',');
fclose(fid);
opts=detectImportOptions(strcat(fpath,files(i).name));
Y=readtable(strcat(fpath,files(i).name),opts);

t=Y{:,1}; ts=t(1); %first column is time variable
ind=find(contains([H{:}],'Resp')); %IP signal is given name Resp by monitor

ip=Y{:,ind};

t=t-ts; %time in ms, subtract the start time as data is exported in milliseconds since 1.1.1970
t=t/1000; %convert from ms to seconds
time_ip=t;
if iscell(ip)
    f=find(cellfun(@isempty,ip));
    ip(f)=[];
    time_ip(f)=[];
    ip=cellfun(@str2num,ip);
else
    f=find(isnan(ip));
    ip(f)=[];
    time_ip(f)=[];

    if isempty(ip)
        time_ip=[];
    end
end


%import annotations of events (i.e. pushbutton triggers from pine box)
fid = fopen(strcat(fpath,files(i).name),'r');
X=textscan(fid,'%s%s%s%s%s','Delimiter',',', 'Headerlines',1);
fclose(fid); 

an=X{5}; %annotations are in column 5
annIndex=[];
for ii=1:length(an)
    if ~isempty(an{ii})
        annIndex=[annIndex,ii];
    end
end
time_annotations=t(annIndex); %time annotations
annotations=an(annIndex);
%find push button events of interest
ind=find(strcmp(annotations,'[PushButton]'));

% remove acciental triggers in study 4
if i==4
    ind([1,3,4,5,6,8])=[];
end
time_trigger=time_annotations(ind);

%epoch events from 1 second before to after the push button events
pre=1; post=1;

fs=62.5; %sampling rate of IP signal

time_epoch=[-pre:1/fs:post];
ip_epoch=zeros(length(time_epoch),length(time_trigger));
for j=1:length(time_trigger)
    ind_epoch=find(time_ip>=time_trigger(j)-pre & time_ip<=time_trigger(j)+post);
    ind_epoch=[ind_epoch(1)-1;ind_epoch];
    if length(ind_epoch)>length(time_epoch)
        ind_epoch(end)=[];
    end
    ip_epoch(:,j)=ip(ind_epoch);

end

ip_epoch_store=[ip_epoch_store,ip_epoch]; %save epoched data
file_no_store=[file_no_store;i*ones(length(time_trigger),1)];

end

%% Produce supplementary Figure for paper and identify time at which artefact occurs in each trial

pre_ind1=find(time_epoch>-0.5,1,'first');
pre_ind2=find(time_epoch<0,1,'last');

figure; count=1;
set(gcf, 'Position',  [100, 100, 1000, 1200])
time_art_store=zeros(size(ip_epoch_store,2),1); %vector to store times at which artefact occurs

for k=1:size(ip_epoch_store,2) %loop through each trial
    
    ip_epoch_store(:,k)=ip_epoch_store(:,k)-mean(ip_epoch_store(pre_ind1:pre_ind2,k)); %baseline correct epoch y subtracting mean
    
    post_min=min(ip_epoch_store(pre_ind2+1:end,k));
    
    
    thres=0.03*post_min; %set threshold as 3% of the minimum signal value in the 1 second after the event mark
    
    if file_no_store(k)==2 || file_no_store(k)==4 % in 2 files the start time was manually adjusted to better fit the data, following visual assessment
        adjust=5;
    else
        adjust=1;
    end
  
    ind_art=find(ip_epoch_store(pre_ind2+adjust:end,k)<thres,1,'first'); %find first point below threshold
    time_art=time_epoch(pre_ind2+ind_art+adjust-1); %time of artefact is the time at this point of crossing threshold
    if ~isempty(time_art)
    time_art_store(k)=time_art;
    end

    %plot epoch
    subplot(7,5,count)
    plot(time_epoch,ip_epoch_store(:,k),'linewidth',1.5)
    box off
    hold on
    plot(time_art,ip_epoch_store(pre_ind2+ind_art+adjust-1,k),'.r','markersize',12)
    plot([0,0],ylim,'k--')
    xlim([-1,1])
    set(gca,'fontsize',11)
    ylim(ylim)
    
    if count==33
        xlabel('Time (s)','fontsize',15)
    end
    if count==16
        ylabel('IP','fontsize',15)
    end
    count=count+1;
    

end
 


%% calculate metrics for paper


min_diff=min(time_art_store*1000)
max_diff=max(time_art_store*1000)
mean_diff=mean(time_art_store*1000)