function [iMatch, iPeakWindows, npks, R, stdWin, stdMat, rmsMat, signal_interp] = peakdelay(signal, minheightSTD, stdWindow, peakWindow, maxDelay, minDelay, sFreqInput, sFreqAnalysis)

% Outline peakdelay function:
% step 0: get size signal & excluded signals
% step 1: interpolate signal to get preferred sample frequency
% step 2: determine STD windows to determine minimum height for selecting peaks
% step 3: find peaks
% step 4: match selected peaks from neighbouring channels within the considered time/sample window
% step 5: compute correlation between matched peaks

[nR,~] = size(signal);
zeroCH = find(~any(signal,1));

%% step 1: interpolate signal

t_signal = double((1:nR)./sFreqInput);
t_interp = linspace(t_signal(1),t_signal(end),(t_signal(end)-t_signal(1))*sFreqAnalysis+1);
signal_interp = interp1(t_signal,signal,t_interp,'spline');

clearvars signal t_interp t_signal % minimize memory burden in case of large data sets 

disp('----------------------')
disp('Interpolation finished')
disp('----------------------')

%% step 2: determine std over windows equal to input value

[nR,nC] = size(signal_interp);
stdWindow = stdWindow*sFreqAnalysis;
NstdWindows = ceil(nR/stdWindow);

stdWin=zeros(NstdWindows,2);stdMat=zeros(NstdWindows,nC);rmsMat=zeros(NstdWindows,nC);
for i = 1:NstdWindows
    if i~=NstdWindows 
        stdWin(i,:) = [i*stdWindow-stdWindow+1,i*stdWindow];
    else
        stdWin(i,:) = [i*stdWindow-stdWindow+1,nR];
    end
        stdMat(i,:) = std(signal_interp(stdWin(i,1):stdWin(i,2),:),[],1);
        rmsMat(i,:) = rms(signal_interp(stdWin(i,1):stdWin(i,2),:),1);
end

tempMinHeight = minheightSTD*min(stdMat,[],1);

stdMat = stdMat.*minheightSTD;

disp('----------------------')
disp('Windowed S.D. finished')
disp('----------------------')

%% step 3: find peaks

loc=cell(2,nC);pk=cell(2,nC);npks=zeros(1,nC);
for c = 1:nC
    if ~ismember(c,zeroCH)
        [pk{1,c},loc{1,c}] = findpeaks(signal_interp(:,c),'MinPeakDistance',sFreqAnalysis/128,'MinPeakHeight',tempMinHeight(1,c));
        [pk{2,c},loc{2,c}] = findpeaks(-signal_interp(:,c),'MinPeakDistance',sFreqAnalysis/128,'MinPeakHeight',tempMinHeight(1,c));

        %remove peaks below windowed std
        for i = 1:size(stdWin,1)
            indpos = logical((loc{1,c} >= stdWin(i,1)) .* (loc{1,c} <= stdWin(i,2)) .* (pk{1,c} < stdMat(i,c)));
            indneg = logical((loc{2,c} >= stdWin(i,1)) .* (loc{2,c} <= stdWin(i,2)) .* (pk{2,c} < stdMat(i,c)));
            loc{1,c}(indpos) = []; pk{1,c}(indpos) = [];
            loc{2,c}(indneg) = []; pk{2,c}(indneg) = [];
        end
        npks(1,c) = length(loc{1,c})+length(loc{2,c});
    end
    disp(['channel analyzed: ',num2str(c)])
end

disp('----------------------')
disp('Finding peaks finished')
disp('----------------------')

%% step 4: peak matching (both directions of travel, ie delay)

iMatch = cell(1,nC-8);
for i=1:nC-8
    if ~ismember(i,unique([zeroCH zeroCH-8]))
        Nit=1000;

        %positive peaks matching
        iMatchpos=[];
        Itwin = 1:Nit:size(loc{1,i},1);Itwin=[Itwin,size(loc{1,i},1)];
        for j=1:length(Itwin)-1
