function qrsM=FecgQRSmDet(Se,fs,cName,graph,dbFlag,saveFig,qrsAf)
% ---------------------------------------------------------------------------------------------
%   Fecg: "Mother" QRS detection
%  - Best mother ECG selection based on a priori information on mother ECG pseudo-periodicity.
%    A raw derivative filter signal obtained as the difference between 
%    the average values on two intervals of 7ms far off 9ms)
%  - The absolute value of the raw derivative signal was filtered by a forward backward 
%    Butterworth bandpass filter (6.3-16.Hz).
%  - mother QRS detection based on this absolute derivative
%
% 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), graph=1; end
if(nargin<5), dbFlag=0; end
if(nargin<6), saveFig=0; end
if(nargin<7), qrsAf=[]; end

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

if(isempty(qrsAf));
    learning=0;   % test set, fetal annotations are not available
else
    learning=1;   % learning set, fetal annotations are available
end

%-------------------------------------------------------------
% recording duration
[ndt, ns]=size(Se);
vtime= [1:ndt]/fs;

if(dbFlag && learning), PlotSgnMrkNc(X, qrsAf*fs, fs, cName); end

for is=1:ns,
    X(:,is)= (Se(:,is)-mean(Se(:,is)))/std(Se(:,is));
end
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,': ECG signals']); end
    end
end
% ----
ecg=X;
ecgf=ecg;

% raw derivative filter coefficients (7ms before and after, 9ms in the
% middle)
nu=ceil(0.0070 * fs); nz=floor(0.0090*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(ecgf(1,:),delay,1);ecgf;repmat(ecgf(end,:),delay,1)];
decgr=filter(B,1,ecgfx);   decgr= decgr(2*delay+1:end,:);

adecg=abs(decgr); % absolute ECG derivative

% -----------    chose the best mother ecg channel
w8=fix(8*fs);              % window of 8 s
w2=fix(2*fs);              % window of 2 s
w02=fix(0.2*fs);           % window of 0.2 s
mD8=zeros(ns,1); mD2=zeros(ns,1); mD02=zeros(ns,1);
for is=1:ns, 
    % compute the average of maximum derivatives on windows of 8s (5% of maxima are discarted)
    mD8(is)=meanMaxSc(adecg(:,is), w8, 0,5);
    % compute the average of maximum derivatives on windows of 2s (5% of maxima are discarted)
    mD2(is)=meanMaxSc(adecg(:,is), w2, 0,5);  
    % compute the average of maximum derivatives on windows of 0.2s (1% of maxima are discarted)
    mD02(is)=meanMaxSc(adecg(:,is), w02, 0,1);  
end
qualFact=mD2./(mD02+mD8);
[dummy,ics]=sort(-qualFact);
% ------------
Ses=Se(:,ics);
ecgs=ecg(:,ics);
decgs=decgr(:,ics);
adecgs=adecg(:,ics);
% ----
if(graph || saveFig)
    figure, set(gcf,'Color','white');
    for is=1:ns,
        subplot(ns,1,is), plot(vtime, Ses(:,is));
        wgmi1= min(Ses(:,is)) -2;
        wgma1= max(Ses(:,is)) +2;
        ylim([wgmi1, wgma1]);
        set(gca,'YTick',[-5 0 5])
        % if(is~=ns), set(gca,'XTickLabel',''); end
        if(is==1), title([cName,': sorted mother ECG']); end
        shg
    end
    if(saveFig), figFmt='png';
        figPath=fullfile('../Figure/',progname);
        if(~exist(figPath,'dir')), mkdir(figPath); end
        figName=fullfile(figPath,[cName,'_mother']);
        print(gcf, ['-d',figFmt],figName);
    end
end

