% PanTomkinsHamilton.m - Calls PanTompkins.dll - the compiled version of Hamilton's
% version placed on Internet designed for 200Hz.
%
%        This software is released under the terms of the GNU General
%        Public License (http://www.gnu.org/copyleft/gpl.html).

function [nbr_heart_beats, r_peak_sample_list] = PanTomkinsHamilton(input_ecg_sig, sampl_freq_ecg_sig)

global UseSignalToolbox

UseSignalToolbox = 1;

if (sampl_freq_ecg_sig ~= 200) % Do Re-Sampling using matlab FCT "resample()" !
    [p,q] = rat((200.0/sampl_freq_ecg_sig),0.0001);
    Resamp=which('upfirdn_rs');
    if(UseSignalToolbox == 1) %| length(Resamp)==0)   % No upfirdn_rs found
        input_ecg_sig_200hz = resample(input_ecg_sig, p, q); % , FilLen,WinBeta); !
    else
        ResampPath = Resamp(1:length(Resamp)-12);
        Coeff = sprintf('%sfircoef_%d_%d.mat',ResampPath, p, q);
        if(exist(Coeff))    % Filter coeff file exists.
            load(Coeff);
        else
            command = sprintf('cd %s',ResampPath);
            dos(command);
            [h,p_aux,q_aux] = ResampFiltDsgn('p','q');
            save(Coeff,'h');
        end;
    end;
    sampl_freq_ecg      = 200; % After ECG-Signal has been ReSampled to 200Hz !
    length_200hz_ecg    = length(input_ecg_sig_200hz);
    resampling_ratio       = sampl_freq_ecg_sig / 200.0; % Ideal target value !
    resampling_ratio_exact = q/p;   % Take care of rounding errors ! 12-04-05 !
    resampling_switch   = 1;
    length_input_ecg = length(input_ecg_sig);  % # Samples with original sampling !

    % -> input_ecg_sig      = input_ecg_sig_200hz;  % Full copy of input_ecg vec. !

else % sampl_freq_ecg_sig == 200 as expected !
    input_ecg_sig_200hz = input_ecg_sig';       % Local copy of input_ecg vec. !
    sampl_freq_ecg      = sampl_freq_ecg_sig;  % ECG-Signal sampled at 200 Hz !
    length_200hz_ecg    = length(input_ecg_sig);   % Idem to length_input_ecg !
    resampling_ratio       = 1.0;
    resampling_ratio_exact = 1.0;
    resampling_switch   = 0;
end % sampl_freq_ecg_sig <> 200 !

aux1 = PanTompkins(input_ecg_sig_200hz);
aux2 = find(aux1>0);
aux3 = (aux2 - aux1(aux2));
r_peak_sample_list = aux3(find(aux3>0));
nbr_heart_beats = length(r_peak_sample_list);

return;