function [fetal_QRSAnn_est,QT_Interval] = physionet2013(tm,ECG)
% Template algorithm for Physionet/CinC competition 2013. This function can
% be used for events 1 and 2. Participants are free to modify any
% components of the code. However the function prototype must stay the
% same:
%
% [fetal_QRSAnn_est,QT_Interval] = physionet2013(tm,ECG) where the inputs and outputs are specified
% below.
%
% inputs:
%   ECG: 4x60000 (4 channels and 1min of signal at 1000Hz) matrix of
%   abdominal ECG channels.
%   tm : Nx1 vector of time in milliseconds
% output:
%   FQRS: FQRS markers in seconds. Each marker indicates the position of one
%   of the FQRS detected by the algorithm.
%   QT_Interval:   1x1 estimated fetal QT duration (enter NaN or 0 if you do wish to calculate)
%
%
% Author: Joachim Behar - IPMG Oxford (joachim.behar@eng.ox.ac.uk)
% Last updated: March 3, 2013 Ikaro Silva


% ---- check size of ECG ----
if size(ECG,2)>size(ECG,1)
    ECG = ECG';
end

fs      = 1000;             % sampling frequency
N       = size(ECG,2);      % number of abdominal channels
debug   = 0;  % enter debug mode?

ECG(find(isnan(ECG)==1)) = 0;
x=ECG(:,4);
x=x'; 
z=wavelet(x);
points=60000;  
level=5;    num_inter=6;   wf='db5';  fs=1000;
clr={'k-','k-','k-','k-','b-'};

[Lo_D,Hi_D,Lo_R,Hi_R]=wfilters(wf);
[swa,swd] = swt(z,level,Lo_D,Hi_D);
if  debug 
   figure(1);
   subplot(level,1,1); plot(real(z)); grid on;axis tight;
   for i=1:level
       subplot(level+1,2,2*(i)+1);
       plot(swa(i,:)); axis tight;grid on;xlabel('time');
       ylabel(strcat('a   ',num2str(i)));
       subplot(level+1,2,2*(i)+2);
       plot(swd(i,:)); axis tight;grid on;
       ylabel(strcat('d   ',num2str(i)));
   end 
end

ddw=zeros(size(swd));
pddw=ddw;
nddw=ddw;
posw=swd.*(swd>0);
pdw=((posw(:,1:points-1)-posw(:,2:points))<0);
pddw(:,2:points-1)=((pdw(:,1:points-2)-pdw(:,2:points-1))>0);
negw=swd.*(swd<0);
ndw=((negw(:,1:points-1)-negw(:,2:points))>0);
nddw(:,2:points-1)=((ndw(:,1:points-2)-ndw(:,2:points-1))>0);
ddw=pddw|nddw;
ddw(:,1)=1;
ddw(:,points)=1;
wpeak=ddw.*swd;
wpeak(:,1)=wpeak(:,1)+1e-10;
wpeak(:,points)=wpeak(:,points)+1e-10;
if debug
figure(2);
subplot(level+1,1,1); plot(real(z)); grid on;axis tight;
   for i=1:level
       subplot(level+1,1,i+1);
       plot(wpeak(i,:)); axis tight;grid on;
   ylabel(strcat('j=   ',num2str(i)));
   end
end

hlevel_wpeak=wpeak(level,:);
selevel_wpeak=wpeak(4,:);
sig=z;
peaks = PeakDetection(hlevel_wpeak,1/fs);
selevel_wpeak2=selevel_wpeak;
 guilinglet=zeros(1,length(selevel_wpeak2));
 for j=1:length(peaks)
     guilinglet(peaks(j))=1;
 end
D4_p=(guilinglet~=0);
D4_p2=(guilinglet~=0);
O_d4=60;%该参数确定在上一级搜索极大值的范围，可以调整。
O_d42=60;
for P_d4=O_d4:(length(guilinglet)-O_d4);
    if D4_p(P_d4)==1; 
        for i=1:O_d4-1;
        D4_p(P_d4-i)=1;
        
        end ;
    end;     
end;
for P_d4=O_d42:(length(guilinglet)-O_d42);
    if D4_p2(P_d4)==1; 
        for i=1:O_d42-1;
        D4_p(P_d4+i)=1;
        
        end ;
    end;     
end;
selevel_wpeak2=selevel_wpeak2.*D4_p;

selevel_wpeak3=selevel_wpeak-selevel_wpeak2;

 fpeaks= PeakDetection(selevel_wpeak3,2/fs);
 if debug
 figure(10);
 plot(selevel_wpeak3);
xlim([0,10000]);
% hold on,grid on;
%  plot(QRSAnn1,sig(QRSAnn1)+(m-1)*2,'x','MarkerSize',7,...
%                          'MarkerEdgeColor','g','LineWidth',2);
      grid on;
   xlim([0,10000]);                
   hold on,plot(fpeaks,sig(fpeaks)+(m-1)*2,'o','MarkerSize',7,...
                         'MarkerEdgeColor','r','LineWidth',2);grid on;   
 end
 
                     
aa=1;
for i=1:length(fpeaks)
     if  fpeaks(i)>0
         refpeaks(aa)=fpeaks(i);
         aa=aa+1;
     end
end

sortlength=20;
sign= abs(max(z))>abs(min(z));
for j=1:length(refpeaks)
     maxorminpost=refpeaks(j);
     if sign 
        for tt=1:sortlength
            if  refpeaks(j)-tt<=0 
                break;
            else  
                if z(refpeaks(j)-tt)<z(maxorminpost)
                   maxorminpost=refpeaks(j)-tt;
                end 
            end
        end
        for tt=1:sortlength 
             if  refpeaks(j)+tt>=points
                 break;
             else
                 if   z(refpeaks(j)+tt)<z(maxorminpost)
                      maxorminpost=refpeaks(j)+tt;
                 end
             end
        end
     else
        for tt=1:sortlength 
            if  refpeaks(j)-tt<=0 
                break;
            else  
                if z(refpeaks(j)-tt)>z(maxorminpost)
                    maxorminpost=refpeaks(j)-tt;
                end 
            end
        end
        for tt=1:sortlength 
             if  refpeaks(j)+tt>=points
                 break;
             else
                 if   z(refpeaks(j)+tt)>z(maxorminpost)
                      maxorminpost=refpeaks(j)+tt;
                 end
             end
        end 
     end
    refpeaks(j) =maxorminpost;
end
if  debug
 figure(13);
 plot(sig + (m-1)*2,clr{m});
                xlim([0,10000])
%                 hold on;grid on
%      plot(QRSAnn,sig(QRSAnn)+(m-1)*2,'x','MarkerSize',7,...
%                          'MarkerEdgeColor','g','LineWidth',2);
      grid on;
   xlim([0,10000]);                
   hold on,plot(refpeaks,sig(refpeaks)+(m-1)*2,'o','MarkerSize',7,...
                         'MarkerEdgeColor','r','LineWidth',2);grid on; 
end
fetal_QRSAnn_est=refpeaks;
QT_Interval=0;

end