function [position,positionaux,messages] = wavedet_3D(sig_all, ext_anot, heasig, messages )
%[position,positionaux,messages] =wavedet_3D(sigdir,headir,matdir,ecgnr,ft,anot,lead,t,qrs_flag,anot_fmt,aname,dirann)
%Significant ECG points detector based on wavelet approach
%
%Input Parameters:
%	sigdir: directory of the original signal (.dat file)
%	headir: directory of the header file
%	matdir: output directory
%	ecgnr: string with the name of the record to analyze
% 	ft: format file:
%           - 0 for MIT header using any recorded leads from same datafile
%           - 1 for LUND header (to be tested)
%           - 2 for data in matlab file (.mat file with a varible fa with sampling frequency and variable sinal with ecg, one lead per column)
%           - 4 for RawData.bin files (one hour) created with Mortara Software SuperECG
%
%          obsolete formats
%           - 3 for MIT header using the Dower xyz (ft=0, leadsynth_flag=1)
%           - 10 for MIT format in PTBDB .dat extension file (ft=0)
%           - 30 for MIT type format in PTBDB .dat extension file use the projected xyz (ft=0,  leadsynth_flag=1)
%           - 20 for MIT type format in PTBDB .xyz extension file (ft=0, leads=leads+12)
%           - 40 for MIT header using 3 PC based all available leads and signal (ft=0, leadsynth_flag=4,flagaux=0)
%           - 41 for MIT header using 3 PC based in 8 leads and sub interval on the beat (ft=40,leadsynth_flag=3,flagaux=1)
%           - 42 for MIT type format in PTBDB .dat extension file use the PC based all signal and 8 leads (ft=0,leadsynth_flag=3,flagaux=0)
%           - 44 for MIT header using for MIT type using lead aVF and 2 PC out of precordial leads and all signal (ft=0,leadsynth_flag=5,flagaux=0)
%           - 123 for RawData.bin files (one hour) created with Mortara Software SuperECG (ft=4)
%           - 124 for RawData.bin files (one hour) created with Mortara Software SuperECG use the projected xyz (ft=4, leadsynth_flag=1)
%           - 50 for MIT type format like in ft=0 file using the PC based all available leads and signal (ft=0, leadsynth_flag=4)
%           - 51 for MIT type format like in ft=0 file using the  PC based in 8 leads and sub interval on the beat (ft=0,leadsynth_flag=3,flagaux=1)
%           - 52 for MIT type format like in ft=0 file use the PC based all signal and 8 leads (ft=40,leadsynth_flag=3,flagaux=0)
%
%	anot: name of the annotation output file
%	lead: vector with the leads number to analyze (1,2 or 3 leads)
%         if lead is a cell array the first element stands for leads number to analyze from the ones read and
%         the second for the leads to be read from the original signal file
%	t=[tbegin tend]: time vector with initial and sample to analyze
%   flags: [qrs_flag leadsynth_flag flagaux]
%           qrs_flag: QRS detection only (2), External (1), internal (0) QRS detector (0 by default)
%           leadsynth_flag: None (0), (default)
%                           Dower matrix derived VCG (1)
%                           Levkov's matrix derived VCG (2)
%                           Principal components based in 8 uncorrelated leads out of 12-lead standard (3)
%                           Principal components based in all available leads (4)
%                           lead aVF and Principal components based on precordial leads (5)
%           flagaux: use all signal (0) (default) or sub interval on the beat (1) for PC calculation, available only for leadsynth_flag>2
%
%   anot_fmt: cell with format in case of external annotation file
%             first field for QRS detector, format of annotation file: MIT (0) or mat file (1)
%             second field for with excluded segments: xls (0), MIT annotations (1)
%   aname: cell with annotation file name with QRS marks in the first field (by default ecgnr)
%               with annotation file name with excluded segments in the second field (by default aname{1})
%   dirann: cell with directory with external annotation files (QRS marks and/or excluded segments, by default sigdir)
%          if diferent annotation files are used, first field refers to QRS marks and second field to excluded segments
%
%Output Parameters:
%	position: struct vector with the detected points locations in samples
%
%   positionaux: struct vector with the detected points using each specified lead
%                   - positionaux.position1: using lead(1)
%                   - positionaux.position2: using lead(2)
%                   - positionaux.position3: using lead(3)
%
%  messages: struct vector with errors, warnings, setup and state
%                   - messages.status=0,1;
%                   - messages.errors
%                   - messages.errors_desc
%                   - messages.warnings
%                   - messages.setup.wavedet.nsamp=nsamp; number of samples to process in each excerpt, consider to calculate some parameters (default 2^16/250*sf)
%                   - messages.setup.wavedet.sig_quality: sig_quality test 1 or 0 (activated or inactivated)
%
% struct description:
%             Pon: P wave onset
%              P: P wave peak
%           Poff: P wave end
%          QRSon: QRS complex onset
%              Q: Q wave peak (in multilead, according to QRSonset best fitted lead)
%      R_inQRSon: R wave peak in QRSonset best fitted lead (multilead only)
%              R: R wave peak (median mark in multilead approach)
%            qrs: QRS complex fiducial mark
%     R_inQRSoff: R wave peak in QRSend best fitted lead (multilead only)
%         Rprima: R' wave peak (in multilead, according to QRSend best fitted lead)
%              S: S wave peak (in multilead, according to QRSend best fitted lead)
%         QRSoff: QRS complex end
%            Ton: T wave onset
%              T: first T wave peak (median mark in multilead approach)
%         Tprima: second T wave peak (in multilead, according to Tend best fitted lead)
%           Toff: T wave end
%          Ttipo: T wave morphology (1:5 corresponding respectively to normal, inverted, upwards only, downwards only, biphasc pos-neg, biphasc neg-pos)
%         Tscale: scale used for T wave detection
%        Ttipoon: T wave morphology, according to Tonset best fitted lead (multilead only)
%       Ttipooff: T wave morphology, according to Tend best fitted lead (multilead only)
%
% warnings description:
%       'warning: there are not 8 uncorrelated leads, all available leads considered instead'
%       leadsynth_flag=3 was used over a file that does not have all precordial leads and at least 2 of the limb leads
%       leadsynth_flag=4 was automatically set
%
%       'warning: lead aVF not available, lead xpto considered instead'
%       leadsynth_flag=5 was used over a file that does not have a lead aVF nor 2 limb leads to calculate it
%       aVF was replaced by lead xpto
%
%       'warning: mV are assumed as units'
%       not units are provided, mV are assumed
%
% The analysis loop has an overlapping structure.  The number of samples
% processed is large and there are an overlap at the beginning and another
% at the end.  The first one is necessary because: a) the first l4-1 samples
% are always incorrectly filtered b) To align the signals, as the filters
% have different lengths, we must discard some samples. c) The algorithms
% must have some possibility of turning back.
% The second overlap is mainly because of the alignment, as the filters have
% different delays.
%
% based on wavedet.m by Juan Pablo Mart�nez Cort�s
% previous versions include wavedetplus.m and wavedet3D.m
% Last update: Rute Almeida  07FEB2012
%
% Designed for MATLAB Version R12; tested with MATLAB Version R13

ft = 2;
str = '1';
positionaux=[];
position=[];
aname_sig_quality_cell{3} = [];
t_initial=0;  % Juan 28/03/2012
lead = 1;

% if (nargin<14) % JB 27/07/2011
%     %     bufferedSignal = [];  % Juan 28/03/2012  No se usan no tienen porque
%     %     estar definidas
%     %     bufferedHeaSig = [];  % Juan 28/03/2012
% else % 03AGO 2011
%     ft=NaN;
%     heasig = bufferedHeaSig;
% %     if ~isfield(heasig,'freq')
% %         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 found on the header information file.'}];
% %         messages.status=0;
% %         return
% %     end
% %     if ~isempty(t) && length(t)==2 %& t(2)>t(1)
% %         if length(bufferedSignal) == t(2)-t(1)+1;
% %             sig_all =bufferedSignal(:,~isnan(bufferedSignal(1,:)));
% %             %             t_initial=t(1);  % Juan 28/03/2012
% %             t=[1 length(sig_all)];
% %         elseif length(bufferedSignal)>t(2)
% %             sig_all = bufferedSignal(:,t(1):t(2));
% %         else
% %             sig_all = bufferedSignal;
% %             t=[1 length(sig_all)];
% %         end
% %     else
% %         sig_all = bufferedSignal;
% %         t=[1 length(sig_all)];
% %     end
% %     
%     
%     
%     if size(sig_all,2)>size(sig_all,1)
%         sig_all=sig_all';
%     end
%     
%     if ~isfield(heasig,'nsamp') || isempty(heasig.nsamp)
%         heasig.nsamp=size(sig_all,1);
%     end
%     if ~isfield(heasig,'nsig')
%         heasig.nsig=size(sig_all,2);
%     end
%     if ~isfield(heasig,'gain')
%         heasig.gain=200.*ones(1,heasig.nsig);  % When heasig.gain = 0 => default 200
%     end
%     if ~isfield(heasig,'spf')
%         %         heasig.spf=ones(1,max(heasig.nsig,leadreading));  % Juan
%         %         28/03/2012  leadreading no esta definido
%         heasig.spf=ones(1,max(heasig.nsig, size(sig_all,1)));
%     end
%     if ~isfield(heasig,'adczero')
%         heasig.adczero=zeros(1,heasig.nsig);
%     end
% end

t=[1 size(sig_all,1)];

%%%%%%%%%%%%%%%%%%%%%%%

if nargin < 4
    messages.warnings=[];
end

if ~isfield(messages,'setup'), messages.setup.wavedet=[]; end
if ~isfield(messages.setup.wavedet,'QRS_detection_only'), messages.setup.wavedet.QRS_detection_only=false; end
if ~isfield(messages.setup.wavedet,'QRS_detection_thr'), messages.setup.wavedet.QRS_detection_thr=ones(5,1); end
if ~isfield(messages,'errors'), messages.errors=[]; end
if ~isfield(messages,'errors_desc'), messages.errors_desc=[]; end
if ~isfield(messages,'warnings'), messages.warnings=[]; end
if isfield(messages,'status')
    if messages.status~=1
        messages.warnings=[messages.warnings {['Initial status=' num2str(messages.status) '. status changed to 1']}];
    end
end
messages.status=1;
messages.setup.wavedet.sig_quality = 2;

% if nargin<7 %21MAR2011
%     messages.errors=[messages.errors {'Fatal error in wavedet_3D: not enough inputs.'}];
%     warning(char(messages.errors(end)))
%     messages.errors_desc=[messages.errors_desc 'Mandatory inputs not defined.'];
%     messages.status=0;
%     return
% elseif nargin<9
%     flags=[];
% end
% if iscell(lead)
%     if size(lead,2)==2
%         leadreading=(lead{2});
%     else
%         leadreading=(lead{1});
%     end
%     lead=lead{1};
% else
%     leadreading=(lead);
% end

% if nargin<=10 %&& anot_fmt==1 %JUL2011
%     aname=ecgnr;
%     dirann=sigdir;
% end%JUL2011
% if iscell(aname)
%     if length(aname)==2
%         aname_sig_quality=aname{2};
%     end
%     aname=aname{1};
% else
%     aname_sig_quality=aname;
% end
% if iscell(dirann)
%     if length(dirann)==2
%         dirann_sig_quality=dirann{2};
%     end
%     dirann=dirann{1};
% else
%     dirann_sig_quality=dirann;
% end

% if isempty(flags)
%     if isfield(messages.setup.wavedet,'qrs_flag')
%         flags(1)=messages.setup.wavedet.qrs_flag;
%     else
%         flags(1)=0;
%     end
% elseif isfield(messages.setup.wavedet,'qrs_flag') &&  flags(1)~=messages.setup.wavedet.qrs_flag;
%     messages.warnings=[messages.warnings {'Obsolete use of flag with inconsistency: flag(1) value differs from messages.setup.wavedet.qrs_flag value. flag(1) ignored.'}];
% else
%     messages.setup.wavedet.qrs_flag=flags(1);
%     messages.warnings=[messages.warnings {['Obsolete use of flag. Use messages.setup.wavedet.qrs_flag=' num2str(flags(1)) 'instead.' ]}];
% end
% if length(flags)<2
%     if isfield(messages.setup.wavedet,'leadsynth_flag')
%         flags(2)=messages.setup.wavedet.leadsynth_flag;
%     else
%         flags(2)=0;
%     end
% elseif isfield(messages.setup.wavedet,'leadsynth_flag') &&  flags(2)~=messages.setup.wavedet.leadsynth_flag;
%     messages.warnings=[messages.warnings {'Obsolete use of flag with inconsistency: flag(2) value differs from messages.setup.wavedet.leadsynth_flag value. flag(2) ignored.'}];
% else
%     messages.setup.wavedet.leadsynth_flag=flags(2);
%     messages.warnings=[messages.warnings {['Obsolete use of flag. Use messages.setup.wavedet.leadsynth_flag=' num2str(flags(2)) 'instead.' ]}];
% end
% if length(flags)<3
%     if isfield(messages.setup.wavedet,'flagaux')
%         flags(3)=messages.setup.wavedet.flagaux;
%     else
%         flags(3)=0;
%     end
% elseif isfield(messages.setup.wavedet,'flagaux') &&  flags(3)~=messages.setup.wavedet.flagaux;
%     messages.warnings=[messages.warnings {'Obsolete use of flag with inconsistency: flag(3) value differs from messages.setup.wavedet.flagaux. flag(3) ignored.'}];
% else
%     messages.setup.wavedet.flagaux=flags(3);
%     messages.warnings=[messages.warnings {['Obsolete use of flag. Use messages.setup.wavedet.flagaux=' num2str(flags(3)) 'instead.' ]}];
% end
% if length(flags)==4
%     if isfield(messages.setup.wavedet,'nsamp') && messages.setup.wavedet.nsamp~=flags(4)
%         messages.warnings=[messages.warnings {'Obsolete use of flag with inconsistency: flag(4) value differs from messages.setup.wavedet.nsamp value. flag(4) ignored.'}];
%     else
%         messages.setup.wavedet.nsamp=flags(4);
%         messages.warnings=[messages.warnings {['Obsolete use of flag. Use messages.setup.wavedet.nsamp=' num2str(flags(4)) 'instead.' ]}];
%     end
% end


