function [Ps,Pm,Pe,Qs,Qm,Qe,Rm,Sm,Se,Ts,Tm,Te]=my_ecg_fun(II,fs)

% Band pass filtering
[bbp,abp]=butter(2,2*[0.5 40]/fs,'bandpass');
fII=filtfilt(bbp,abp,II);

% Differentiaitiun and power
dII=diff([fII(1); fII]);
dII_2=dII.^2;

Tma=round(0.05*fs); % 50 ms MA filter
b=ones(Tma,1)/Tma; a=1;
madII_2=filtfilt(b,a,dII_2); % energic envelope

th=0.18*max(madII_2); % 10% of MA envelope

R_suspect=madII_2>th; % suspect section with R peak
start=find(diff([0;R_suspect])>0); % zeros before, >0 
stop=find(diff([R_suspect;0])<0); % zeros after, <0

% local maxima detection in suspect R-segment ------------------------------
Rm=zeros(length(start),1);
for i=1:length(start)
    II_segment=fII(start(i):stop(i)); % ECG 1-40 Hz
    [~,pos]=max(II_segment);
    Rm(i)=start(i)+pos-1;
end



% waves parametrization in R-R intervals ---------------------------------
for i=1:length(Rm)-1
    % IN R-R INTERVAL
    rr=Rm(i):Rm(i+1)-1; % idexes of |R-R| intervals
    seg=fII(rr); % % R-R segment II 0.5-40 Hz

        
    % Q, S waves
%     df=diff([seg(1); seg]); % 1. derivation (local extreme=0 for continuous signal)
%     sgn=sign(df); % decrease to increase change (local minimum for digital signals)
%     S=find(diff(sgn)>0,1); % 1st local minimum
%     Q=find(diff(sgn)>0,1,'last'); % last local minimum
    
    % T, P waves
%     [~,T]=max(seg(2*S:floor(end/2))); % maximum between S-minimum and |R-R|/2
%     T=T+2*S-1;
%     [~,P]=max(seg(floor(end/2):Q)); % maximum between |R-R|/2 and Q-minimum
%     P=P+floor(length(seg)/2)-1;

    
    % start/end of T, P
%     Tstart=triangle_detector(seg(2*S:T)); % double of R-S interval (2*S) and T-maximum
%     Tstart=Tstart+2*S-1;
%     Tend=triangle_detector(seg(T:floor(end/2))); % triangle between T and |R-R|/2
%     Tend=Tend+T-1;
    
%     Pstart=triangle_detector(seg(floor(end/2):P)); % trinagle between |R-R|/2 and P-maximum
%     Pstart=Pstart+floor(length(seg)/2)-1;
%     QR=length(seg)-Q; % QR interval
% 
%     % exception: if 2*|Q-R| is earlyer than P, use only 1*|Q-R|
%     if ((length(seg)-2*QR)-P)<1
%         Pend=triangle_detector(seg(P:end-QR)); % triangle between P-maximum and R-1*|QR|
%         Pend=P+Pend-1;
%     else
%         Pend=triangle_detector(seg(P:end-2*QR)); % % triangle between P-maximum and R-2*|QR|
%         Pend=P+Pend-1;
%     end
    
    % start/end of Q, S
%     Qstart=triangle_detector(seg(Pend:Q)); % triangle between P-end and Q-minimum
%     Qstart=Qstart+Pend-1;
%     
%     Qend=find(seg(Q:end)>=seg(Qstart),1); % Q-end is first corssing of Q-start isoline 
%     Qend=Qend+Q-1;
%     
%     Send=triangle_detector(seg(S:Tstart)); % triangle between S-minimum and T-start
%     Send=S+Send-1;
    
    % IN SIGNAL INDEXING -------------------------------------------- 
%     Ps(i,1)=Pstart+Rm(i)-1;
%     Pm(i,1)=P+Rm(i)-1;
%     Pe(i,1)=Pend+Rm(i)-1;
%     Qs(i,1)=Qstart+Rm(i)-1;
%     Qm(i,1)=Q+Rm(i)-1;
%     Qe(i,1)=Qend+Rm(i)-1;
%     Sm(i,1)=S+Rm(i)-1;
%     Se(i,1)=Send+Rm(i)-1;
%     Ts(i,1)=Tstart+Rm(i)-1;
%     Tm(i,1)=T+Rm(i)-1;
%     Te(i,1)=Tend+Rm(i)-1;
    Ps(i,1)=0;
    Pm(i,1)=0;
    Pe(i,1)=0;
    Qs(i,1)=0;
    Qm(i,1)=0;
    Qe(i,1)=0;
    Sm(i,1)=0;
    Se(i,1)=0;
    Ts(i,1)=0;
    Tm(i,1)=0;
    Te(i,1)=0;
end
end




function [bx_max,S]=triangle_detector(segment)
 
ax=1; % x-position of vertex A
ay=segment(ax); % y-position of vertex A
cx=length(segment); % x-position of vertex C
cy=segment(cx); % y-position of vertex C
% triangle area S=f(B)
S=zeros(length(segment),1);
for k=1:length(segment) % from A to C
    bx=k; % x-position of vertex B
    by=segment(k); % y-position of vertex B
    S(k)=0.5*abs((cx-ax)*(by-ay)-(cy-ay)*(bx-ax)); % area
end
[~,bx_max]=max(S); % x-position of vertex B correspond to the maximal area
end