if(dbFlag && learning), PlotSgnMrkNc(Ses', qrsAf*fs, fs, [cName,' - eS']); end
% ----
ecg1=ecgs(:,1);
decg1=decgs(:,1);
adecg1=adecgs(:,1);

% The absolute value of the raw derivative signal was filtered by 
% a forward-backward  Butterworth bandpass filter (6.3-16.Hz).
fmind=5; fmaxd=20;
Wn = [fmind, fmaxd]/(fs/2);  % normalized cut-off frequency
[b,a]= butter(1,Wn);
adecg1=filtfilt(b,a,adecg1);


npx=fix(1*fs);
decg1x=[zeros(npx,1);decg1;zeros(npx,1)];
adecg1x=[zeros(npx,1);adecg1;zeros(npx,1)];

% --- QRS detection -----
% qrsM(1,:) QRSref    (reference point, max signed derivative )
% qrsM(2,:) QRSonset  (derivative overcomes half threshold)
% qrsM(3,:) supThDer  (derivative overcomes threshold)
% qrsM(4,:) maxAbsDer (max absolute derivative)
% qrsM(5,:) infThDer  (derivative becomes lower than threshold)
% qrsM(5,:) QRSoffset (derivative decreses below half threshold)
pth=0.45; qrsM=QRSdetectorM(adecg1x,decg1x,fs,pth,0.85,1);

qrsM=qrsM-npx;
qrsRef=qrsM(1,:);
nQRS=length(qrsRef);
fprintf('Number of detected mother QRSs= %d\n', nQRS);

RRc= diff(qrsRef);  RRs= RRc/fs;

if(learning), RRafs= diff(qrsAf); end    % annotated fetal qrs
[nb,vp]=hist(RRs,[0.4:0.05:1.9]);
if(graph)
    figure, bar(vp,nb,1); title([cName,': detected RR hist']);
end
% [nbm,ibm]=max(nb);  RRst=vp(ibm);   % most frequent value
% iqrsiT=find(RRst-0.05<RRs & RRs<RRst+0.05);

RRmean=meansc(RRs,4,4);
RRstd= std(RRs);
fprintf('RR mean= %5.3f,  stdev=%5.3f\n', RRmean, RRstd);
if(graph && dbFlag)
    figure
    subplot(2,1,1),plot(adecg1,'b');
    %    DrawVertMarker(qrsM(2,:)','c',':','none');
    %    DrawVertMarker(qrsM(end,:)','c',':','none');
    DrawVertMarker(qrsM(1,:)','r',':','none');           % reference point
    subplot(2,1,2),plot(ecg1,'b');
    %    DrawVertMarker(qrsM(2,:)','c',':','none');
    %    DrawVertMarker(qrsM(end,:)','c',':','none');
    DrawVertMarker(qrsM(1,:)','r',':','none');           % reference point
end

if(graph || saveFig)
    figure, set(gcf,'Color','white'); hold on;
    if(learning), plot(qrsAf, [RRafs(1);RRafs],'r.-'); end
    plot(qrsM(1,:)/fs, [RRs(1),RRs],'b.-');
    ylim([0.3, 2]);
    if(learning), title([cName,': mother (b) & annotated fetal(r) RR series']);
    else title([cName,': mother (b) RR series']); end
    figResize(0, 1, 1, .35);
    if(saveFig), figFmt='png';
        figPath=fullfile('../Figure/',progname);
        if(~exist(figPath,'dir')), mkdir(figPath); end
        figName=fullfile(figPath,[cName,'_RRmother']);
        print(gcf, ['-d',figFmt],figName);
    end
end
if(graph)
    figure, set(gcf,'Color','white');
    for is=1:ns,
        subplot(ns,1,is), plot(vtime, Ses(:,is));
        wgmi1= min(Ses(:,is)) -2;
        wgma1= max(Ses(:,is)) +2;
        ylim([wgmi1, wgma1]);
        %   DrawVertLines(vtime', qrsM(1,:)/fs, 'r', ylim);
        DrawVertMarker(qrsM(1,:)'/fs,'r',':','none');
        %     DrawVertMarker(qrsM(2,:)'/fs,'b',':','none');
        %     DrawVertMarker(qrsM(end,:)'/fs,'r',':','none');
        set(gca,'YTick',[-5 0 5])
        % if(is~=ns), set(gca,'XTickLabel',''); end
        if(is==1), title([cName,': ECGs with mother QRS marker ']); end
        shg
    end
end
if(graph || saveFig)
    PlotSgnMrkNc(ecgs(:,1), qrsM(1,:), fs, cName);
    shg
    if(saveFig), figFmt='png';
        figPath=fullfile('../Figure/',progname);
        if(~exist(figPath,'dir')), mkdir(figPath); end
        figName=fullfile(figPath,[cName,'_mother']);
        print(gcf, ['-d',figFmt],figName);
    end
end

return
end %== function ================================================================
%