wavedet.m

function position = wavedet(sigdir,headir,matdir,ecgnr,ft,anot,lead,t,qrs_flag,anot_fmt,aname,dirann)
%function position = mainwav(sigdir, headir, matdir, ecgnr, ft, anot, lead, t)
%Significant ECG points detector based on wavelet approach %Syntaxis:
%Input Parameters:

%	sigdir: directory of the original signal
%	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 1 for LUND header 2 for a mat file
% 	without header.
%	anot: name of the annotation output file
%	lead: scalar with the lead number to analyze (1,2, ...) (only one lead)
%	t=[tbegin tend]: time vector with initial and sample to analyze
%    qrs_flag: External (1) or internal (0) QRS detector
%    anot_fmt: in case of external QRS detector, format of annotation file:
%    aname: annotation file name with QRS marks
%    dirann: directory with external annotation file
%Output Parameters:
%	position: struct vector with the detected points


% 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.
% Juan Pablo Martínez Cortés
% 1/3/2006 option ft=2 was added from wavedetplus.
%%%%%alteraçoes para estudar a importancia do limiar de regularidade %%%%%%%
%%%%%% Rute 04/09/02
global regularity
regularity.time=[];
regularity.amp2=[];
regularity.amp3=[];
regularity.amp4=[];
regularity.alphalinha=[];
regularity.alpha1=[];
regularity.alpha2=[];
regularity.alpha3=[];
regularity.alpha4=[];  %%%%% para incluir a escala 5 %% Rute 22/05/02
regularity.thalpha1=[];
regularity.thalpha2=[];
regularity.thalpha3=[];
regularity.th2=[];
regularity.th3=[];
regularity.th4=[];

regularity.thlinha=[];
regularity.escolhidos=[];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

warning off
ultimo_anot=1;
timeoffset=0;
leadsfile = [];
nsamp = 2^16;   % Number of samples per excerpt for the WT

if ft==0   % MIT format files
    if isunix,          %%%%%%% Rute 13/03/02
        sep = '/';
    else sep = '\';     %%%%%%% Rute 13/03/02
    end                 %%%%%%% Rute 13/03/02

    heasig = readheader([headir sep ecgnr '.hea']);

    %%%%%%%%%%%%%%  ANA  PARA BIOPAC %%%%%%%%%%%%%%%%%%%%%%
    if isfield(heasig,'spf')
        heasig.freq=heasig.freq*heasig.spf(lead);
        heasig.nsamp=heasig.nsamp*heasig.spf(lead);
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if (heasig.fmt(lead)==16) || (heasig.fmt(lead)==212) || (heasig.fmt(lead)==61)
        formato = num2str(heasig.fmt(lead));
    else
        error('This format is not supported by the program')
    end

    for kk=1:heasig.nsig
        leadsfile(kk)=strcmp(heasig.fname(lead,:),heasig.fname(kk,:));
    end
    leadsfile = cumsum(leadsfile);
    leadinfile = leadsfile(lead);  % this is necessary when there are leads in different files.   e.g. ptbdb database  JP 2006
    nsiginfile = leadsfile(end);

    if strcmp(formato,'16')||strcmp(formato,'61'),
        % fid = fopen([sigdir ecgnr '.dat'],'rb');
        % if fid == -1,
        fid = fopen([sigdir heasig.fname(lead,:)],'rb');
        % end
        fseek(fid,0,-1);  % Rewind the file
        if strcmp(heasig.fname(lead,1),'_'),  % Siemens card recordings with MIT-type header
            timeoffset=512;
            fseek (fid, timeoffset,-1); % offset
        end
    end
elseif ft==1  % LUND format files
    hdsig=gethdsig(sigdir,ecgnr); % Reading header information
    heasig = hdsig2heasig(hdsig);       
    freq = heasig.freq;
elseif ft==2  % data in a Matlab file
    aux=[sigdir ecgnr];
    load ([aux '.mat'])
    heasig.nsamp=length(sinal);
    heasig.freq=fa;
    freq = heasig.freq;
    if exist('gain','var')
        heasig.gain=gain;
    else
        heasig.gain=200;  % When heasig.gain = 0 => default 200
    end
else
    error('Format ft should be 0, 1 or 2')
end
heasig.gain(heasig.gain==0)=200;  % When heasig.gain = 0 => default 200

if nargin >=8,
    t=[max(t(1),1) min(heasig.nsamp,t(2))];
else
    qrs_flag=0;% default value was missing 01/04/02 Rute
    t=[1 heasig.nsamp];
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 ecgnr '.' aname],'file')
                s=readannot([dirann ecgnr '.' 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 '.mat'],'file')
                s=load ([dirann aname]);
                kk=getfield(s);
                eval(['ext_anot=s.' kk ';']);
            else error('QRS annotation file not found');
            end
        otherwise
            error('Bad annotation format for external QRS annotation file');
    end
