function [new,messages] = searchbk (w, i, thres,freq,thres_alpha,messages)
%performs the searching back when no qrs has been detected for
% 1.5 times the RR before. i is the beginning sample. thres
% is the threshold for scale 4;
% Last update: Rute Almeida 27Jan2012
if nargin<6 || (isempty(freq) && (~isfield(messages.setup.wavedet,'freq') || ~(messages.setup.wavedet.freq>0)))
messages.errors=[messages.errors {'Fatal error in wavedet_3D: no sampling frequency variable found.'}];
warning(char(messages.errors(end)))
messages.errors_desc=[messages.errors_desc {'No sampling frequency variable give to seachbk from fiducialf.m.'}];
new=0;
return
elseif ~isempty(freq) && isfield(messages.setup.wavedet,'freq') && messages.setup.wavedet.freq>0 && freq~=messages.setup.wavedet.freq
messages.warnings=[messages.warnings {'sampling frequency used in searchbk is different from messages.setup.wavedet.freq.'}];
elseif isempty(freq)
freq=messages.setup.wavedet.freq;
end
if ~isfield(messages.setup.wavedet,'refrper') || ~isfield(messages.setup.wavedet,'peakcriteria') || ~isfield(messages.setup.wavedet,'timelapthr') || ~isfield(messages.setup.wavedet,'timelapthr') || ~isfield(messages.setup.wavedet,'intvlthr1') || ~isfield(messages.setup.wavedet,'intvlthr2') || ~isfield(messages.setup.wavedet,'nghbhd')
messages.errors=[messages.errors {'Fatal error in wavedet_3D: missing parameters on seachbk.'}];
warning(char(messages.errors(end)))
messages.errors_desc=[messages.errors_desc {'Setup value is missing on seachbk from fiducialf.m.'}];
new=0;
return
end
if ~isfield(messages.setup.wavedet,'intvlthr1_2_sbk')
messages.setup.wavedet.intvlthr1_2_sbk=2;%reduction on threshold interval for considering redundancy if it subsists
end
if ~isfield(messages.setup.wavedet,'thfraction_sbk')
if isfield(messages.setup.wavedet,'thfraction')
messages.setup.wavedet.thfraction_sbk=messages.setup.wavedet.thfraction;
else
messages.errors=[messages.errors {'Fatal error in wavedet_3D: missing parameters on seachbk.'}];
warning(char(messages.errors(end)))
messages.errors_desc=[messages.errors_desc {'Setup value thfraction is missing on seachbk from fiducialf.m.'}];
new=0;
return
end
end
if ~isfield(messages.setup.wavedet,'nghbhd')
messages.setup.wavedet.nghbhd=0.025; % neighbourwood fro maximum across scales in sec
end
thfraction_sbk=messages.setup.wavedet.thfraction_sbk;
peakcriteria=messages.setup.wavedet.peakcriteria;
timelap = ceil(messages.setup.wavedet.timelapthr*freq);
thres = thres/2; % When searching back we take half
thres(4) = thres(4)/2; % so new thres(4) = old thres (4)/4
n = modmax(w(:,4),2,thres(4),0); % Search maximum moduli at scale 4
neighb = ceil(freq*messages.setup.wavedet.nghbhd); % Neighbourhood = 25 ms
m = zeros(4,length(n));
signo = sign(w(n,4))';
if ~isempty(n),
m(4,:) = n';
end
new = []; %#ok<NASGU>
for k = 1:size(m,2),
window = [max(n(k)-neighb,1), min(n(k)+neighb,size(w,1))];
n3 = modmax(w(window(1):window(2),3),2,thres(3),signo(k));
n3 =window(1)-1+n3;
num = length(n3);
if num>0,
if num==1,
m(3,k)= n3;
elseif num>1
if length(find(max(abs(w(n3,3)))./abs(w(n3,3))<peakcriteria))==1,
[aux1,ind]=max(abs(w(n3,3))); %#ok<ASGLU> % greatest modulus
m(3,k)= n3(ind);
else % 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);
m(:,ind) = [];
signo(ind) = [];
% Search for maximum moduli in the neighborhood at scale 2
for k = 1:size(m,2),
window = [max(m(3,k)-neighb,1), min(m(3,k)+neighb,size(w,1))];
n2 =modmax(w(window(1):window(2),2),2,thres(2),signo(k));
n2 =window(1)-1+n2;
num = length(n2);
if num>0,
if num==1,
m(2,k)= n2;
elseif num>1
if length(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);
m(:,ind) = [];
signo(ind) = [];
for k = 1:size(m,2),
window = [max(m(2,k)-neighb,1), min(m(2,k)+neighb,size(w,1))];
n1 =modmax(w(window(1):window(2),1),2,thres(1),signo(k));
n1 =window(1)-1+n1;
num = length(n1);
if num>0,
if num==1,
m(1,k)= n1;
elseif num>1
if length(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
m(:,ind) = []; % associated maximum at scale 1
signo(ind) = [];
% Regularity Exponent Validation
% alpha proportional to log(a3(nk3))-log(a1(nk1))
alpha = log(abs(w(m(3,:),3))); %%!!!!
% alpha = log(abs(w(m(3,:),3))) - log (abs(w(m(1,:),1))); !!!
ind = find(alpha <= thres_alpha-thfraction_sbk); %%%% !!!
m(:,ind) = [];
signo(ind) = [];
thresinterval = ceil(messages.setup.wavedet.intvlthr2 * freq); % 120 ms. (Li)
if size(m,2)>2,
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 (m(1,2)-m(1,1))>thresinterval,
ind = [1 ind];
end
if (m(1,end)-m(1,end-1))>thresinterval,
ind = [ind size(m,2)];
end
m(:,ind) =[]; % Discard isolated maximum lines
signo(ind) = [];
elseif size(m,2)==2, % If only two lines
if (m(1,2)-m(1,1))>thresinterval, % discard them if too separated
m(:,1:2)=[];
signo(1:2)=[];
end
elseif size(m,2)==1, % If only one, discard it
m(:,1)=[];
signo(1)=[];
end
% Threshold interval for considering redundancy
redundant= [];
thresinterval = ceil(messages.setup.wavedet.intvlthr1*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_sbk*thresinterval)&(m(3,:)<m(3,l)+messages.setup.wavedet.intvlthr1_2_sbk*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
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
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)=[];
redundant = [redundant ind]; %#ok<AGROW> % RULE 1
end
end
end
m(:,redundant)=[];
signo(redundant)=[];
%%%%%%%%%%%%%% isolated maximum lines resulting from Discarded
%%%%%%%%%%%%%% redundant OUT1011
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
m(:,ind) = []; %Discard isolated maximum lines
signo(ind)=[];
if size(m,2)<2
m=[];
signo=[];
end
eliminar=[];
%%%%extra protection% OUT2011
for ii=1:size(m,2)
pa = picant(w(max(1,m(2,ii)-round(messages.setup.wavedet.pictime*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*freq)),2),m(2,ii));
if isempty(pa) || isempty(pp)
eliminar=[eliminar ii]; %#ok<AGROW>
end
end
m(:,eliminar)=[];
signo(eliminar)=[];
% QRS peak detection / wavelet Zero cross detection
time = [];
aux=[]; %Rute26Jun09
if length(signo)>1,
for l = find(signo>0),
if (l==1)
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> %Rute26Jun09
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> %Rute26Jun09
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> %Rute26Jun09
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> %Rute 26Jun09
end
end
end
% Refractary period after a QRS detection (200 ms).
rr = (time(2:end)-time(1:end-1));
ind = find(rr<ceil(messages.setup.wavedet.refrper(end)*freq));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%RUTE 26Jun09
%time(ind+1)=[]; % always the second one is eliminated!!! to be changed
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
time(ind)=[]; %RUTE 27Jun11
aux(ind)=[]; %#ok<NASGU>
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
new = time +i -1;