function [fetal_QRSAnn_est,QT_Interval,best_index] = 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

warning off;

% ---- code constants and parameters----

fs = 1000;             % sampling frequency
blocksize = 2000;
sample = load('samplewaves.mat');

% ---- preprocessing and maternal peak detection----

[nleads,nsamples]=size(ECG);

if nleads>nsamples % make sure ECG is a row matrix
    ECG = ECG';
    [nleads,nsamples]=size(ECG);
end


I=MissedValCheck(ECG);

Y=preprocessing(I,fs,0.05); % filtered ECG

[Y,accY,peaks,goodint] = MECGPeakDetect(Y,blocksize,sample.m,fs);

ind1=find(peaks==1);
if length(ind1)<20 || median(abs(accY(ind1)))<=50
    Y=preprocessing(I,fs,0.2);
    [Y,accY,peaks,goodint] = MECGPeakDetect(Y,blocksize,sample.m,fs);
end

%----- finding principal componenets-----

[pcM, Ucomp] = findPrinComp(Y,blocksize);

% ---- MECG cancellation ----
L=1;
mcount=1;
[Mout,Fout,corrcoef]=FECG_extract_multichan(Y,fs,peaks,L);
while (max(corrcoef)>=0.2 && mcount<4)
    [Mout,Fout,corrcoef]=FECG_extract_multichan(Fout,fs,peaks,L);
    mcount=mcount+1;
end
% ---- Post processing and channel selection ----

[FECG,post_order] = postFilterSelect(Fout,blocksize);

% ---- FQRS detection ----
avg_power=median(abs(FECG'))';
FECG_norm=FECG./repmat(avg_power,1,nsamples);

[factorcheat,best_choice]=SigMerge(FECG_norm,sample.f,goodint,fs);

if any(factorcheat)
    [fpr1,Ftest_round1,fsample,findnc1,findc1]=Fetal_Peakdetection2(FECG_norm,sample.f,fs,factorcheat,0);
    
    
    if sum(sample.f.*fsample)<0
        [fpr1,Ftest_round1,fsample,findnc1,findc1]=Fetal_Peakdetection2(FECG_norm,sample.f,fs,-factorcheat,0);
    end

else
    findc1=0;
    fsample=zeros(size(sample.f));
end

% comparing the detected peaks with those detected from principal
% components

best_index = 1;
if (length(findc1)<85)

    fkurt=zeros(1,3);
    fkurt(1)=3; % this is to avoid canceling fpr1 in the competition at this level. 
    
    p_order = SetSelect(pcM(1:3,:),fs);
        
    [Ftest,fpr2,findc2,psample2] = PCMPeakDetect(pcM(p_order(1),:),peaks,sample.f,goodint,fs,blocksize);
    fkurt(2) = kurtosis(psample2);
    
    u_order = SetSelect(Ucomp,fs);
    
%     [Ftest,fpr3,findc3,psample3] = PCMPeakDetect(Ucomp(u_order(1),:),peaks,sample.f,goodint,fs,blocksize);
%     fkurt(3) = kurtosis(psample3);
    [Ftest,fpr3,findc3,psample3] = PCMPeakDetect(Ucomp(u_order(2),:),peaks,sample.f,goodint,fs,blocksize);
    fkurt(3) = kurtosis(psample3);
    ll=[length(findc1) length(findc2) length(findc3)];
    
    % selecting the best answer
    [fkurt_max, best_index] = max(ll.*(ll<140).*(fkurt>2.8));
end

switch best_index
    case 1
        fpr=fpr1;
        findc=findc1;
    case 2
        fpr=fpr2;
        findc=findc2;
    case 3
        fpr=fpr3;
        findc=findc3;
end

    


fetal_QRSAnn_est = fpr-1;
QT_Interval = 0;

end