function [features] = ECG_features_158_old(ecg, fs, P_index, Q_index, R_index, S_index, T_index)

% 
% Copyright (C) 2017 
% Shreyasi Datta
% Chetanya Puri
% Ayan Mukherjee
% Rohan Banerjee
% Anirban Dutta Choudhury
% Arijit Ukil
% Soma Bandyopadhyay
% Rituraj Singh
% Arpan Pal
% Sundeep Khandelwal
% 
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.
% 
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
% 
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.


features = [];
noisefeat = [];
RR = [];
pde = [];
ecg_f = [];
QRS = [];
ploc = [];

if length(ecg)<300 || length(R_index) < 6 || isempty(Q_index) || isempty(S_index) || isempty(T_index) 
    features = 0;
    return
end

P_to_R_ratio = length(P_index)/length(R_index);

%% Invert ECG if necessary
amp = ecg(R_index);
if mean(amp) < 0 && median(amp) < 0
    ecg = -ecg;
end

%% Filter ECG
ecg_f = filter_ecg( ecg );

%% Find RR Interval and oulier Removal
QRS = R_index;
RR = diff(R_index')/fs;
%[RR, ind]  = outlier_removal(RR,3);

%% Plotting
% plot(ecg_f); hold on; 
% plot(P_index,ecg_f(P_index),'k^');
% plot(Q_index,ecg_f(Q_index),'c^');
% plot(R_index,ecg_f(R_index),'r*');
% plot(S_index,ecg_f(S_index),'m^');
% plot(T_index,ecg_f(T_index),'g^');
% legend('ECG','P','Q','R','S','T');

%% Feature Extraction

%% SoA Features from RR Interval
% AF score according to SoA
[AFEv, OriginCount, IrrEv, PACEv, DensityEv, AniEv] = comput_AFEv(RR);

% Coefficient of variation of RR and delta(RR)
CVrr = std(RR)/mean(RR);
CVdrr = std(diff(RR))/mean(RR);

% Poincare Features
poincare = poincare_features(RR);

% Autocorrelation Features
% [SimilarityMatrix, degree] = noise_autocorr(ecg_f, fs);
% degree1 = std(degree)/mean(degree);
% [count,d] = time_delay_features(ecg_f, fs);
       
% Sample Entropy Features
% se_feat = sample_entropy_features(ecg_f);

feat1 = [AFEv OriginCount IrrEv PACEv DensityEv AniEv CVrr CVdrr poincare];

%% Statistcal charateristics of RR Interval
mean_RR = mean(RR);
median_RR = median(RR);
min_RR = min(RR);
max_RR = max(RR);
sk_RR = skewness(RR);
kurt_RR = kurtosis(RR);
range_RR = max(RR) - min(RR);
var_RR = var(RR);

feat2 = [mean_RR median_RR min_RR max_RR sk_RR kurt_RR range_RR var_RR];

%% Wavelet Entropy
E = wentropy(RR,'shannon');
%[E,unique] = ShannonEn(RR,1,5)

%% Hjorth Parameters
[ACTIVITY, MOBILITY, COMPLEXITY] = hjorth(RR,0);

%% Features from the Probability Density Estimate of RR
% pts = 0.33 : 0.05 : 1.5;
[f,xi] = ksdensity(RR);
%f = mapminmax(f,0,1);
pde = [f' xi'];
kurt_f = kurtosis(f);
skew_f = skewness(f);
[pks,loc] = findpeaks(f);
num_pks = length(pks);
maxdist = 0;
mindist = 0;
if num_pks > 1
    [dompeak dompeakloc] = max(pks);
    pks1 = pks; loc1 = loc;
    pks1(dompeakloc) = [];
    loc1(dompeakloc) = [];
    for ii = 1 : length(pks1)
       distance(ii) = abs(xi(loc(dompeakloc))-xi(loc(ii)));
    end
    maxdist = max(distance);
    mindist = min(distance);
end
    
%% Features from the Probability Density Estimate of delta(RR)
dRR = diff(RR);
[f1,xi1] = ksdensity(dRR);
%f1 = mapminmax(f1,0,1);
pde1 = [f1' xi1'];
kurt_f1 = kurtosis(f1);
skew_f1 = skewness(f1);

feat3 = [kurt_f skew_f num_pks maxdist mindist kurt_f1 skew_f1];

%% Amplitude based features for noise
% Variance in amplitude of QRS 
% amplitude = ecg_f(QRS);
amplitude = ecg(QRS);
var_amp = std(amplitude)/mean(amplitude);

% Statistical features on amplitude of QRS
[f3,xi3] = ksdensity(amplitude);
pde = [f3' xi3'];
kurt_f3 = kurtosis(f3);
skew_f3 = skewness(f3);
[pks3,loc3] = findpeaks(f3);
num_pks3 = length(pks3);

% Variance of Energy between R peaks
for jj = 1 : length(QRS)-1
   temp(jj) = mean(ecg_f(QRS(jj)+1:QRS(jj+1)-1));
end
var_en = std(temp)/mean(amplitude);

feat4 = [kurt_f3 skew_f3 num_pks3 var_amp var_en];

%% HRV Features
diffNN = diff(RR);
RMSSD = sqrt(mean(diffNN.^2));
SDSD  = std(diffNN);
NN50 = length(find(abs(diffNN)>0.5));
pNN50 = NN50/length(RR);
NN20 = length(find(abs(diffNN)>0.2));
pNN20 = NN20/length(RR);

feat_HRV = [pNN20 pNN50 SDSD];
%% Cluster based features by ADC
% if length(dRR) < 3
%     fclust = 0;
% else
%     c = clusterdata(dRR,'maxclust',3,'linkage','median');
%     c1 = find(c==1);
%     c2 = find(c==2);
%     c3 = find(c==3);
%     f4 =  min([ std(dRR(c1))/median(c1)  std(dRR(c2))/median(c2)  std(dRR(c3))/median(c3)] );
%     m =  [median(dRR(c1)) median(dRR(c2)) median(dRR(c3))];
%     fclust = std(m);
% end

%% Depth of S (RS amplitude difference)
ecg_f(R_index(find(ecg_f(R_index) == 0))) = 0.00001;
rsfeat = abs(ecg_f(S_index)./ecg_f(R_index));
%rsfeat = rs./ecg_f(R_index);
[rsfeat_clean, r_ind] = outlier_removal(rsfeat,1);
rsfeat(r_ind) = [];
rs1 = std(rsfeat)/mean(rsfeat);
rs2 = max(rsfeat) - min(rsfeat);
rs3 = median(rsfeat);

%% Features from ST segment
% Slope of ST

for i = 1 : length(S_index)
   slope_ST(i) = (ecg_f(T_index(i)) - ecg_f(S_index(i)))/(T_index(i) - S_index(i));
end
slope_ST(find(isnan(slope_ST))) = [];

if ~isempty(slope_ST)
    slope_ST = outlier_removal(slope_ST,2);
end

if ~isempty(slope_ST)
    med_st = median(slope_ST);
    var_st = std(slope_ST)/mean(slope_ST);
    num_neg_st = length(find(slope_ST<0))/length(QRS);
else
    med_st = 0;
    var_st = 0;
    num_neg_st = 0;
end

% Crossing of ST segment w.r.t. baseline
deep_s = 0;
inflec_dist_s = 0;
mod_T_prev = [];
mod_S = [];
mod_T = [];
mod_Q = [];
deep_s1 = find(ecg_f(S_index)<ecg_f(Q_index));

if(~isempty(deep_s1))
    deep_s = length(deep_s1)/length(S_index);
    mod_S = S_index(deep_s1);
    mod_T = T_index(deep_s1);
    mod_Q = Q_index(deep_s1);
    
    for jj = 1 : length(deep_s1)
        if deep_s1(jj) == 1
            mod_T_prev(jj) = 1;
        else
            mod_T_prev(jj) = T_index(deep_s1(jj)-1);
        end
    end
    
    for z = 1 : length(mod_S)
        temp = ecg_f(mod_S(z):1:mod_T(z));
        if z == 1
            ind = find(temp>ecg_f(mod_Q(z)));
        else
            ind = find(temp>median(ecg_f(mod_T_prev(z):mod_Q(z))));
        end
        if isempty(ind)
            inflec_dist_s(z) = 0;
        else
            inflec_dist_s(z) = ind(1);
        end
    end
end
inflec_dist_s = median(inflec_dist_s);

%% TR Amplitude
tr_amp = ecg_f(T_index)./ecg_f(R_index);
% ecg_f(find(tr_amp<0)) = [];
if ~isempty(tr_amp)
    tr_amp = median(tr_amp);
else
    tr_amp = 0;
end

%% QS Width and QR Width
QRS_width = (S_index - Q_index)./fs;
QR_width = (R_index - Q_index)./fs;
Q1 = median(QRS_width);
Q2 = std(QRS_width)/mean(QRS_width);
Q3 = median(QR_width);
Q4 = std(QR_width)/mean(QR_width);

%% Height of R or Depth of Q
rqfeat = abs(ecg_f(Q_index)./ecg_f(R_index));
[rqfeat_clean, r_ind1] = outlier_removal(rqfeat,1);
rqfeat(r_ind1) = [];
rq1 = std(rqfeat)/mean(rqfeat);
rq2 = max(rqfeat) - min(rqfeat);
rq3 = median(rqfeat);

%% Ratio of Depth of S to Height of R (w.r.t. Q)
S_depth = ecg_f(Q_index) - ecg_f(S_index);
R_height = ecg_f(R_index) - ecg_f(Q_index);
SR_ratio = S_depth./R_height;
SR_ratio(find(isinf(SR_ratio))) = [];
SR_ratio(find(isnan(SR_ratio))) = [];
SR1 = median(SR_ratio);
SR2 = std(SR_ratio)/mean(SR_ratio);

%% Slope of QR / RS / Sx
% Slope of QR
for i = 1 : length(Q_index)
   slope_QR(i) = (ecg_f(R_index(i)) - ecg_f(Q_index(i)))/(R_index(i) - Q_index(i));
   slope_RS(i) = (ecg_f(S_index(i)) - ecg_f(R_index(i)))/(S_index(i) - R_index(i));
end

for i = 1 : length(S_index) - 1
    x = S_index(i) + 15;
    slope_Sx(i) = (ecg_f(x) - ecg_f(S_index(i)))/(x - S_index(i));
end

slope_QR(find(isnan(slope_QR))) = [];
slope_RS(find(isnan(slope_RS))) = [];
slope_Sx(find(isnan(slope_Sx))) = [];

if ~isempty(slope_QR)
    slope_QR = outlier_removal(slope_QR,2);
end
if ~isempty(slope_RS)
    slope_RS = outlier_removal(slope_RS,2);
end
if ~isempty(slope_Sx)
    slope_Sx = outlier_removal(slope_Sx,2);
end

if ~isempty(slope_QR)
    med_qr = median(slope_QR);
    var_qr = std(slope_QR)/mean(slope_QR);
else
    med_qr = 0;
    var_qr = 0;
end

if ~isempty(slope_RS)
    med_rs = median(slope_RS);
    var_rs = std(slope_RS)/mean(slope_RS);
else
    med_rs = 0;
    var_rs = 0;
end

if ~isempty(slope_Sx)
    med_sx = median(slope_Sx);
    var_sx = std(slope_Sx)/mean(slope_Sx);
else
    med_sx = 0;
    var_sx = 0;
end

%% QTc related features
QT = T_index - Q_index;
QT = QT(1:end-1);
sqrtRR = sqrt(RR);
cubrtRR = nthroot(RR,3);
QTcB = QT./sqrtRR'; % Bazett's Formula
QTcF = QT./cubrtRR'; % Fridericia's Formula
QTlc = QT + 0.154.*(1-RR'); % Sagie's Formula

QT1 = median(QTcB); QT2 = std(QTcB)/mean(QTcB);
QT3 = median(QTcF); QT4 = std(QTcF)/mean(QTcF);
QT5 = median(QTlc); QT6 = std(QTlc)/mean(QTlc);

feat5 = [rs1 rs2 rs3 med_st var_st num_neg_st deep_s inflec_dist_s tr_amp P_to_R_ratio ...
    Q1 Q2 Q3 Q4 rq1 rq2 rq3 SR1 SR2 med_qr var_qr med_rs var_rs med_sx var_sx QT1 QT2 QT3 QT4 QT5 QT6];

%%
features = [feat1 feat2 E MOBILITY COMPLEXITY feat3 feat4 feat5 feat_HRV];
% features = [AFEv OriginCount IrrEv PACEv DensityEv AniEv mean_RR median_RR min_RR max_RR sk_RR kurt_RR range_RR var_RR kurt_f skew_f num_pks CVrr CVdrr kurt_f1 skew_f1 poincare E MOBILITY COMPLEXITY maxdist mindist kurt_f3 skew_f3 num_pks3 var_amp var_en out rs1 rs2 rs3 med_st var_st num_neg_st inflec_dist_s tr_amp];


end

