function onset_data = SMIDIBT_percussionpad_onsets(data_file,params)
% [onset_data] = SMIDIBT_percussionpad_onsets(data_file,params)
%   Extracts the send MIDI onset and offset times and audio onset times
%
%   INPUTS:
%   data_file       The name of a ".wav" file ending in "_trigs" that
%                   contains the MIDI onset and offset triggers. This file
%                   MUST have a complementary audio file ending in "_audio"
%                   to extract the audio onsets.
%
%   params          A file containing the parameters for onset extraction
%                   (e.g., amplitude and timing thresholds). The parameters
%                   used in Schultz (2017, Experiment 3) are set as default.
%
%   OUTPUT:
%   onset_data      The times (in milliseconds) of the FSR onsets (first 
%                   column), FSR offsets (second column), MIDI received
%                   onsets from the read MIDI (third column), and the audio
%                   onsets (fourth column).
%
%   This version has dependencies appended at the bottom
%
%   2017-09-28 ben.schultz@maastrichtuniversity.nl
%   Copyright (c) 2017, Benjamin Schultz, Maastricht University.

%   This script is described in more detail in the publication:
%   Schultz, B. G. (submitted). The Schultz MIDI Benchmarking toolbox for
%   MIDI interfaces, percussion pads, and sound cards. Behavior Research 
%   Methods.
%
%   The SMIDIB Toolbox is distributed in the hope that it will be useful, but
%   WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
%   General Public License for more details.
% 
%   You should have received a copy of the GNU General Public License
%   along with the SMIDIB Toolbox; if not, write to the Free Software
%   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
%   02110-1301 USA
% 
%   See the file "COPYING" for the text of the license.


%% Check inputs
% check input parameters and set defaults
params = check_params(params);

% check file names
fsr_ind = strfind(data_file,'_fsr');
audio_file = strrep(data_file,'_fsr.wav','_audio.wav');
midi_file = strrep(data_file,'_fsr.wav','_trigs.wav');

if isempty(fsr_ind)
    error('Filename does not have the correct format. Should end in ''_trigs''');
end

if ~exist(data_file,'file') || ~exist(audio_file,'file') || ~exist(midi_file,'file')
    error('One or more files do not exist in specified location');
end

% load fsr signal, midi triggers, and audio
[fsr_data,fsr_Fs]=audioread(data_file);
[trig_data,midi_Fs]=audioread(midi_file);
[aud_data,aud_Fs]=audioread(audio_file);

if fsr_Fs~=aud_Fs || fsr_Fs~=midi_Fs
     error('Sample rates do not match. Files may not be synchronous.');
end

% get parameter information
moving_max_wind_samps = round(params.moving_max_wind_ms*(aud_Fs/1000));
time_thresh_samps = round(params.time_thresh_ms*(fsr_Fs/1000));
off_time_thresh_samps = round(params.off_time_thresh_ms*(fsr_Fs/1000));
fsr_off_time_thresh_samps = round(params.fsr_off_time_thresh_ms*(fsr_Fs/1000));
trig_time_thresh_samps = round(params.trig_time_thresh_ms*(midi_Fs/1000));

% perform any signal tranformations
fsr_data = fsr_data/max(abs(fsr_data));
fsr_data = fsr_data*-1;

trig_data = abs(trig_data)/max(abs(trig_data));
trig_data = (trig_data-min(trig_data))/(max(trig_data)-min(trig_data));

aud_data = abs(aud_data);
aud_data = getMovingMax(aud_data,moving_max_wind_samps);
aud_data = (aud_data-min(aud_data))/(max(aud_data)-min(aud_data));

% get the onsets
prev_offset = 1;
cur_pos = 1;
cur_onset_temp = find(fsr_data(prev_offset:end)>params.fsr_on_thresh,1,'first')+prev_offset-1;
onset_data = ones(params.n_trigs,4)*NaN;