if( nargin < 2 || isempty(ext_anot) )
    qrs_flag=0;
else
    qrs_flag=1;
end

leadsynth_flag=0;
flagaux=0;

flags = [ qrs_flag leadsynth_flag flagaux];

messages.setup.wavedet.qrs_flag=qrs_flag;
messages.setup.wavedet.leadsynth_flag=leadsynth_flag;
messages.setup.wavedet.flagaux=flagaux;

% if ~isfield(messages.setup.wavedet,'sig_quality') || messages.setup.wavedet.sig_quality==0
%     messages.setup.wavedet.sig_quality=0;
%     messages.warnings=[messages.warnings {'No signal quality test used: messages.setup.wavedet.sig_quality=0'}];
% else
%     if isfield(messages.setup.wavedet,'dirann_sig_quality')
%         dirann_sig_quality=messages.setup.wavedet.dirann_sig_quality;
%     elseif ~exist('dirann_sig_quality','var') || isempty(dirann_sig_quality)
%         dirann_sig_quality=sigdir;
%     end
%     if isfield(messages.setup.wavedet,'aname_sig_quality')
%         aname_sig_quality=messages.setup.wavedet.aname_sig_quality;
%     elseif ~exist('aname_sig_quality','var') || isempty(aname_sig_quality)
%         aname_sig_quality=ecgnr;
%     end
%     messages.setup.wavedet.sig_quality=1;
%     messages.warnings=[messages.warnings {['Signal quality information used: messages.setup.wavedet.sig_quality='  num2str(messages.setup.wavedet.sig_quality)]}];
% end

if length(flags)==3
    if  leadsynth_flag<3
        flagaux=0;
        messages.warnings=[messages.warnings {'flagaux>0 only available for leadsynth_flag>2. flagaux=0 used'}];
    end
end

%obsolete formats
% if ft==3;
%     ft=0; leadsynth_flag=1;
% end %obsolete formats
% if ft==10
%     ft=0;
% end %obsolete formats
% if ft==20
%     ft=0;
%     lead=lead+12;
%     leadreading=leadreading+12;
% end %obsolete format
% if ft==30
%     ft=0;
%     leadsynth_flag=1;
%     leadreading=1:12;
% end%obsolete formats
% if ft==41
%     ft=0;
%     leadsynth_flag=3;
%     flagaux=1;
% end %obsolete format
% if ft==42||ft==52
%     ft=0;
%     leadsynth_flag=3;
%     flagaux=0;
% end%obsolete formats
% if ft==50
%     ft=0;
%     leadsynth_flag=4;
%     flagaux=0;
% end%obsolete formats
% if ft==51
%     ft=0;
%     leadsynth_flag=4;
%     flagaux=1;
% end%obsolete formats
% if ft==123
%     ft=4;
% end% change format notation
% if ft==124
%     ft=4;
%     leadsynth_flag=1;
% end% change format notation
% if ft==44
%     ft=0;
%     leadsynth_flag=5;
%     flagaux=0;
% end%obsolete formats
% if ft==40
%     ft=0;
%     leadsynth_flag=4;
%     flagaux=0;
% end%obsolete formats
% 
% %number of leads from 1 to 3
% if length(lead)>3
%     lead(4:end)=[];
%     messages.warnings=[messages.warnings {['Too many leads to process, just first 3 leads are considered: lead=[' num2str(lead) ']']}];
% elseif isempty(lead)
%     lead=1;
%     leadreading=lead;
%     messages.warnings=[messages.warnings {'No leads indicated to process: first lead in file is considered'}];
% end
% 
% if leadsynth_flag==1
%     leadstring='synt';
% end
% if leadsynth_flag==2
%     leadstring='levk';
% end
% if leadsynth_flag==1 || leadsynth_flag==2
%     %%%%%%%%%%%%%%%%%%%%%%%
%     %verify leads to process
%     invalid_lead= find(lead>3); %#ok<EFIND>
%     if ~isempty(invalid_lead)
%         [C,IA,IB]=intersect(leadreading,lead); %#ok<NASGU> % Juan 28/03/2012
%         if length(C)==length(lead)
%             lead=IA;
%             invalid_lead= find(lead>3); % Juan 28/03/2012
%             if isempty(invalid_lead)
%                 messages.warnings=[messages.warnings {['the lead(s) chosen to be processed are the lead(s) ' num2str(leadreading(lead))]}];
%             else
%                 alead=lead;
%                 lead(invalid_lead)=[];
%                 if isempty(lead)
%                     lead=1:min(length(lead),3);
%                 end
%                 messages.warnings=[messages.warnings {[ 'Invalid lead for VCG analysis (' num2str(alead) '): transformed lead(s) ' num2str(1:min(length(lead),3)) ' considered instead.' ]}];
%             end
%         elseif ~isempty(C)
%             messages.warnings=[messages.warnings {[ 'Invalid lead for VCG analysis (' num2str(lead) '): transformed lead(s) ' num2str(leadreading(IA)) ' considered instead.' ]}];
%             lead=IA;
%         else
%             messages.warnings=[messages.warnings {[ 'Invalid lead for VCG analysis (' num2str(lead) '): transformed lead(s) ' num2str(1:min(length(lead),3)) ' considered instead.' ]}];
%             lead=1:min(length(lead),3);
%         end
%     end
%     %%%%%%%%%%%%%%%%
% elseif leadsynth_flag==3
%     invalid_lead=find(lead>8);
%     lead(invalid_lead)=[];
%     if ~isempty(invalid_lead)
%         messages.warnings=[messages.warnings {'Invalid lead after the chosen PC transformation: 8 PC are produced.'}];
%     end
% elseif leadsynth_flag==4
%     invalid_lead= find(lead>max(leadreading)); % Juan 28/03/2012
%     lead(invalid_lead)=[];
%     if ~isempty(invalid_lead)
%         messages.warnings=[messages.warnings {['Invalid lead after the chosen PC transformation: ' num2str(leadreading) ' PC are produced. Lead(s) ' num2str(lead)  ' considered instead']}];
%     end
%     if isempty(lead)
%         lead=1;
%         messages.warnings=[messages.warnings {['Invalid lead after the chosen PC transformation: lead=[' num2str(lead) '] considered instead.']}];
%     end
%     
% elseif leadsynth_flag==5
%     if isempty(find(lead==1,1))
%         messages.warnings=[messages.warnings { 'Processing lead 1, standing for aVF, is mandatory in the chosen PC transformation.'}];
%         lead=[1 lead];
%     end
%     invalid_lead=find(lead>7);
%     lead(invalid_lead)=[];
%     if ~isempty(invalid_lead)
%         messages.warnings=[messages.warnings {'Invalid lead after the chosen PC transformation: aVF plus 6 PC are available.'}];
%     end
%     if isempty(lead)
%         lead=1:3;
%         messages.warnings=[messages.warnings {['Invalid lead after the chosen PC transformation: lead=[' num2str(lead) '] considered instead.']}];
%     elseif length(lead)>3
%         lead(4:end)=[];
%         messages.warnings=[messages.warnings {['Too many leads to process, just first 3 leads are considered: lead=[' num2str(lead) ']']}];
%     elseif length(lead)==1
%         messages.warnings=[messages.warnings { 'Single lead over lead aVF selected.' }];
%     end
% end

% global OPT
allfig=0;figuresoff=1; %#ok<NASGU>

% OPT=optimset(@fminbnd);
% OPT.MaxFunEvals=1000;
% OPT.Display='off';
% warning off %#ok<WNOFF>
ultimo_anot=1;
ultimo_anottyp=0;
% timeoffset=0;   % Juan 28/03/2012  Solo se define para el formato lund si
% es necesario dentro
count=0;

%errores=struct('errores',[]);

sep = filesep;

% if isunix,          %%%%%%% Rute 13/03/02
%     sep = '/';
% else sep = '\';     %%%%%%% Rute 13/03/02
% end                 %%%%%%% Rute 13/03/02

PC_aux_ind=[];
% aux=find(ecgnr==sep);
% headir=[headir ecgnr(1:aux)];
% sigdir=[sigdir ecgnr(1:aux)];
% matdir=[matdir ecgnr(1:aux)];
% if ~isempty(aux)
%     ecgnr=ecgnr(aux+1:end);
% end

