function cMQRS = adjust_mqrs_location(ecg,MQRS,fs,debug)
% the point of this function is to adjust the MQRS location since
% MQRS are not at the exact same location from a lead to another
% due to conduction delay and/or quantification limitation.
% This can cause the template substraction technique to be
% imprecise (a small lack of alignement can cause poor results) if not
% performed. This function showed to improve the results of the PCinC2013.
%
% inputs
% ecg: the ecg signals
% MQRS: MQRS location (in ms)
% fs: sampling frequency
%
% outputs
% cMQRS: corrected MQRS location
% (if nargout is empty then plot the initial and corrected time series)
%
%
% FECG extraction toolbox, version 1.0, Sept 2013
% Released under the GNU General Public License
%
% Copyright (C) 2013 Joachim Behar
% Oxford university, Intelligent Patient Monitoring Group - Oxford 2013
% joachim.behar@eng.ox.ac.uk
%
% Last updated : 28-08-2013
%
% 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.
% == managing inputs
if nargin<2; error('adjust_mqrs_location: wrong number of input arguments \n'); end;
if nargin<3; fs=1000; end;
if nargin<4; debug=0; end;
if size(ecg,2)>size(ecg,1); ecg=ecg'; end;
% == general
cMQRS = zeros(length(MQRS),1); % corrected MQRS vector
NB_MQRS = length(MQRS);
WINDOW = ceil(fs*0.030); % allow a maximum of 30ms shift
NB_SAMPLES = length(ecg);
signVec = zeros(30,1);
% == core function
try
% == estimates if we are looking for local max or min on the first few
% peaks
for i=1:30
startLooking = MQRS(i)-WINDOW;
stopLooking = MQRS(i)+WINDOW;
if startLooking<=0; startLooking=1; end;
if stopLooking>=NB_SAMPLES; stopLooking=NB_SAMPLES; end;
[~,indm] = max(abs(ecg(startLooking:stopLooking)));
if startLooking==1; signVec(i) = MQRS(i)+indm-1;
else signVec(i) = MQRS(i)+indm-WINDOW; end;
end
SIGN = sign(median(ecg(signVec)));
% == local refinement to correct for sampling freq resolution
for i=1:NB_MQRS
if MQRS(i)>WINDOW && MQRS(i)+WINDOW<NB_SAMPLES
if SIGN>0
[~,indm] = max(ecg(MQRS(i)-WINDOW:MQRS(i)+WINDOW));
cMQRS(i) = MQRS(i)+indm-WINDOW-1;
else
[~,indm] = min(ecg(MQRS(i)-WINDOW:MQRS(i)+WINDOW));
cMQRS(i) = MQRS(i)+indm-WINDOW-1;
end
elseif MQRS(i)<WINDOW
% managing left boder
if SIGN>0
[~,indm] = max(ecg(1:MQRS(i)+WINDOW));
cMQRS(i) = MQRS(i)+indm-1;
else
[~,indm] = min(ecg(1:MQRS(i)+WINDOW));
cMQRS(i) = MQRS(i)+indm-1;
end
elseif MQRS(i)+WINDOW>NB_SAMPLES
% managing right border
if SIGN>0
[~,indm] = max(ecg(MQRS(i)-WINDOW:end));
cMQRS(i) = MQRS(i)+indm-WINDOW-1;
else
[~,indm] = min(ecg(MQRS(i)-WINDOW:end));
cMQRS(i) = MQRS(i)+indm-WINDOW-1;
end
else
cMQRS(i) = MQRS(i);
end
end
catch ME
for enb=1:length(ME.stack); disp(ME.stack(enb)); end;
cMQRS = MQRS;
end
% == plots
if debug || isempty(nargout)
plot(ecg); hold on, plot(MQRS, ecg(MQRS),'ok');
hold on, plot(cMQRS, ecg(cMQRS),'or');
legend('ecg','MQRS initial','MQRS adjusted');
set(findall(gcf,'type','text'),'fontSize',14,'fontWeight','bold');
end
end