function Xr=FecgQRSmCanc(X,qrsM,fs,cName,graph,dbFlag,saveFig,qrsAfs)
% ---------------------------------------------------------------------------------------------
% Fecg: "Mother" ECG canceling using  PQRST approximation obtained by 
% Singular Value Decomposition (SVD).
% A trapezoidal window is used to select and weight the signal around each detected mother QRS.
%
% 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<4), cName=''; end
if(nargin<5), graph=1; end
if(nargin<6), dbFlag=1; end
if(nargin<7), saveFig=0; end
if(nargin<8), qrsAfs=[]; end

grafEigD=0;

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

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

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

if(dbFlag && learning), PlotSgnMrkNc(X, qrsAfs*fs, fs, [cName, ' - fetal QRS ann.']); end

vtime= [1:ndt]/fs;
% -----
% 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)

qrsR=qrsM(1,:);        % the reference point is the point of max signed derivative
if(dbFlag), PlotSgnMrkNc(X, qrsR, fs, [cName, ' - mother QRS det.']); end
nQRS=length(qrsR);
fprintf('Number of QRSs= %d\n', nQRS);

RRc= diff(qrsR);
RRs= RRc/fs;
RRmean=meansc(RRs,4,4);
%RRmean= mean(RRs);
RRstd= std(RRs);
fprintf('RR mean= %d,  stdev=%d\n', RRmean, RRstd);

% the reference point is the point of max signed derivative
npp=fix(0.2*fs);               % number of samples before the qrs reference
npd=fix(min(0.5, 0.8*(RRmean-0.1))*fs);  % number of samples after the qrs reference
npt=1+npp+npd;

% extend the signals in order to manage the first and the last QRS
Xx=X;
npqp=fix(0.12*fs);
if(qrsR(1)-npqp < 1), qi=2; 
    npxp=0;
    Xx(1:npqp,:)=repmat(Xx(npqp+1,:),npqp,1);
else qi=1;
    npxp=max(0,npp+1-qrsR(1));  % number of samples added to left of each signal
    Xx=[repmat(X(1,:),npxp,1);X];
end
npqd=fix(0.14*fs);
qf=nQRS;
if(qrsR(end)+npqd > size(X,1)), qf=qf-1;
    Xx(end-npqd+1:end,:)=repmat(Xx(end-npqd-1,:),npqd,1);
end
if(qrsR(qf)+0.85*fs*median(RRs(end-4:end)) < size(X,1))
    npep=fix(max(0.1*fs, 0.15*fs*mean(RRs(end-4:end))));
    Xx(end-npep+1:end,:)=repmat(Xx(end-npep-1,:),npep,1);
end
npxd=max(0,qrsR(qf)+npd-ndt);  % number of samples added to right of each signal
Xx=[Xx;repmat(X(end,:),npxd,1)];

ndtx=size(Xx,1);
nqe= qf-qi+1;
vtimex= [1-npxp:ndt+npxd]/fs;


% built vectors of indexes (start and end) of QRS window
iqw=qrsR(qi:qf);
iw=npxp+iqw-npp;      %  start of QRS window
fw=npxp+iqw+npd;      %  end of QRS window
A=zeros(npt,nqe);
wwg=weightFun2(npp,npd,fs);
Xc=zeros(ndtx,ns);
for is=1:ns
    for iq=1:nqe
        iwq=iw(iq); fwq=fw(iq);
        A(:,iq)=Xx(iwq:fwq,is).*wwg;
    end
    % Singular value decomposition
    % S = diagonal matrix containing the singular values
    % U and V = unitary matrices so that X = U*S*V'.
%    tic;
    [U,S,V] = svd(A,0);
%    fprintf('SVD (%d,%d) computation time = %3.2f\n', size(A), toc);
    
%    fprintf('%f,  ', diag(S));   fprintf('\n');
    if(grafEigD), figure; plot(diag(S),'r-*');
        title([num2str(is),'- Singular values']); drawnow;
    end
    
    mt=nqe;
        
    nds=3;     % select the number of singular values  <===
    
    sv=diag(S);
    pownds(is)=sum(sv(1:nds).^2)/sum(sv(1:nqe).^2);
    if(dbFlag), fprintf('is=%d, pownds=%f\n', is, pownds(is)); end
    
%    tic
    Sr=S; for k=1:mt-nds-1, Sr(mt-k,mt-k)=0; end
    Ar=U*Sr*V';
%    fprintf('Ar (%d,%d) rebuilding time = %3.2f\n', size(Ar), toc);
    
    for iq=1:nqe
        iwq=iw(iq); fwq=fw(iq);
%        Xc(iwq:fwq,is)=Ar(:,iq);
        Xc(iwq:fwq,is)=Ar(:,iq)./wwg;   % .*windLR;
    end
    
    % smoothing connections between successive windows
    for iq=1:nqe-1
        fwq=fw(iq); iwqs=iw(iq+1);
        if(iwqs>fwq)
            dv=Xc(iwqs,is)-Xc(fwq,is);
            pv=dv/(iwqs-fwq);
            Xc(fwq+1:iwqs-1,is)=Xc(fwq,is)+pv*(1:iwqs-fwq-1)';
        end
    end