%%%%%%%%%%%%%%%%%%%%%%% file format %%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if ft==0   % MIT format files
%     name=[headir sep ecgnr '.hea'];
%     heasig = readheader(name);
%     [auxleads] = ecgleads(heasig.desc); %V1 V2 V3 V4 V5 V6 I II III aVR aVL aVF X Y Z VX VY VZ
%     
%     %verify leads to read
%     if sum(leadreading>heasig.nsig)>0
%         messages.warnings=[messages.warnings {[ 'The selected file contains only ' num2str(heasig.nsig) ' signals: unable to read leads [' num2str(leadreading) '], leads out of the range will be ignored' ]}];
%         leadreading(leadreading>heasig.nsig)=[];
%     end
%     
%     if leadsynth_flag==0 % should be leadreading=lead
%         if length(leadreading)~=length(lead) | (length(leadreading)==length(lead) & sort(leadreading)~=sort(lead)) %#ok<AND2,OR2>
%             [C,IA,IB]=intersect(leadreading,lead); %#ok<NASGU>
%             if length(C)==length(lead) % extra leads on leadreading will be deleted
%                 leadreading=IA;
%                 messages.warnings=[messages.warnings {['the lead(s) chosen to be processed are the lead(s) ' num2str(leadreading)]}];
%             else % extra leads on lead or different leads
%                 if ~isempty(lead>heasig.nsig)
%                     messages.warnings=[messages.warnings {[ 'The selected file contains only ' num2str(heasig.nsig) ' signals: unable to read leads [' num2str(lead) '], leads out of the range will be ignored' ]}];
%                     lead=lead(lead<=heasig.nsig);
%                 end
%                 if length(leadreading)~=length(lead) | (length(leadreading)==length(lead) & sort(leadreading)~=sort(lead)) %#ok<AND2,OR2>
%                     messages.warnings=[messages.warnings {[ 'unable to interpret the lead(s) number to process (' num2str(leadreading) '): lead(s) ' num2str(lead(lead<=heasig.nsig)) ' considered instead.' ]}];
%                     leadreading=lead;
%                 end
%             end
%         else
%             messages.warnings=[messages.warnings {['the lead(s) chosen to be processed are the lead(s) ' num2str(leadreading)]}];
%         end
%     elseif leadsynth_flag==1 || leadsynth_flag==2 || leadsynth_flag==3% Leads V1 to V6 and 2 limb leads required
%         if size(leadreading)<8 %wrong number of leads
%             if ~isnan(auxleads(1:8))
%                 leadreading=(auxleads(1:8));
%             elseif ~isnan(auxleads([1:7 9]))
%                 leadreading=(auxleads([1:7 9]));
%             elseif ~isnan(auxleads([1:6 8 9]))
%                 leadreading=(auxleads([1:6 8 9]));
%             else
%                 messages.errors=[messages.errors {'Fatal error in wavedet_3D: not enough leads in file to apply the chosen transformation.'}];
%                 warning(char(messages.errors(end)))
%                 messages.errors_desc=[messages.errors_desc 'Leads V1 to V6 and 2 limb leads are required.'];
%                 messages.status=0;
%                 return
%             end
%             messages.warnings=[messages.warnings {['Not enough leads to apply the chosen transformation in given leadreading: leadreading=['  num2str(leadreading) '] considered instead' ]}];
%         else
%             if sum(isnan(auxleads(1:6)))>0 || sum(isnan(auxleads(7:9)))>1 % heasig.desc does not include Leads V1 to V6 and 2 limb leads
%                 if heasig.nsig>=8
%                     messages.warnings=[messages.warnings {'Leads V1 to V6 and 2 limb leads required to apply the chosen transformation but cannot be recognized. It is assumed that the first 8 leads of leadreading are V1 V2 V3 V4 V5 V6 I II'}];
%                     leadreading=leadreading(1:8);
%                     auxleads=NaN*ones(size(auxleads));
%                     auxleads(1:8)=1:8;
%                 else
%                     messages.errors=[messages.errors {'Fatal error in wavedet_3D: not enough leads in file to apply the chosen transformation.'}];
%                     warning(char(messages.errors(end)))
%                     messages.errors_desc=[messages.errors_desc 'Leads V1 to V6 and 2 limb leads are required.'];
%                     messages.status=0;
%                     return
%                 end
%             else % heasig.desc includes Leads V1 to V6 and 2 limb leads
%                 leadreading=leadreading(auxleads(auxleads<=length(leadreading)));% sort leadreding as V1 V2 V3 V4 V5 V6 I II III aVR aVL aVF X Y Z VX VY VZ
%             end
%         end
%     elseif leadsynth_flag==4 && size(leadreading)==1
%         leadsynth_flag=0;
%         messages.warnings=[messages.warnings {'To aply PC transformation more than one lead is requires in leadreading: no trasnformation applied (leadsynth_flag=0)'}];
%     elseif leadsynth_flag==5 && length(lead)~=1 %  aVL and leads V1 to V6 required
%         if size(leadreading)<7 %wrong number of leads
%             if ~isnan(auxleads([12 1:6])) %lead aVF and V1 to V6 available
%                 leadreading=auxleads([1:6 12]);
%             elseif ~isnan(auxleads(1:8))
%                 leadreading=auxleads(1:8);
%             elseif ~isnan(auxleads([1:7 9]))
%                 leadreading=auxleads([1:7 9]);
%             elseif ~isnan(auxleads([1:6 8 9]))
%                 leadreading=auxleads([1:6 8 9]);
%             else
%                 messages.errors=[messages.errors {'Fatal error in wavedet_3D: not enough leads in file to apply the chosen transformation.'}];
%                 warning(char(messages.errors(end)))
%                 messages.errors_desc=[messages.errors_desc 'Leads aVF and V1 to V6 are required. Sugestion: use leadsynth_flag=4 instead'];
%                 messages.status=0;
%                 return
%             end
%         else
%             if sum(isnan(auxleads(1:6)))>0 || (isnan(auxleads(12)) && sum(isnan(auxleads(7:9)))>1) % heasig.desc does not include Leads V1 to V6 and aVF cannot be calculated
%                 if size(leadreading)==7 &&  heasig.nsig>=7
%                     messages.warnings=[messages.warnings {'Leads V1 to V6 and aVF limb leads required to apply the chosen transformation but cannot be recognized. It is assumed that the 7 leads of leadreading are V1 V2 V3 V4 V5 V6 aVF' }];
%                     leadreading=leadreading(1:7);
%                     auxleads=NaN*ones(size(auxleads));
%                     auxleads([1:6 12])=1:7;
%                 elseif heasig.nsig>=8
%                     messages.warnings=[messages.warnings {'Leads V1 to V6 and aVF limb leads required to apply the chosen transformation but cannot be recognized. It is assumed that the first 8 leads of leadreading are V1 V2 V3 V4 V5 V6 aVF' }];
%                     leadreading=leadreading(1:8);
%                     auxleads=NaN*ones(size(auxleads));
%                     auxleads(1:8)=1:8;
%                 else
%                     messages.errors=[messages.errors {'Fatal error in wavedet_3D: not enough leads in file to apply the chosen transformation.'}];
%                     warning(char(messages.errors(end)))
%                     messages.errors_desc=[messages.errors_desc 'Leads aVF and V1 to V6 are required. Sugestion: use leadsynth_flag=4 instead'];
%                     messages.status=0;
%                     return
%                 end
%             else % heasig.desc includes Leads V1 to V6 and 2 limb leads
%                 leadreading=leadreading(auxleads(auxleads<=length(leadreading)));% sort leadreding as %V1 V2 V3 V4 V5 V6 I II III aVR aVL aVF X Y Z VX VY VZ
%             end
%         end
%     elseif leadsynth_flag==5 && length(lead)==1
%         if ~isnan(auxleads(12))
%             leadreading=auxleads(12);
%         elseif sum(~isnan(auxleads(7:9)))>1
%             leadreading=auxleads(~isnan(auxleads(7:9)));
%         elseif length(leadreading)==1
%             messages.warnings=[messages.warnings {'Lead aVF lead is required but cannot be recognized. It is assumed that the first lead of leadreading is aVF'}];
%             leadreading=1;
%             auxleads=NaN*ones(size(auxleads));
%             auxleads(12)=1;
%         else
%             messages.warnings=[messages.warnings {'Lead aVF lead is required but cannot be recognized. It is assumed that the first 2 lead of leadreading are I and II'}];
%             auxleads=NaN*ones(size(auxleads));
%             auxleads(7:8)=1:2;
%         end
%     end
%     
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% this is necessary when there are leads in different files.   e.g. ptbdb database  JP 2006
%     leadsfile=NaN*ones(heasig.nsig,1);
%     leadinfile=NaN*ones(heasig.nsig,1);
%     nsiginfile=NaN*ones(heasig.nsig,1);
%     for jj=1:heasig.nsig
%         for kk=1:heasig.nsig
%             leadsfile(kk)=strcmp(heasig.fname(jj,:),heasig.fname(kk,:));
%         end
%         leadsfile = cumsum(leadsfile);
%         leadinfile(jj) = leadsfile(jj);  % lead indexes in the file
%         nsiginfile(jj) = leadsfile(end); % number of the leads in th file nsiginfile<=heasig.nsig
%     end
%     
%     if(min(lead)>nsiginfile(1))
%         auxleads(auxleads<=nsiginfile(1))=NaN;
%     else
%         auxleads(auxleads>nsiginfile(1))=NaN;
%     end
%     n_auxleads=find(auxleads>0);
%     auxleads(n_auxleads)=leadinfile(auxleads(n_auxleads));
%     %     nsiginfile=nsiginfile(lead(1)); % Juan 28/03/2012
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     %%%%%%%%%%%%%%%%%%%%%%%
%     if (heasig.fmt(lead)==16) | (heasig.fmt(lead)==212) | (heasig.fmt(lead)==61), %#ok<OR2>
%         %%%%%%% multilead%%% Rute 02.Dec.04
%         if size(lead,2)==1
%             formato = num2str(heasig.fmt(lead));
%         elseif   diff(heasig.fmt(lead))==0
%             formato = num2str(heasig.fmt(lead(1)));
%         end
%         %%%%%%%%%%%%%
%     else
%         messages.errors=[messages.errors {'This format is not supported by the program.'}];
%         warning(char(messages.errors(end)))
%         messages.errors_desc=[messages.errors_desc ['Format ' num2str(heasig.fmt(lead)) ' data. For MIT data, formats 16, 212 and 61 are supported.'] ];
%         messages.status=0;
%         return
%     end
%     if strcmp(formato,'16')||strcmp(formato,'61'),
%         fid = fopen([sigdir heasig.fname(lead(1),:)],'rb');
%         fseek(fid,0,-1);  % Rewind the file
%         if strcmp(heasig.fname(1,1),'_'),  % Siemens card recordings with MIT-type header
%             timeoffset=512;
%             fseek (fid, timeoffset,-1); % offset
%         end
%     end
%     
% elseif ft==1  % LUND format files
%     try
%         hdsig=gethdsig(sigdir,ecgnr); % Reading header information
%         heasig = hdsig2heasig(hdsig);
%         leadinfile=lead;
%         
%         [auxleads] = ecgleads(heasig.desc);
%     catch me
%         messages.errors=[messages.errors {'Fatal error in wavedet_3D: unable to read header signal on Lund format.'}];
%         warning(char(messages.errors(end)))
%         messages.errors_desc=[messages.errors_desc {me.message}]; % Juan 28/03/2012
%         messages.status=0;
%         return
%     end
% elseif ft==2  % data in a Matlab file
%     try
%         [sig_all,messages1,heasig] = readsignal(sigdir,ecgnr,t,ft,[],leadreading,messages);  % Juan 28/03/2012
%         if messages1.status==0
%             error('Fatal error in wavedet_3D: unable to read signal on mat format.');  % Juan 28/03/2012
%         end
%     catch me
%         messages.errors=[messages.errors {messages1.errors}];  % Juan 28/03/2012
%         warning(char(messages.errors(end)))
%         messages.errors_desc=[messages.errors_desc {me.message}];  % Juan 28/03/2012
%         messages.status=0;
%         return
%     end
%     if size(sig_all,2)>size(sig_all,1)
%         sig_all=sig_all';
%     end
%     %should not be needed;   I will remove that part, all of that is done in
%     %readsignal  Juan 28/03/2012
%     if ~exist('heasig','var')
%         heasig.nsamp=size(sig_all,1);
%         try
%             if exist('fa','var')
%                 heasig.freq=fa; %#ok<NODEF>
%             else
%                 heasig.freq=messages.setup.wavedet.freq;
%             end
%         catch me
%             messages.errors=[messages.errors {'Fatal error in wavedet_3D: no sampling frequency variable found.'}];
%             warning(char(messages.errors(end)))
%             me.message = 'No sampling frequency variable found on the mat file.';  % Juan 28/03/2012
%             messages.errors_desc=[messages.errors_desc {me.message}];  % Juan 28/03/2012
%             messages.status=0;
%             return
%         end
%         heasig.nsig=size(sig_all,2);
%         heasig.gain=200.*ones(1,heasig.nsig);  % When heasig.gain = 0 => default 200
%         heasig.adczero=zeros(1,heasig.nsig);
%     end
%     %%%%%%%%%%%%%%%%%%%%%%%
%     
%     leadinfile=lead;
%     %     nsiginfile=heasig.nsig;  %Juan 28/03/2012
%     heasig.spf_ecg=1;
%     
%     if isfield(heasig,'desc')
%         [auxleads] = ecgleads(heasig.desc);
%     else
%         auxleads = 1:heasig.nsig; % assuming that all are ECG
%     end
%     
% elseif ft==4  % Ana e Tiago 11-04-2006 (copy from wavedetplus)
%     fa=1000;
%     %     aux=[sigdir ecgnr '.bin'];  Juan 28/03/2012
%     if nargin >=8,
%         if isempty(t)
%             t=[1 fa*60*60];
%             %         else
%             %             t=[t(1) min(t(1)+fa*60*60,t(2))];
%         end
%     else
%         t=[1 fa*60*60];
%     end
%     %     first_sample=t(1); %ver com Sonia se '1 -1 pq  first_sample=0 por
%     %     defeito!!!!  Juan 28/03/2012
%     %     n_sample_per_lead=t(end);Juan 28/03/2012
%     %     lead_number=lead;Juan 28/03/2012
%     heasig.nsamp=t(2)-t(1);
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%
%     heasig.freq=fa;
%     heasig.gain=200;
%     heasig.nsig=12;
%     heasig.units=char(ones(heasig.nsig,1)*'microV');
%     heasig.desc=['I  ';'II ';'V1 ';'V2 ';'V3 ';'V4 ';'V5 ';'V6 ';'III';'aVR';'aVL';'aVF'];
%     [auxleads] = ecgleads(heasig.desc);
%     leadinfile=lead;
%     %     nsiginfile=heasig.nsig;
%     
% elseif isnan(ft) % data in buffer % 03AGO 2011
%     leadinfile=lead;
%     %     nsiginfile=heasig.nsig;           Juan 28/03/2012
%     heasig.spf_ecg=1;
%     if ~isfield(heasig,'nsamp') || heasig.nsamp<=0
%         heasig.nsamp=size(sig_all,1);
%     end
%     if isfield(heasig,'desc')
%         [auxleads] = ecgleads(heasig.desc);
%     else
%         auxleads = 1:heasig.nsig; % assuming that all are ECG
%     end
% else
%     try
%         [sig_all,messages1,heasig] = readsignal(sigdir,ecgnr,t,ft,[],leadreading,messages);
%         if messages1.status==0
%             error('Fatal error in wavedet_3D: no admissible format.');  % Juan 28/03/2012
%         end
%     catch me
%         messages.errors=[messages.errors {messages1.errors}];
%         me.message = 'Format ft should be one of {0,1,2,4}. See help wavedet3D for details.'; % Juan 28/03/2012
%         warning(char(messages.errors(end)))
%         messages.errors_desc=[messages.errors_desc {me.message}];  % Juan 28/03/2012
%         messages.status=0;
%         return
%     end
%     if size(sig_all,2)>size(sig_all,1)
%         sig_all=sig_all';
%     end
% end
% leadinfile=lead;  % Juan 28/03/2012  definido dentro de cada formato
% nsiginfile=heasig.nsig;  % Juan 28/03/2012 parece que no se usa m�s no es
% necesario definir
if ~isfield(heasig,'spf')
    heasig.spf=ones(1,heasig.nsig);
end
% try
%     heasig.spf_ecg=heasig.spf(leadinfile(1));
% catch me
%     me.message = 'It is taken the first lead in file.';
%     messages.warnings = [messages.warnings {me.message}];
%     heasig.spf_ecg=heasig.spf(1);
% end
heasig.spf_ecg = 1;
messages.setup.wavedet.freq=heasig.spf_ecg*heasig.freq; %13ENE2010
% if isfield(heasig,'desc')
%     [auxleads] = ecgleads(heasig.desc);
% else
%     auxleads = 1:heasig.nsig; % assuming that all are ECG
% end
% at this point heasig must exist
if isfield(messages.setup.wavedet,'nsamp')
    nsamp = messages.setup.wavedet.nsamp; % Number of samples per excerpt for the WT as input
else
    nsamp=0;
end
if nsamp==0
    nsamp = 2^16/250*messages.setup.wavedet.freq; % Number of samples per excerpt for the WT corresponding to 2^16 samples at sf=250
end
messages.setup.wavedet.nsamp=nsamp;
% if nargin >=8,
%     if isempty(t)
%         t=[1 heasig.nsamp]; % pode falhar em formato lund...
%     else
%         t(1)=max(t(1),1);
%         if ~isempty(heasig.nsamp)
%             t(2)=min(t(1)+heasig.nsamp-1,t(2));
%         end
%     end
% else
%     t=[1 heasig.nsamp];
% end
t=[1 heasig.nsamp];

