function [hrt, hr] = myecg2hr2(ecgt, ecg, fs, t1, t2)

% algorithm:
%   1. simple bandpass filter ECG
%   2. select pulse
%   3. match filter
%   4. take absolute value of match filter output
%   5. generate initial threshold curve, detect beats
%   6. determine bpm
%   7. generate comb filter bank
%   8. weight each time point according to estimated bpm
%   9. detect new beats
%   10. possibly low pass filter this function
%   11. go back to 6, using tighter controls

% optional method to detect beats:
%   take difference between threshold function and ecg
%   lowpass filter difference

%   find local maxima
lowpass = 60;  % hertz
highpass = 1;   % hertz

hrlowpass = 1; % hertz

freqperiods = 4;        % number of periods to avrage over for frequency

max_bpm  = 150;
min_bpm = 30;

% comb filter bank parameters
comb_coeff = 1.4;       % amount to scale the comb bank filters
comb_mean = .85;        % deptch comb cuts
center_strength = 1; % amount to overemphasize centerpoint
bank_size = 40;         % number of filter banks
filter_length = 500;  % length of comb filter
comb_cycles = 1;
filter_window = 2;     % size of pulse to detect
min_period = .066;     % half minimum period of heart rate (s)
follow_peak_scale = .5;
min_period = floor(min_period * fs);

initial_filter_length = 500;
initial_filter_mean = 0.95;
initial_filter_coeff = 5/(comb_coeff*initial_filter_mean);

disp_filter_bank = 0;

comb_offset = zeros(1,length(ecg));

if(nargin < 3) 
    error "Error, not enough arguments\n";
elseif(nargin < 4) 
    plot(ecgt,ecg);
    t1 = input('Enter early time of good beat ==> ');
    t2 = input('Enter late time of good beat ==> ');
elseif(nargin < 5)
    t1 = [];
    t2 = [];
end

% initial bandpass filtering
hpf = fir1(5000,highpass/fs,'high');
lpf = fir1(5000,lowpass/fs);
ecgfilt = real(ifft(fft(ecg,length(ecg)+length(lpf)-1) .* fft(hpf,length(ecg)+length(lpf)-1) .* fft(lpf,length(ecg)+length(lpf)-1)));
ecgfilt = ecgfilt((length(lpf)):end);

if(~isempty(t1) & ~isempty(t2)) 
    times = find((ecgt <= t2) & (ecgt > t1));
    disp('Match filtering');
    [ecgfilt, ecgt] = centermatch(ecgfilt, ecgfilt(times), ecgt);
end

% take abs for improved combing
ecgfilt = (ecgfilt) .* (ecgfilt > 0);
% ecgfilt = (ecgfilt .* (ecgfilt > 0)).^2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Get initial HR estimate

