% -------------------------------------------------------------------------------------------------
% 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