function out = run_sqrs(ecg,HRVparams,rs)
%out = run_sqrs(ecg,s,rs)
%	OVERVIEW:
%       QRS onset detector
%   INPUT:
%       ecg = single channel of ECG in digital values
%       s   = settings struct
%       rs  = resampling option; 1 = resample original signal, 0 = don't
%   OUTPUT:
%       out = sample points of onset of each QRS complex
%   DEPENDENCIES & LIBRARIES:
%   REFERENCE: 
%	REPO:       
%       https://github.com/cliffordlab/hrv_toolbox
%   ORIGINAL SOURCE AND AUTHORS:     
%       Converted from sqrs.c at:
%       http://www.physionet.org/physiotools/wfdb/app/sqrs.c
%       to  sqrs.m (Matlab)     by J. Perry,    July        2006
%           sqrs.c	(C)         by	G.B. Moody  27 October  1990
%
%	Last revised by John:  25 February 2006
%	Last revised by Gari:  07 February 2007
%   
%   03-06-2017
%   Edited by Adriana Vest
%   Now requires settings struct HRVparams, which defines debug mode and 
%   sampling frequency, and resampling option rs.
%   LICENSE AND COPYRIGHT:    
%       See Below
% 
% if debug>0, then plot final data (default; debug=0)
% 
% This algorithm is sensitive to high frequency noise, so you 
% might want to use an FIR band-pass filter first.
% 
% Avaialble under the GPL (see below)
%-------------------------------------------------------------------------------
% sqrs: Single-channel QRS detector
% Copyright (C) 1990-2006 George B. Moody
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License as
% published by the Free Software Foundation; either version 2 of the
% License, or (at your option) any later version.
%
% This program 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 this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
% 02111-1307, USA.
%
% You may contact the author by e-mail (george@mit.edu) or postal
% mail (MIT Room E25-505A, Cambridge, MA 02139 USA).  For updates to
% this software, please visit PhysioNet (http://www.physionet.org/).
% _________________________________________________________________
%
% The detector algorithm is based on example 10 in the WFDB
% Programmer's Guide, which in turn is based on a Pascal program
% written by W.A.H. Engelse and C. Zeelenberg, "A single scan
% algorithm for QRS-detection and feature extraction", Computers in
% Cardiology 6:37-42 (1979).  `sqrs' does not include the feature
% extraction capability of the Pascal program.  The output of `sqrs'
% is an annotation file (with annotator name `qrs') in which all
% detected beats are labelled normal;  the annotation file may also
% contain `artifact' annotations at locations which `sqrs' believes
% are noise-corrupted.
%
% 'sqrs' has been optimized for adult human ECGs. For other ECGs,
% it may be necessary to experiment with the input sampling
% frequency and the time constants indicated below. 
%
% This program is provided as an example only, and is not intended
% for any clinical application.  At the time the algorithm was
% originally published, its performance was typical of
% state-of-the-art QRS detectors.  Recent designs, particularly
% those that can analyze two or more input signals, may exhibit
% significantly better performance.
%%

if nargin<2
    error('Incorrect number of input arguments: must provide ECG signal and HRV paramters struct')
end
if nargin < 3
    rs = 1;
end

%debug = s.debug;
debug = 0;
fs = HRVparams.Fs;

if ~rs
    %% 1. Use Fs of Input
    time = 0;
    now = 10;

    %freq = 256;
    freq = fs;
    ms160 = ceil(0.16*freq);
    ms200 = ceil(0.2*freq);
    s2 = ceil(2*freq);
    scmin = 500;
    % number of ADC units corresponding to scmin microvolts
    % scmin = muvadu(signal, scmin); 
    scmax = 10 * scmin;
    slopecrit = 10 * scmin;
    maxslope = 0;
    nslope = 0;
    out = [];

    while (now < length(ecg)) %   && (to == 0)
        filter = [1 4 6 4 1 -1 -4 -6 -4 -1] * ecg((now-9):now);
        if (mod(time, s2) == 0)
	    % Adjust slope 
            if (nslope == 0)
	      slopecrit = max(slopecrit - slopecrit/16, scmin);
            elseif (nslope >= 5)
	      slopecrit = min(slopecrit + slopecrit/16, scmax);
            end
        end
        if (nslope == 0  && abs(filter) > slopecrit)
            nslope = nslope + 1; 
	    maxtime = ms160;
	    if (filter > 0) 
	      sign = 1;
	    else
	      sign = -1;
	    end
            qtime = time;
        end
        if (nslope ~= 0)
            if (filter * sign < -slopecrit)
                sign = -sign;
                nslope = nslope + 1;
                if (nslope > 4) 
                  maxtime = ms200;
                else
                  maxtime = ms160;
                end
            elseif (filter * sign > slopecrit &&  abs(filter) > maxslope)
                maxslope = abs(filter);
            end
            
            if (maxtime < 0)
                if (2 <= nslope && nslope <= 4)
                    slopecrit = slopecrit + ((maxslope/4) - slopecrit)/8;
                    if (slopecrit < scmin)
                        slopecrit = scmin;
                    elseif (slopecrit > scmax) 
                        slopecrit = scmax;
                    end
                    out = [out; now - (time - qtime) - 4];
                    %annot.anntyp = NORMAL; 
                    time = 0;
                elseif (nslope >= 5)
                    out = [out; now - (time - qtime) - 4];
                    %annot.anntyp = ARFCT; 
                end
                nslope = 0;
            end
            maxtime = maxtime - 1;
        end
	time = time + 1;
	now = now + 1;
    end

    out=out-1; % adjust for 1 sample offset problem.

    if debug > 0
       plot(ecg,'b');
       hold on;
       plot(out,ecg(out),'m*'); 
    end
elseif rs
    %% 2. Resample at 256 Hz
    % These time constants may need adjustment for pediatric or
    % small mammal ECGs.
    ecg_rs = resample(ecg,256,fs);
    
    time = 0;
    now = 10;

    freq = 256;
    ms160 = ceil(0.16*freq);
    ms200 = ceil(0.2*freq);
    s2 = ceil(2*freq);
    %scmin = 500;
    % ANV changed this to work with real units
    scmin = 500;
    % number of ADC units corresponding to scmin microvolts
    % scmin = muvadu(ecg_rs, scmin); 
    scmax = 10 * scmin;
    slopecrit = 10 * scmin;
    maxslope = 0;
    nslope = 0;
    out = [];

    while (now < length(ecg_rs)) %   && (to == 0)
        filter = [1 4 6 4 1 -1 -4 -6 -4 -1] * ecg_rs((now-9):now);
        if (mod(time, s2) == 0)
            % Adjust slope 
            if (nslope == 0)
                slopecrit = max(slopecrit - slopecrit/16, scmin);
            elseif (nslope >= 5)
                slopecrit = min(slopecrit + slopecrit/16, scmax);
            end
        end
        if (nslope == 0  && abs(filter) > slopecrit)
            nslope = nslope + 1; 
            maxtime = ms160;
            if (filter > 0) 
                sign = 1;
            else
                sign = -1;
            end
            qtime = time;
        end
        if (nslope ~= 0)
            if (filter * sign < -slopecrit)
                sign = -sign;
                nslope = nslope + 1;
                if (nslope > 4) 
                    maxtime = ms200;
                else
                  maxtime = ms160;
                end
            elseif (filter * sign > slopecrit &&  abs(filter) > maxslope)
                maxslope = abs(filter);
            end
            if (maxtime < 0)
                if (2 <= nslope && nslope <= 4)
                    slopecrit = slopecrit + ((maxslope/4) - slopecrit)/8;
                    if (slopecrit < scmin)
                        slopecrit = scmin;
                    elseif (slopecrit > scmax) 
                        slopecrit = scmax;
                    end
                    out = [out; now - (time - qtime) - 4];
                    %annot.anntyp = NORMAL; 
                    time = 0;
                elseif (nslope >= 5)
                    out = [out; now - (time - qtime) - 4];
                    %annot.anntyp = ARFCT; 
                end
                nslope = 0;
            end
            maxtime = maxtime - 1;
        end
	time = time + 1;
	now = now + 1;
    end

    out=out-1; % adjust for 1 sample offset problem.

    if debug > 0
       plot(ecg_rs,'b');
       hold on;
       plot(out,ecg_rs(out),'m*'); 
    end

end