--- a
+++ b/data preprocessing_Matlab/qsPeaks.m
@@ -0,0 +1,210 @@
+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)
\ No newline at end of file