end
if(dbFlag)
    for is=1:ns,
        figure
        set(gcf,'Color','white');
        subplot(2,1,1), plot(vtimex,[Xx(:,is),Xc(:,is)]);
        wgmi1= min(Xx(:,is)) -2;
        wgma1= max(Xx(:,is)) +2;
        ylim([wgmi1, wgma1]);
        set(gca,'YTick',[-5 0 5])
        DrawVertMarker(iw'/fs,'r',':','none');
        DrawVertMarker(qrsM(1,:)'/fs,'g',':','none');
        DrawVertMarker(fw'/fs,'m',':','none');
        title([cName,' - ',num2str(is),': ecg & ecgM']);
        subplot(2,1,2), plot(vtimex,[Xx(:,is)-Xc(:,is)]);
        ylim([wgmi1, wgma1]);
        set(gca,'YTick',[-5 0 5])
        DrawVertMarker(iw'/fs,'r',':','none');
        DrawVertMarker(qrsM(1,:)'/fs,'g',':','none');
        DrawVertMarker(fw'/fs,'m',':','none');
        title([cName,' - ',num2str(is),': ecg - ecgM']);
        shg
    end
end
Xc= Xc(npxp+1:end-npxd,:);
Xx= Xx(npxp+1:end-npxd,:);
Xe=Xx-Xc;
iw=iw-npxp;  fw=fw-npxp;
if(graph)
    for is=1:ns,
        figure,  set(gcf,'Color','white');
        subplot(2,1,1), plot(vtime,[X(:,is),Xc(:,is)]);
        wgmi1= min(X(:,is)) -2;
        wgma1= max(X(:,is)) +2;
        ylim([wgmi1, wgma1]);
        set(gca,'YTick',[-5 0 5])
        DrawVertMarker(iw'/fs,'r',':','none');
        DrawVertMarker(qrsM(1,:)'/fs,'g',':','none');
        DrawVertMarker(fw'/fs,'m',':','none');
        title([cName,' - ',num2str(is),': ecg & ecgM']);
        subplot(2,1,2), plot(vtime,Xe(:,is));
        ylim([wgmi1, wgma1]);
        set(gca,'YTick',[-5 0 5])
        DrawVertMarker(iw'/fs,'r',':','none');
        DrawVertMarker(qrsM(1,:)'/fs,'g',':','none');
        DrawVertMarker(fw'/fs,'m',':','none');
        title([cName,' - ',num2str(is),': ecg - ecgM']);
        shg
    end
end
%--------------------------------------------------------------------------------------------
if(graph || saveFig)
    figure, set(gcf,'Color','white');
    for is=1:ns,
        subplot(ns,1,is), plot(vtime,Xe(:,is));
        wgmi1= min(Xe(:,is)) -2;
        wgma1= max(Xe(:,is)) +2;
        ylim([wgmi1, wgma1]);
        DrawVertMarker(iw'/fs,'r',':','none');
        DrawVertMarker(qrsM(1,:)'/fs,'g',':','none');
        DrawVertMarker(fw'/fs,'m',':','none');
        set(gca,'YTick',[-5 0 5])
        % if(is~=ns), set(gca,'XTickLabel',''); end
        if(is==1), title([cName,': mixed foetal ECG']); end
    end
    shg
    if(saveFig), figFmt='png';
        figPath=fullfile('../Figure/',progname);
        if(~exist(figPath,'dir')), mkdir(figPath); end
        figName=fullfile(figPath,[cName,'_mCanc']);
        print(gcf, ['-d',figFmt],figName);
        %    saveas(gcf, figName,'fig');
    end
end
%-------------------------------------------------------------

Xr=Xe;

return
end % = function ===============================================================
% ------------------------------------------------------------------
function wwg=weightFun2(npp,npd,fs)
nppc1=fix(0.06*fs); npdc1=fix(0.06*fs);
nppc2=fix(0.08*fs); npdc2=min(fix(0.2*fs),npd-npdc1);
ii1=1; ie1=npp-nppc1-nppc2;
ii2=ie1+1; ie2=ie1+nppc2;
ii3=ie2+1; ie3=ie2+nppc1+npdc1+1;
ii4=ie3+1; ie4=ie3+npdc2;
ii5=ie4+1; ie5=npp+npd+1;
wwg(ii1:ie1)=0.20;
wwg(ii2:ie2)=0.20+0.8*(1:nppc2)/nppc2;
wwg(ii3:ie3)=1;
wwg(ii4:ie4)=1-0.8*(1:npdc2)/npdc2;
wwg(ii5:ie5)=0.20;
wwg=wwg';
% figure, plot(wwg);
return
end % = function ===============================================================
%