% if flagaux==1
%     if isfield(messages.setup.wavedet,'leadsynth')
%         leadref= messages.setup.wavedet.leadsynth;
%         messages.warnings=[messages.warnings {'Lead' num2str(leadref) ' will be considered for PCA time interval.'}];
%     else
%         if ~isnan(auxleads(8)); %lead II
%             leadref=auxleads(8);
%         else
%             messages.warnings=[messages.warnings {'Lead II not available. First selected lead would be considered instead.'}];
%             warning(char(messages.warnings(end)));
%             leadref=lead(1);
%         end
%     end
%     if ~exist([matdir ecgnr '.s' num2str(leadref)],'file')
%         wavedet_3D(sigdir,headir,matdir,ecgnr,ft,['s' num2str(leadref)],leadref);
%     end
%     marks=readannot([matdir ecgnr '.s' num2str(leadref)]);
%     if exist('marks','var')
%         marks=posmat2(marks);
%         marks(marks==0)=NaN;
%     end
%     PC_aux_ind= [];
%     tol1=0;tol2=0; tol3=200;
%     
%     mark_on=4;mark_off=10;mark_on2=5; mark_off2=6;   % QRSonset-tend
%     %         mark_on=7;mark_off=10;mark_on2=6;mark_off2=10; % ton-tend
%     %         mark_on=8;mark_off=10;mark_on2=6;mark_off2=10;% tpeak-tend
%     %         mark_on=4;mark_off=6;mark_on2=4;mark_off2=6;% QRS based
%     
%     for tt=1:size(marks,1)
%         if ~isnan(marks(tt,mark_off)) & ~isnan(marks(tt,mark_on))%#ok<AND2> %
%             PC_aux_ind=[PC_aux_ind round((marks(tt,mark_on)-tol1)):round((marks(tt,mark_off)+tol2))]; %#ok<AGROW>
%         elseif ~isnan(marks(tt,mark_on))& ~isnan(marks(tt,mark_off2))%#ok<AND2> %
%             PC_aux_ind=[PC_aux_ind round((marks(tt,mark_on)-tol1)):(marks(tt,mark_off2)+tol3)]; %#ok<AGROW>
%         elseif ~isnan(marks(tt,mark_on2))& ~isnan(marks(tt,mark_off))%#ok<AND2> %
%             PC_aux_ind=[PC_aux_ind round((marks(tt,mark_on2)-tol3)):(marks(tt,mark_off)+tol2)]; %#ok<AGROW>
%         else
%             PC_aux_ind=[PC_aux_ind round((marks(tt,5)-tol3)):(marks(tt,5)+tol3)]; %#ok<AGROW>
%         end
%     end
% else
%     PC_aux_indaux= 1:(t(2)-t(1)+1);   %EXARLE OJO %PC_aux_indaux=[];%alternativa
% end
% try
%     heasig.gain(heasig.gain==0)=200;  % When heasig.gain = 0 => default 200
% catch me
%     me.message = 'Signal gain is set to 200'; % Juan 28/03/2012
%     messages.warnings = [messages.warnings {me.message}]; % Juan 28/03/2012
% end


% Reading of external QRS annotation file if qrs_flag
% if qrs_flag==1     % The QRS fiducial point is read from external annotator
%     switch (anot_fmt);
%         case 0        % MIT annotation file
%             if exist([dirann filesep aname],'file')%Jul2011 % Jbolea 01/12/11
%                 s=readannot([dirann aname],t);
%                 s=isqrs(s,heasig,t);
%                 ext_anot=s.time';
%             else error('QRS annotation file not found');
%             end
%         case 1        % mat file
%             if exist([dirann aname],'file')%Jul2011
% %                 if exist([dirann aname '.mat'],'file')%Jul2011
%                 
%                 try
%                     s= readannot_mat(dirann, aname,'tm_qrs');%OCT2012
%                     s.tm_qrs.pos = round(s.tm_qrs.pos*heasig.freq/1000);
%                     ext_anot=s.tm_qrs.pos(s.tm_qrs.pos >= t(1)&s.tm_qrs.pos <= t(2));%OCT2012
%                 catch
%                     
%                     s= readannot_mat(dirann, aname,'qrs');%DEC2011
%                     s.qrs.pos = round(s.qrs.pos*heasig.freq/1000);
%                     ext_anot=s.qrs.pos(s.qrs.pos >= t(1)&s.qrs.pos <= t(2));%Jul2011 % Jbolea 01/12/11
%                 end
%                 %                 aux_t=t/heasig.freq/60/60; %%%???? ver heasig.freq % Juan
%                 %                 28/03/2012
%                 
%                 %s=load ([dirann aname]);
%                 %kk=getfield(s);%JUL2011
%                 % eval(['ext_anot=s.qrs;']);
%             else error('QRS annotation file not found');
%             end
%         otherwise
%             messages.warnings=[messages.warnings {'Bad annotation format for external QRS annotation file; internal annotator would be considered instead.'}];
%             warning(char(messages.warnings(end)));
%             qrs_flag=0;
%     end
% end
% 
messages.setup.wavedet.qrs_flag=qrs_flag;
messages.setup.wavedet.leadsynth_flag=leadsynth_flag;
messages.setup.wavedet.flagaux=flagaux;

%%%%%%%%%%%%%%%%%%%%%%% Initiliazation %%%%%%%%%%%%%%%%%%%%%%%%%%%%
inisamp = t(1);
endsamp = t(1)-1; %15MAR09
timeqrs1 = [];timeqrs2 = [];timeqrs3 = [];
lastqrs1 = 1;lastqrs2 = 1;lastqrs3 = 1;

if qrs_flag==0 || qrs_flag==2 % estimated maximum number of beats
    maxlength = round((t(2)-t(1))/messages.setup.wavedet.freq*3);
else maxlength=length(ext_anot);
end

if qrs_flag==1,
%     sel=find(ext_anot<messages.setup.wavedet.freq);
%     if ~isempty(sel)
%         ext_anot(sel)=[];
        maxlength= length(ext_anot);
%     end
end

nanvec = nan(1,maxlength);
position=struct('Pon',nanvec,'P',nanvec,'Poff',nanvec,...
    'Pprima',nanvec,'Pscale',nanvec,'Ptipo',nanvec,...
    'QRSon',nanvec,'Q',nanvec,'R',nanvec,'Rprima',nanvec,'S',nanvec,'QRSoff',nanvec,...
    'qrs',zeros(1,maxlength),...
    'R_inQRSoff',nanvec,'R_inQRSon',nanvec,...
    'Ton',nanvec,'T',nanvec,'Tprima', nanvec,'Toff',nanvec,'Ttipo',nanvec,...
    'Tscale',nanvec,'Ttipoon',nanvec,'Ttipooff',nanvec,'contadorToff',nanvec,...
    'QRSonsetcriteria',nanvec,'QRSoffcriteria',nanvec,...
    'QRSpa',nanvec,'QRSpp',nanvec,...
    'QRSmainpos',nanvec,'QRSmaininv',nanvec);


position1=position;
position2=position;
position3=position;
position0=position;
%messages.errors=[messages.errors errores.errores];


if qrs_flag==1, position.qrs=ext_anot;  end
% numlatdet = 0;      % Number of detected beats until now  % Juan
% 28/03/2012  no se usa despues
numlatdet1 = 0;
numlatdet2 = 0;
numlatdet3 = 0;

if ~isfield(messages.setup.wavedet,'filter_bank_design') || messages.setup.wavedet.filter_bank_design~=1 % Rute 24/11/11
    messages.setup.wavedet.filter_bank_design=0; % Rute 24/11/11
    [q1,q2,q3,q4,q5,messages] = qspfilt5(messages.setup.wavedet.freq,messages);  % WT equivalent filters % Rute 22/05/02
    if messages.status==1
        l1=length(q1);l2=length(q2);l3=length(q3);l4=length(q4);
        d1=floor((l1-1)/2);d2=floor((l2-1)/2);
        d3=floor((l3-1)/2);d4=floor((l4-1)/2);
        l5=length(q5);d5=floor((l5-1)/2); % Rute 22/05/02
        begoverlap = l5+2* messages.setup.wavedet.freq;    % l5 + 2 sec.
        endoverlap = d5;
%         messages.warnings=[messages.warnings {'Default filter banks are used.'}];
%         warning(char(messages.warnings(end)));
    else
        messages.setup.wavedet.filter_bank_design=1; % Rute 24/11/11
        messages.status=1;
    end
end
%%%% % Rute 24/11/11
if  messages.setup.wavedet.filter_bank_design==1;
    MaxScales=5;
    filters_cache_filename = ['wt_filters_' num2str(MaxScales) ' scales_' num2str(messages.setup.wavedet.freq) ' Hz.mat' ];
    % check it form setup
    if( exist(filters_cache_filename, 'file') )
        
        db_stat = dbstatus();
        bCaughtErrors = false;
        for ii = 1:length(db_stat)
            if( strcmpi(db_stat(ii).cond, 'caught error') )
                bCaughtErrors = true;
                dbclear if caught error
                break
            end
        end
        
        load( filters_cache_filename );
        
        if(bCaughtErrors)
            dbstop if caught error
        end
    else
        q_filters = qs_filter_design(1:MaxScales, messages.setup.wavedet.freq); %ver spf
        try
            save([fileparts(which('wavedet_3D')) sep filters_cache_filename],'q_filters')
        catch me
            me.messages = ['Fatal error saving ' filters_cache_filename ' on ' fileparts(which('wavedet_3D')) '; please check folder permissions'];
            messages.warnings = [messages.warnings {me.messages}];
        end
    end
    begoverlap=messages.setup.wavedet.freq;  % Juan, Julia 15/06/2012
    endoverlap=messages.setup.wavedet.freq;  % Juan, Julia 15/06/2012
    messages.warnings=[messages.warnings {'arbitrary sampling frequency designed filter banks are used.'}];
%     warning(char(messages.warnings(end)));
end
%%%%%%%%%%

% The two seconds overlap makes sure that the last/first beat
% are completely within the excerpt processed

samp=1; %heasig.spf_ecg*heasig.freq;    %Initialization of samp (1 sec)
intervaloall=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
heasigreading=heasig; %11JUL2011
%%%%%%%%%%%%%%%%%%%%%%%%%%%% excerpt processing %%%%%%%%%%%%%%%%%%%%%%%%%%%
while ((endsamp+1) < t(2))
    fileevolution = [(endsamp+1) t(2)];
    
    firstnewsamp = samp(end);   % last sample analyzed
    endsamp = min(inisamp + nsamp -1,t(2));
    yy=1; %#ok<NASGU>
    %%%%% Protection 25/08/2010
    inisamp = round(inisamp);
    endsamp = round(endsamp);
%     display(['File Evolution (analyzing): ' timestr(round(inisamp/messages.setup.wavedet.freq*1000))...
%         ' to ' timestr(round(endsamp/messages.setup.wavedet.freq*1000))]);
    nsamp = round(nsamp);
    
    %%%%%%%%%%%%%%%%%%%%%%% loading ECG signal %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if isnan(ft)||ft==0||ft==1||ft==2||ft==4||ft>4
%         if ft==0||ft==4 || ft==1 %MIT , Lund, Mortara %21MAR2011
%             [sig,messages1,heasig] = readsignal(sigdir,ecgnr,[inisamp endsamp],ft,heasigreading,leadreading,messages); % read leadreading
%             if messages1.status==0
%                 messages.errors=[messages.errors messages1.errors];
%                 warning(char(messages.errors(end)))
%                 messages.errors_desc=[messages.errors_desc messages1.errors_desc];
%                 messages.status=0;
%                 return;
%             else
%                 messages.warnings=[messages.warnings messages1.warnings ];%Jul2011
%                 if ~isfield(messages.setup,'readsignal') || ~isfield(messages.setup.readsignal,'file')
%                     messages.setup.readsignal.file=messages1.setup.file;
%                 end
%                 if ~isfield(messages.setup.readsignal,'t')
%                     messages.setup.readsignal.t=messages1.setup.t;
%                 else
%                     messages.setup.readsignal.t(2)=messages1.setup.t(2);
%                 end
%                 if ~isfield(messages.setup.readsignal,'nsig'), messages.setup.readsignal.nsig=messages1.setup.nsig; end
%             end
%         else
%             sig=sig_all(inisamp-t(1)+1:min(endsamp-t(1)+1,size(sig_all,1)),:);
%         end
        
        sig=sig_all(inisamp-t(1)+1:min(endsamp-t(1)+1,size(sig_all,1)));
        
        
%         if size(sig,1)==length(leadreading) % Jbolea 21/03/2011
%             sig=sig';
%         end
% 
%         n_auxleads=find(auxleads>0); %#ok<NASGU>

