% -------------------------------------------------------------------------------------------------
% QrsAvgCancelTestClPcan.m: QRST complex canceling
% 
% This is an experimental routine. 
% I think many things should be changed to improve QRST canceling.
%
% Input parameter:
% mainPath_in = main path of data files ('..\..\www.physionet.org\physiobank\database\aftdb')
% casiDir     = subdirectory with data files (learning-set, test-set-a, test-set-b)
% cname       = name of the case to process
% mainPath    = main path for output data files ('..\..\')
%
% input file:
%    mainPath_in \ casiDir \ cname .dat
%    (example: '..\..\www.physionet.org\physiobank\database\aftdb\learning-set\n01.dat'):
% output files (directory: mainPath \DatiQc):
%    signal file (cname .dat),  containing: 4 channels, the last of them obtained canceling QRST from ecg signals.
%    parameter file (cname .mat),  containing: 'RRs', 'RRmean', 'RRstd', 'nc', 'crrf', 'rp1fn', 'rp2fn'
%    (example: '..\..\datiQc\learning-set\n01.dat'):
%
% Copyright (C) 2004 Maurizio Varanini, Clinical Physiology Institute, CNR, Pisa, Italy
%
% 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.
%
% You should have received a copy of the GNU General Public License along with this program;
% if not, write to the Free Software Foundation, Inc.,
% 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
%
% For any comment or bug report, please send e-mail to: maurizio.varanini@ifc.cnr.it
% -------------------------------------------------------------------------------------------------

function QrsAvgCancelTestClpCan(mainPath_in, casiDir, cname, mainPath)
%nargin=0
echo off
close all
debug=0;
debugQRS=0;
graf=0;
savedat=1;

if nargin<1
    debug=1;
    debugQRS=0;
    graf=1;
    mainPath_in='..\..\www.physionet.org\physiobank\database\aftdb';
    mainPath='..\..\';

    savedat=1;
%   savedat=minput('Overwrite previous results? (yes=1)','%d',savedat);
%    casiDir='learning-set';
%    cname='n01';
    casiDir='test-set-a';
    cname='a10';
%    cname='a13';
%    cname='a19';
%    cname='a29';
%    casiDir='test-set-b';
%    cname='b20';

end 
    
dir_out='DatiQc';
if(debug) dir_out=[dir_out, 'D']; end
mainPath_out=fullfile(mainPath, dir_out);
if ~exist(mainPath_out,'dir') mkdir(mainPath, dir_out); end
filepath_out= fullfile(mainPath_out, casiDir);
if ~exist(filepath_out,'dir') mkdir(mainPath_out, casiDir); end
fpname_out= fullfile(filepath_out, [cname, '.dat']);

dataType='short';
ns=2;
extension='.dat';
freq=128;

[progpath, progname] = fileparts(which(mfilename));

fprintf('\n --------------------------------------------------------- \n');
fprintf('Case name: %s\n', cname);

fpname_in = fullfile(mainPath_in, casiDir, [cname extension]);
fprintf('Reading signal from file: %s\n', fpname_in);
fid  = fopen(fpname_in,'r');
[data, count] = fread(fid,[ns,inf],dataType);
fclose(fid);
%-------------------------------------------------------------

% recording duration 
[ns ndt]=size(data);
lendts=ndt/freq;
fprintf('Sampling frequency= %g\n', freq);
fprintf('Total recording length = %g (samples) = %g (seconds)\n', ndt, lendts);

if(graf)
    PlotSgnMrkN([data(1,:);data(2,:)], [], freq, cname);
end

x1=detrend(data(1,:)');
x2=detrend(data(2,:)');

fmind=0.2;
fmaxd=40;
Wn = [fmind fmaxd]/(freq/2);  % The cut-off frequency Wn must be between 0 < Wn < 1.0, with 1.0 corresponding to half the sample rate
delay=1000;
ncf=delay*2+1;
B = FIR1(ncf,Wn); % default -> Hamming
x1 = filter(B,1,[x1; zeros(delay,1)]);   x1= x1(delay+1:end);
x2 = filter(B,1,[x2; zeros(delay,1)]);   x2=x2(delay+1:end);

if(graf)
    PlotSgnMrkN([x1,x2], [], freq, [cname, '  filtered']);
end

% interpolation
fs=1024;
x1 = InterpFFT(x1, freq, fs);
x2 = InterpFFT(x2, freq, fs);
freq=fs;
xtime=0:length(x1)-1; xtime=xtime/freq;
B=[1;1;1;1;1;1;1-1; 0.5; 0;0;0; -0.5; -1;-1;-1;-1;-1; ];
delay=floor(length(B)/2);
% compute the derivative signals
dx1=filter(B,1,[x1; zeros(delay,1)]);   dx1=dx1(delay+1:end);
dx2=filter(B,1,[x2; zeros(delay,1)]);   dx2=dx2(delay+1:end);
adx1=abs(dx1);
adx2=abs(dx2);
vadx=adx1+adx2;
if(graf & debugQRS)
    % plot derivative signals and spatial velocity
    figure;
    subplot(3,1,1), plot(xtime,x1,'b',xtime,dx1,'g');
    subplot(3,1,2), plot(xtime,x2,'b',xtime,dx2,'g');
    subplot(3,1,3), plot(xtime,vadx,'g');
end

w02=fix(0.2*freq);
mD02=meanMaxSc(vadx, w02, 1,1);

w2=fix(2*freq);
mD2=meanMaxSc(vadx, w2, 1,1);

pth=0.35;
%mvadx=madx(:,2);
mvadx=vadx;
% QRS detection
qrsM=QRSdetector(mvadx,freq,pth);

nQRS=size(qrsM,2);
fprintf('Number of QRSs= %d\n', nQRS);

RRc= diff(qrsM(1,:));
RRs= RRc/freq;
RRmean=meansc(RRs,4,4);
%RRmean= mean(RRs);
RRstd= std(RRs);
fprintf('RR mean= %d,  stdev=%d\n', RRmean, RRstd);

if(graf & debugQRS)
    figure
    plot(RRs);
    figure
    plot(vadx,'b');
    DrawVertMarker(qrsM(1,:)','r',':','none');
    DrawVertMarker(qrsM(2,:)','b',':','none');
    DrawVertMarker(qrsM(3,:)','r',':','none');
end    
if(graf) PlotSgnMrkN([x1,x2], qrsM(1,:), freq, cname);  end

x1=(x1-mean(x1))/std(x1);
QRSw=meansc(qrsM(3,:)-qrsM(1,:),2,2);
m=QRSw+2;
npp=fix(0.1*freq); 
npd=fix(min(0.5, RRmean-0.1)*freq);
npt=1+npp+npd;
wigp=10;
wigd=10;
% the reference point is the first time (derivate > threshold)
if(qrsM(1,1)-npp -wigp < 1) qi=2; else qi=1; end 
if(qrsM(1,nQRS)+npd +wigd > length(x1)) qf=nQRS-1; else qf=nQRS; end 
nq= qf-qi+1;

% built matrix of QRS templates
iw=qrsM(1,qi:qf)-npp;
fw=qrsM(1,qi:qf)+npd;
ifw=find(iw>0 & fw<size(x1,1));
iw=iw(ifw);
fw=fw(ifw);
nqe=size(iw,2);
% for i=1:nqe;
%     x1M(:,i)= x1(iw(i):fw(i));
%     x2M(:,i)= x2(iw(i):fw(i));
% end
% if(graf)
%     figure
%     plot(x1M);
%     title('ecg 1: all QRS template')
%     figure
%     plot(x2M);
%     title('ecg 2: all QRS template')
% end

beatGroup1= beatGroupingSC1(x1, iw, npt);
vq{1}=beatGroup1.vq;
vij{1}=beatGroup1.vij;
vih{1}=beatGroup1.vih;
vqcc{1}=beatGroup1.vqcc;
avgC1Q=beatGroup1.avgC1Q;
nc=size(avgC1Q,2);
fprintf(' number of classes for channel 1= %d\n', nc);

beatGroup2= beatGroupingSC1(x2, iw, npt);
vq{2}=beatGroup2.vq;
vij{2}=beatGroup2.vij;
vih{2}=beatGroup2.vih;
vqcc{2}=beatGroup2.vqcc;
avgC2Q=beatGroup2.avgC1Q;
nc=size(avgC2Q,2);
fprintf(' number of classes for channel 2= %d\n', nc);

%fprintf('min, max shift of template= %f,  %f\n', min(iw-vij), max(iw-vij))

% Plot the QRS template of each class
if(graf)
    figure, hold on
    subplot(2,1,1), plot(avgC1Q);
    title('ecg 1: QRS template classes')
    subplot(2,1,2), plot(avgC2Q);
    title('ecg 2: QRS template classes')
    hold off
end
% For each class, plot all QRS belonging to it
% if(graf)
%     figure, hold on
%     for ic=1:nc
%         subplot(2,nc,ic), plot(x1M(:,vq{1}{ic}));
%         title(['ecg 1, class=', num2str(ic), ': QRS '])
%         subplot(2,nc,nc+ic), plot(x2M(:,vq{2}{ic}));
%         title(['ecg 2, class=', num2str(ic), ': QRS '])
%     end
%     hold off
% end
% building signal by the class centroid
sq1=zeros(length(x1),1);
sq2=sq1;
vp1=0; vp2=0;
jp1=1; jp2=1;
for iq=1:nqe
    ic1=vqcc{1}(iq);  
    ij1=vij{1}(iq);  
    fj1=ij1+npt-1;   
    ih1=vih{1}(iq);
    ic2=vqcc{2}(iq);  ij2=vij{2}(iq);  fj2=ij2+npt-1;   ih2=vih{2}(iq);
    va1=avgC1Q(1,ic1)+ih1;  dv1=va1-vp1;  
    va2=avgC2Q(1,ic2)+ih2;  dv2=va2-vp2;  
    if(ij1>jp1)
        pv1=dv1/(ij1-jp1);
        sq1(jp1+1:ij1-1)=vp1+pv1*(1:ij1-jp1-1);
    end
    if(ij2>jp2)
        pv2=dv2/(ij2-jp2);
        sq2(jp2+1:ij2-1)=vp2+pv2*(1:ij2-jp2-1);
    end

    sq1(ij1:fj1)= avgC1Q(:,ic1)+ih1;
    sq2(ij2:fj2)= avgC2Q(:,ic2)+ih2;
    vp1= avgC1Q(end,ic1)+ih1;  
    vp2= avgC2Q(end,ic2)+ih2;
    jp1=fj1; 
    jp2=fj2;
end
ij1=qrsM(1,end)-npp;
ij2=qrsM(2,end)-npp;
va1=x1(ij1);  dv1=va1-vp1;  pv1=dv1/(ij1-jp1);
sq1(jp1+1:ij1-1)=vp1+pv1:pv1:va1-pv1;
sq1(ij1:end)=x1(ij1:end);
va2=x2(ij2);  dv2=va2-vp2;  pv2=dv2/(ij2-jp2);
sq2(jp2+1:ij2-1)=vp2+pv2:pv2:va2-pv2;
sq2(ij2:end)=x2(ij2:end);

e1=x1-sq1;
e2=x2-sq2;

if(graf)
    % channel 1
    PlotSgnMrkN([x1,sq1,e1], qrsM(1,:), freq, [cname '_c1: x,y,e']);
    % channel 2
    PlotSgnMrkN([x2,sq2,e2], qrsM(1,:), freq, [cname '_c2: x,y,e']);
end

% Begin code for channel selection
w02=fix(0.2*freq);  % averaged value of maxima on 0.2s wide windows
mD02=meanMaxSc(e1, w02, 1,1);
w2=fix(2*freq);   % averaged value of maxima on 2.0s wide windows
mD2=meanMaxSc(e1, w2, 1,1);
rq1=mD2/mD02;

w02=fix(0.2*freq);  % averaged value of maxima on 0.2s wide windows
mD02=meanMaxSc(e2, w02, 1,1);
w2=fix(2*freq);   % averaged value of maxima on 2.0s wide windows
mD2=meanMaxSc(e2, w2, 1,1);
rq2=mD2/mD02;
ic= 1+ (rq2<rq1);
fprintf('rq1=%f; rq2=%f; ic=%d\n', rq1, rq2, ic);

% [miy, may] = mimaxsc(e1,5,5);
[miy, may] = mimaxsc(e1,5,5);
e1n=e1/(may-miy);
[miy, may] = mimaxsc(e2,5,5);
e2n=e2/(may-miy);
if(graf)
    PlotSgnMrk(e1n, qrsM(1,:), freq, [cname, ':  e1n']);
    PlotSgnMrk(e2n, qrsM(1,:), freq, [cname, ':  e2n']);
end
% clipping
mima=[-2, 2];
e1nb=boundV(e1n,mima);
e2nb=boundV(e2n,mima);

pn1b=e1nb'*e1nb;
pn2b=e2nb'*e2nb;
crb=e1nb'*e2nb;
crrb=crb/sqrt(pn1b*pn2b);
fprintf('crrb=%f\n', crrb);
if(graf)
    PlotSgnMrk([e1n,e2n], qrsM(1,:), freq, [cname, ':  e1nb, e2nb']);
end

% clipping for computing cr
mima=[-0.8, 0.8];
e1nbcr=boundV(e1n,mima);
e2nbcr=boundV(e2n,mima);
fmind=3;
fmaxd=8;
Wn = [fmind fmaxd]/(freq/2);  % The cut-off frequency Wn must be between 0 < Wn < 1.0, with 1.0 
                              % corresponding to half the sample rate
B = FIR1(401,Wn); % default -> Hamming
e1f = filter(B,1,[e1nbcr; zeros(delay,1)]);    e1f=e1nbcr(delay+1:end);
e2f = filter(B,1,[e2nbcr; zeros(delay,1)]);    e2f=e2nbcr(delay+1:end);
if(graf)
    PlotSgnMrk([e1f,e2f], qrsM(1,:), freq, cname);
end

crf=e1f'*e2f; % correlation coefficient of filtered residual signal
pf1=e1f'*e1f;
pf2=e2f'*e2f;
crrf=crf/sqrt(pf1*pf2);
fprintf('crf=%f, pf1=%f, pf2=%f, crrf=%f\n', crf, pf1, pf2, crrf);
rp1fn=pf1/pn1b;
rp2fn=pf2/pn2b;
fprintf('rp1fn=%f, rp2fn=%f\n', rp1fn, rp2fn);
%
yd=e1nb-e2nb;
ys=e1nb+e2nb;
if(graf)
    PlotSgnMrk(yd, qrsM(1,:), freq, [cname, ':  e1nb-e2nb']);
    PlotSgnMrk(ys, qrsM(1,:), freq, [cname, ':  e1nb+e2nb']);
end
if(abs(crrf) < 0.1) 
    if(rq1 > rq2) 
        e1s=e1nb; e2s=e2nb; 
        fprintf(' --> channel 2\n');
    else 
        e1s=e2nb; e2s=e1nb; 
        fprintf(' --> channel 1\n');
    end
else
    if(crrf > 0) 
        e1s= yd; e2s=ys; 
        fprintf(' --> sum (1+2)\n');
    else 
        e1s= ys; e2s=yd;  
        fprintf(' --> difference (1-2)\n');
    end
end
if savedat
    sgn_out=[x1, x2, e1s, e2s];
    fprintf('Writing signals on file: %s\n', fpname_out);
    fid  = fopen(fpname_out,'wb');
    count = fwrite(fid, sgn_out', 'double');
    fclose(fid);
    
    fpnamePar_out= fullfile(filepath_out, [cname,'.mat']);
    fprintf('Writing parameters on file: %s\n', fpnamePar_out);
    save(fpnamePar_out, 'RRs', 'RRmean', 'RRstd', 'nc', 'crrf', 'rp1fn', 'rp2fn');
end