function [sQRS,chNb] = select_channel(ecg,ecgType,varargin)
% this function automates the channel selection onto
% which to perform sQRS detection for whether to get MQRS or FQRS. 
% In the case of FQRS extraction, the key ideas are to to use the regularity 
% of the extracted HR to decide what FECG channel to choose and to use bxb 
% (proportion of matching sQRS btw two sQRS time series) to check that we 
% are not picking up the MECG.
%
% inputs
%   ecg:                the preprocessed FECG matrix
%   ecgType:            'FECG' or 'MECG'
%   varargin
%       MQRS:           MQRS position (in ms)
%       paramStruct:    structure containing challenge parameters
% outputs
%   sQRS:               selected QRS location (in ms)
%   chNb:               selected ABD channel
%
% FECG extraction 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.

% == managing inputs
MQRS = [];
paramStruct.detThres = 0.5;
paramStruct.fs = 1000;
paramStruct.qrsdet_wind = 15;
debug = 0;

switch nargin
    case 3
        MQRS = varargin{1};
    case 4
        MQRS = varargin{1};
        paramStruct = varargin{2};
    case 5
        MQRS = varargin{1};
        paramStruct = varargin{2};
        debug = varargin{3};
    otherwise
        error('select_channel: wrong number of input arguments \n');
end

if size(ecg,2)>size(ecg,1)
    ecg = ecg';
end

NB_AECG = size(ecg,2); % nb of ABD ecg
QRSce = cell(NB_AECG,1); % sQRS cell
sQRS = cell(1,1);
SMI = zeros(NB_AECG,1); % smoothness indices
bxb = zeros(NB_AECG,1); % matching btw MQRS and another QRS time series

try
    % == core function
    if strcmp(ecgType,'FECG')
        % == FECG channel selection
        for ch=1:NB_AECG
            % Normalize residuals to avoid big bumps that crashed the FQRS detector (this showed improvement)
            [ecgn,~] = normalise_ecg(ecg(:,ch),1,5*paramStruct.fs);

            % sQRS detection
            QRSce{ch} = run_qrsdet_by_seg(ecgn,paramStruct.fs,paramStruct.qrsdet_wind,paramStruct.detThres,ecgType);
            
            % smooth the time series (physiological)
            if paramStruct.smoothRR
                QRSce{ch} = smooth_rr(QRSce{ch});
            end
            
            % compute smoothness indice
            [SMI(ch),~] = assess_regularity(QRSce{ch},1);

            % check that the FQRSs are not matching the MQRSs
            bxb(ch) = bsqi(MQRS/paramStruct.fs,QRSce{ch}/paramStruct.fs);
        end

        % if time series too similar to the one from the mother then we are
        % probably picking up the mother residuals. Then put SMI to inf to
        % not select this channel
        SMI(bxb>0.4) = Inf;

        % in order to select the ABD channel onto which the FQRS will be
        % identified we look at SMI (lower SMI -> most likely to be a good detection)
        [~,ind_std_trans] = sort(SMI);
        chNb = ind_std_trans(1);
        sQRS = QRSce{chNb};

        % == plots
        if debug || isempty(nargout)
            tm = 1/paramStruct.fs:1/fs:length(ecg)/paramStruct.fs;
            [m,n] = size(ecg);
            ax = zeros(m,1);
            for pp=1:n
                ax(pp) = subplot(ceil(n/2),2,pp); plot(tm,ecg(:,pp)); 
                legend(strcat('SMI:',num2str(SMI(pp)),' bxb:',num2str(bxb(pp))));
                linkaxes(ax,'x');
            end
        end

        
        
        elseif strcmp(ecgType,'MECG')
             % == MECG channel selection
             % this step is performed to get more robust MQRS detection.
             % This is aimed at finding the MQRS with the best accuracy but
             % ALSO to distinguish between MQRS and FQRS time series since
             % in some instances the FECG is of higher or similar amplitude 
             % that the MECG. In practice a chest electrode should be used
             % to avoid this to happen as the chest signal will be of
             % better quality and will ensure that there is no confusion
             % between MQRS and FQRS time series.
             [~,icaECG] = jade(ecg); 
             ecg = [ecg';icaECG]';
             NB_AECG = 2*NB_AECG; 
             lengthQRS  = zeros(4,1);

             for i=1:NB_AECG
                % run sQRS detector
                sQRS{i} = run_qrsdet_by_seg(ecg(:,i),paramStruct.fs,10,0.6,ecgType);
                lengthQRS(i) = length(sQRS{i});
                
                % compute SMI of MHR at 96% CI
                [SMI(i),~] = assess_regularity(sQRS{i});
             end

             if isempty(find(SMI~=inf,1))
                [~,chNb] = min(lengthQRS);
                sQRS = sQRS{chNb};
             else
                % this set of rules is specific to the PCinC2013
                % challenge and should not be kept for practical
                % applications
                indLowSMI = find(SMI<10);
                if length(indLowSMI)>1
                    [mi,indMin] = min(lengthQRS(indLowSMI));
                    [ma,indMax] = max(lengthQRS(indLowSMI));
                    if ma-mi>7
                        % if there is no big gap in length btw the 'good' 
                        % time series then take the most regular
                        chNb = indLowSMI(indMin);
                    else
                        % if there is the gap then we might be detecting
                        % both the MQRS and FQRS in which case it is safer
                        % to assume that the MQRS is the shortest time
                        % series so take the longuest for the foetus
                        chNb = indLowSMI(indMax);
                    end
                else
                    [~,chNb] = min(SMI);
                end
                sQRS = sQRS{chNb};
             end  
    else          
            disp('Wrong entry \n');
    end

catch ME
    for enb=1:length(ME.stack); disp(ME.stack(enb)); end;
    sQRS = [20 400 600 1000]; chNb=1;
end

end


% alternative to bSQI is to use the Physionet bxb function but this requires converting 
% .mat annotations to WFDB format. This is unfortunately taking a lot of computational time!!
% (0.116627sec)
%
%         mat2wfdbanns('qrs','fqrs',QRSce{ch},paramStruct.fs);
%         mat2wfdbanns('qrs','mqrs',MQRS,paramStruct.fs);
%         bxb('qrs','mqrs','fqrs','bxbReport.txt','00:00:00','00:01:00');
%         testFil = importdata('bxbReport.txt');
%         delete('bxbReport.txt');
%         bxb(ch) = testFil.data(1,1)/(testFil.data(1,1)+testFil.data(6,1));