%         if isfield(heasig,'units')
%             for auxlead=1:length(leadreading) %
%                 if strcmp(deblank(upper(heasig.units(leadreading(auxlead),:))),'MV')||strcmp(deblank(upper(heasig.units(leadreading(auxlead),:))),'MILLIVOLTS')||...
%                         strcmp(deblank(upper(heasig.units(leadreading(auxlead),:))),'MILIVOLTIOS')
%                     sig(leadreading(auxlead),:) = sig(leadreading(auxlead),:)*1e3; % conversion to microV
%                     messages.warnings=[messages.warnings {'signal units changed: converted from milivolts to microvolts'}];
%                     warning(char(messages.warnings(end)))
%                     heasig.units=(char(ones(heasig.nsig,1)*'microvolts')); % Jbolea 21/3/2011
%                 elseif strcmp(deblank(upper(heasig.units(leadreading(auxlead),:))),'nV')||strcmp(deblank(upper(heasig.units(leadreading(auxlead),:))),'NANOVOLTS')||...
%                     strcmp(deblank(upper(heasig.units(leadreading(auxlead),:))),'NANOVOLTIOS')
%                     sig(:,leadreading(auxlead)) = (sig(:,leadreading(auxlead)))/1e3; % conversion to microV
%                     messages.warnings=[messages.warnings {'signal units changed: converted from volts to microvolts'}];
%                     warning(char(messages.warnings(end)))
%                     heasig.units=(char(ones(heasig.nsig,1)*'microvolts'));
%                 elseif strcmp(deblank(upper(heasig.units(leadreading(auxlead),:))),'V')||strcmp(deblank(upper(heasig.units(leadreading(auxlead),:))),'VOLTS')||...
%                         strcmp(deblank(upper(heasig.units(leadreading(auxlead),:))),'VOLTIOS')
%                     sig(leadreading(auxlead),:) = (sig(leadreading(auxlead),:))*1e6; % conversion to microV
%                     messages.warnings=[messages.warnings {'signal units changed: converted from volts to microvolts'}];
%                     warning(char(messages.warnings(end)))
%                     heasig.units=(char(ones(heasig.nsig,1)*'microvolts'));
%                 end
%             end
%         else
%             for auxlead=1:length(leadreading) %
%                 sig(leadreading(auxlead),:) = sig(leadreading(auxlead),:).* 1e3;  %% by default
%             end
%             messages.warnings=[messages.warnings {'no signal units given: milivolts assumed'}];
%             warning(char(messages.warnings(end)))
%             messages.warnings=[messages.warnings {'signal units changed: converted from milivolts to microvolts'}];
%             warning(char(messages.warnings(end)))
%             heasig.units=(char(ones(heasig.nsig,1)*'microvolts'));
%         end
        