while ~isempty(cur_onset_temp)
    
    % estimate FSR onset    
    cur_onset = cur_onset_temp;
    
    % go back to onset
    while fsr_data(cur_onset-1)<fsr_data(cur_onset)
        cur_onset = cur_onset-1;
    end
    cur_onset = cur_onset+1;
    
    cur_offset = find(fsr_data(cur_onset_temp+fsr_off_time_thresh_samps:end)<params.fsr_off_thresh,1,'first')+cur_onset_temp+off_time_thresh_samps;
    
    % get MIDI onset
    cur_trig_onset_temp = find(trig_data(cur_onset_temp:cur_onset_temp+time_thresh_samps)>params.trig_on_thresh,1)+cur_onset_temp;
    
    if isempty(cur_trig_onset_temp)
        cur_trig_onset = NaN;
        cur_max = NaN;
    else
        
        cur_trig_onset = cur_trig_onset_temp;
        
        % go back to start of trigger
        while trig_data(cur_trig_onset-1)<trig_data(cur_trig_onset)
            cur_trig_onset = cur_trig_onset-1;
        end
        
        cur_trig_onset = cur_trig_onset+1;
        
        % get trigger offset
        cur_trig_offset = find(trig_data(cur_trig_onset_temp+trig_time_thresh_samps:end)<params.trig_off_thresh,1,'first')+cur_trig_onset_temp+trig_time_thresh_samps-1;
        
        % go back to start of offset
        cur_max = cur_trig_offset;
        while trig_data(cur_max-1)> trig_data(cur_max)
            cur_max = cur_max-1;
        end
        cur_max = cur_max+1;
        
    end
    
    % get audio onset
    cur_aud_onset_temp = find(aud_data(cur_onset_temp:cur_onset_temp+time_thresh_samps)>params.aud_on_thresh,1)+cur_onset_temp;
    
    if isempty(cur_aud_onset_temp)
        cur_aud_onset = NaN;
    else
        cur_aud_onset = cur_aud_onset_temp;
    end
    
    % set time threshold for next onset
    prev_offset = cur_offset+off_time_thresh_samps;
        
    % add onsets to the output
    onset_data(cur_pos,1:4) = [cur_onset,cur_aud_onset,cur_trig_onset,cur_max];
    
    % move to the next row
    cur_pos = cur_pos+1;
    
    %find next onset
    cur_onset_temp = find(fsr_data(prev_offset:end)>params.fsr_on_thresh,1,'first')+prev_offset-1;
    
end

% remove excess NaNs and turn into milliseconds
onset_data = onset_data(~isnan(onset_data(:,1)),:)/(fsr_Fs/1000);

%% Save onsets (recommended for large datasets)
if isfield(params,'out_dir')
    [~,data_filename,data_ext] = fileparts(data_file);
    out_filename = strrep([data_filename,data_ext],'_fsr.wav','_onsets.mat');
    save(fullfile(params.out_dir,out_filename),'onset_data','fsr_Fs');
end

%% Print plot
if isfield(params,'plot_dir')
    t = (1:length(trig_data))/fsr_Fs;
    cur_fig = figure;
    plot(t,fsr_data,'k'); 
    hold on;
    plot(t,trig_data,'b');
    plot(t,aud_data,'r');
    hold off;
    [~,data_filename,~] = fileparts(data_file);
    saveas(cur_fig,fullfile(params.plot_dir,strrep(data_filename,'_fsr','')),'fig');
    close(cur_fig);    
end

%% FUNCTIONS
function out = getMovingMax(data,wind)
% Gets moving maximum

out = ones(size(data))*NaN;

for i = 1:length(data)
    
    if i <= wind
        cur_data = data(1:i);
    else
        cur_data = data(i-wind:i);
    end
    
    out(i)=max(cur_data);
    
end

function params = check_params(params)

% check parameters and thresholds (set defaults)
if ~isfield(params,'time_thresh_ms');
    params.time_thresh_ms = 20; % time after trig to look for audio onset
end
if ~isfield(params,'off_time_thresh_ms');
    params.off_time_thresh_ms = 80; % time after trig to look for audio onset
end
if ~isfield(params,'fsr_on_thresh');
    params.fsr_on_thresh = 0.1; % threshold to detect MIDI onset
end
if ~isfield(params,'fsr_off_thresh');
    params.fsr_off_thresh = -0.1; % threshold to detect MIDI offset
end
if ~isfield(params,'trig_time_thresh_ms');
    params.trig_time_thresh_ms = 0.5; % time after trig to look for offset
end
if ~isfield(params,'trig_on_thresh');
    params.trig_on_thresh = 0.1; % threshold to detect MIDI onset
end
if ~isfield(params,'trig_off_thresh');
    params.trig_off_thresh = 0.1; % threshold to detect MIDI offset
end
if ~isfield(params,'aud_on_thresh');
    params.aud_on_thresh = 0.1; % threshold to detect audio onset
end
if ~isfield(params,'moving_max_wind_ms');
    params.moving_max_wind_ms = 1; % window for moving maximum
end
if ~isfield(params,'off_time_thresh_ms');
    params.off_time_thresh_ms = 80; % time to wait for next MIDI onset
end
if ~isfield(params,'n_trigs');
    params.n_trigs = 201; % expected number of onsets
end
if ~isfield(params,'fsr_off_time_thresh_ms');
    params.fsr_off_time_thresh_ms = 60; % time to wait for next FSR onset
end



