function qrsFes=FecgQRSfDet(X,fs,cName,QRSm,graph,dbFlag,saveFig,saveFigRRf,qrsAfs)
%---------------------------------------------------------------------------------------------
%   Fecg: Fetal QRS detection
%   - Fetal QRS pre-detection on each signal
%   - Best fetal ECG channel selection based on a priori information on fECG pseudo-periodicity
%   - Identification of a good interval
%   - Fetal QRS detection in forward direction starting from the end
%     and backward starting from the beginning of such interval
% --------------------------------------------------------------------------------------------
% Author: Maurizio Varanini, Clinical Physiology Institute, CNR, Pisa, Italy
% For any comment or bug report, please send e-mail to: maurizio.varanini@ifc.cnr.it
% --------------------------------------------------------------------------------------------
% 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 "as is" and "as available" in the hope that it will be useful,
% but WITHOUT ANY WARRANTY of any kind; without even the implied warranty of MERCHANTABILITY
% or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
% --------------------------------------------------------------------------------------------
if(nargin<3), cName=''; end
if(nargin<4), QRSm=[]; end
if(nargin<5), graph=1; end
if(nargin<6), dbFlag=0; end
if(nargin<7), saveFig=0; end
if(nargin<8), saveFigRRf=0; end
if(nargin<9), qrsAfs=0; end

fprintf('\n --------------------------------------------------------- \n');
[progpath, progname] = fileparts(which(mfilename));
fprintf('Program: %s,   record name: %s\n', progname, cName);

if(saveFig), figFmt='png';
    figPath=fullfile('../Figure/',progname);
    if(~exist(figPath,'dir')), mkdir(figPath); end
end

if(~isempty(qrsAfs));
    learning=1;   % learning set, annotations are available
else
    learning=0;  %test set
end

%-------------------------------------------------------------

% recording duration
[ndt, ns]=size(X);

time= [1:ndt]/fs;
if(dbFlag && learning), PlotSgnMrkNc(X, qrsAfs*fs, fs, [cName,': faQRS']); end

% -----
qrsM=QRSm(1,:)';
if(dbFlag), PlotSgnMrkNc(X, qrsM, fs, [cName,': mQRS']); end
qrsMs=qrsM/fs;
%--------------------------------------------------------------------------------------------
for is=1:ns,
    X(:,is)= (X(:,is)-mean(X(:,is)))/std(X(:,is));
