function hbMark=dopRawHBmarks(leftCh,varargin)
%% dopRawHBmarks
%
% hbMark=dopRawHBmarks(leftCh,[sampleRate],[displayPlot]);
%
% find peaks of heart beat and return a channel with markers at these
% peaks.
% this can later be used for heart cycle correction
%
% accepts an optional sample rate so as to adjust the range to explore the
% heart cycles within - e.g., higher sampling rate requires greater
% distance between markers.
%
% notes: wondering about looking for heart beats in both channels and
% taking an average to be more accurate.  Not entirely happy with the way
% this operates on atypical heart cycles
%

% -------------------------------------------------------------------------
%     'dopOSCCI' summarises functional transcranial Doppler ultrasonography
%     (fTCD) data.  The primary function of the software is to assess the
%     hemispheric lateralization of cognitive function.
%
%     Copyright (C) 2011 the Chancellor, Masters and Scholars of the
%     University of Oxford.
% 
%     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 3 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, see <http://www.gnu.org/licenses/> 
%     or write to the Free Software Foundation, Inc., 51 Franklin Street,
%     Fifth Floor, Boston, MA 02110-1301, USA.
%
%     Authors: Nic Badcock, Georgina Holt, Anneka Holden, & Dorothy Bishop.
%       contact: nicholas.badcock@gmail.com
%
%     Type 'dopOSCCIlicense' for license details or see dopOSCCIlicense.txt
%
try
    sampleRate=25; % default is 25 Hz
    plotRange=[];
    if nargin>1
        sampleRate=varargin{1};
    end
    if nargin>2
        plotRange=varargin{2};
    end
    
    %% >> find peaks
    % peaks=1; % treat first point as peak, even though it might not be
    peaks=[];
    windowSize=3; % formerly 10, July 2010...some people's heart beats are weird...
    % filtCh=leftCh;
    filtCh=filter(ones(1,windowSize)/windowSize,1,leftCh);
    
    r=round(sampleRate/3.7879); % search range/windows
    
    r2=r*2;
    for ii=2:r2+1
        if filtCh(ii-1)>=filtCh(ii-ii+1:ii)
            if filtCh(ii-1)>=filtCh(ii:ii+r)
                peaks=[peaks ii-1]; %peaks holds times of peaks
            end
        end
    end
    for i=1+r:length(leftCh)-1-r
        if isempty(peaks)
            if and(filtCh(i)>=filtCh(i-r:i-1),filtCh(i)>=filtCh(i+1:i+r))
                peaks=[peaks i]; %peaks holds times of peaks
            end
        else
            if and(and(filtCh(i)>=filtCh(i-r:i-1),filtCh(i)>=filtCh(i+1:i+r)),(i-max(peaks))>r) %find a peak
                peaks=[peaks i]; %peaks holds times of peaks
            end
        end
    end
    peaks=[peaks length(leftCh)]; % last one needs to be one as well
    hbMark=zeros(length(leftCh),1);
    hbMark(peaks)=1;
    
    if ~isempty(plotRange)
        h=figure; hold;range=plotRange(1):plotRange(2); plot(leftCh(range),'g'); ...
            plot(filtCh(range));  plot(hbMark(range)*500,'r');
        pause; close(h);
    end
catch err
    %% catch dopOSCCI error
    save(dopOSCCIdebug(mfilename));rethrow(err); 
end