function [hbCor isnan_count] = dopRawHBcor(inData,hbMark)
% dopOSCCI Mark II: dopRawHBcor~ 11-jan-12 (last edit)
%
% hbCor=dopRawHBcor(inData,hbMark)
%
% correction for heart beat - average activation between hbMarkers
%
%

% -------------------------------------------------------------------------
%     '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
    marksOnly=find(hbMark);
    lastMark=1;
    hbCor=inData;
    isnan_count = 0;
    % multiple event markers may use 3 dimensional matrix
    evChs=1;
    ds=size(inData); if length(ds)==3; evChs=ds(3); end
    for i = 1 : numel(marksOnly) % for each hb maker
        range=lastMark:marksOnly(i);
        for ii=1:evChs
            for j=1:2 % left right channels
                if sum(isnan(hbCor(range,j,ii)))
                    % replace nan with 0 = still rejected but can use the
                    % data
                    tmp = hbCor(range,j,ii); filt = isnan(hbCor(range,j,ii));
                    tmp(filt) = 0;
                    hbCor(range,j,ii) = tmp;
                    fprintf('*** %u NaN (not a number) samples replaced with zero during Heart Cycle correction (function dopRawHBcor)\n',sum(filt));
                    isnan_count = isnan_count + sum(filt);
                end
                
                if min(hbCor(range,j,ii))==0 % if dropout - reduce influence
                    hbCor(range,j,ii)=median(hbCor(range,j,ii));
                else
                    hbCor(range,j,ii)=mean(hbCor(range,j,ii));
                end
                lastMark=marksOnly(i);
            end
        end
    end
    clear inData hbMark marksOnly lastMark evChs ds range
catch err
    %% catch dopOSCCI error
    save(dopOSCCIdebug(mfilename));rethrow(err);
end