function [FQRSout,QTout,outExtra] = physionet2013extract(tm,ecg,paramStruct)
% 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
% [FQRSout,QTout] = physionet2013(tm,ecg) where the inputs and outputs are specified
% below.
%
% inputs
%   tm :         Nx1 vector of time in milliseconds
%   ecg:         4x60000 (4 channels and 1min of signal at 1000Hz) matrix of
%                abdominal ecg channels.
%   paramStruct: structure with the key parameters for the different
%                algorithms. if empty then the defaults parameters are
%                used. 
% 
% output
%   FQRSout:    FQRS markers in seconds. Each marker indicates the position of one
%               of the FQRS detected by the algorithm.
%   QTout:      1x1 estimated fetal QTout duration
%   CHout:      the channel nb that was selected
%
%
%
% Safe Foetus Monitoring Toolbox, version 1.0, Sept 2013
% Released under the GNU General Public License
%
% Copyright (C) 2013 Joachim Behar
% Oxford university, Intelligent Patient Monitoring Group - Oxford 2013
% joachim.behar@eng.ox.ac.uk
%
% Last updated : 02-08-2013
%
% This program is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by the
% Free Software Foundation; either version 2 of the License, or (at your
% option) any later version.
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
% Public License for more details.

% == check the inputs
if nargin<2; error('physionet2013extract: wrong number of input arguments \n'); end;
if size(ecg,2)>size(ecg,1); ecg = ecg'; end;


