% use adaptive threshold (delta + lambda * local maximum) to find peaks in
% the detection function
% df is detection function;
% dfTime is corresponding time vector;
% delta, lambda, M is in the adaptive threshold model above, M is the
% adjacent frame length
%
% Author: Yuchao Song
% 2008/09/08

function [onsettime, ot_samp, peakval, rectdiff] = adpthresholding(dfTime, df, delta, lambda, M)
onsettime = [];
ot_samp = [];   % onset time in samples
peakval = [];
difference = zeros(1, length(df));
rectdiff = zeros(1, length(df));

dt = dfTime(2) - dfTime(1);
for n = 1:length(df)
    % tempthres = delta + lambda * median(abs(df(  max(1, n-M):min(length(df), n+M)  )));
    tempthres = delta + lambda * max(abs(df(  max(1, n-M):min(length(df), n+M)  )));
    difference(n) = df(n) - tempthres;
    if(difference(n) > 0)
        rectdiff(n) = difference(n);
    end
end

for n = 2:(length(difference)-1)
    if (difference(n)>difference(n-1)) && (difference(n)>difference(n+1)) && (difference(n)>0)
        % with parabolic interpolation
        %         p = 0.5 * (difference(n-1) - difference(n+1))/(difference(n-1)-2*difference(n)+difference(n+1));
        %         onsettime = [onsettime, dfTime(n) + p * dt]; % offset the group delay caused by the filter
        %         peakval = [peakval, (df(n) - 0.25*(df(n-1)-df(n+1))*p)];

        % without parabolic interpolation
        onsettime = [onsettime, dfTime(n)]; % offset the group delay caused by the filter
        ot_samp = [ot_samp, n];
        peakval = [peakval, df(n)];
    end
end