function [position,timeqrs,lastqrs,intervalo,numlatdet,time,messages]= fiducialf(position,firstnewsamp,samp,heasig,w,eps,nlead,lastqrs,timeqrs,numlatdet,ultimo_anot,messages)
% [position,timeqrs,lastqrs,intervalo,numlatdet,time,messages]= fiducialf(position,firstnewsamp,samp,heasig,w,eps,nlead,lastqrs,timeqrs,numlatdet,messages)
%
% This script performs the QRS fiducial point detection of ECG
% In this version the thresholds eps(x) are not beat to beat actualized.
%
%Input Parameters:
% position: struct vector with the detected points
% firstnewsamp: last sample analyzed in the previous excerpt
% samp: samples included in the current excerpt (borders excluded)
% heasig: struct vector with header information
% w: matrix with WT scales 1 to 5
% eps: thresholds
% nlead: lead to be delineated
% lastqrs: last detecter QRS time (shoud be timeqrs(end))
% timeqrs: QRS times
% numlatdet: Number of detected beats until now
% ultimo_anot: last annotation in the previous segment
%Output Parameters:
% intervalo: numeration of the beats processed in this segment
% time: beats processed in this segment
% Actualized parameters: position,timeqrs,lastqrs,numlatdet,timeqrs
%
% Last update: Rute Almeida 07FEB2012
%
% Designed for MATLAB Version R12; tested with MATLAB Version R13
if nargin<12 || ~isfield(messages.setup,'wavedet')
messages.setup.wavedet=[];
end
if isfield(messages.setup.wavedet,'auxff') && messages.setup.wavedet.auxff==1
ff=1;
if isfield(messages.setup.wavedet,'aux_aa')
aa=messages.setup.wavedet.aux_aa;
else
aa=[1 1500];
end
else
ff=0;
end
if isempty(heasig)
end
if ff==1
figure
plot(w(:,1),'k')
fig1=gca;
set(gca,'xlim',aa)
hold on
ylabel('W_{2^1}')
figure
plot(messages.sig,'k')
fig2=gca;
hold on
set(gca,'xlim',aa)
ylabel('ECG')
figure
plot(w(:,3),'k')
fig3=gca;
set(gca,'xlim',aa)
hold on
ylabel('W_{2^3}')
figure
subplot(4,1,1)
plot(w(:,4),'k')
ylabel('W_{2^4}')
hold on
figWT4=gca;
set(gca,'xlim',aa)
subplot(4,1,2)
plot(w(:,3),'k')
ylabel('W_{2^3}')
hold on
figWT3=gca;
set(gca,'xlim',aa)
subplot(4,1,3)
plot(w(:,2),'k')
hold on
ylabel('W_{2^2}')
figWT2=gca;
set(gca,'xlim',aa)
subplot(4,1,4)
plot(w(:,1),'k')
hold on
ylabel('W_{2^1}')
figWT1=gca;
set(gca,'xlim',aa)
figure
subplot(4,1,1)
plot(messages.sig,'k')
ylabel('ECG')
hold on
% figECG=gca;
set(gca,'xlim',aa)
subplot(4,1,2)
plot(w(:,3),'k')
ylabel('W_{2^3}')
hold on
figECGWT3=gca;
set(gca,'xlim',aa)
subplot(4,1,3)
plot(w(:,2),'k')
hold on
ylabel('W_{2^2}')
figECGWT2=gca;
set(gca,'xlim',aa)
subplot(4,1,4)
plot(w(:,1),'k')
hold on
ylabel('W_{2^1}')
figECGWT1=gca;
set(gca,'xlim',aa)
end
if~isfield(messages.setup.wavedet,'firstmin')
messages.setup.wavedet.firstmin=0.050; % minimum time in sec to discart at the beggining
end
if ~isfield(messages.setup.wavedet,'nghbhd')
messages.setup.wavedet.nghbhd=0.025; % neighbourwood for maximum across scales in sec
end
if ~isfield(messages.setup.wavedet,'peakcriteria')
messages.setup.wavedet.peakcriteria=1.2; % If a modulus is more than peakcriteria times greater than others in the sane scale the greatest is chosen
end
if ~isfield(messages.setup.wavedet,'thmax')
messages.setup.wavedet.thmax=5+log(5); %maximum value for log(abs(W_3)) in a maximun line
end
if ~isfield(messages.setup.wavedet,'th_schb')
messages.setup.wavedet.th_schb=5+log(5); %maximum value for log(abs(W_3)) in a maximun line
end
if ~isfield(messages.setup.wavedet,'thfraction')
messages.setup.wavedet.thfraction=0.5; % tolerance below the average log(abs(W_3)) in a maximun line
end
if ~isfield(messages.setup.wavedet,'intvlthr2')
messages.setup.wavedet.intvlthr2=0.15;%threshold interval to consider isolated maximum lines
end
if ~isfield(messages.setup.wavedet,'intvlthr1')
messages.setup.wavedet.intvlthr1=0.12;%threshold interval for considering redundancy
end
if ~isfield(messages.setup.wavedet,'intvlthr1_2')
messages.setup.wavedet.intvlthr1_2=1.5;%reduction on threshold interval for considering redundancy if it subsists
end
if ~isfield(messages.setup.wavedet,'pictime')
messages.setup.wavedet.pictime=0.1;%length of the search window for associated maximumlines in extra verification step
end
if ~isfield(messages.setup.wavedet,'timelapthr')
messages.setup.wavedet.timelapthr=0.125;% maximum interval between two assiciated maximumlines in sec
end
if ~isfield(messages.setup.wavedet,'refrper')
messages.setup.wavedet.refrper=0.275;% Refractary period after a QRS detection in sec
end
if ~isfield(messages.setup.wavedet,'refrper_reduc')% 22NOV2011
messages.setup.wavedet.refrper_reduc=0;% 22NOV2011
end
if ~isfield(messages.setup.wavedet,'rrmaxmaxlim')
messages.setup.wavedet.rrmaxmaxlim=0.3; % search back criteria: searchback if RR>wavedet.rrmaxmaxlim sec
end
if ~isfield(messages.setup.wavedet,'rrmaxlim')
messages.setup.wavedet.rrmaxlim=1.5;% search back criteria: searchback if RR>wavedet.rrmaxlim * median(RR)
end
if ~isfield(messages.setup.wavedet,'rrmaxlimsec')
messages.setup.wavedet.rrmaxlimsec=1.5;% search back criteria: searchback if RR>wavedet.rrmaxlimsec
end
if ~isfield(messages.setup.wavedet,'rrmax_rbt')
messages.setup.wavedet.rrmax_rbt=0.1; % search back criteria only if >wavedet.rrmax_rbt sec o
end
if ~isfield(messages.setup.wavedet,'gapthrs')
messages.setup.wavedet.gapthrs=2; % fraction of mean(rr) to decide if there are a final "gap"
end
if ~isfield(messages.setup.wavedet,'pwaveout')
messages.setup.wavedet.pwaveout=0.2;% mimimum time in sec before each beat in order to ensure a complete beat
end
if ~isfield(messages.setup.wavedet,'twaveout')
messages.setup.wavedet.twaveout=0.45; % mimimum time in sec after each beat in order to ensure a complete beat
end
if ~isfield(messages.setup.wavedet,'refrperdef')
messages.setup.wavedet.refrperdef=messages.setup.wavedet.refrper;%
end
if messages.setup.wavedet.refrper_reduc>=1 % 22NOV2011
messages.setup.wavedet.refrper_reduc=0.9;% 22NOV2011
end
refrper=messages.setup.wavedet.refrper(end);
peakcriteria=messages.setup.wavedet.peakcriteria;
%global regularity
intervalo=[];
first = firstnewsamp-round(samp(1))+1;
first = max(ceil(messages.setup.wavedet.freq*messages.setup.wavedet.firstmin), first - ceil(1*messages.setup.wavedet.freq));
n = modmax(w(:,4),first,eps(4,nlead),0); % Maximum moduli at scale 4 % Rute 02.Dec.04
if ~isempty(n),
signo = sign(w(n,4))'; % Maximum or minimum?
neighb = ceil(messages.setup.wavedet.freq*messages.setup.wavedet.nghbhd); % Definition of neighborhood 25 ms.
%it makes sence to verify if all m(4,:) diffre at least neighb....
m = zeros(4,length(n));
m(4,:) = n'; % Positions of maximum moduli
for k = 1:size(m,2), % For each maximum at scale 4
n3 = modmax(w(max(1,n(k)-neighb):min(n(k)+neighb,size(w,1)),3),2,eps(3,nlead),signo(k));% Rute 02.Dec.04
n3 =n(k)-neighb-1+n3; % Check the neighborhood at scale 3
num = length(n3);
if num>0,
if num==1, % If only one neighbour maximum
m(3,k)= n3; % at scale 3
elseif num>1 % If more than one
if ~isempty(find(max(abs(w(n3,3))./(abs(w(n3,3)))<peakcriteria))==1),
% If a modulus is more than peakcriteria times greater than others
[aux1,ind]=max(abs(w(n3,3))); %#ok<ASGLU> % the greatest is chosen
m(3,k)= n3(ind);
else % Choose the one with minimum distance
[aux1,ind]=min(abs(m(4,k)-n3)); %#ok<ASGLU>
m(3,k)= n3(ind);
end
end
end
end
ind = find(m(3,:)==0); % Discard all maximum lines with no
if ff==1
axes(figWT4);
plot(figWT4,m(4,:),w(m(4,:),4),'ob')
plot(figWT4,[aa(1) aa(2)],[eps(4,nlead) eps(4,nlead)],':b')
plot(figWT4,[aa(1) aa(2)],-[eps(4,nlead) eps(4,nlead)],':b')
if ~isempty(ind)
plot(figWT4,m(4,ind),w(m(4,ind),4),'xr')
plot(figWT4,[m(4,ind(1))-neighb m(4,ind(1))-neighb],[-5000 5000],':r')
plot(figWT4,[m(4,ind(1))+neighb m(4,ind(1))+neighb],[-5000 5000],':r')
axes(figWT3);
plot(figWT3,[m(4,ind(1))-neighb m(4,ind(1))-neighb],[-5000 5000],':r')
plot(figWT3,[m(4,ind(1))+neighb m(3,ind(1))+neighb],[-5000 5000],':r')
end
plot(figWT3,m(3,m(3,:)~=0),w(m(3,m(3,:)~=0),3),'ob')
plot(figWT3,[aa(1) aa(2)],[eps(3,nlead) eps(3,nlead)],':b')
plot(figWT3,[aa(1) aa(2)],-[eps(3,nlead) eps(3,nlead)],':b')
end
m(:,ind) = []; % associated maximum at scale 3
signo(ind)= [];
% Search for maximum modula in the neighborhood at scale 2
for k = 1:size(m,2),
n2 =modmax(w(max(1,m(3,k)-neighb):min(m(3,k)+neighb,size(w,1)),2),2,eps(2,nlead),signo(k));% Rute 02.Dec.04
n2 =m(3,k)-neighb-1+n2;
num = length(n2);
if num>0,
if num==1, % If only one neighbour maximum
m(2,k)= n2; % at scale 2
elseif num>1 % If more than one
if ~isempty(find(max(abs(w(n2,2))./(abs(w(n2,2)))<peakcriteria))==1),
[aux1,ind]=max(abs(w(n2,2))); %#ok<ASGLU> % greatest modulus
m(2,k)= n2(ind);
else % shortest distance
[aux1,ind]=min(abs(m(3,k)-n2)); %#ok<ASGLU>
m(2,k)= n2(ind);
end
end
end
end
ind = find(m(2,:)==0); % Discard all maximum lines with no
if ff==1
% axes(figWT3);
if ~isempty(ind)
plot(figWT3,[m(3,ind(1))-neighb m(3,ind(1))-neighb],[-5000 5000],':r')
plot(figWT3,[m(3,ind(1))+neighb m(3,ind(1))+neighb],[-5000 5000],':r')
plot(figWT3,m(3,ind),w(m(3,ind),3),'xr')
axes(figWT2);
plot(figWT2,[m(3,ind(1))-neighb m(3,ind(1))-neighb],[-5000 5000],':r')
plot(figWT2,[m(3,ind(1))+neighb m(3,ind(1))+neighb],[-5000 5000],':r')
end
plot(figWT2,m(2,m(2,:)~=0),w(m(2,m(2,:)~=0),2),'ob')
plot(figWT2,[aa(1) aa(2)],[eps(2,nlead) eps(2,nlead)],':b')
plot(figWT2,[aa(1) aa(2)],-[eps(2,nlead) eps(2,nlead)],':b')
end
m(:,ind) = []; % associated maximum at scale 2
signo(ind)=[];
for k = 1:size(m,2),
n1 = modmax(w(max(1,m(2,k)-neighb):min(m(2,k)+neighb,size(w,1)),1),2,eps(1,nlead),signo(k)); % Rute 02.Dec.04
n1 =m(2,k)-neighb-1+n1;
n1 = n1(n1 > 0 & n1 <= size(w,1) );
num = length(n1);
if num>0,
if num==1, % If only one neighbour maximum
m(1,k)= n1; % at scale 1
elseif num>1 % If more than one
if ~isempty(find(max(abs(w(n1,1))./(abs(w(n1,1)))<peakcriteria))==1),
[aux1,ind]=max(abs(w(n1,1))); %#ok<ASGLU> % greatest modulus
m(1,k)= n1(ind);
else % shortest distance
[aux1,ind]=min(abs(m(2,k)-n1)); %#ok<ASGLU>
m(1,k)= n1(ind);
end
end
end
end
ind = find(m(1,:)==0); %Discard all maximum lines with no
if ff==1
% axes(figWT2);
if ~isempty(ind)
plot(figWT2,m(2,ind),w(m(2,ind),2),'xr')
end
plot(figWT1,m(1,m(1,:)~=0),w(m(1,m(1,:)~=0),1),'ob')
plot(figWT1,[aa(1) aa(2)],[eps(1,nlead) eps(1,nlead)],':b')
plot(figWT1,[aa(1) aa(2)],-[eps(1,nlead) eps(1,nlead)],':b')
end
m(:,ind) = []; % associated maximum at scale 1
signo(ind)=[];
%%% Regularity Exponent Validation %%%%%%%%%
% alpha is proportional to log(a3(nk3))-log(a1(nk1))
alpha = log(abs(w(m(3,:),3))); % !!!!
% alphalinha = log(abs(w(m(3,:),3))) - log (abs(w(m(1,:),1))); % !!! (1) % cometario na versao original%%%%%% Rute 24/04/02
%
% %%%% SO
% alpha1 = log(abs(w(m(2,:),2))./abs(w(m(1,:),1)));
% alpha2 = log(abs(w(m(3,:),3))./abs(w(m(2,:),2)));
% alpha3 = log(abs(w(m(4,:),4))./abs(w(m(3,:),3)));
% %%%% SO: you are right alpha=alpha1+alpha2 is calculated acording to (1)
% %%%% SO: por ahorro computacional se puede eliminar el log, y en el
% %%%% SO: valor del threshold si queremos.
%
% th was not being using at all in many cases, because thmax < min(alpha)
% th=min(messages.setup.wavedet.thmax,mean(alpha)-messages.setup.wavedet.thfraction);
th=mean(alpha)-messages.setup.wavedet.thfraction;
% thlinha=min(5+log(5),mean(alphalinha)-.5); %%%%%% Rute 24/04/02
ind = find(alpha<=th); % always empty?????!????1
% ind=find(alpha3>=1 || alpha1<=0.9 || alpha2<=0.9);% elimination of large T waves
if ff==1
plot(fig1,m(1,:),w(m(1,:)),'ob')
legend(fig1,{'WT_1','mm lines'},'Location','NorthEastOutside')
plot(m(1,ind),w(m(1,ind)),'xr')
if ~isempty(ind)
ll={'mm lines','rejected mm lines'};
else
ll={'mm lines'};
end
plot(figECGWT1,m(1,:),w(m(1,:),1),'ob')
plot(figECGWT1,m(1,ind),w(m(1,ind),1),'xr')
plot(figECGWT2,m(2,:),w(m(2,:),2),'ob')
plot(figECGWT2,m(2,ind),w(m(2,ind),2),'xr')
plot(figECGWT3,m(3,:),w(m(3,:),3),'ob')
plot(figECGWT3,m(3,ind),w(m(3,ind),3),'xr')
plot(fig3,m(3,:),w(m(3,:),3),'ob')
legend(fig3,{'WT_3','mm lines'},'Location','NorthEastOutside')
plot(m(3,ind),w(m(3,ind),3),'xr')
% axes(fig1);
legend(fig1,['WT_1',ll],'Location','NorthEastOutside')
% axes(fig3);
legend(fig3,['WT_3',ll],'Location','NorthEastOutside')
figure;
hs = subplot(3,1,1);
plot(hs,w(:,3),'k')
ylabel(hs,'W_{2^3}')
hold on
plot(hs,m(3,:),w(m(3,:),3),'ob')
plot(hs,m(3,ind),w(m(3,ind),3),'xr')
set(hs,'xlim',aa)
hs3 = subplot(3,1,3);
plot(hs3,messages.sig,'k')
hold on
set(hs3,'xlim',aa)
ylabel(hs3,'ECG')
set(hs3,'xlim',aa)
hs2 = subplot(3,1,2);
plot(hs2,m(3,:),alpha,'k')
hold on
plot(hs2,m(3,:),alpha,'ob')
plot(hs2,[0 261770],[mean(alpha)-messages.setup.wavedet.thfraction mean(alpha)-messages.setup.wavedet.thfraction],':g')
plot(hs2,[0 261770],[th th],':r'),
plot(hs2,m(3,ind),alpha(ind),'xr')
set(hs2,'xlim',aa)
ylabel(hs2,'\alpha')
%legend('alpha','alpha','mean alpha-thfraction','th','Location','NorthEastOutside')
end
if ~isempty(ind) % and high frequency noise peaks
m(:,ind) = [];
signo(ind)=[];
end
thresinterval = ceil(messages.setup.wavedet.intvlthr2* messages.setup.wavedet.freq);
% threshold interval to consider isolated maximum lines
ind = find( ((m(1,2:end-1)-m(1,1:end-2))>thresinterval) ...
& (m(1,3:end)- m(1,2:end-1))>thresinterval)+1;
if size(m,2)<2 && ~isempty(ind), ind=1; end % when only one maximum %16DEZ08
if ff==1
% axes(fig1);
plot(fig1,m(1,ind),w(m(1,ind)),'+r')
if ~isempty(ind)
ll=[ll {'isolated mm lines'}];
end
% axes(fig1);
legend(fig1,['WT_1',ll],'Location','NorthEastOutside')
% axes(fig3);
plot(fig3,m(3,ind),w(m(3,ind),3),'+r')
legend(['WT_3',ll],'Location','NorthEastOutside')
try
% axes(figECGWT1);
plot(figECGWT1,m(1,ind),w(m(1,ind),1),'+r')
plot(figECGWT1,[m(1,ind(1))-thresinterval m(1,ind(1))-thresinterval],[-5000 5000],':r')
plot(figECGWT1,[m(1,ind(1))+thresinterval m(1,ind(1))+thresinterval],[-5000 5000],':r')
% axes(figECGWT2);
plot(figECGWT2,m(2,ind),w(m(2,ind),2),'+r')
% axes(figECGWT3);
plot(figECGWT3,m(3,ind),w(m(3,ind),3),'+r')
catch me
me.message = 'Error plotting WT signal';
end
end
m(:,ind) = []; %Discard isolated maximum lines
signo(ind)=[];
% Threshold interval for considering redundancy
redundant= [];
thresinterval = ceil(messages.setup.wavedet.intvlthr1*messages.setup.wavedet.freq); % 120 ms. (li)
for l = find(signo>0), % For each positive maximum line
if ~any(redundant ==l), % If it has not been declared redundant yet
ind=find((m(3,:)>m(3,l)-messages.setup.wavedet.intvlthr1_2*thresinterval)&(m(3,:)<m(3,l)+messages.setup.wavedet.intvlthr1_2*thresinterval)&signo>0);
% index of positive lines near the present one (including it)
if length(ind)>1, % If more than one --> redundancy
[mx,ind2]= max(abs(w(m(3,ind),3))); %#ok<ASGLU>
ind(ind2)=[];
redundant = [redundant ind]; %#ok<AGROW> % All but the greatest are redundant
end
end
end
if ff==1
% axes(fig1);
plot(fig1,m(1,redundant),w(m(1,redundant)),'sm')
% axes(fig3);
plot(fig3,m(3,redundant),w(m(3,redundant),3),'sm')
if ~isempty(redundant)
ll=[ll {'redundant mm lines'}];
end
% axes(figECGWT1);
plot(figECGWT1,m(1,redundant),w(m(1,redundant)),'sm')
% axes(figECGWT2);
plot(figECGWT2,m(2,redundant),w(m(2,redundant),2),'sm')
% axes(figECGWT3);
plot(figECGWT3,m(3,redundant),w(m(3,redundant),3),'sm')
plot(figECGWT3,[m(3,206)-messages.setup.wavedet.intvlthr1_2*thresinterval m(3,206)-messages.setup.wavedet.intvlthr1_2*thresinterval],[-5000 5000],':r')
plot(figECGWT3,[m(3,206)+messages.setup.wavedet.intvlthr1_2*thresinterval m(3,206)+messages.setup.wavedet.intvlthr1_2*thresinterval],[-5000 5000],':r')
end
m(:,redundant)=[]; % Discard redundant lines
signo(redundant)=[];
redundant =[];
for l = find(signo>0), % For each remaining positive maximum line
ind=find((m(3,:)>m(3,l)-thresinterval)&(m(3,:)<m(3,l)+thresinterval)&signo<0);
% Search for negative minima near it
if length(ind)>1, % If more than one ---> redundancy
aux = abs(w(m(3,ind),3)./(m(3,ind)'-m(3,l)));
% auxiliary variable: height over distance
[mx,indmx]=max(aux);
ind2=ind;
aux(indmx)=[]; ind2(indmx)=[];
if all((mx./aux)>peakcriteria), % RULE 2 (see PFC or Li's paper)
redundant = [redundant ind2]; %#ok<AGROW>
else
[aux,aux2]=min(abs(m(3,l)-m(3,ind))); %#ok<ASGLU>
ind(aux2)=[]; % RULE 1 (see PFC or Li's paper)
redundant = [redundant ind]; %#ok<AGROW>
end
end
end
if ff==1
% axes(fig1);
plot(fig1,m(1,redundant),w(m(1,redundant)),'*r')
% axes(fig3);
plot(fig3,m(3,redundant),w(m(3,redundant),3),'*r')
if ~isempty(redundant)
ll=[ll {'redundant mm lines'}];
end
try
% axes(figECGWT1);
plot(figECGWT1,m(1,redundant(end)),w(m(1,redundant(end))),'sm')
% axes(figECGWT2);
plot(figECGWT2,m(2,redundant(end)),w(m(2,redundant(end)),2),'sm')
% axes(figECGWT3);
plot(figECGWT3,m(3,redundant(end)),w(m(3,redundant),3),'sm')
plot(figECGWT3,[m(3,182)-messages.setup.wavedet.intvlthr1_2*thresinterval m(3,182)-messages.setup.wavedet.intvlthr1_2*thresinterval],[-5000 5000],':r')
plot(figECGWT3,[m(3,182)+messages.setup.wavedet.intvlthr1_2*thresinterval m(3,182)+messages.setup.wavedet.intvlthr1_2*thresinterval],[-5000 5000],':r')
catch me
me.message = 'Error plotting WT signal';
end
end
m(:,redundant)=[]; % Discard redundant lines
signo(redundant)=[];
redundant =[];
for l = find(signo<0), % For each remaining negative minimum line
ind = find((m(3,:)>m(3,l)-thresinterval)&(m(3,:)<m(3,l)+thresinterval)&signo>0);
% Search for positive maxima near it
if length(ind)>1, % If more than one ----> redundancy
aux = abs(w(m(3,ind),3)./(m(3,ind)'-m(3,l)));
% auxiliary variable: height over distance
[mx,indmx]=max(aux);
ind2=ind;
aux(indmx)=[]; ind2(indmx)=[];
if all((mx./aux)>peakcriteria),
redundant = [redundant ind2]; %#ok<AGROW> % RULE 2
else
[aux,aux2]=min(abs(m(3,l)-m(3,ind))); %#ok<ASGLU>
ind(aux2)=[]; % RULE 1
redundant = [redundant ind]; %#ok<AGROW>
end
end
end
if ff==1
% axes(fig1);
plot(fig1,m(1,redundant),w(m(1,redundant)),'*r')
% axes(fig3);
plot(fig3,m(3,redundant),w(m(3,redundant),3),'*r')
if ~isempty(redundant)
ll=[ll {'redundant mm lines'}];
end
end
m(:,redundant)=[]; % Discard redundant lines
signo(redundant)=[];
%%%% isolated maximum lines resulting from Discarded redundant
ind = find( ((m(1,2:end-1)-m(1,1:end-2))>thresinterval) ...
& (m(1,3:end)- m(1,2:end-1))>thresinterval)+1;
if size(m,2)<2 && ~isempty(ind), ind=1; end % when only one maximum %16DEZ08
if ff==1
% axes(fig1);
legend(fig1,['WT_1',ll],'Location','NorthEastOutside')
% axes(fig3);
legend(fig3,['WT_3',ll],'Location','NorthEastOutside')
% axes(fig1);
plot(fig1,m(1,ind),w(m(1,ind)),'+m')
% axes(fig3);
plot(fig3,m(3,ind),w(m(3,ind),3),'+m')
if ~isempty(ind)
ll=[ll {'isolated mm lines'}];
end
% axes(fig1);
legend(fig1,['WT_1',ll],'Location','NorthEastOutside')
% axes(fig3);
legend(fig3,['WT_3',ll],'Location','NorthEastOutside')
end
m(:,ind) = []; %Discard isolated maximum lines
signo(ind)=[];
if size(m,2)<2
m=[];
signo=[];
end
eliminar=[];
%%%%extra protection% 23MAR09
% discard if no peak before and after messages.setup.wavedet.pictime ms
for ii=1:size(m,2)
pa = picant(w(max(1,m(2,ii)-round(messages.setup.wavedet.pictime*messages.setup.wavedet.freq)):m(2,ii),2),m(2,ii));
% first peak before detected qrs position at scale 2
pp = picpost(w(m(2,ii):min(size(w,1),m(2,ii)+round(messages.setup.wavedet.pictime*messages.setup.wavedet.freq)),2),m(2,ii));
%if isempty(pa) || isempty(pp)
if isempty(pa) && isempty(pp) %NOV2011
eliminar=[eliminar ii]; %#ok<AGROW>
end
end
if ff==1
% axes(fig1);
plot(fig1,m(1,eliminar),w(m(1,eliminar)),'.r')
% axes(fig3);
plot(fig3,m(3,eliminar),w(m(3,eliminar),3),'.r')
if ~isempty(eliminar)
ll=[ll {'extra protection'}];
end
% axes(fig1);
legend(fig1,['WT_1',ll],'Location','NorthEastOutside')
% axes(fig3);
legend(fig3,['WT_3',ll],'Location','NorthEastOutside')
end
m(:,eliminar)=[];
signo(eliminar)=[];
% QRS peak detection / wavelet Zero cross detection
timelap = ceil(messages.setup.wavedet.timelapthr*messages.setup.wavedet.freq); % the two maximum lines defining a QRS
% should not be separated more than timelapthr
time = [];
aux=[];
auxm=[];
if size(m,2)>1 % 18 OUT2011
for l = find(signo>0), % For each positive maximum line
if (l==1) % Special case: first line
if (signo(2)<0)&&(m(1,2)-m(1,1)<thresinterval),
ind = zerocros(w(m(1,l):min(m(1,l)+timelap,size(w,1)),1));
% Zero crossing at scale 1
time = [time ind+m(1,l)-1]; %#ok<AGROW>
aux=[aux abs(w(m(1,l))-w(m(1,l+1)))]; %#ok<AGROW>
auxm=[auxm; [m(1,l) m(1,l+1)]]; %#ok<AGROW>
end
elseif (l==size(m,2)) % Special case: last line
if (signo(l-1)<0)&&(m(1,end)-m(1,end-1)<thresinterval),
ind = zerocros(w(m(1,l-1):min(m(1,l-1)+timelap,size(w,1)),1));
time = [time ind+m(1,l-1)-1]; %#ok<AGROW>
aux=[aux abs(w(m(1,l-1))-w(m(1,l)))]; %#ok<AGROW>
auxm=[auxm ; [m(1,l-1) m(1,l)]]; %#ok<AGROW>
end
elseif signo(l+1)<0 && ((signo(l-1)>0) ...
|| ((m(1,l+1)-m(1,l))<(m(1,l)-m(1,l-1))))
ind = zerocros(w(m(1,l):min(m(1,l)+timelap,size(w,1)),1));
time = [time ind+m(1,l)-1]; %#ok<AGROW>
aux=[aux abs(w(m(1,l))-w(m(1,l+1)))]; %#ok<AGROW>
auxm=[auxm;[ m(1,l) m(1,l+1)]]; %#ok<AGROW>
elseif signo(l-1)<0,
ind = zerocros(w(m(1,l-1):min(m(1,l-1)+timelap,size(w,1)),1));
time = [time ind+m(1,l-1)-1]; %#ok<AGROW>
aux=[aux abs(w(m(1,l-1))-w(m(1,l)))]; %#ok<AGROW>
auxm=[auxm;[ m(1,l-1) m(1,l)]]; %#ok<AGROW>
end
end
end
if ff==1
plot(fig2,time,messages.sig(time),'ob')
ll_ECG={'ECG','QRS candidates'};
legend(fig2,ll_ECG,'Location','NorthEastOutside')
end
% Refractary period after a QRS detection. Retain only the waves with
% more energy within the refractary period.
if ~isempty(time),
rr = (time(2:end)-time(1:end-1));
ind = find(rr<ceil(refrper*messages.setup.wavedet.freq));
for auxi=1:length(ind) %RUTE 27Jun11
[M,ii]=min([aux(ind(auxi)) aux(ind(auxi)+1)]); %#ok<ASGLU> % 18JUL2011
if ii==2, ind(auxi)=ind(auxi)+1; end
end
% ind=ind+1; % olde version % always the second one is eliminated!!!
if ff==1
% axes(fig2);
plot(fig2,time(ind),messages.sig(time(ind)),'xr')
% axes(fig1);
plot(fig1,m(1,ind),w(m(1,ind)),'xr')
% axes(fig3);
plot(fig3,m(3,ind),w(m(3,ind),3),'xr')
if ~isempty(ind)
ll=[ll {'rejected by the refractary period'}];
ll_ECG=[ll_ECG {'rejected by the refractary period'}];
end
% axes(fig2);
legend(fig2,ll_ECG,'Location','NorthEastOutside')
% axes(fig1);
legend(fig1,['WT_1',ll],'Location','NorthEastOutside')
% axes(fig3);
legend(fig3,['WT_3',ll],'Location','NorthEastOutside')
end
time(ind)=[]; %RUTE 27Jun11
aux(ind)=[];
if abs(samp(time(1)) - lastqrs)< refrper*messages.setup.wavedet.freq, % Refractory period for
time(1)=[];
aux(1)=[];%#ok<NASGU> % last beat of previous
end % fragment
end
% search back if no detection in messages.setup.wavedet.rrmaxlim RR median
if ~isempty(time),
%%%%%%%%%%%%%%%% caso em que rr tem dimensao inferior ao necessario
if length(time)>1 %%%%%%%%%%%%%%%%introduzido a 17/05/02 Rute
rr = ([max(samp(time(1))-lastqrs,time(1)) time(2:end)-time(1:end-1)]); %rr = ([samp(time(1))-lastqrs, time(2:end)-time(1:end-1)]);
rrmed = [rr(1) 0.5*(rr(1)+rr(2)) median([rr(3:end);rr(2:end-1);rr(1:end-2)])];
ind = find((rr(2:end)>messages.setup.wavedet.rrmaxlim*rrmed(1:end-1)) & (rr(1:end-1)>messages.setup.wavedet.rrmaxmaxlim*messages.setup.wavedet.freq) | rr(2:end)>messages.setup.wavedet.rrmaxlimsec*messages.setup.wavedet.freq);
%%%% 22NOV2011
if messages.setup.wavedet.refrper_reduc>0 && refrper>0.250 % 22NOV2011
while sum(rrmed<refrper*messages.setup.wavedet.freq*4/3)/length(rrmed)>0.6
refrper=refrper*messages.setup.wavedet.refrper_reduc;
end % 22NOV2011
if messages.setup.wavedet.refrper(end)~=refrper
messages.setup.wavedet.refrper(end+1)=refrper;
messages.setup.wavedet.refrperdef(end+1)=samp(time(1));
end
end
if ff==1
if ~isempty(ind)
ll=[ll {'searchback'}];
ll_ECG=[ll_ECG {'searchback'}];
% axes(fig2);
plot(fig2,time(ind),messages.sig(time(ind))+50,'vc')
legend(ll_ECG,'Location','NorthEastOutside')
% axes(fig1);
plot(fig1,time(ind),w(time(ind),1)+200,'vc')
% axes(fig3);
plot(fig3,time(ind),w(time(ind),3)+200,'vc')
% axes(fig1);
legend(fig1,'WT_1',ll,'Location','NorthEastOutside')
% axes(fig3);
legend(fig3,'WT_3',ll,'Location','NorthEastOutside')
end
end
for l = ind, % For all rr intervals greater than 1.5RR
interv = (time(l)+ceil(refrper*messages.setup.wavedet.freq)) : (time(l+1)-ceil(refrper*messages.setup.wavedet.freq));
if length(interv)>messages.setup.wavedet.rrmax_rbt*messages.setup.wavedet.freq, % for robustness
% [nuevo, messages] = searchbk(w(interv,:),interv(1),eps(:,nlead)',messages.setup.wavedet.freq,messages.setup.wavedet.th_schb,messages);
[nuevo, messages] = searchbk(w(interv,:),interv(1),eps(:,nlead)',messages.setup.wavedet.freq,th,messages); %FEB16_2012
% searchbk returns the new time of occurrence of a QRS
time = [time nuevo]; %#ok<AGROW> % New times are added
% time = sort(time); % and sorted
if ff==1
% axes(fig2);
plot(fig2,nuevo,messages.sig(nuevo),'oc')
if ~isempty(nuevo)
ll_ECG=[ll_ECG {'new candidates by searchback'}]; %#ok<AGROW>
end
legend(fig2,ll_ECG,'Location','NorthEastOutside')
end
end
end %%%%%%%%%%%%%%%%introduzido a 17/05/02 Rute
else % only one detection DEZ2011
if time > max(messages.setup.wavedet.gapthrs,ceil(refrper*messages.setup.wavedet.freq))%initial gap % 24ABRIL2010
% if time > messages.setup.wavedet.gapthrs %initial gap
if (~exist('th','var') || isnan(th)), th = messages.setup.wavedet.thmax ; end
[nuevo, messages]= searchbk(w(1:(time-ceil(refrper*messages.setup.wavedet.freq))),1,eps(:,nlead)',messages.setup.wavedet.freq,th,messages);
time = [time nuevo];
end
end
else % in case there is no detection, searchback anyway
if (~exist('th','var') || isnan(th)), th = messages.setup.wavedet.thmax ; end
%nuevo = searchbk(w,1,eps,messages.setup.wavedet.freq, messages.setup.wavedet.th_schb); %Rute 02.Dec.04
[nuevo, messages]= searchbk(w,1,eps(:,nlead)',messages.setup.wavedet.freq, th,messages);
time = nuevo;
if ff==1
% axes(fig2);
plot(fig2,nuevo,messages.sig(nuevo),'oc')
if ~isempty(nuevo)
ll_ECG=[ll_ECG {'new candidates by searchback'}];
end
legend(fig2,ll_ECG,'Location','NorthEastOutside')
end
end
if ~isempty(time),
time = sort(time); % Rute 14DEz
if length(time)>1
rr = ([samp(time(1))-lastqrs, time(2:end)-time(1:end-1)]);
rrmed = [rr(1) 0.5*(rr(1)+rr(2)) median([rr(3:end);rr(2:end-1);rr(1:end-2)])]; %RUTE OUT2011
else
rrmed =1;
end
if length(w)-time(end) > messages.setup.wavedet.gapthrs*rrmed, % If there is a final "gap"
interv= (time(end)+ceil(refrper*messages.setup.wavedet.freq) : length(w)); %RUTE OUT2011
if length(interv)>messages.setup.wavedet.rrmax_rbt*messages.setup.wavedet.freq,
[nuevo, messages] = searchbk(w(interv,:),interv(1),eps(:,nlead)',messages.setup.wavedet.freq, th,messages);
time = [time nuevo];
if ff==1
% axes(fig2);
plot(fig2,nuevo,messages.sig(nuevo),'oc')
if ~isempty(nuevo)
ll_ECG=[ll_ECG {'new candidates by searchback'}];
end
legend(fig2,ll_ECG,'Location','NorthEastOutside')
end
end
end
%%% The first and last beats must be complete, to detect P-wave and T-wave %%%
%%% We use an overlap that assures that the beat will be complete always in
%%% one of the segments.
if (time(1) <= messages.setup.wavedet.pwaveout*messages.setup.wavedet.freq),
time(1)=[];
end
% If the P-wave is not in the present segment, discard the beat, because the
% the beat should have been detected in the last segment.
if ~isempty(lastqrs)
time(samp(time)<= lastqrs + refrper*messages.setup.wavedet.freq)=[];
end
if ~isempty(time) %%%%%%% Rute 06/05/02 caso time(indrep)=time
% if beat before ultimo_anot in the previous segment %6AGO2007
if samp(time(1))<=ultimo_anot
time(1)=[];
end
if ~isempty(time) && (length(w)-time(end) < messages.setup.wavedet.twaveout*messages.setup.wavedet.freq), % si no ye a onda T drento
time(end)=[]; % de sig
end
% If last beat's T-wave is not in 'sig', the beat is detected in next one
if ~isempty(time) %Rute 10.04.05
timeqrs = [timeqrs samp(time)]; % QRS times
lastqrs = samp(time(end));
end
end
intervalo = [numlatdet+1 numlatdet+length(time)];
% numeration of the beats processed in this segment
numlatdet = intervalo(2);
position.qrs(intervalo(1):intervalo(2)) = samp(time); % Fill QRS position
end
if ff==1
% axes(fig2);
plot(fig2,time,messages.sig(time),'.g')
ll_ECG=[ll_ECG {'final SL marks'}];
legend(fig2,ll_ECG,'Location','NorthEastOutside')
% axes(fig1);
legend(fig1,['WT_1',ll],'Location','NorthEastOutside')
% axes(fig3);
legend(fig3,['WT_3',ll],'Location','NorthEastOutside')
end
else
intervalo=[];time=[]; %Rute 13Set06
end