% == entering core of the function
try      
    % == allocate space    
    NB_ABD_CHANNELS = size(ecg,2);
    MQRS = cell(1,NB_ABD_CHANNELS);
    %LogFileID = fopen('log.txt','a+');
    
    % == preprocessing
    ecgf = preprocessing([flipdim(ecg(1:paramStruct.fs*1-1,:),1)' ecg'],paramStruct);
    ecgf = ecgf(paramStruct.fs*1:end,:); % for border filtering reasons..
    
    if ~paramStruct.extract_qt
        % ===========================================================
        % ================= FECG extraction block ===================
        % ===========================================================
        % == MQRS detection
        [rMQRS,~] = select_channel(ecgf,'MECG',[],paramStruct);
        
        tic;
        for cc=1:NB_ABD_CHANNELS
            %fprintf(LogFileID,'ch %d: ',cc);
            % == adjust MQRS location
            MQRS{cc} = adjust_mqrs_location(ecgf(:,cc),rMQRS,paramStruct.fs,0); %MQRS{cc} = rMQRS;
            % == now remove MECG
            if find(strcmp(paramStruct.method,{'TS','TS-SUZANNA','TS-CERUTTI','TS-LP','TS-PCA'}))
                % TS techniques
                fecg.m1(cc,:) = mecg_cancellation(MQRS{cc},ecgf(:,cc)',paramStruct.method,paramStruct.TSnbCycles,paramStruct.NbPC);
            elseif strcmp(paramStruct.method,'TS-EKF')
                % EKF
                [fecg.m1(cc,:),errorInit,~] = run_kalman(ecgf(:,cc),MQRS{cc},paramStruct);
                %fprintf(LogFileID,' EKF errorInit=%d',round(100*errorInit)/100);
            elseif strcmp(paramStruct.method,'FUSE')
                % FUSE
                fecg.m2(cc,:) = mecg_cancellation(MQRS{cc},ecgf(:,cc)','TS',paramStruct.TSnbCycles);    % can also be TS/PCA
                % fecg.m3(cc,:) = mecg_cancellation(MQRS{cc},ecgf(:,cc)','TS',paramStruct.TSnbCycles);
                % [fecg.m4(cc,:),~,~] = run_kalman(ecgf(:,cc),MQRS{cc},paramStruct);
                % plot(ecgf(:,cc)); hold on, plot(fecg.m2(cc,:),'r')
            end        
            %fprintf(LogFileID,' \n'); 
        end

        if strcmp(paramStruct.method,'PCA')
            % PCA
            [~,fecg.m1] = princomp(ecgf); fecg.m1=fecg.m1';
        elseif strcmp(paramStruct.method,'ICA')
            % ICA
            [~,fecg.m1] = jade(ecgf);
            %fecg.m1 = fastica(ecgf');
        elseif strcmp(paramStruct.method,'DEFLATION')
            % ICA as preprocessing
            [~,fecg.m1] = jade(ecgf);        
        elseif strcmp(paramStruct.method,'FUSE')
            %[~,fecg.m4] = princomp(ecgf); fecg.m4=fecg.m4';
            [~,fecg.m1] = jade(ecgf);
        end
        %fprintf(LogFileID,paramStruct.method);

        % == BSS step
        if find(strcmp(paramStruct.method,{'TS','TS-SUZANNA','TS-CERUTTI','TS-LP','TS-PCA','TS-EKF'})) ...
                & strcmp(paramStruct.methodComplexity,'COMPLEX')
            % TS/PCA/EKF: add a BSS step on the residuals
            ecgout  = bss_step(fecg.m1,ecgf,rMQRS,paramStruct);
            fecg.m1 = [ecgout;fecg.m1];
        elseif strcmp(paramStruct.method,'DEFLATION')
            ecgout  = bss_step([],fecg.m1,rMQRS,paramStruct);
            fecg.m1 = ecgout;    
        elseif strcmp(paramStruct.method,'FUSE')
            ecgout1 = bss_step([],fecg.m1,rMQRS,paramStruct);
            paramStruct.method = 'TS';
            ecgout2 = bss_step(fecg.m2,fecg.m1,rMQRS,paramStruct);        
            fecg.m1 = [ecgout1;ecgout2;fecg.m1;fecg.m2];       % [ICA-TS/ICA-TS-ICA] [TS-ICA] [ICA] [TS]
        end
        timeExtrac = toc; fprintf('Execution time FECG extraction: %f \n',timeExtrac);
        
        % == channel selection & FQRS detection & post processing
        [FQRStemp,CHout] = select_channel(fecg.m1,'FECG',rMQRS,paramStruct); 
        FQRS = FQRStemp;

        fprintf('Selected channel number is: %d \n',CHout);
        %fprintf(LogFileID,'SelectFECGCha: %d ',CHout);

        % ===========================================================
        % ================= Challenge twicks ========================
        % ===========================================================
        % NOTE: Consider removing parts of that for an actual application!
        outExtra.rFQRS = FQRS;
        outExtra.rMQRS = rMQRS;
        
        if paramStruct.cinc_match
            [STD,~] = assess_regularity(FQRS,1,paramStruct.fs,1);  
            fprintf('STD: %f \n',STD);
            if length(FQRS)<85 || length(FQRS)>200
                % identify rubbish! In the case the time
                % series is rubbish then just create one at 134bpm
                disp('abnormal length(FQRS)<85 || length(FQRS)>200');
                FQRSout = 0.1:0.42:60;
                FQRSout = round(FQRSout*1000);
                QTout = 0;
            elseif STD>paramStruct.STDcut
                disp('abnormal STD');
                % if time series that we get is too irregular then better to look 
                % for predominant HR mode and guess the time series!
                [FQRSmode,m1,m2] = guess_fqrs(FQRS,rMQRS,paramStruct.fs);
                %m = mode(60./diff(FQRS./paramStruct.fs)); 
                if m1>125 && m1<170 && m2>125 && m2<175
                    disp('- using MODE');
                    FQRSout = FQRSmode;
                else
                    disp('- just guessing');
                    FQRSout = 0.1:0.42:60; % cc.e. 143bpm
                end
                FQRSout = round(FQRSout*1000);
                QTout = 250;
            else
                FQRSout = FQRS;
                QTout = 250;

                % then add the magic smoothing step while possibly remove first
                % and last peaks
                FQRSout = smooth_rr_lots(FQRSout(FQRSout>paramStruct.fs*0.200 & ...
                FQRSout<paramStruct.fs*59.9)/paramStruct.fs,paramStruct.fs);
            end
        else
            FQRSout = FQRS;
            QTout = 250;
        end
        % ===========================================================
        
        
    else
        % ===========================================================
        % ================= QT extraction ===========================
        % ===========================================================

        % in previous step need to output rMQRS and rFQRS
        MED_RR = median(diff(paramStruct.rFQRS));
        NB_ITER = 5;
        Qtrans = zeros(NB_ABD_CHANNELS*NB_ITER,1);
        Ttrans = zeros(NB_ABD_CHANNELS*NB_ITER,1);
        compt = 0;
        tic
        for cc=1:NB_ABD_CHANNELS
            fprintf('QT processing channel: %f \n',cc);
            % adjust MQRS
            MQRS{cc} = adjust_mqrs_location(ecgf(:,cc),paramStruct.rMQRS,paramStruct.fs,0);
            % source separation
            fecg.m1(cc,:) = mecg_cancellation(MQRS{cc},ecgf(:,cc)',paramStruct.method,paramStruct.TSnbCycles,paramStruct.NbPC);
            % template construction
            relevantMode = fecg_template_construction(paramStruct.rFQRS,fecg.m1(cc,:));
            for kk=1:NB_ITER
                compt = compt+1;
                [~,Qtrans(compt),Ttrans(compt),~] = measure_qt(relevantMode.cycleMean,MED_RR);
            end
        end
        
        Qtrans2 = Qtrans;
        Ttrans2 = Ttrans; 
        
        Qtrans2(Qtrans==0 | Qtrans==1)=[]; % remove when it failed
        Ttrans2(Ttrans==0 | Ttrans==1)=[];
        
        Q = median(Qtrans2);
        T = median(Ttrans2);
        QTout = T-Q;
                
        if QTout== 0|| QTout>320 || QTout<125; QTout=250; end;

        CHout = 1; % how to select a good channnel for plot??
        FQRSout=[]; outExtra=[];
        toc
        % ===========================================================
    end
    
    
    % ==================== plots for analysis ===================
    if paramStruct.debug==1 && paramStruct.extract_qt==0
        figure(1);
        title(strcat('Channel selected is: ',CHout));
        one_ecgres = fecg.m1(CHout,:);
        hrv = 60./(diff(FQRS)/paramStruct.fs);
        hrvss = 60./(diff(FQRSout)/paramStruct.fs);
        ax(1) = subplot(311); 
            plot(tm,ecg(:,1),tm(MQRS{1}),ecg(MQRS{1},1),'+r','LineWidth',2);
            plot(tm,ecgf(:,1),tm(MQRS{1}),ecgf(MQRS{1},1),'+r','LineWidth',2);
            legend('One ABD ecg','MQRS');
        ax(2) = subplot(312);
            hold on, plot(tm,one_ecgres,'k','LineWidth',2); 
            hold on, plot(tm(FQRS),one_ecgres(FQRS),'+r');
            legend('fecg (residual)','FQRS detected');
        ax(3) = subplot(313); plot(FQRS(1:end-1)/paramStruct.fs,hrv,'-ro');
            hold on, plot(FQRSout(1:end-1)/paramStruct.fs,hrvss,'--ko');
            legend('FHRc','FHRcc'); xlabel('Time (sec)'); ylabel('Amplitude (NU)');
        linkaxes(ax,'x'); set(findall(gcf,'type','text'),'fontSize',14,'fontWeight','bold');
    elseif paramStruct.debug==1 && paramStruct.extract_qt==1
        % extract the template ecg closest from T point identified (rather than plotting a random one..)
        k = dsearchn(Ttrans-Qtrans,QTout);
        ch = k+1-floor(k/NB_ABD_CHANNELS)*NB_ABD_CHANNELS;
        relevantMode = fecg_template_construction(paramStruct.rFQRS,fecg.m1(ch,:));

        figure(2);
        fprintf('QT: Number of cycles: %f \n',relevantMode.NbCycles);
        one_ecgres = fecg.m1(CHout,:);
        title('QT analysis');
        subplot(211);
            hold on, plot(tm,one_ecgres,'k','LineWidth',2); 
            hold on, plot(tm(paramStruct.rFQRS),one_ecgres(paramStruct.rFQRS),'+r');
            legend('fecg (residual)','FQRS detected');
        subplot(212); 
            plot(relevantMode.cycleMean,'LineWidth',2);
            QTpos = round([Q T]*250/MED_RR);
            if QTpos(1)>0 && QTpos(2)<250
                hold on, plot(QTpos,relevantMode.cycleMean(QTpos),'+r','LineWidth',3);
            else
                hold on, plot(QTpos,0,'+r','LineWidth',3);
            end
            xlabel('Phase Wrap'); ylabel('Amplitude [NU]');
    end
    % ===========================================================
    
    
%fclose(LogFileID);    
catch ME
    for enb=1:length(ME.stack); disp(ME.stack(enb)); end;   
    FQRSout = 0.1:0.42:60;
    QTout = 250;
end

end