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)