%-------------------------------------------------------------------------
% QRS Detection - Caceres method
%--------------------------------------------------------------------------
fp=fopen('G:\A1001D1.DAT','r'); % Open file and create file pointer , Fs=250Hz, V=5uV/bit
RecLen=3000;
d=fread(fp,RecLen,'int16'); % Read 3000 discrets from the file
d=d*0.005; % convert digits to millivolts
%Change sample frequency from 250Hz to 1000Hz
d500=zeros(size(d)*2);
for i=1:size(d)-1
d500((i-1)*2+1)=d(i);
d500((i-1)*2+2)=(d(i)+d(i+1))/2;
end
d1000=zeros(size(d500)*2);
for i=1:size(d500)-1
d1000((i-1)*2+1)=d500(i);
d1000((i-1)*2+2)=(d500(i)+d500(i+1))/2;
end
Time=0.001:0.001:RecLen*0.001*4; % Time vector
ECG=d1000;
ECG_Delta=d1000(11:length(d1000)) - d1000(1:length(d1000)-10);
Time_Delta=0.001:0.001:(RecLen)*0.001*4-0.01;
figure(1);
subplot(5,1,1); plot(Time,ECG); axis([1 4 -5 5]); title('ECG Record'); ylabel('[mV]')
subplot(5,1,2); plot(Time_Delta,ECG_Delta); axis([1 4 -5 5]); title('ECG Delta'); ylabel('[mV]')
subplot(5,1,3); plot(Time,ECG); title('ECG Full Record'); ylabel('[mV]')
% 370 ms periods and find Dmin and Dmax for any one
i=1;
Num_Intervals=round(length(ECG_Delta)/370+1);
Minimums=zeros(Num_Intervals);
Maximums=zeros(Num_Intervals);
while (i<Num_Intervals)
ECG370=ECG_Delta((i-1)*370+1:i*370+1);
Minimums(i)=min(ECG370);
Maximums(i)=max(ECG370);
i=i+1;
end
% Find the maximums and the minimums in all intervals
Dmin=0;
Dmax=0;
for i=1:Num_Intervals
if (Dmin>Minimums(i))
Dmin=Minimums(i);
end
if (Dmax<Maximums(i))
Dmax=Maximums(i);
end
end
% Determination of the polarity
Min_Max=0;
if (abs(Dmax)>abs(Dmin))
K=5/8*Dmax; %work with maximums
Min_Max=1;
else
K=5/8*Dmin; %work with minimums
Min_Max=-1;
end
%find minimums or maximums by 5/8 criteria
QRS=zeros(size(ECG_Delta));
for i=1:length(ECG_Delta)
if ((Min_Max>0) && (ECG_Delta(i)>K))
QRS(i)=ECG_Delta(i);
end
if ((Min_Max<0) && (ECG_Delta(i)<K))
QRS(i)=ECG_Delta(i);
end
end
subplot(5,1,4); plot(Time_Delta,QRS); title('QRS');
%clear minimums or maximums if they are too close (dist<184mS)
QRS_Clear=QRS;
for i=92:length(QRS)-93
for j=-91:92
if (abs(QRS(i))<abs(QRS(j+i)))
QRS_Clear(i)=0;
end
end
end
subplot(5,1,5); plot(Time_Delta,QRS_Clear); title('Final QRS');
%--------------------------------------------------------------------------
% QRS Detector
%--------------------------------------------------------------------------
fp=fopen('D:\AHA_250Hz\A1001D1.DAT','r'); % Open file and initialize file pointer , Fs=250Hz, V=5uV/bit
RecLen=3000;
d=fread(fp,RecLen,'int16'); % Read 3000 discrets from the file
d=d*5; % convert digits to uV
Time=0.004:0.004:RecLen*0.004; % Time vector
% Param's Initialise
MAP=100; % Mean Amplitude Positive = 100uV
MAN=-100; % Mean Amplitude Negative = -100uV
MSL=0; % Mean Slope = 0
MAP_Arr=[100 100 100 100];
MAN_Arr=[100 100 100 100];
MSL_Arr=[0 0 0 0];
MAP_Idx=1;
MAN_Idx=1;
MSL_Idx=1;
ATP=100; % Amplitude Treshold Positive = 100uV
ATN=-100; % Amplitude Treshold Negative = -100uV
ATP_Arr=zeros(size(d)); %ATP Array
ATN_Arr=zeros(size(d)); %ATN Array
SLT=5; %Slope Treshold = 5uV
SLT_Arr=zeros(size(d)); %Slopes Treshold Array
TempS=0; %Current Slope
TempS_Arr=zeros(size(d)); %Slopes Array
Mute=0; %60mS in samples
R_Peak=false;
Peaks=zeros(size(d));
for i=11:RecLen-11
if (Mute>0) % 240mS (Mute=60) pause after QRS detection
ATP_Arr(i)=ATP_Arr(i-1);
ATN_Arr(i)=ATP_Arr(i-1);
SLT_Arr(i)=SLT_Arr(i-1);
Mute=Mute-1;
continue;
end
if (ATP>100) % Modify ATP
ATP=ATP-5;
ATP_Arr(i)=ATP;
end
if (ATN<-100) % Modify ATN
ATN=ATN+5;
ATN_Arr(i)=ATN;
end
if (SLT>1) % Modify SLT
SLT=SLT-SLT/128;
SLT_Arr(i)=SLT;
end
% Calculate current slope
TempS_Arr(i)=abs((d(i)-d(i-1)) + (d(i)-d(i-2)) + (d(i)-d(i-3)) + (d(i)-d(i-4)) + (d(i)-d(i-5)) + (d(i)-d(i+5)) +(d(i)-d(i+4)) + (d(i)-d(i+3)) + (d(i)-d(i+2)) +(d(i)-d(i+1)));
% Is a Peak?
if (((d(i)-d(i-5))*(d(i)-d(i+5))>0) && ((d(i)-d(i-4))*(d(i)-d(i+4))>0) && ((d(i)-d(i-3))*(d(i)-d(i+3))>0) && ((d(i)-d(i-2))*(d(i)-d(i+2))>0) && ((d(i)-d(i-1))*(d(i)-d(i+1))>=0))
TempS=TempS_Arr(i);
MAP=sum(MAP_Arr)/4; % Calculate current moving averages
MAN=sum(MAN_Arr)/4;
MSL=sum(MSL_Arr)/4;
if ( (TempS>SLT) && ((d(i)>ATP) || (d(i)<ATN)) ) % I-st QRS detection rules
R_Peak=true;
end
if ( (TempS>0.75*SLT) && ((d(i)>1.2*ATP) || (d(i)<1.2*ATN)) ) % II-nd QRS detection rules
R_Peak=true;
end
if (R_Peak==true) % True QRS. Set new tresholds values
if (d(i)>0) % Update positive amplitude tresholds
if (ATP*1.5<MAP)
ATP=d(i);
else
ATP=MAP;
end
ATP_Arr(i)=ATP;
MAP_Arr(MAP_Idx)=d(i);
MAP_Idx=MAP_Idx+1;
if (MAP_Idx>=5)
MAP_Idx=1;
end
end
if (d(i)<0) % Update negative amplitude tresholds
if (ATN*1.5>MAN)
ATN=d(i);
else
ATN=MAN;
end
ATN_Arr(i)=ATN;
MAN_Arr(MAN_Idx)=d(i);
MAN_Idx=MAN_Idx+1;
if (MAN_Idx>=5)
MAN_Idx=1;
end
end
MSL_Arr(MSL_Idx)=TempS; %Update slope tresholds
SLT=TempS;
SLT_Arr(i)=SLT;
MSL_Idx=MSL_Idx+1;
if (MSL_Idx>=5)
MSL_Idx=1;
end
Peaks(i)=1;
R_Peak=false;
Mute=60;
end
end
end
figure(1)
subplot(4,1,1), plot(Time, d), title('ECG Record'), xlabel('[Sec]'), ylabel('[uV]'); % Original ECG
subplot(4,1,2), plot(Time, ATP_Arr, Time, d), title('Amplitude Treshold'), xlabel('[Sec]'), ylabel('[uV]'); % Ampitude and amplitude treshold
subplot(4,1,3), plot(Time, SLT_Arr, Time, TempS_Arr), title('Slope Treshold'), xlabel('[Sec]'), ylabel('[uV]'); % Current slope and slope treshold
subplot(4,1,4), plot(Time, Peaks), title('QRS Markers'), xlabel('[Sec]'), ylabel('QRS'); % QRS Marker