end

inisamp = t(1);
endsamp = 0;
timeqrs = [];
lastqrs = 1;


if qrs_flag==0 % estimated maximum number of beats
    maxlength = (t(2)-t(1))/heasig.freq*2;
else maxlength=length(ext_anot);
end

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

nanvec = nan*ones(1,maxlength);
position=struct('Pon',nanvec,'P',nanvec,'Poff',nanvec,'QRSon',nanvec,...
    'Q',nanvec,'R',nanvec,'Fiducial',nanvec,'qrs',zeros(1,maxlength),...
    'Rprima',nanvec,'S',nanvec,'QRSoff',nanvec,'Ton',nanvec,'T',nanvec,...
    'Tprima', nanvec,'Toff',nanvec,'Ttipo',nanvec);

if qrs_flag==1, position.qrs=ext_anot;  end


numlatdet = 0;      % Number of detected beats until now

[q1,q2,q3,q4,q5] = qspfilt5(heasig.freq);  % WT equivalent filters %%%%% para incluir a escala 5 %% Rute 22/05/02
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); %%%%% para incluir a escala 5 %% Rute 22/05/02
begoverlap = l5+2*heasig.freq;    % l5 + 2 sec. %%%%% passa a ser em relaçao a escala 5 %% Rute 22/05/02
endoverlap = d5;                                %%%%% passa a ser em relaçao a escala 5 %% Rute 22/05/02
% The two seconds overlap makes sure that the last/first beat
% are completely within the excerpt processed

samp=1; %freq;    %Initialization of samp (1 sec)