myinput = 'a';
initialthreshold = centermatch(ecgfilt, makecomb(0.0001, initial_filter_length, filter_window, fs, initial_filter_mean), ecgt);
filter_coeff = initial_filter_coeff;
for n = 1:2
    threshold = comb_coeff.*filter_coeff.*initialthreshold+comb_offset;
    threshold = threshold .* (threshold >= 0);
    thresheld = ecgfilt-threshold;
    beats = maxcenter(thresheld, min_period);
    rr = [diff(ecgt(beats))]; % resample the heart rate
    [hrt, hr] = hrtach(ecgt(beats(1)), rr*1000, ones(size(rr)), fs, ecgt);
    averagerr = round(mean(rr)*fs);
    filtlen = round(averagerr*freqperiods);
    [hr(1).*ones(1,filtlen) hr' hr(end)*ones(1,filtlen)];
    %hr = centermatch([hr(1).*ones(1,filtlen) hr' hr(end)*ones(1,filtlen)],ones(1,filtlen));
    %hr = hr((1+filtlen):(end-filtlen));
    % Display
    disp_ecg2hr(ecgt,ecg,ecgfilt,rr,hr,threshold,thresheld,beats);
    
    noisethreshold = centermatch(ecgfilt, makecomb(0.0001, initial_filter_length, 0, fs, 0.999), ecgt);
    mean_peak = spline(ecgt(beats), ecgfilt(beats), ecgt)';
    mean_peak = mean(mean_peak) + (mean_peak - mean(mean_peak))*follow_peak_scale;
    filter_coeff = mean_peak./noisethreshold;
    quality = mean(thresheld(beats))
    std = mean((thresheld(beats)-quality).^2)
    a = input('Would you like to continue? [Y/n] ', 's');
    if(~isempty(a) & (a == 'n' | a == 'N'))
        [hrt, hr] = hrtach(ecgt(beats(1)), rr*1000, ones(size(rr)), 64);
        return
    end
end
    %filter_coeff = mean_peak/mean(noisethreshold);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Perform iterative convergence

%low_bpm = max(min(hr), min_bpm);
%high_bpm = min(max(hr), max_bpm);

% calculate the filter banks
%bpm = linspace(low_bpm, high_bpm, bank_size);
%for i = 1:bank_size;
%    filter_bank(i,:) = makecomb(bpm(i), filter_length, filter_window, fs, comb_mean, center_strength);
%    filter_ecg(i,:) = centermatch(ecgfilt, filter_bank(i,:),ecgt);
    %ampl_ecg(1,:) = centermatch(ecgfilt, ones(size(filter_bank(i,:))), ecgt);
    %end

%if(disp_filter_bank == 1)
%    for i = 1:bank_size;
%        subplot(211);
%        plot(0:1/fs:((length(filter_bank)-1)/fs),filter_bank(i,:));
%        title(bpm(i));
%        subplot(212);
%        plot(ecgt,comb_coeff.*filter_ecg(i,:),ecgt,ecgfilt);
%        title('Filtered ecg');
%        pause
%    end
%end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Repeat

% comb_coeff = comb_coeff * 4;
for(n = 1:30)
    if(mod(n,5) == 1)
        low_bpm = max([min(hr), min_bpm, mean(hr)-2*sqrt(mean((hr-mean(hr)).^2))])
        high_bpm = min([max(hr), max_bpm, mean(hr)+2*sqrt(mean((hr-mean(hr)).^2))])
        % calculate the filter banks
        bpm = linspace(low_bpm, high_bpm, bank_size);
        for i = 1:bank_size;
            filter_bank(i,:) = makecomb(bpm(i), filter_length, filter_window, fs, comb_mean, center_strength);
            filter_ecg(i,:) = centermatch(ecgfilt, filter_bank(i,:),ecgt);
            %ampl_ecg(1,:) = centermatch(ecgfilt, ones(size(filter_bank(i,:))), ecgt);
        end

        center_strength_decay = .8;
        center_strength = center_strength_decay * (1 + center_strength) - 1;
        
        if(disp_filter_bank == 1)
            for i = 1:bank_size;
                subplot(211);
                plot(0:1/fs:((length(filter_bank)-1)/fs),filter_bank(i,:));
                title(bpm(i));
                subplot(212);
                plot(ecgt,comb_coeff.*filter_ecg(i,:),ecgt,ecgfilt);
                title('Filtered ecg');
                pause
            end
        end
    end
       
    hr = repmat(hr', bank_size, 1);
    l = length(filter_ecg(1,:));
    bpm2 = repmat(bpm', 1, l);
    [diverge, bank_sel] = min(abs(bpm2-hr));

    % weight the filter banks
    bank_sel_coef = zeros(bank_size,l);
    for(i = 1:l)
        bank_sel_coef(bank_sel(i),i) = 1;
        %[trash,index] = min(filter_ecg(:,i));
        %bank_sel_coef(index,i) = 1;
    end

    threshold = comb_coeff.*filter_coeff.*sum(filter_ecg.*bank_sel_coef)/comb_mean+comb_offset;
    thresheld = ecgfilt-threshold;
    beats = maxcenter(thresheld, min_period);   % pick out the beats
    rr = [diff(ecgt(beats))];                   % resample the heart rate
    [hrt, hr] = hrtach(ecgt(beats(1)), rr*1000, ones(size(rr)), fs, ecgt);

    % Display
    disp_ecg2hr(ecgt,ecg,ecgfilt,rr,hr,threshold,thresheld,beats);

    quality = mean(thresheld(beats))
    std = mean((thresheld(beats)-quality).^2)
    a = input('Would you like to continue? [Y/n] ', 's');
    if(~isempty(a) & (a == 'n' | a == 'N')) 
        [hrt, hr] = hrtach(ecgt(beats(1)), rr*1000, ones(size(rr)), 64);
        return
    end
end