function [ECGpeaks] = QSpeaks( ECG,Rposition,fs )
%Q,S peaks detection
% point to point time duration is determined by sampling frequency: 1/fs
% the duration of QRS complex varies from 0.1s to 0.2s
% complex--fs*0.2
aveHB = length(ECG)/length(Rposition);
fid_pks = zeros(length(Rposition),7);
% P QRSon Q R S QRSoff
%% set up searching windows
windowS = round(fs*0.1); windowQ = round(fs*0.05);
windowP = round(aveHB/3); windowT = round(aveHB*2/3);
windowOF = round(fs*0.04);
% initialization
for i = 1:length(Rposition)
thisR = Rposition(i);
if i==1
fid_pks(i,4) = thisR;
fid_pks(i,6) = thisR+windowS;
elseif i==length(Rposition)
fid_pks(i,4) = thisR;
%(thisR+windowT) < length(ECG) && (thisR - windowP) >=1
fid_pks(i,2) = thisR-windowQ;
else
if (thisR+windowT) < length(ECG) && (thisR - windowP) >=1
% Q S peaks
fid_pks(i,4) = thisR;
[Sv,Sp] = min(ECG(thisR:thisR+windowS));
thisS = Sp + thisR-1;
fid_pks(i,5) = thisS;
[Qv,Qp] = min(ECG(thisR-windowQ:thisR));
thisQ = thisR-(windowQ+1) + Qp;
fid_pks(i,3)=thisQ;
% onset and offset detection
interval_q = ECG(thisQ-windowOF:thisQ);
[ ind ] = onoffset(interval_q,'on' );
thisON = thisQ - (windowOF+1) + ind;
interval_s = ECG(thisS:thisS+windowOF);
[ ind ] = onoffset( interval_s,'off' );
thisOFF = thisS + ind-1;
fid_pks(i,2) = thisON;
fid_pks(i,6) = thisOFF;
% % T and P waves detection
% lastOFF = fid_pks(i-1,6);
% nextON =
end
end
end
%%
% P,T detection
% Detection T waves first and distinguish the type of T waves
for i = 2:length(Rposition)-1
lastOFF = fid_pks(i-1,6);
thisON = fid_pks(i,2);
thisOFF = fid_pks(i,6);
nextON = fid_pks(i+1,2);
if thisON>lastOFF && thisOFF<nextON
Tzone = (thisOFF:(nextON-round((nextON-thisOFF)/3)));
Pzone = ((lastOFF+round(2*(thisON-lastOFF)/3)):thisON);
[Tpks,thisT] = max(ECG(Tzone));
[Ppks,thisP] = max(ECG(Pzone));
fid_pks(i,1) = (lastOFF+round(2*(thisON-lastOFF)/3)) + thisP -1;
fid_pks(i,7) = thisOFF + thisT -1;
end
end
ECGpeaks = [];
Ind = zeros(length(Rposition),1);
for i = 1:length(Rposition)
Ind(i) = prod(fid_pks(i,:));
if Ind(i) ~= 0
ECGpeaks = [ECGpeaks;fid_pks(i,:)];
end
end
end
% Tpks = [];
% Tposition =[];
% Ttype = [];
% Ppks = [];
% Pposition =[];
% if length(QRS_OFF) <= length(QRS_ON)
% for i = 1:length(QRS_OFF)-1
% lengthTzone = int64((QRS_ON(i+1)-QRS_OFF(i))*2/3);
% if lengthTzone <= 0
% lengthTzone = abs(lengthTzone);
% disp('abnormal heart beat');
% end
% if QRS_OFF(i)+lengthTzone-1 > length(ECG)
% Tzone = ECG(QRS_OFF(i):length(ECG));
% else
% Tzone = ECG(QRS_OFF(i):QRS_OFF(i)+lengthTzone-1);
% end
% % if length(max(Tzone))~=1
% % deb = 1;
% % end
% [Tpks(end+1),Tind] = max(Tzone);
% Tposition(end+1) = QRS_OFF(i)+Tind-1;
% lengthPzone = int64((-QRS_OFF(i)+QRS_ON(i+1))/3);
% if lengthPzone <= 0
% lengthPzone = abs(lengthPzone);
% disp('abnormal heart beat');
% end
% Pzone = ECG(QRS_ON(i+1)-lengthPzone-1:QRS_ON(i+1));
% [Ppks(end+1),Pind] = max(Pzone);
% Pposition(end+1) = QRS_ON(i+1)+Pind-lengthPzone;
% end
% else
% for i = 1:length(QRS_ON)-1
% lengthTzone = int64((QRS_ON(i+1)-QRS_OFF(i))*2/3);
% if QRS_OFF(i)+lengthTzone-1 > length(ECG)
% Tzone = ECG(QRS_OFF(i):length(ECG));
% else
% Tzone = ECG(QRS_OFF(i):QRS_OFF(i)+lengthTzone-1);
% end
% [Tpks(end+1),Tind] = max(abs(Tzone));
% Tposition(end+1) = QRS_OFF(i)+Tind-1;
% lengthPzone = int64((-QRS_OFF(i)+QRS_ON(i+1))/3);
% Pzone = ECG(QRS_ON(i+1)-lengthPzone-1:QRS_ON(i+1));
% [Ppks(end+1),Pind] = max(Pzone);
% Pposition(end+1) = QRS_ON(i+1)+Pind-lengthPzone;
% end
% end
% for i = 1:length(Tpks)
% if Tpks(i)>0
% display('T Upright');
% else
% display('T Inverted');
% end
% end
% hold on;plot(Tposition,Tpks,'o');
% hold on;plot(Pposition,Ppks,'*');
%%
% [pks,locs] = findpeaks(A5);hold on; plot(locs,pks,'*');
%
% for i = 1:length(S)
% SearchArea = [];
% for j = 1:length(locs)
% if locs(j) > S(i)
% SearchArea(end+1) = locs(j);
% end
% end
% su = [];
% for k = 1:length(SearchArea)
% su(end+1) = SearchArea(k)-S(i);
% end
% if ~isempty(su)
% Tpks(end+1) = S(i) + min(su);
% end
% end
%
% % P detection
%
%
% for i = 1:length(Q)
% SearchArea = [];
% for j = 1:length(locs)
% if locs(j) < Q(i)
% SearchArea(end+1) = locs(j);
% end
% end
% su = [];
% for k = 1:length(SearchArea)
% su(end+1) = Q(i)-SearchArea(k);
% end
% if ~isempty(su)
% Ppks(end+1) = Q(i) - min(su);
% end
% end
% % check the result
% for i = 1:length(Tpks)
% T(i) = A5(Tpks(i));
% end
% for i = 1:length(Ppks)
% P(i) = A5(Ppks(i));
% end
% figure(5);plot(A5);
% hold on;plot(Tpks,T,'*');
% hold on;plot(Ppks,P,'+');
%%
%Oops after interpolation, we slightly deviate from the true value
%Thus we need to search for the local minimum in the whole signal again
% Qposition = int64(Q.*Scale);
% Sposition = int64(S.*Scale);
%
% for i = 1:length(R)
% [Sv,Sp] = min(ECG(Sposition(i):Sposition(i)+20));
% Sposition(i) = Sp + Sposition(i)-1;
% [Qv,Qp] = min(ECG(Qposition(i)-20:Qposition(i)));
% Qposition(i) = Qposition(i)-21 + Qp;
% end
%
% Pposition = int64(Ppks.*Scale);
% Tposition = int64(Tpks.*Scale);
%
% for i = length(Pposition)
% [Pv,Pp] = max(ECG(Pposition(i)-20:Pposition(i)));
% Pposition(i) = Pposition(i)-21 + Pp;
% end
%
% for i = length(Tposition)
% [Tv,Tp] = max(ECG(Tposition(i):Tposition(i)+20));
% Tposition(i) = Tp + Tposition(i)-1;
% end
% ecgS = [];
% for i = 1:2271
% ecgS(end+1) = issorted(ECGpeaks(1,:));
% end
% prod(ecgS)
Datasets
Models