while ((endsamp+1) < t(2))
    firstnewsamp = samp(end);   % last sample analyzed
    endsamp = min(inisamp + nsamp -1,t(2));

    if ft==0 % MIT format
        if strcmp(formato,'212'),  % !!!!! Different formats (heasig)
            sig = rdsign212([sigdir ecgnr '.dat'],nsiginfile,inisamp, endsamp);  % I changed heasig.nsig by nsiginfile  ---> see above JP2006
        elseif strcmp(formato,'16')|strcmp(formato,'61'),

            %%%%%%%%%%%%%%  ANA  PARA BIOPAC %%%%%%%%%%%%%%%%%%%%%%
            if isfield(heasig,'spf')
                if(~isempty(heasig.spf) & find(heasig.spf ~= 1) )
                    agrupacion = sum(heasig.spf);
                    %         ecg = zeros(heasig.nsig,heasig.nsamp*(max(heasig.spf)));
                end
                fid = fopen([sigdir ecgnr '.ecg'],'rb');
                %      fid = fopen([prmBrowse.ecgdir prmBrowse.ecgnr],'rb');
                if(length(find(heasig.fmt == 16)) == heasig.nsig)
                    %fseek(fid,agrupacion*(inisamp-1)*2,'bof');
                    fseek(fid,agrupacion*round((inisamp-1)/heasig.spf(lead))*2,'bof');
                    x = fread(fid,[agrupacion (endsamp-inisamp+1)/heasig.spf(lead)], 'int16');
                    ini_pos_lead = sum(heasig.spf(1:lead))-heasig.spf(lead)+1;
                    end_pos_lead = sum(heasig.spf(1:lead));
                    sig = reshape(x(ini_pos_lead:end_pos_lead,:),[],1);
                    % Para no tener que cambiar el codigo
                    sig=[zeros(length(sig),leadinfile-1) sig];
                    clear('x')
                end
                fclose('all');
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

            else
                fseek(fid, timeoffset+2*(inisamp-1)*nsiginfile, -1); % Locate the pointer
                sig = fread(fid,[nsiginfile endsamp-inisamp+1] ,'int16')';
            end
            if (strcmp(computer,'SOL2')&strcmp(formato,'16'))|(~strcmp(computer,'SOL2')&strcmp(formato,'61')),
                sig=swap16(sig);
            end
        end
        sig = sig(:,leadinfile);   % Only selected lead.  --> JP 2006  I changed lead by leainfile... see at the top

        if exist('heasig','var') && isfield(heasig,'units')
            if upper(heasig.units(lead))=='MV'
                sig = (sig - heasig.adczero(lead))/heasig.gain(lead)*1e3; % conversion to microV
            elseif upper(heasig.units(lead))=='V'
                sig = (sig - heasig.adczero(lead))/heasig.gain(lead)*1e6; % conversion to microV
            else
                sig = (sig - heasig.adczero(lead))/heasig.gain(lead);
            end
        else
            sig = (sig - heasig.adczero(lead))/heasig.gain(lead)*1e3; % conversion to microV

        end

        %         sig = (sig - heasig.adczero(lead))/heasig.gain(lead)*1e3;  % conversion to microV
    elseif ft==1  % LUND format
        sig=getsig(sigdir,ecgnr,[inisamp,endsamp],lead)';
        %sig=sig/5;   % Mejor en uV JP
    elseif ft==2  % data in a Matlab file
        sig=sinal(lead,inisamp:endsamp)';
    end

    wt = wavt5(sig',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
    w = zeros(size(wt,1)-d5,5);
    w(:,5) = wt(d5+1:end,5);
    w(:,4) = wt(d4+1:end+d4-d5,4);
    w(:,3) = wt(d3+1:end+d3-d5,3);
    w(:,2) = wt(d2+1:end+d2-d5,2);
    w(:,1) = wt(d1+1:end+d1-d5,1);
    sig=sig(l5:end-d5); % piece of signal processed in this iteration

    clear wt;
    % Initialization of the thresholds for QRS detection (eps)
    eps= 0.5*sqrt(mean(w.^2));
    eps(4) = eps(4)*2;

    if qrs_flag==0  % The QRS fiducial point is calculated with wavedet
        %fiducial_altera_novo;  % uso de Criterios alternativos para selecao de candidatos a QRS % Rute 4/9/02
        fiducial; % Detection of fiducial point % 22/05/02 Rute including scale 5 % 07/06/02 Rute exclui escala5
    else           % The QRS fiducial point is read from external annotator
        first = max(firstnewsamp-ceil(heasig.freq), samp(1)-1+ceil(heasig.freq*0.050));
        % first is calculated according with the first lines of fiducial JP

        sel=find(position.qrs>=first & position.qrs<=samp(end));
        time=position.qrs(sel)-samp(1)+1;


        % If the P-wave is not in the present segment, discard the beat, because the
        % the beat should have been detected in the last segment.
        if abs(samp(time(1)) - lastqrs)< 0.275*heasig.freq,
            time(1)=[];  sel(1)=[];              % In order not to repeat beats
        end

        if (length(sig)-time(end) < 0.45*heasig.freq),
            time(end)=[]; sel(end)=[];
        end
        % If last beat's T-wave is not in 'sig', the beat is detected in next one

        if time(1)-0.35*heasig.freq<=0, time(1)=[]; sel(1)=[]; end

        timeqrs = [timeqrs samp(time)];  % QRS times

        lastqrs = samp(time(end));

        intervalo = [sel(1) sel(end)];       %%%%%%
        % numeration of the beats processed in this segment
        %keyboard
    end
    if ~isempty(time)
        %% Individual wave detection  %%
        qrswave;
        %%%%%%%%%%%%%%%% caso em que rr tem dimensao inferior ao necessario
        if length(time)>1  %%%%%%%%%%%%%%%%introduzido a 17/05/02 Rute
            twave,  %%%%%%%% 07/06/0 Rute including scale 5
        end %%%%%%%%%%%%%%%%introduzido a 17/05/02 Rute
        pwave;
        % Last annotated position
        if ~isnan(position.Toff(intervalo(2))),
            ultimo_anot = position.Toff(intervalo(2));
        elseif ~isnan(position.Tprima(intervalo(2))),
            ultimo_anot = position.Tprima(intervalo(2));
        elseif ~isnan(position.T(intervalo(2))),
            ultimo_anot = position.T(intervalo(2));
        elseif ~isnan(position.QRSoff(intervalo(2))),
            ultimo_anot = position.QRSoff(intervalo(2));
        else
            ultimo_anot = position.qrs(intervalo(2));
        end
    end

    inisamp = endsamp+2-endoverlap-begoverlap;
end


% Remove void annotations at the end
aux = find (position.qrs == 0);
position.Pon(aux)=[];
position.P(aux)=[];
position.Poff(aux)=[];
position.QRSon(aux)=[];
position.Q(aux)=[];
position.R(aux)=[];
position.Fiducial(aux)=[];
position.qrs(aux)=[];
position.Rprima(aux)=[];
position.S(aux)=[];
position.QRSoff(aux)=[];
position.Ton(aux)=[];
position.T(aux)=[];
position.Tprima(aux)=[];
position.Toff(aux)=[];
position.Ttipo(aux)=[];


annstruct = pos2ann(position);
if ~isempty(position.qrs)
    writeannot ([matdir ecgnr '.' anot 's' num2str(lead)], annstruct);   % mexfile, faster!
end

%escribeanot(fann,position,anot,lead, ecgnr, freq, matdir, 'anot');

fclose('all');
clear global regularity