%         if leadsynth_flag==1||leadsynth_flag==2
%             
%             if isnan(auxleads(8)) %lead II is not available
%                 sinal= leadcalc([sig(:,auxleads(1:7))  sig(:,auxleads(7))+sig(:,auxleads(9)) ],leadstring );  %Juan 23/06/2011
%             else
%                 sinal= leadcalc(sig,leadstring);
%             end
%             sig=sinal(:,~isnan(sinal(1,:)))';
%         elseif leadsynth_flag==3||leadsynth_flag==4||leadsynth_flag==5 % PC based delineation
%             if leadsynth_flag==4
%                 PCleads=1:size(sig,2); % considering all available leads
%             elseif leadsynth_flag==3 % considering 8 uncorrelated leads
%                 if ~isnan(auxleads(1:7)) & (~isnan(auxleads(8))||~isnan(auxleads(9))) %#ok<AND2>
%                     if isnan(auxleads(8)) %lead II is not available
%                         PCleads=[1:7 9]; %leads defined in leadreading, plus constructed lead II
%                         sig(:,end+1)=sig(:,7)+sig(:,8); %#ok<AGROW>
%                     else
%                         PCleads=1:8; %leads are already defined in leadreading
%                     end
%                 else
%                     messages.warnings=[messages.warnings {'There are not 8 uncorrelated leads, all available leads considered instead'}];
%                     warning(char(position.warningss(end)));
%                     PCleads=1:size(sig,2);
%                     leadsynth_flag=4;
%                 end
%             elseif leadsynth_flag==5
%                 if length(lead)>1 % considering precordial uncorrelated leads
%                     PCleads=2:7; %leads defined in leadreading, plus aVF
%                 else
%                     PCleads=[];
%                 end
%             end
%             if ~exist('PC_aux_ind','var') || isempty(PC_aux_ind)
%                 PC_aux_indaux= 1:length(sig);
%             else
%                 if flagaux == 1                    
%                     PC_aux_indaux=PC_aux_ind(PC_aux_ind>=firstnewsamp & PC_aux_ind<=endsamp)-firstnewsamp+1;
%                 end
%             end
%             %no NaN allowed in sig 29JUN2011
%             indx = find(isnan(sig(:,1)) == 1);
%             if ~isempty(indx)
%                 sig(indx,:) = sig(indx-1,:);
%             end
%             
%             [eigvectors,eigvalues]=eig(cov(sig(PC_aux_indaux,PCleads)));
%             [s_eigvalues,order]=sort(diag(eigvalues)); %#ok<ASGLU>
%             %pc=(eigvectors(:, order(end:-1:(end-length(lead)+1))))'*sig(:,PCleads)';
%             pc=(eigvectors(:, order(end:-1:(end-2))))'*sig(:,PCleads)'; %11JUL2011
%             if leadsynth_flag==5
%                 if isnan(auxleads(12)) %lead aVF is not available
%                     if ~isnan(auxleads(7)) && ~isnan(auxleads(8))
%                         sig=[sig(:,auxleads(8))-1/2*sig(:,auxleads(7)) pc(1:min(2,(length(lead)-1)),:)'];
%                     elseif ~isnan(auxleads(7)) && ~isnan(auxleads(9))
%                         sig=[sig(:,auxleads(9))+1/2*sig(:,auxleads(7)) pc(1:min(2,(length(lead)-1)),:)'];
%                     elseif ~isnan(auxleads(8)) && ~isnan(auxleads(9))
%                         sig=[1/2*(sig(:,auxleads(8))+sig(:,auxleads(9))) pc(1:min(2,(length(lead)-1)),:)'];
%                     else
%                         lead=auxleads(find(~isnan(auxleads([7:end 1:6])))); %#ok<FNDSB>
%                         messages.errors=[messages.errors {'lead aVF not available, PC based on precordial leads used only'}];
%                         warning(char(messages.errors(end)));
%                         sig= pc(1:min(3,(length(lead)-1)),:)';
%                         leadsynth_flag=4;
%                     end
%                 elseif length(lead)>1
%                     sig=[sig(:,7) pc'];
%                     sig=sig(:,lead);
%                 end
%             else
%                 sig=pc(lead,:)';
%             end
%         else %leadsynth_flag==0: should be leadsreading=lead
%             if lead<=size(sig,2)
%                 sig=sig(:,lead); %leadinfile=leadinfile(lead);
%                 try
%                     messages.warnings=[messages.warnings {[ 'the lead(s) chosen to be processed are the lead(s)' num2str(leadreading(lead))]}];
%                 catch me
%                     me.message = 'Error including warning information';
%                 end
%             else % if leads are in leadreading
%                 [C,IA,IB]=intersect(leadreading,lead); %#ok<NASGU>
%                 if length(C)==length(lead)
%                     sig=sig(:,IA);
%                     messages.warnings=[messages.warnings {['the lead(s) chosen to be processed are the lead(s) ' num2str(leadreading(IA))]}];
%                     lead=leadreading(IA); %7ABRIL2011
%                 else  %should not append
%                     lead=leadreading(1:min(3,end));
%                     messages.warnings=[messages.warnings {[ 'unable to interpret the lead(s) number to process (' num2str(lead) '): lead(s) ' num2str(leadreading) ' considered instead.' ]}];
%                 end
%             end
%         end
        
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Modificado Juan 9/03/11 Rute 9May2011
    if size(sig,2)>1
        str = '';
    else
        str = '1';
    end
    %%%%%%%%%%%%%%%%%%%%%%% WT construction and parameters setting %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if  messages.setup.wavedet.filter_bank_design==0;% RUTE 24NOV2011
        wt = wavt5(sig(:,1)',q1,q2,q3,q4,q5);  % WT equivalent filtering % 22/05/02 Rute including scale 5
        wt=wt(l5:end,1:5);            % Remove "incorrect" samples
        samp = (inisamp + l5 -1):endsamp-d5;
        % First l5-1 samples are not correctly filtered (border effect)
        % Last d5 samples are discarded in order to alineate all the
        % filtered signals, taking into acount the filter delays
        
        % Synchronizing filtered signals at different scales
        w1 = zeros(size(wt,1)-d5,5);
        w1(:,5) = wt(d5+1:end,5);
        w1(:,4) = wt(d4+1:end+d4-d5,4);
        w1(:,3) = wt(d3+1:end+d3-d5,3);
        w1(:,2) = wt(d2+1:end+d2-d5,2);
        w1(:,1) = wt(d1+1:end+d1-d5,1);
        
    else % RUTE 24NOV2011
        wtECG = qs_wt(sig,1:5, messages.setup.wavedet.freq, q_filters);% RUTE 24NOV2011 ver pf
        clear w1
        w1(1:size(wtECG,1),1:5) = wtECG(:,1,:);% RUTE 24NOV2011
        samp =inisamp:endsamp;% RUTE 24NOV2011
        l5=1;d5=0;
    end
    interval_to_processed=[inisamp+l5+1 endsamp-d5];
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%check quality of the signal
    swt=w1.^2;
    % Initialization of the thresholds for QRS detection (eps)
    %step= messages.setup.wavedet.freq; %internal
%     if messages.setup.wavedet.sig_quality==2 %internal
%         step=heasig.spf_ecg*heasig.freq; %internal
%         indexes = signaltest(w1(:,1),step,1,sqrt(mean(swt(:,1))));
%         indexes(indexes>length(w1))=[];
%         
%         if ~isempty(indexes)
%             segments= indexes([find(diff(indexes)~=1) 1+find(diff(indexes)~=1)])';
%             segments=[indexes(1) segments(:)' indexes(end)];
%             segments=reshape(segments,2,length(segments)/2)'-d5-l5+1+t(1)+fileevolution(1)-2;
%             nsegments=size(segments,1);
%             resultados2xls(char(ones(nsegments,1)*ecgnr),[ones(nsegments,1)*[ft lead(1)] segments],[dirann_sig_quality aname_sig_quality '_sig_quality'],date);
%         end
%     elseif messages.setup.wavedet.sig_quality==1 %external
%         if ~iscell(aname_sig_quality)
%             aname_sig_quality_cell{1}=aname_sig_quality(:,1);
%         end
%         anots = readannot_mat(dirann_sig_quality,aname_sig_quality_cell{1},'no_quality_segments');
%         data = [anots.no_quality_segments.pos anots.no_quality_segments.info{1}];% timing and duration? units?
%         if  strcmpi(anots.no_quality_segments.units,'sec') ||   strcmpi(anots.no_quality_segments.units,'seconds')
%             data=messages.setup.wavedet.freq*data;
%         elseif  strcmpi(anots.no_quality_segments.units,'msec') || strcmpi(anots.no_quality_segments.units,'miliseconds')
%             data=messages.setup.wavedet.freq*data/1000;
%         else
%             messages.warnings=[messages.warnings {'no quality segments assumed to be in samples'}];
%         end
%         % for data in samples
%         indexes=data(data(:,1)>=interval_to_processed(1) & data(:,2)<=interval_to_processed(2),:);
%         indexes = intervalmatrix(indexes);
%         indexes(indexes>length(w1))=[];
%     else %no quality check
%         indexes=[];
%     end
%     if ~isempty(indexes)
%         w1(indexes,1)=0;
%         swt(indexes,1)=0;
%     end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     eps= 0.5*sqrt(nanmean(swt));
    eps= rowvec(0.5*sqrt(nanmedian(swt))) .* rowvec(messages.setup.wavedet.QRS_detection_thr);
    eps(4) = eps(4)*2;
    
    %%%%%%% multilead%%% Rute 02.Dec.04
%     if size(sig,2)>1
%         % lead 2
%         if  messages.setup.wavedet.filter_bank_design==0;% RUTE 24NOV2011
%             wt = wavt5(sig(:,2)',q1,q2,q3,q4,q5);  % WT equivalent filtering % 22/05/02 Rute including scale 5
%             wt=wt(l5:end,1:5);            % Remove "incorrect" samples
%             %samp = (inisamp + l5 -1):endsamp-d5;
%             % First l5-1 samples are not correctly filtered (border effect)
%             % Last d5 samples are discarded in order to alineate all the
%             % filtered signals, taking into acount the filter delays
%             % Synchronizing filtered signals at different scales
%             w2 = zeros(size(wt,1)-d5,5);
%             w2(:,5) = wt(d5+1:end,5);
%             w2(:,4) = wt(d4+1:end+d4-d5,4);
%             w2(:,3) = wt(d3+1:end+d3-d5,3);
%             w2(:,2) = wt(d2+1:end+d2-d5,2);
%             w2(:,1) = wt(d1+1:end+d1-d5,1);
%         else
%             clear w2
%             w2(1:size(wtECG,1),1:5) = wtECG(:,2,:);% RUTE 24NOV2011
%         end
%         swt=w2.^2;
%         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%check quality of the signal
%         if messages.setup.wavedet.sig_quality==1 %external
%             if iscell(aname_sig_quality) && size(aname_sig_quality,2)>1
%                 aname_sig_quality_cell{2}=aname_sig_quality(:,2);
%                 anots = readannot_mat(dirann_sig_quality,aname_sig_quality_cell{2},'no_quality_segments');
%                 data = [anots.no_quality_segments.pos anots.no_quality_segments.info{1}];% timing and duration? units?
%                 if  strcmpi(anots.no_quality_segments.units,'sec') ||   strcmpi(anots.no_quality_segments.units,'seconds')
%                     data=messages.setup.wavedet.freq*data;
%                 elseif  strcmpi(anots.no_quality_segments.units,'msec') || strcmpi(anots.no_quality_segments.units,'miliseconds')
%                     data=messages.setup.wavedet.freq*data/1000;
%                 else
%                     messages.warnings=[messages.warnings {'no quality segments assumed to be in samples'}];
%                 end                % for data in samples
%                 indexes=data(data(:,1)>=interval_to_processed(1) & data(:,2)<=interval_to_processed(2),:);
%                 indexes = intervalmatrix(indexes);
%             end
%             indexes(indexes>length(w2))=[];
%         else %no quality check
%             indexes=[];
%         end
%         if ~isempty(indexes)
%             w2(indexes,1)=0;
%             swt(indexes,1)=0;
%         end
%         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%         
%         % Initialization of the thresholds for QRS detection (eps)
%         eps(2,:)= 0.5*sqrt(mean(swt));
%         eps(2,4) = eps(2,4)*2;
%         
%         if size(sig,2)>2
%             % lead 3
%             if  messages.setup.wavedet.filter_bank_design==0;% RUTE 24NOV2011
%                 wt = wavt5(sig(:,3)',q1,q2,q3,q4,q5);  % WT equivalent filtering % 22/05/02 Rute including scale 5
%                 wt=wt(l5:end,1:5);            % Remove "incorrect" samples
%                 %samp = (inisamp + l5 -1):endsamp-d5;
%                 % First l5-1 samples are not correctly filtered (border effect)
%                 % Last d5 samples are discarded in order to alineate all the
%                 % filtered signals, taking into acount the filter delays
%                 
%                 % Synchronizing filtered signals at different scales
%                 w3 = zeros(size(wt,1)-d5,5);
%                 w3(:,5) = wt(d5+1:end,5);
%                 w3(:,4) = wt(d4+1:end+d4-d5,4);
%                 w3(:,3) = wt(d3+1:end+d3-d5,3);
%                 w3(:,2) = wt(d2+1:end+d2-d5,2);
%                 w3(:,1) = wt(d1+1:end+d1-d5,1);
%             else
%                 clear w3
%                 w3(1:size(wtECG,1),1:5) = wtECG(:,3,:);% RUTE 24NOV2011
%             end
%             swt=w3.^2;
%             %%%%%%%%%%%%%%%%%%%%%%%%%%%%%check quality of the signal
%             if messages.setup.wavedet.sig_quality==1 %external
%                 if iscell(aname_sig_quality) && size(aname_sig_quality,2)>2
%                     aname_sig_quality_cell{3}=aname_sig_quality(:,3);
%                     anots = readannot_mat(dirann_sig_quality,aname_sig_quality_cell{3},'no_quality_segments');
%                     data = [anots.no_quality_segments.pos anots.no_quality_segments.info{1}];% timing and duration? units?
%                     if  strcmpi(anots.no_quality_segments.units,'sec') ||   strcmpi(anots.no_quality_segments.units,'seconds')
%                         data=messages.setup.wavedet.freq*data;
%                     elseif  strcmpi(anots.no_quality_segments.units,'msec') || strcmpi(anots.no_quality_segments.units,'miliseconds')
%                         data=messages.setup.wavedet.freq*data/1000;
%                     else
%                         messages.warnings=[messages.warnings {'no quality segments assumed to be in samples'}];
%                     end
%                     % for data in samples
%                     indexes=data(data(:,1)>=interval_to_processed(1) & data(:,2)<=interval_to_processed(2),:);
%                     indexes = intervalmatrix(indexes);
%                 end
%                 indexes(indexes>length(w3))=[];
%             else %no quality check
%                 indexes=[];
%             end
%             if ~isempty(indexes)
%                 w3(indexes,1)=0;
%                 swt(indexes,1)=0;
%             end
%             %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%             % Initialization of the thresholds for QRS detection (eps)
%             eps(3,:)= 0.5*sqrt(mean(swt));
%             eps(3,4) = eps(3,4)*2;
%         end
%     end
%     
    sig=sig(l5:end-d5,:); % piece of signal processed in this iteration
    %messages.lixo=sig;
    eps=eps';
    clear wt swt;
    
%%    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%% beat detection %%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    if qrs_flag==0 || qrs_flag==2   % The QRS fiducial point is calculated with wavedet
        messages.sig=sig;
        [position1,timeqrs1,lastqrs1,intervalo1,numlatdet1,time1,messages]= fiducialf(position1,firstnewsamp,samp,heasig,w1,eps,1,lastqrs1,timeqrs1,numlatdet1,ultimo_anot,messages);
        timeqrs=timeqrs1; % comentar isto pode dar erro no single lead
        indexes=1:length(timeqrs1);
        
%         if size(sig,2)>1
%             [position2,timeqrs2,lastqrs2,intervalo2,numlatdet2,time2,messages]= fiducialf(position2,firstnewsamp,samp,heasig,w2,eps,2,lastqrs2,timeqrs2,numlatdet2,ultimo_anot,messages);
%             if size(sig,2)>2
%                 [position3,timeqrs3,lastqrs3,intervalo3,numlatdet3,time3,messages]= fiducialf(position3,firstnewsamp,samp,heasig,w3,eps,3,lastqrs3,timeqrs3,numlatdet3,ultimo_anot,messages);
%             else
%                 timeqrs3=[];
%             end
%             [timeqrs,indexes]=qrscandidatesnew(timeqrs1,timeqrs2,timeqrs3,heasig,messages.setup.wavedet.refrper(end)/2* messages.setup.wavedet.freq);
%             
%             %2JUN08 %RUTE no candidates at all do not mark
%             if ~isempty(timeqrs)
%                 lastqrs1=timeqrs(1,~isnan(timeqrs(1,:)));
%                 if ~isempty(lastqrs1)
%                     lastqrs1=lastqrs1(end);
%                 end
%                 
%                 if sum(~isnan(timeqrs(2,:)))>0
%                     lastqrs2=timeqrs(2,~isnan(timeqrs(2,:)));
%                 end
%                 if ~isempty(lastqrs2)
%                     lastqrs2=lastqrs2(end);
%                 end
%                 if size(sig,2)>2
%                     if sum(~isnan(timeqrs(3,:)))>0
%                         lastqrs3=timeqrs(3,~isnan(timeqrs(3,:)));
%                     end
%                     if ~isempty(lastqrs3)
%                         lastqrs3=lastqrs3(end);
%                     end
%                 end
%             end
%         end
        
        intervaloall=[intervaloall(end)+1 size(timeqrs,2)];
        intervalo=intervaloall;
        time= timeqrs(:,intervaloall(1):intervaloall(end))-samp(1)+1;
        
        if ~isempty(time) && size(sig,2)>1 %16OUT08
            position0.qrs(intervalo(1):intervalo(end))=round(nanmedian(time)+samp(1)-1); % 09AGO2011
            position.qrs(intervalo(1):intervalo(end))=round(nanmedian(time)+samp(1)-1);  % 09AGO2011
        end
    else           % The QRS fiducial point is read from external annotator % multilead need to be checked
        
        
        first = max(firstnewsamp-ceil(messages.setup.wavedet.freq), samp(1)-1+ceil( messages.setup.wavedet.freq*0.050));
        % first is calculated according with the first lines of fiducial JP
        sel=find(ext_anot>=first & ext_anot<=samp(end)); %
%         position1.qrs(sel)=ext_anot(sel); %4ABRIL2011
        position1.qrs = rowvec(ext_anot);
        timeqrs(1,sel)=position1.qrs(sel)-samp(1)+1; %global
        time=timeqrs(1,sel);%4ABRIL2011 ?????%local
        indexes= 1:length(timeqrs(1,:));%4ABRIL2011
        lastqrs1=timeqrs(1,~isnan(timeqrs(1,:)));
        if ~isempty(lastqrs1)
            lastqrs1=lastqrs1(end);
        end
        if ~isempty(sel)
            intervaloall=[sel(1) sel(end)];
            % time=[timeqrs(:,intervaloall(1):intervaloall(end))-samp(1)+1];
        else
            intervaloall=[];
            %time=[];
        end
        intervalo=intervaloall;
        intervalo1=intervalo;%4ABRIL2011
        time1=time;%4ABRIL2011
        
        if size(sig,2)>1
            position2.qrs(sel)=ext_anot(sel);%4ABRIL2011
            timeqrs(2,sel)=position2.qrs(sel)-samp(1)+1;
            time(2,:)= time(1,:);%4ABRIL2011
            indexes=[indexes; indexes];%#ok<AGROW>%4ABRIL2011
            if sum(~isnan(timeqrs(2,:)))>0
                lastqrs2=timeqrs(2,~isnan(timeqrs(2,:)));
            end
            if ~isempty(lastqrs2)
                lastqrs2=lastqrs2(end);
            end
            intervalo2=intervalo;%4ABRIL2011
            time2=time;
            
            if size(sig,2)>2
                position3.qrs(sel)=ext_anot(sel);%4ABRIL2011
                timeqrs(3,sel)=position3.qrs(sel)-samp(1)+1;
                time=[time;time(1,:)];%#ok<AGROW> %4ABRIL2011
                indexes=[indexes; indexes];%#ok<AGROW>%4ABRIL2011
                if sum(~isnan(timeqrs(3,:)))>0
                    lastqrs3=timeqrs(3,~isnan(timeqrs(3,:)));
                end
                if ~isempty(lastqrs3)
                    lastqrs3=lastqrs3(end);
                end
                
                intervalo3=intervalo;%4ABRIL2011
                time3=time;
            else
                indexes=[indexes;NaN*(1:length(timeqrs(1,:)))];%#ok<AGROW>%9MAYO2011
            end
        end
        
        if ~isempty(time) &&  size(sig,2)>2 %Rute 13Set06 %by default position is the median mark
            position0.qrs(intervalo(1):intervalo(end))=round(nanmedian(time)+samp(1)-1); % 09AGO2011
            position.qrs(intervalo(1):intervalo(end))=round(nanmedian(time)+samp(1)-1); % 09AGO2011
        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.
        
        %%%%%%%%%%%%%%%%% Rute 03.Dec.04?
        %         if abs(samp(time(1)) - lastqrs)< 0.275*heasig.spf_ecg*heasig.freq,
        %             time(1)=[];  sel(1)=[];              % In order not to repeat beats
        %         end
        %         if (length(sig)-time(end) < 0.45*heasig.spf_ecg*heasig.freq),
        %             time(end)=[]; sel(end)=[];
        %         end
        % If last beat's T-wave is not in 'sig', the beat is detected in next one
        %%% timeqrs = [timeqrs samp(time)];  % QRS times ??????????????? Rute 03.Dec.04?
        %lastqrs = samp(time(end));Rute 03.Dec.04?
        % numeration of the beats processed in this segment
    end
    if size(sig,2)>1
        timenew=round(nanmedian(time)); %RUTE 09/07/2007
    else
        timenew=time;
    end
    
    
    if( ~messages.setup.wavedet.QRS_detection_only )
    
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%% QRS delineation %%%%%%%%%%%%%%%%%%%%%%%%%%%
        %if ~isempty(time(1,:)) %2JUN08 RUTE
        if ~isempty(time) %2JUN08 RUTE
            QRSon=indexes(:,intervalo(1):intervalo(end));
            qrspiconsetall=NaN*ones(3,intervalo(end)-intervalo(1)+1);
            QRSoff=indexes(:,intervalo(1):intervalo(end));
            qrspicoffsetall=NaN*ones(3,intervalo(end)-intervalo(1)+1);
            if ~isempty(intervalo1)
                [position1,qrspiconset1,qrspicoffset1,messages]= qrswavef(heasig,samp,time1,position1,w1,intervalo1,messages);
                auxindexes=indexes(1,intervalo(1):intervalo(end));
                %NOTE: CHANGED 10.SET.05
                qrspiconsetall(1,(~isnan(indexes(1,intervalo(1):intervalo(end)))&(QRSon(1,:)>=intervalo1(1)&QRSon(1,:)<=intervalo1(2))))=qrspiconset1(QRSon(1,~isnan(QRSon(1,:))&(QRSon(1,:)>=intervalo1(1)&QRSon(1,:)<=intervalo1(2)))-intervalo1(1)+1); %17Mar06
                QRSon(1,~isnan(QRSon(1,:)))=position1.QRSon(QRSon(1,~isnan(QRSon(1,:))));

                qrspicoffsetall(1,(~isnan(indexes(1,intervalo(1):intervalo(end)))&(QRSoff(1,:)>=intervalo1(1)&QRSoff(1,:)<=intervalo1(2))))=qrspicoffset1(QRSoff(1,~isnan(QRSoff(1,:))&(QRSoff(1,:)>=intervalo1(1)&QRSoff(1,:)<=intervalo1(2)))-intervalo1(1)+1);
                QRSoff(1,~isnan(QRSoff(1,:)))=position1.QRSoff(QRSoff(1,~isnan(QRSoff(1,:))));
            end
            if size(sig,2)>1
                [position2,qrspiconset2,qrspicoffset2,messages]= qrswavef(heasig,samp,time2,position2,w2,intervalo2,messages);
                if ~isempty(intervalo2)
                    auxindexes=indexes(2,intervalo(1):intervalo(end));
                    qrspiconsetall(2,(~isnan(indexes(2,intervalo(1):intervalo(end)))&(QRSon(2,:)<=intervalo2(2)&QRSon(2,:)>=intervalo2(1))))=qrspiconset2(QRSon(2,~isnan(QRSon(2,:))&(QRSon(2,:)<=intervalo2(2)&QRSon(2,:)>=intervalo2(1)))-intervalo2(1)+1); %Rute 17Mar06
                    QRSon(2,~isnan(QRSon(2,:)))=position2.QRSon(QRSon(2,~isnan(QRSon(2,:))));
                    qrspicoffsetall(2,(~isnan(indexes(2,intervalo(1):intervalo(end)))&(QRSoff(2,:)<=intervalo2(2)&QRSoff(2,:)>=intervalo2(1))))=qrspicoffset2(QRSoff(2,~isnan(QRSoff(2,:))&(QRSoff(2,:)<=intervalo2(2)&QRSoff(2,:)>=intervalo2(1)))-intervalo2(1)+1); %Rute 17Mar06
                    QRSoff(2,~isnan(QRSoff(2,:)))=position2.QRSoff(QRSoff(2,~isnan(QRSoff(2,:))));
                end
                if size(sig,2)>2 && ~isempty(intervalo3)
                    [position3,qrspiconset3,qrspicoffset3,messages]= qrswavef(heasig,samp,time3,position3,w3,intervalo3,messages); %#ok<ASGLU>
                    auxindexes=indexes(3,intervalo(1):intervalo(end));
                    if max(QRSon(3,~isnan(QRSon(3,:))&(QRSon(3,:)>=intervalo3(1)&QRSon(3,:)<=intervalo3(2)))-intervalo3(1)+1)> length(qrspiconset3)%17.Mar.06
                        qrspiconsetall(3,(QRSon(3,~isnan(QRSon(3,:))&(QRSon(3,:)<=intervalo3(2)&QRSon(3,:)>=intervalo3(1)))-intervalo3(1)+1)<=length(qrspiconsetall(3,(~isnan(indexes(3,intervalo(1):intervalo(end)))&(QRSon(3,:)<=intervalo3(2)&QRSon(3,:)>=intervalo3(1)))))&(~isnan(indexes(3,intervalo(1):intervalo(end)))&(QRSon(3,:)>=intervalo3(1))))=qrspiconset3((QRSon(3,~isnan(QRSon(3,:))&(QRSon(3,:)>=intervalo3(1)))-intervalo3(1)+1)<=length(qrspiconsetall(3,(~isnan(indexes(3,intervalo(1):intervalo(end)))&(QRSon(3,:)>=intervalo3(1)))))& QRSon(3,~isnan(QRSon(3,:))&(QRSon(3,:)>=intervalo3(1)))-intervalo3(1)+1);%17Mar06
                    else%16.Fev.06
                        qrspiconsetall(3,(~isnan(indexes(3,intervalo(1):intervalo(end)))&(QRSon(3,:)>=intervalo3(1))))=qrspiconset3(QRSon(3,~isnan(QRSon(3,:))&(QRSon(3,:)>=intervalo3(1)))-intervalo3(1)+1);
                    end %16.Fev.06
                    QRSon(3,~isnan(QRSon(3,:)))=position3.QRSon(QRSon(3,~isnan(QRSon(3,:))));

                    if max(QRSoff(3,~isnan(QRSoff(3,:))&(QRSoff(3,:)>=intervalo3(1)&QRSoff(3,:)<=intervalo3(2)))-intervalo3(1)+1)> length(qrspiconset3)%17.Mar.06
                        qrspicoffsetall(3,(QRSoff(3,~isnan(QRSoff(3,:))&(QRSoff(3,:)<=intervalo3(2)&QRSoff(3,:)>=intervalo3(1)))-intervalo3(1)+1)<=length(qrspicoffsetall(3,(~isnan(indexes(3,intervalo(1):intervalo(end)))&(QRSoff(3,:)<=intervalo3(2)&QRSoff(3,:)>=intervalo3(1)))))&(~isnan(indexes(3,intervalo(1):intervalo(end)))&(QRSoff(3,:)>=intervalo3(1))))=qrspiconset3((QRSoff(3,~isnan(QRSoff(3,:))&(QRSoff(3,:)>=intervalo3(1)))-intervalo3(1)+1)<=length(qrspicoffsetall(3,(~isnan(indexes(3,intervalo(1):intervalo(end)))&(QRSoff(3,:)>=intervalo3(1)))))& QRSoff(3,~isnan(QRSoff(3,:))&(QRSoff(3,:)>=intervalo3(1)))-intervalo3(1)+1);%17Mar06
                    else%16.Fev.06
                        qrspicoffsetall(3,(~isnan(indexes(3,intervalo(1):intervalo(end)))&(QRSoff(3,:)>=intervalo3(1))))=qrspiconset3(QRSoff(3,~isnan(QRSoff(3,:))&(QRSoff(3,:)>=intervalo3(1)))-intervalo3(1)+1);
                    end %16.Fev.06
                    QRSoff(3,~isnan(QRSoff(3,:)))=position3.QRSoff(QRSoff(3,~isnan(QRSoff(3,:))));
                end
            end
            auxqrspiconset=nanmin(qrspiconsetall);
            auxqrspicoffset=nanmax(qrspicoffsetall);
            aux=find(isnan(auxqrspiconset));
            auxaux=intervalo(1):intervalo(end);
            %
            if ~isempty(aux)  % one beat is from the former interval  1.AGO.07
%                 if aux(1)>1
%                     messages.errors=[messages.errors {['QRS onset not found at beat(s): ' num2str(auxaux(aux ))]}];
%                     warning(char(messages.errors(end)));
%                 end
                auxqrspiconset(aux)=[];
                QRSon(:,aux)=[];  %31.07.07
                auxqrspicoffset(aux)=[];
                QRSoff(:,aux)=[]; %31.07.07
                intervalo(end)=intervalo(end)-aux(end); %16MAR09
                time(:,aux)=[];%31.07.07
                timenew(:,aux)=[];%31.07.07
                if length(auxaux)>size(QRSon,2)
                    auxaux(aux)=[];
                end
            end

            if size(sig,2)>1 %Rute 16Out08
                %by default QRS onset in the earliest of the marks
                %by default QRS end in the latest of the marks
                position0.QRSon(auxaux)=min(QRSon);
                position.QRSon(auxaux)=min(QRSon);
                position0.QRSoff(auxaux)=max(QRSoff);
                position.QRSoff(auxaux)=max(QRSoff);
            end
            if size(sig,2)>1
                if ~exist('w3','var')
                    w3=[];
                end
                timeaux=time; %time is sorted %18.07.07
                auxaux=NaN*ones(3,intervalo(end)-intervalo(1)+1);
                auxindexes=indexes(:,intervalo(1):intervalo(end));
                auxaux(1,~isnan(indexes(1,intervalo(1):intervalo(end))))=position1.R(auxindexes(1,~isnan(indexes(1,intervalo(1):intervalo(end)))));
                auxaux(2,~isnan(indexes(2,intervalo(1):intervalo(end))))=position2.R(auxindexes(2,~isnan(indexes(2,intervalo(1):intervalo(end)))));
                auxaux(3,~isnan(indexes(3,intervalo(1):intervalo(end))))=position3.R(auxindexes(3,~isnan(indexes(3,intervalo(1):intervalo(end)))));
                position.R(intervalo(1):intervalo(end))=nanmin(auxaux);

                %timenew is the median, use min instead?
                [position0,position,messages]=multiqrson(heasig,w1,w2,w3,auxqrspiconset,QRSon,timeaux,timenew,position0,position,intervalo,samp,messages);
                timeaux=time;

                position.R(intervalo(1):intervalo(end))=round(nanmedian(auxaux)); %RUTE 09AGO2011 %MEDIAN %by default position is the median mark of all R
                %timenew is the median, use max instead?
                [position0,position,messages]=multiqrsend(heasig,w1,w2,w3,auxqrspicoffset,QRSoff,timeaux,timenew,position0,position, intervalo,samp,messages);
                position.R(intervalo(1):intervalo(end))=round(nanmedian(auxaux));%RUTE 09AGO2011
                %            squaresignal=(sig(:,1)/max(sig(:,1))).^2+(sig(:,2)/max(sig(:,2))).^2+(sig(:,3)/max(sig(:,3))).^2;
                %             for ii=intervalo(1):intervalo(end)
                %             if ~isempty(position.QRSon(ii):position.QRSoff(ii))
                %                 [M,position.qrsnew(ii)]=max(squaresignal(position.QRSon(ii):position.QRSoff(ii)));
                %                 position.qrsnew(ii)=position.qrsnew(ii)+position.QRSon(ii)-1;
                %             end
                %             end
            end
            %         lastqrs1=timeqrs(1,end);
            %%%%%%%%%%%%%%%% caso em que rr tem dimensao inferior ao necessario
            %%%%%%%%%%%%%%%% protection with respect to ultimo_anot

            while eval(['position' str '.QRSon(intervalo(1))<ultimo_anot'])
                if ultimo_anottyp==10; % Toff
                    if eval(['(position' str '.qrs(intervalo(1))-position' str '.QRSon(intervalo(1))) < (position' str '.Toff(intervalo(1)-1)-position' str '.T(intervalo(1)-1))'])
                        eval(['position' str '.Toff(intervalo(1)-1)=NaN;'])
                        if eval(['isnan(position' str '.Tprima(intervalo(1)-1))'])
                            eval(['ultimo_anot=position' str '.T(intervalo(1)-1);'])
                            ultimo_anottyp=8;
                        else
                            eval(['ultimo_anot=position' str '.Tprima(intervalo(1)-1);'])
                            ultimo_anottyp=9;
                        end
                    else
                        eval(['position' str '.QRSon(intervalo(1))=NaN;'])
                    end
                elseif ultimo_anottyp==9 || ultimo_anottyp==8 || ultimo_anottyp==5; % Tprima | Tpeak | qrs
                    eval(['position' str '.QRSon(intervalo(1))=NaN;'])
                else  % ultimo_anottyp==6 % QRSoff
                    if eval(['(position' str '.qrs(intervalo(1))-position' str '.QRSon(intervalo(1))) < (position' str '.QRSoff(intervalo(1)-1)-position' str '.qrs(intervalo(1)-1))'])
                        eval(['position' str '.qrsoff(intervalo(1)-1)=NaN;'])
                        eval(['ultimo_anot=position' str '.qrs(intervalo(1)-1);'])
                        ultimo_anottyp=5;
                    else
                        eval(['position' str '.QRSon(intervalo(1))=NaN;'])
                    end
                end
            end
            % Modificado Juan 9/03/11
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            if qrs_flag~=2
                %%%%%%%%%%%%%%%%%%%%%%%%%%%% T wave
                if ~isempty(time1)
                    [position1,janelas1,messages]= twavef(heasig,samp,time1,position1,w1,intervalo1,messages); %21AGO09
                end
                %end
                if size(sig,2)>1 && ~isempty(time2)
                    if length(time2(1,:))>1%23SET08
                        [position2,janelas2,messages]= twavef(heasig,samp,time2,position2,w2,intervalo2,messages);%21AGO09
                    end
                    if size(sig,2)>2  && ~isempty(intervalo3)
                        if length(time3(1,:))>1%23SET08
                            [position3,janelas3,messages]= twavef(heasig,samp,time3,position3,w3,intervalo3,messages);%21AGO09
                        end
                    end
                end
                if size(sig,2)>1
                    subintervalo=indexes(:,intervalo(1):intervalo(end));
                    %NOTE are excluded the ones from the previous seg
                    if ~isempty(intervalo1)
                        aux1=subintervalo(1,~isnan(subintervalo(1,:)))-intervalo1(1)+1;
                    end
                    if size(sig,2)>1 & ~isempty(intervalo2) %#ok<AND2>
                        aux2=subintervalo(2,~isnan(subintervalo(2,:)))-intervalo2(1)+1;
                    end
                    if size(sig,2)>2 & ~isempty(intervalo3) %#ok<AND2>
                        aux3=subintervalo(3,~isnan(subintervalo(3,:)))-intervalo3(1)+1;
                    end
                    %26.04.05
                    janelas=NaN*ones(size(subintervalo,2),9);
                    janelas(:,1)=(intervalo(1):intervalo(end))';
                    janelas(:,4)=(intervalo(1):intervalo(end))';
                    janelas(:,7)=(intervalo(1):intervalo(end))';
                    if ~isempty(intervalo1)
                        janelas(~isnan(subintervalo(1,:)-intervalo1(1)+1),2)=aux1';
                    end
                    if ~isempty(intervalo2)
                        janelas(~isnan(subintervalo(2,:)-intervalo2(1)+1),5)=aux2';
                    end
                    if ~isempty(intervalo1)
                        aux1=aux1(aux1>0);
                    end
                    aux2=aux2(aux2>0);
                    if size(sig,2)>2 & ~isempty(intervalo3) %#ok<AND2>
                        janelas(~isnan(subintervalo(3,:)-intervalo3(1)+1),8)=aux3';
                        aux3=aux3(aux3>0);
                    end
                    janelas(janelas<=0)=NaN;
                    if exist('janelas1','var')&& sum(~isnan(janelas(:,2)))~=0
                        janelas(~isnan(janelas(:,2)) & janelas(:,2)>0,2:3)=janelas1(aux1,2:3);
                    end
                    if exist('janelas2','var')& ~isempty(janelas(~isnan(janelas(:,5)) & janelas(:,5)>0 ,5:6)) %#ok<AND2>
                        if max(aux2) > length(janelas2)%17.Mar.06
                            janelas(~isnan(janelas(:,5)) & janelas(:,5)>0 ,5:6)=[janelas2((aux2) <= length(janelas2),2:3); NaN*ones(sum(aux2 > length(janelas2)),2)];%14.Mar.06
                        else%16.Fev.06
                            janelas(~isnan(janelas(:,5)) & janelas(:,5)>0 ,5:6)=janelas2(aux2,2:3);
                        end%16.Fev.06
                    end
                    if size(sig,2)>2 & ~isempty(intervalo3) %#ok<AND2>
                        %ver altera�ao por erro de indexes em irec5 de
                        %politec!!!!%16.Fev.06
                        if exist('janelas3','var')
                            if max(aux3) > length(janelas3)%16.Fev.06
                                %janelas(~isnan(janelas(:,8)) & janelas(:,8)>0 ,8:9)=[janelas3(aux3(1:(end-1)),2:3);NaN NaN];%16.Fev.06
                                janelas(~isnan(janelas(:,8)) & janelas(:,8)>0 ,8:9)=[janelas3((aux3) <= length(janelas3),2:3); NaN*ones(sum(aux3 > length(janelas3)),2)];%14.Mar.06
                            else%16.Fev.06
                                janelas(~isnan(janelas(:,8)) & janelas(:,8)>0 ,8:9)=janelas3(aux3,2:3);
                            end%16.Fev.06
                        end
                    end
                    %cases were begwin>endwin!!!
                    janelas((janelas(:,3))<=janelas(:,2)|isnan(janelas(:,3)),2)=NaN;
                    janelas(janelas(:,3)<=janelas(:,2)|isnan(janelas(:,2)),3)=NaN;
                    janelas(janelas(:,6)<=janelas(:,5)|isnan(janelas(:,6)),5)=NaN;
                    janelas(janelas(:,6)<=janelas(:,5)|isnan(janelas(:,5)),6)=NaN;
                    janelas(janelas(:,9)<=janelas(:,8)|isnan(janelas(:,9)),8)=NaN;
                    janelas(janelas(:,9)<=janelas(:,8)|isnan(janelas(:,8)),9)=NaN;
                    begwin= min(janelas(:,[2 5 8])'); %#ok<UDIM>
                    endwin= max(janelas(:,[3 6 9])'); %#ok<UDIM>
                    intreg=[begwin' endwin'];

                    count=count+1;
                    if ~exist('w3','var')
                        w3=[];
                    end
                    [position0,position,messages]= delineate3D(w1,w2,w3,intreg,timenew,position0,position,intervalo,position1,position2,position3,indexes,samp,[],messages); %   05MAY2011
                end
                %end

                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 24AGO2011
                if ~exist('str','var'), str=''; end
                eval(['position' str '.QRSon(position' str '.R<=(position' str '.QRSon+1))=NaN;']);
                eval(['position' str '.QRSoff(position' str '.QRSoff>=(position' str '.Ton-1))=NaN;']);
                eval(['position' str '.QRSoff(position' str '.QRSoff>=(position' str '.T-1))=NaN;']);
                eval(['position' str '.QRSoff(position' str '.QRSoff<=(position' str '.qrs+1))=NaN;']);
                eval(['position' str '.QRSoff(position' str '.R>=(position' str '.QRSoff-1))=NaN;']);
                eval(['position' str '.Toff(position' str '.Toff<=(position' str '.T+1))=NaN;']);
                eval(['position' str '.Ton(position' str '.Ton>=(position' str '.T-1))=NaN;']);
                eval(['ii=1:sum(position' str '.qrs>0);']);
                eval(['position' str '.QRSoff(position' str '.QRSoff(ii(1:end-1))>=(position' str '.R(ii(2:end))))=NaN;']);
                eval(['position' str '.QRSoff(position' str '.QRSoff(ii(1:end-1))>=(position' str '.qrs(ii(2:end))))=NaN;']);
                eval(['position' str '.Toff(position' str '.Toff(ii(1:end-1))>=position' str '.qrs(ii(2:end)))=NaN;']);

                %%%%%%%%%%%%%%%%%%%%%%%%%%%% P wave %%%%%%%%%%%%%%%%%%%%%%%%%%%
                %11.04.05 P wave delineation inside the condition of existing some
                [position1,picon_all,picoff_all,messages]= pwavef(heasig,samp,time1,position1,w1,intervalo1,ultimo_anot,messages); %#ok<ASGLU>
                if size(sig,2)>1
                    [position2,picon_all,picoff_all,messages]= pwavef(heasig,samp,time2,position2,w2,intervalo2,ultimo_anot,messages); %#ok<ASGLU>
                    if size(sig,2)>2 & ~isempty(intervalo3) %#ok<AND2>
                        [position3,picon_all,picoff_all,messages]= pwavef(heasig,samp,time3,position3,w3,intervalo3,ultimo_anot,messages); %#ok<ASGLU>
                        auxaux(1,~isnan(indexes(1,intervalo(1):intervalo(end))))=position1.P(auxindexes(1,~isnan(indexes(1,intervalo(1):intervalo(end)))));
                        auxaux(2,~isnan(indexes(2,intervalo(1):intervalo(end))))=position2.P(auxindexes(2,~isnan(indexes(2,intervalo(1):intervalo(end)))));
                        auxaux(3,~isnan(indexes(3,intervalo(1):intervalo(end))))=position3.P(auxindexes(3,~isnan(indexes(3,intervalo(1):intervalo(end)))));
                        position0.P(intervalo(1):intervalo(end))=round(nanmedian(auxaux)); %RUTE 09AGO2011; %MEDIAN %by default position is the median mark of all R
                        [position0,position,messages]= delineateP3D(heasig,w1,w2,w3,timenew,position0,position,intervalo,samp,ultimo_anot,messages);
                    end
                end
                %check the position structure size
                eval(['position' str '.Pon((length(position' str '.Pon)+1):length(position' str '.qrs))=NaN;']); %6May2011
                eval(['position' str '.P((length(position' str '.P)+1):length(position' str '.qrs))=NaN;']); %6May2011
                eval(['position' str '.Poff((length(position' str '.Poff)+1):length(position' str '.qrs))=NaN;']); %6May2011
                eval(['position' str '.Poff((length(position' str '.Poff)+1):length(position' str '.qrs))=NaN;']); %6May2011
                if size(sig,2)>2 %6May2011
                    position0.Pon((length(position0.Pon)+1):length(position.qrs))=NaN;
                    position0.P((length(position0.P)+1):length(position.qrs))=NaN;
                    position0.Poff((length(position0.Poff)+1):length(position.qrs))=NaN;
                end
            end
            %%%%%%%%%%%%%%%%%%%%%%%%%%%% finishing excerpt %%%%%%%%%%%%%%%%%%%%%%%%%%%
            % Last annotated position
            if intervalo(2)>0
                if size(sig,2)>1
                    if ~isnan(position.Toff(min(intervalo(2),length(position.Toff)))),
                        ultimo_anot = position.Toff(min(intervalo(2),length(position.Toff)));
                        ultimo_anottyp=10;
                    elseif ~isnan(position.Tprima(min(intervalo(2),length(position.Tprima)))),
                        ultimo_anot = position.Tprima(min(intervalo(2),length(position.Tprima)));
                        ultimo_anottyp=9;
                    elseif ~isnan(position.T(min(intervalo(2),length(position.T)))),
                        ultimo_anot = position.T(min(intervalo(2),length(position.T)));
                        ultimo_anottyp=8;
                    elseif ~isnan(position.QRSoff(min(intervalo(2),length(position.QRSoff)))),
                        ultimo_anot = position.QRSoff(min(intervalo(2),length(position.QRSoff)));
                        ultimo_anottyp=6;
                    else
                        ultimo_anot = ceil(position.qrs(min(intervalo(2),length(position.qrs))));
                        ultimo_anottyp=5;
                    end
                else
                    if ~isnan(position1.Toff(min(intervalo(2),length(position1.Toff)))),
                        ultimo_anot = position1.Toff(min(intervalo(2),length(position1.Toff)));
                        ultimo_anottyp=10;
                    elseif ~isnan(position1.Tprima(min(intervalo(2),length(position1.Tprima)))),
                        ultimo_anot = position1.Tprima(min(intervalo(2),length(position1.Tprima)));
                        ultimo_anottyp=9;
                    elseif ~isnan(position1.T(min(intervalo(2),length(position1.T)))),
                        ultimo_anot = position1.T(min(intervalo(2),length(position1.T)));
                        ultimo_anottyp=8;
                    elseif ~isnan(position1.QRSoff(min(intervalo(2),length(position1.QRSoff)))),
                        ultimo_anot = position1.QRSoff(min(intervalo(2),length(position1.QRSoff)));
                        ultimo_anottyp=6;
                    else
                        ultimo_anot = ceil(position1.qrs(min(intervalo(2),length(position1.qrs))));
                        ultimo_anottyp=5;
                    end
                end
            end
        end %%%%%%%%%%%%%%%%introduzido a 17/05/02 Rute

    end
    
    
    inisamp = endsamp +2-endoverlap-begoverlap;
end

if( ~messages.setup.wavedet.QRS_detection_only )

    eval(['position' str '.Pon((length(position' str '.Pon)+1):length(position' str '.qrs))=NaN;']);
    eval(['position' str '.P((length(position' str '.P)+1):length(position' str '.qrs))=NaN;']);
    eval(['position' str '.Pprima((length(position' str '.Pprima)+1):length(position' str '.qrs))=NaN;']);
    eval(['position' str '.Poff((length(position' str '.Poff)+1):length(position' str '.qrs))=NaN;']);
    eval(['position' str '.Ton((length(position' str '.Ton)+1):length(position' str '.qrs))=NaN;']);
    eval(['position' str '.T((length(position' str '.T)+1):length(position' str '.qrs))=NaN;']);
    eval(['position' str '.Tprima((length(position' str '.Tprima)+1):length(position' str '.qrs))=NaN;']);
    eval(['position' str '.Toff((length(position' str '.Toff)+1):length(position' str '.qrs))=NaN;']);
    eval(['position' str '.Ttipo((length(position' str '.Ttipo)+1):length(position' str '.qrs))=NaN;']);
    eval(['position' str '.Tscale((length(position' str '.Tscale)+1):length(position' str '.qrs))=NaN;']);
    eval(['position' str '.Pon(position' str '.Pon>=(position' str '.qrs))=NaN;']);
    eval(['position' str '.P(position' str '.P>=(position' str '.qrs))=NaN;']);
    eval(['position' str '.Pprima(position' str '.Pprima>=(position' str '.qrs))=NaN;']);
    eval(['position' str '.Poff(position' str '.Poff>=(position' str '.qrs))=NaN;']);
    eval(['position' str '.Pon(position' str '.Pon>=(position' str '.QRSon))=NaN;']);
    eval(['position' str '.P(position' str '.P>=(position' str '.QRSon))=NaN;']);
    eval(['position' str '.Pprima(position' str '.Pprima>=(position' str '.QRSon))=NaN;']);
    eval(['position' str '.Poff(position' str '.Poff>=(position' str '.QRSon))=NaN;']);
    eval(['position' str '.Pon(position' str '.Pon>=(position' str '.P))=NaN;']);
    eval(['position' str '.Poff(position' str '.Poff<=(position' str '.P))=NaN;']);
    eval(['position' str '.Ton(position' str '.Ton>=(position' str '.T))=NaN;']);
    eval(['position' str '.Toff(position' str '.Toff<=(position' str '.T))=NaN;']);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%% position assigment %%%%%%%%%%%%%
% Remove void annotations at the end
%strf = {'Pon';'P';'Pprima';'Ptipo';'Poff';'QRSon';'Q';'R';'qrs';'Rprima';'S';'QRSoff';...
%    'QRSpa';'QRSpp';'Ton';'T';'Tprima';'Toff';'Ttipo';'Tscale';'QRSmainpos';'QRSmaininv'};
strf = [];
for i = 1:3
    eval(['strf =fieldnames(position' num2str(i) ');']); % RUTE 4 AGO2011
    eval(['aux = find(position' num2str(i) '.qrs == 0| isnan(position' num2str(i) '.qrs));']);
    
    for j = 1:length(strf)
        eval(['aux_l=aux(aux<=size(position' num2str(i) '.' strf{j} ',2));'])
        eval(['position' num2str(i) '.' strf{j} '(aux_l) = [];']);
    end
end

if size(lead,2)==1
    position=position1; % 28FEB2012
    if ~isempty(position.QRSmainpos) && ~isempty(position.QRSmaininv)
        a = length(find(position.QRSmainpos == position.qrs));
        b = length(find(position.QRSmaininv == position.qrs));
        if a >= b
            position.qrs_hrv = '+';
        else
            position.qrs_hrv = '-';
        end
    end
else
    %     strf = {'Pon';'P';'Poff';'QRSon';'Q';'R';'qrs';'Rprima';'S';'QRSoff';...
    %         'Ton';'T';'Tprima';'Toff';'Ttipo';'Tscale'};
    %     strf2 = {'Ttipoon';'Ttipooff'};
    %     strf3 = {'QRSonsetcriteria';'QRSoffcriteria';'contadorToff';...
    %         'R_inQRSon';'R_inQRSoff'};
    strn = {'0';''};
    for i = 1:2
        strM = 'MM = max([';
        eval(['aux = find(position' strn{i} '.qrs == 0 | isnan(position' strn{i} '.qrs));']);
        eval(['strf =fieldnames(position' strn{i} ');']); % RUTE 4 AGO2011
        strd = ',';
        for j = 1:length(strf)
            eval(['aux_l=aux(aux<=size(position' strn{i} '.' strf{j} ',2));'])
            eval(['position' strn{i} '.' strf{j} '(aux_l) = [];']);
            if j == length(strf)
                strd = ']);';
            end
            strM = [strM 'size(position' strn{i} '.' strf{j} ')' strd]; %#ok<AGROW>
        end
        eval(strM);
        for j = 1:length(strf)
            eval(['position' strn{i} '.' strf{j} '(:,end+1:MM)=NaN*ones(size(position' strn{i} '.' strf{j} ',1),MM-size(position' strn{i} '.' strf{j} ',2));'])
        end
    end
end
% try
%     annstruct = pos2ann(position);
%     if ~isempty(annstruct.time)
%         writeannot([matdir ecgnr '.' anot], annstruct);   % mexfile, faster!
%     end
% catch me
%     me.message = 'No MIT format annotation file saved.';
%     messages.warnings=[messages.warnings {me.message}];
% end
for j = 1:4
    eval(['positionaux.position' num2str(j-1) ' = position' num2str(j-1) ';']); % RUTE 09.07.2007
end
% clear global OPT ecgnr
if exist('fid','var') && fid > 2
    fclose(fid);  % JBOLEA 11MAY2010
end