end
% X(:,ns+1)= mean(X')';  % if ICA fails S/N improvment could be obtained
% ns=ns+1;

if(graph)
    figure, set(gcf,'Color','white');
    for is=1:ns,
        subplot(ns,1,is), plot(vtime,X(:,is));
        wgmi1= min(X(:,is)) -2;
        wgma1= max(X(:,is)) +2;
        ylim([wgmi1, wgma1]);
        set(gca,'YTick',[-5 0 5])
        % if(is~=ns), set(gca,'XTickLabel',''); end
        if(is==1), title([cName,': ICA fetal ECG']); end
    end
end

ecg=X;

% raw derivative filter coefficients
% 5ms before and after, 3ms in the middle
nu=ceil(0.005 * fs); nz=floor(0.0030*fs /2)*2 +1;  % nz=nearest odd value
B=[ones(nu,1);zeros(nz,1);-ones(nu,1)];
delay=floor(length(B)/2);
% ---   compute the derivative signal
ecgfx=[repmat(ecg(1,:),delay,1);ecg;repmat(ecg(end,:),delay,1)];
decgr=filter(B,1,ecgfx);   decgr= decgr(2*delay+1:end,:);
if(dbFlag)
    figure; freqz(B,1,1000,fs);
end
% ----- Butterworth forward and backward bandpass filtered (1-6Hz)
fmind=0.7; fmaxd=8;
Wn = [fmind, fmaxd]/(fs/2);  % normalized cut-off frequency (0,1)
[b,a]= butter(1,Wn);
decgBp=filtfilt(b,a,decgr);
% -----   absolute derivative ------
adecg=abs(decgr);
if(graph)
    figure, set(gcf,'Color','white');
    wgmima= mimaxscG(adecg,0,0,0.1);
    for is=1:ns,
        subplot(ns,1,is), plot(vtime,adecg(:,is));
        ylim(wgmima);
        % if(is~=ns), set(gca,'XTickLabel',''); end
        if(is==1), title([cName,': ECG absolute derivative']); end
    end
end
% ----- Butterworth forward and backward bandpass filtered (1-6Hz)
fmind=0.7; fmaxd=8;
Wn = [fmind, fmaxd]/(fs/2);  % normalized cut-off frequency (0,1)
[b,a]= butter(1,Wn);
adecgBp=filtfilt(b,a,adecg);

if(graph)
    figure, set(gcf,'Color','white');
    wgmima= mimaxscG(adecgBp(fs:end-fs),0,0,0.1);
    for is=1:ns,
        subplot(ns,1,is), plot(vtime,adecgBp(:,is));
        ylim(wgmima);
        % if(is~=ns), set(gca,'XTickLabel',''); end
        if(is==1), title([cName,': ECG filtered absolute derivative']); end
    end
end

% ----
if(learning),  nQRSaf=length(qrsAfs);
    fprintf('Number of annotated fetal QRSs= %d\n', nQRSaf);
    RRafs= diff(qrsAfs);          % annotated fetal RR
end
RRmc= diff(qrsM); RRms= RRmc/fs;  % detected mother RR
%
% -------------------------------------------------------------------------------
qrsFs=cell(1,4);
RRfs=cell(1,4);
madRR = zeros(ns,1);
rmssd = zeros(ns,1);

for ic=1:ns,
    
    if(dbFlag && learning), PlotSgnMrkNc(xs, qrsAfs*fs, fs, [cName,' - eS']); end
    % ----
    xic=X(:,ic);
    ecgic=ecg(:,ic);
    decgic=decgr(:,ic);
    decgBpic=decgBp(:,ic);
    adecgic=adecg(:,ic);
    adecgBpic=adecgBp(:,ic);
    
    % ----  QRS detection   (first step) ---------------------------------------------
    %        qrsM=QRSdetectorQeA2T(vadx,vdx,fs,pth,RRts,pmQT)
    pth=0.35; qrsFp=QRSdetectorF1(adecgBpic,decgic,fs,pth,0.45,0.8);
    %
    qrsFp1=qrsFp(1,:)';
    if(dbFlag) fprintf('Number of pre-detected fetal QRSs= %d\n', length(qrsFp1)); end
    RRfc1= diff(qrsFp1); RRfs1= RRfc1/fs;  % predetected fetal RR
    
    if(graph)
        figure, set(gcf,'Color','white'); hold on;
        plot(qrsMs, [RRms(1);RRms],'m.-');
        if(learning),   plot(qrsAfs, [RRafs(1);RRafs],'b.-');  end
        plot(qrsFp1/fs, [RRfs1(1);RRfs1],'r.-'); ylim([0.1, 1]);
        if(learning), title([cName,': detected mother(m), annotated (b) & pre-detected(r) fetal RR series']);
        else title([cName,': detected mother (m) & pre-detected fetal(r) RR series']); end
        figResize(0, 1, 1, .35);
    end
    
    % [nb,vp]=hist(RRfs,[0.25:0.05:0.9]);
    RRfsc1=sort(RRfs1); RRfsc1=RRfsc1(round(0.05*length(RRfsc1)):round(0.95*length(RRfsc1)));
    [nb,vp]=hist(RRfsc1);
    if(graph)
        figure, bar(vp,nb,1); title([cName,': pre-detected RR hist']);
    end
    [nbm,ibm]=max(nb);  RRst=vp(ibm);   % most frequent value
    
    adRR=abs(diff(RRfs1));
    adRRx=[repmat(adRR(1),2,1);adRR;repmat(adRR(end),2,1)];
    %adRRx=medfilt1(adRRx,5);
    mvadRR=filter(ones(1,5),1,adRRx); mvadRR= mvadRR(2*2+1:end,:);
    if(graph), figure, plot(qrsFp(1,3:end)'/fs,mvadRR); end
    thad=0.4;
    RRoklen=0;
    while(~RRoklen)
        for j=1:4
            isRRok=[]; ieRRok=[];
            jint=1;  RRokFlag=0;
            for i=1:length(mvadRR)
                if(RRst-0.1<RRfs1(i) && RRfs1(i)<RRst+0.1 && mvadRR(i)<thad)
                    if(~RRokFlag), RRokFlag=1; isRRok(jint)=i; else ieRRok(jint)=i; end
                else
                    if(RRokFlag), ieRRok(jint)=i; jint=jint+1; RRokFlag=0; end
                end
            end
            if(RRokFlag), ieRRok(jint)=i+1; RRokFlag=0; end
            if(isempty(isRRok)), RRoklen=0;
            else [RRoklen,im]=max(ieRRok-isRRok);
            end
            if(RRoklen>5), break; end
            thad=thad+0.1;
        end
    end
    RRsmi=mean(RRfs1(isRRok(im):ieRRok(im)));
    iqis=qrsFp(1,isRRok(im)); iqie=qrsFp(1,ieRRok(im));
    sgnOklens=(iqie-iqis)/fs;
    pth=0.2; qrsFb=QRSdetectorF2(adecgBpic(iqie:-1:1),decgic(iqie:-1:1),fs,pth,RRst,RRsmi,sgnOklens,1);  % backward
    qrsFb=iqie-qrsFb(end:-1:1);
    pth=0.2; qrsFf=QRSdetectorF2(adecgBpic(iqis:end),decgic(iqis:end),fs,pth,RRst,RRsmi,sgnOklens,0);    % forward
    qrsFf=iqis+qrsFf;
    if(isempty(qrsFb))
        qrsFic=qrsFf';
    elseif(isempty(qrsFf))
        qrsFic=qrsFb';
    else
        ibe=min(find(qrsFb>qrsFf(1)-0.5*RRst*fs));
        ifs=min(find(qrsFf>qrsFb(ibe)+0.5*RRst*fs));
        qrsFic=[qrsFb(1:ibe), qrsFf(ifs:end)]';
    end
    nQRSdf=length(qrsFic);
    if(dbFlag), fprintf('Channel=%d, number of detected fetal QRSs= %d\n', ic, nQRSdf); end
    qrsFsic=qrsFic/fs;
    RRfsic= diff(qrsFsic);   % detected fetal RR
    [nb,vp]=hist(RRfsic,[0.25:0.05:0.9]);
    if(graph)
        figure, bar(vp,nb,1);  title([cName,': detected RR hist']);
    end
    [nbm,ibm]=max(nb);  RRst=vp(ibm);
    
    RRfmean=meansc(RRfsic,4,4);
    %RfRmean= mean(RRfsic);
    RRfstd= std(RRfsic);
    if(dbFlag), fprintf('fetal RR mean= %d,  stdev=%d\n', RRfmean, RRfstd); end
    debugQRS=1;
    if(dbFlag && debugQRS)
        figure
        subplot(2,1,1),plot(adecgic,'b'); title([cName,': filter absolute derivative']);
        DrawVertMarker(qrsFic,'r',':','none');
        %     DrawVertMarker(qrsF(2,:)','b',':','none');
        %     DrawVertMarker(qrsF(end,:)','r',':','none');
        subplot(2,1,2),plot(ecgic(:,1),'b'); title('ecg signal');
        DrawVertMarker(qrsFic,'r',':','none');
        %     DrawVertMarker(qrsF(2,:)','b',':','none');
        %     DrawVertMarker(qrsF(end,:)','r',':','none');
    end
    
    if(graph || saveFig)
        if(graph), figure; else figure('visible','off'); end
        set(gcf,'Color','white'); hold on;
        plot(qrsMs, [RRms(1);RRms],'m.-');
        if(learning),  plot(qrsAfs, [RRafs(1);RRafs],'b.-');  end
        plot(qrsFsic, [RRfsic(1);RRfsic],'r.-');
        ylim([0.1, 1]);
        if(learning), title([cName,': ic=',num2str(ic),', detected mother(m), annotated (b) & detected(r) fetal RR series']);
        else title([cName,': ic=',num2str(ic),', detected mother (m) & fetal(r) RR series']); end
        figResize(0, 1, 1, .35);
        if(saveFigRRf), figFmtRR='png';
            figPathRR='../Figure/';
            if(~exist(figPathRR,'dir')), mkdir(figPathRR); end
            figName=fullfile(figPathRR,[cName,'_',num2str(ic),'_RRfetal']);
            print(gcf, ['-d',figFmtRR],figName);
        end
    end
    
    qrsFs{ic}=qrsFsic;
    RRfs{ic}=RRfsic;
    
end

% ----------------- choosing the best fetal RR series -----------------------
for ic=1:ns,
    meRR(ic) = median(RRfs{ic});
    if(meRR(ic)<0.35), cfact(ic)=0.3*(0.35-meRR(ic));
    elseif(meRR(ic)>0.5), cfact(ic)=0.05*(meRR(ic)-0.5);
    else cfact(ic)=0;
    end
    dRRfs=diff(RRfs{ic}); ddRRfs=diff(dRRfs);
    madRR(ic) = meansc(abs(dRRfs),0,5);
    maddRR(ic) = meansc(abs(ddRRfs),0,5);
    %rmssd(ic) = std(diff(RRfs{ic}));
    %dRRfsz=dRRfs-mean(dRRfs); rmssd(ic) = sqrt(meansc(dRRfsz.*dRRfsz,0,8));
    cmpRes=QRSdet_ann_cmp(qrsFs{ic},qrsMs, 0.1, 0);
    %FMsim(ic)=cmpRes.nTP/(cmpRes.nFP+cmpRes.nFN);
    cFMsim(ic)=0.05*cmpRes.nTP/(0.05*cmpRes.nTP+cmpRes.nFP+cmpRes.nFN);
    fprintf('channel=%1d,  meRR=%7.4f,  madRR=%7.4f,  maddRR=%7.4f, cfact=%7.4f, cFMsim=%7.4f\n', ...
        ic, meRR(ic), madRR(ic),maddRR(ic),cfact(ic),cFMsim(ic));
    chanIqf(ic)=madRR(ic)+maddRR(ic)+cfact(ic)+cFMsim(ic);
end
for ic=1:ns,
    fprintf('channel=%1d,  chanIqf=%7.4f\n', ic, chanIqf(ic));
end
[dummy,icsV]=sort(chanIqf);  % select the best channel
ics=icsV(1);
fprintf('Selected channel=%1d\n', ics);

qrsFes=qrsFs{ics};
RRfes=RRfs{ics};

xs=X(:,ics);
ecgs=ecg(:,ics);
% decgs=decgr(:,ics);
% adecgs=adecg(:,ics);

if(graph)
    figure, set(gcf,'Color','white');
    for is=1:ns,
        subplot(ns,1,is), plot(vtime, X(:,is));
        wgmi1= min(X(:,is)) -2;
        wgma1= max(X(:,is)) +2;
        ylim([wgmi1, wgma1]);
        set(gca,'YTick',[-5 0 5])
        % if(is~=ns), set(gca,'XTickLabel',''); end
        if(is==1), title([cName,': sorted fetal ECG']); end
        shg
    end
end
% -----------------------------------------------------------------------------
% if(dbFlag) PlotSgnMrkNc(ecgs, qrsF(1,:), fs, cName);  end

if(learning && graph), PlotSgnMrkNc(ecgs(:,1), {qrsAfs'*fs,qrsFe*fs}, fs, [cName,': fecg, ann.(m) & det.(b) qrs']); end
if(saveFig)
    figName=fullfile(figPath,[cName,'_ECGfetal']);
    print(gcf, ['-d',figFmt],figName);
end

if(graph || saveFig)
    figure, set(gcf,'Color','white');
    for is=1:ns,
        subplot(ns,1,is), plot(vtime, X(:,is));
        wgmi1= min(X(:,is)) -2;
        wgma1= max(X(:,is)) +2;
        ylim([wgmi1, wgma1]);
        %   DrawVertLines(vtime', qrsF(1,:)/fs, 'r', ylim);
        if(learning), DrawVertMarker(qrsAfs,'k',':','+'); end
        DrawVertMarker(qrsFes,'r',':','none');
        set(gca,'YTick',[-5 0 5])
        % if(is~=ns), set(gca,'XTickLabel',''); end
        if(is==1)
            if(learning), title([cName,': annotated (k) detected (r) fetal qrs']);
            else title([cName,': detected (r) fetal qrs']); end
        end
    end
    shg
    if(saveFig)
        figName=fullfile(figPath,[cName,'_fetal']);
        print(gcf, ['-d',figFmt],figName);
    end
end
%-------------------------------------------------------------

if(graph || saveFigRRf)
    if(graph), figure; else figure('visible','off'); end
    set(gcf,'Color','white'); hold on;
    plot(qrsMs, [RRms(1);RRms],'m.-');
    if(learning),  plot(qrsAfs, [RRafs(1);RRafs],'b.-');  end
    plot(qrsFes, [RRfes(1);RRfes],'r.-');
    ylim([0.1, 1]);
    if(learning), title([cName,': detected mother(m), annotated (b) & detected(r) fetal RR series']);
    else title([cName,': detected mother (m) & fetal(r) RR series']); end
    figResize(0, 1, 1, .35);
    if(saveFigRRf), figFmtRR='png';
        figPathRR='../FigureRRf/';
        if(~exist(figPathRR,'dir')), mkdir(figPathRR); end
        figName=fullfile(figPathRR,[cName,'_RRfetal']);
        print(gcf, ['-d',figFmtRR],figName);
    end
end


return
end %== function ================================================================
%
% -------------------------------------------------------------------------------------------------
%   Routine to compare QRS detection annotation
%
% Author: Maurizio Varanini, Clinical Physiology Institute, CNR, Pisa, Italy
% For any comment or bug report, please send e-mail to: maurizio.varanini@ifc.cnr.it
% -------------------------------------------------------------------------------------------------

function cmpRes=QRSdet_ann_cmp(qrsD,qrsA, diffMax, debug)
if(nargin<4) debug=0; end
nqrsA=length(qrsA);
nqrsD=length(qrsD);
qrsA_TP=zeros(1,nqrsA);
qrsA_diff=zeros(1,nqrsA);
nTP = 0;    % number of matching (TP)
nFN = 0;    % number of False Positive
for i = 1:nqrsA
    TP=0;
    qrsAi=qrsA(i);
    iQs = min(find(qrsD > qrsAi));
    iQp = max(find(qrsD <= qrsAi));
    if(~isempty(iQs))
        qrsDis = qrsD(iQs);
        diffQPs= qrsDis - qrsAi;
        if(diffQPs <= diffMax)
            TP=1;
            qrsA_TP(i) = qrsDis;
            qrsA_diff(i) = diffQPs;
        end
    end
    if(~isempty(iQp))
        qrsDip = qrsD(iQp);
        diffQPp= qrsAi - qrsDip;
        if(diffQPp <= diffMax && diffQPp<=diffQPs)
            TP=1;
            qrsA_TP(i) = qrsDip;
            qrsA_diff(i) = -diffQPp;
        end
    end
    if (TP)
        nTP = nTP + 1;
        qrsA_T(i) = qrsAi;
    else
        nFN = nFN + 1;
    end
end

nFP = max(0, nqrsD - nTP);

indTP=find(qrsA_TP>0);
if(~isempty(indTP))
    qrsA_diffTP=qrsA_diff(indTP);
    meanDiff=mean(qrsA_diffTP);
else
    meanDiff=0;   % patch
end
if(debug)
    fprintf('Number of mother QRSs= %d\n', nqrsA);
    fprintf('Number of fetal  QRSs= %d\n', nqrsD);
    
    fprintf('Number of true positives= %d\n', nTP);
    fprintf('Number of false positives= %d\n', nFP);
    fprintf('Number of false negatives= %d\n', nFN);
end

cmpRes.nTP=nTP;
cmpRes.nFP=nFP;
cmpRes.nFN=nFN;

cmpRes.meanDiff=meanDiff;
return
end     %== function ================================================================
%