%--------------------------------------------------------------------------
% Filter DxN + QRS Detector + V-beats recognition
%--------------------------------------------------------------------------
fp=fopen('G:\A1001D1.DAT','r'); % Open file and create file pointer , Fs=250Hz, V=5uV/bit
RecLen=10000;
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
% Filter 10x17
D=10;
N=17;
ECG2=zeros(size(d));
for i=1+D*round(N/2):RecLen-D*round(N/2)
m=0;
for j=-round(N/2)+1:round(N/2)-1
m=m+d(i+j*D);
end
m=m/N;
ECG2(i)=d(i)-m;
end
d=ECG2;
% QRS detector - 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));
PPP=zeros(size(d));
% Ventricular beats detector matrix initialise
M=zeros([64 64]);
row=zeros(1,'int8');
Last_QRS_Pos=1;
New_QRS_Pos=1;
Last5_QRS=ones(5);
QRS5_Cntr=1;
QRS_Delta=0;
RatePM_Arr=zeros(size(d));
QRS_Cntr=1;
DifRR=zeros(1);
DifRR_Arr=zeros(size(d));
Mean_QRS=zeros(1);
V_Beats=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)=ATN_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)));
TempS_Arr(i)=TempS_Arr(i)^2;
d1=(d(i)-d(i-1))*(d(i)-d(i+1));
d2=(d(i)-d(i-2))*(d(i)-d(i+2));
d3=(d(i)-d(i-3))*(d(i)-d(i+3));
d4=(d(i)-d(i-4))*(d(i)-d(i+4));
d5=(d(i)-d(i-5))*(d(i)-d(i+5));
% Peak?
if (d1>=0) && (d2>0) && (d3>0) && (d4>0) && (d5>0) %
PPP(i)=4000;
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) % QRS detected. Set the new threshold 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;
% Check QRS position
New_QRS_Pos=i-Last_QRS_Pos;
Mean_QRS=sum(Last5_QRS)/5;
DifRR=floor(100*(New_QRS_Pos-Mean_QRS(1))/New_QRS_Pos);
DifRR_Arr(i)=DifRR;
Last5_QRS(QRS5_Cntr)=New_QRS_Pos;
QRS5_Cntr=QRS5_Cntr+1;
if (QRS5_Cntr>=6)
QRS5_Cntr=1;
end
Last_QRS_Pos=i;
% Fill and check Matrix
k=1;
for j=i-16:i+47
y=cast(d(j)/5,'int16'); %Convert amplitude to discrets +/- 10 bits accuracy
y=(y+1024)/32; % Remove sign and convert to 64 bit accuracy
if (M(y,k)>3)
RatePM_Arr(i)=RatePM_Arr(i)+1;
end
if M(y,k)<10
M(y,k)=M(y,k)+1;
M(y+1,k)=M(y+1,k)+1;
M(y-1,k)=M(y-1,k)+1;
end
k=k+1;
end
% Partially erase of the Matrix (after every 10-th QRS)
QRS_Cntr=QRS_Cntr+1;
if (QRS_Cntr>=10)
QRS_Cntr=0;
for i=1:64
for j=1:64
if M(i,j)>1
M(i,j)=M(i,j)-1;
end
end
end
end
% Desigion rules for V-beat
if (RatePM_Arr(i)<30)
V_Beats(i)=1;
end
if (RatePM_Arr(i)>=30) && (RatePM_Arr(i)<40) && (DifRR<0)
V_Beats(i)=1;
end
if (RatePM_Arr(i)>=40) && (RatePM_Arr(i)<50) && (DifRR<=-10)
V_Beats(i)=1;
end
if (RatePM_Arr(i)>=50) && (RatePM_Arr(i)<60) && (DifRR<=-15)
V_Beats(i)=1;
end
end
end
end
for i=1:64
for j=1:64
M(i,j)=M(i,j)/10;
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, DifRR_Arr), title('DifRR'), xlabel('[Sec]'), ylabel('[%]'); % Ampitude and amplitude treshold
subplot(4,1,3), plot(Time, RatePM_Arr), title('RatePM'), xlabel('[Sec]'), ylabel('[Num]'); % Ampitude and amplitude treshold
subplot(4,1,4), plot(Time, V_Beats), title('V Beats'), xlabel('[Sec]'); % Current slope and slope treshold
figure(2)
pcolor (M);
colormap(flipud(gray));