%             disp(j)
            iMatchpos_temp=[];
            if j~=1
                olap=10;
            else 
                olap=0;
            end
            [iPr,iPc] = find( logical( (abs(loc{1,i}(Itwin(j)-olap:Itwin(j+1))-loc{1,i+8}') <= maxDelay) .* (abs(loc{1,i}(Itwin(j)-olap:Itwin(j+1))-loc{1,i+8}') >= minDelay) ));
            iPr_corr = iPr+Itwin(j)-olap-1;
            iMatchpos_temp(:,1) = single(loc{1,i}(iPr_corr));
            iMatchpos_temp(:,2) = single(loc{1,i+8}(iPc));
            iMatchpos_temp(:,3) = single(loc{1,i+8}(iPc)-loc{1,i}(iPr_corr));
            iMatchpos_temp(:,4) = single(pk{1,i}(iPr_corr));
            iMatchpos_temp(:,5) = single(pk{1,i+8}(iPc));
            iMatchpos = [iMatchpos;iMatchpos_temp];
        end
        iMatchpos = unique(iMatchpos, 'stable', 'rows');

        %negative peaks matching
        iMatchneg=[];
        Itwin = 1:Nit:size(loc{2,i},1);Itwin=[Itwin,size(loc{2,i},1)];
        for j=1:length(Itwin)-1
            iMatchneg_temp=[];
            if j~=1
                olap=10;
            else 
                olap=0;
            end
            [iPr,iPc] = find( logical( (abs(loc{2,i}(Itwin(j)-olap:Itwin(j+1))-loc{2,i+8}') <= maxDelay) .* (abs(loc{2,i}(Itwin(j)-olap:Itwin(j+1))-loc{2,i+8}') >= minDelay) ));
            iPr_corr = iPr+Itwin(j)-olap-1;
            iMatchneg_temp(:,1) = single(loc{2,i}(iPr_corr));
            iMatchneg_temp(:,2) = single(loc{2,i+8}(iPc));
            iMatchneg_temp(:,3) = single(loc{2,i+8}(iPc)-loc{2,i}(iPr_corr));
            iMatchneg_temp(:,4) = single(pk{2,i}(iPr_corr));
            iMatchneg_temp(:,5) = single(pk{2,i+8}(iPc));
            iMatchneg = [iMatchneg;iMatchneg_temp];
        end
        iMatchneg = unique(iMatchneg, 'stable', 'rows');
    else
        iMatchpos = zeros(1,5);iMatchneg = zeros(1,5);
    end
        
    % collect matching positive and negative peaks
    iMatchTemp = [iMatchpos(~isnan(iMatchpos(:,1)),1), iMatchpos(~isnan(iMatchpos(:,2)),2), iMatchpos(~isnan(iMatchpos(:,3)),3), iMatchpos(~isnan(iMatchpos(:,4)),4), iMatchpos(~isnan(iMatchpos(:,5)),5);...
                  iMatchneg(~isnan(iMatchneg(:,1)),1), iMatchneg(~isnan(iMatchneg(:,2)),2), iMatchneg(~isnan(iMatchneg(:,3)),3), iMatchneg(~isnan(iMatchneg(:,4)),4), iMatchneg(~isnan(iMatchneg(:,5)),5)];
    [iMatch{1,i}(:,1), isort] = sort(iMatchTemp(:,1),1);
    iMatch{1,i}(:,2) = iMatchTemp(isort,2); iMatch{1,i}(:,3) = iMatchTemp(isort,3);
    iMatch{1,i}(:,4) = iMatchTemp(isort,4); iMatch{1,i}(:,5) = iMatchTemp(isort,5);
    
    disp(['Peak matching CH',num2str(i),' done'])
end

disp('----------------------')
disp('Peak matching finished')
disp('----------------------')


%% step 5: compute R between matched peaks

peakWindow = ceil(peakWindow*sFreqAnalysis); %preferebly peakWindow = 2^N

% assuming a window equal to 2^N is used as input, the window is created as
% follows: N-1 - iPeak - N
% peak windows exceeding signal_interp dimensions are ignored.
iPeakWindows=cell(1,nC-8);
for i=1:nC-8
    if ~ismember(i,unique([zeroCH zeroCH-8]))
        idim = logical( (iMatch{1,i}(:,1)-(peakWindow/2-1) >= 1) .* (iMatch{1,i}(:,1)+(peakWindow/2) <= nR) .* (iMatch{1,i}(:,2)-(peakWindow/2-1) >= 1) .* (iMatch{1,i}(:,2)+(peakWindow/2) <= nR) );
        iPeakWindows{1,i} = [iMatch{1,i}(idim,1)-(peakWindow/2-1), iMatch{1,i}(idim,1)+(peakWindow/2), iMatch{1,i}(idim,2)-(peakWindow/2-1), iMatch{1,i}(idim,2)+(peakWindow/2)];
        iMatch{1,i} = iMatch{1,i}(idim,:);
    end
end

R=cell(1,nC-8);
%     nPeaks = cellfun(@size,iPeakWindows,'UniformOutput',false);
for ch=1:nC-8
    if ~ismember(ch,unique([zeroCH zeroCH-8]))
        for nP = 1:size(iPeakWindows{1,ch},1)
%             disp(num2str(nP))
            Rmat = corrcoef(signal_interp(iPeakWindows{1,ch}(nP,1):iPeakWindows{1,ch}(nP,2),ch), signal_interp(iPeakWindows{1,ch}(nP,3):iPeakWindows{1,ch}(nP,4),ch+8));
            R{1,ch}(nP,1) = Rmat(1,2);
        end
    end
    disp(['channel analyzed: ',num2str(ch)])
end
disp('----------------------')
disp('Peak2Peak R calculated')
disp('----------------------')

end
