%% This code is for the Physionet 2017 challenge to distinguish Normal, Atrial fibrillation, other and noisy ECG rhtyhms.

 
close all
clear all
clc
 
load('allrawdataPhysionet2017.mat') ;

%% Resample and filter the data for outlier detection
sfreq=200;
df=designfilt('lowpassfir','FilterOrder',35,'PassbandFrequency',1,'StopbandFrequency',35,'SampleRate',sfreq);

[trash idx]=sort([allrawdata{:,2}],'ascend');
new=allrawdata(idx,:);

myx1={};
% myx={};
for i=1:size(allrawdata,1)          
    
    % Resample data to 200 Hz
    myx=resample(new{i,3},sfreq,300);
    
    % Filter data
    myx1(i,1) = {filtfilt(df,myx)};
end
 
%% Check if the ECG signals are reversed
 
% Check for the outlier data and remove the outliers
cop_myx1=myx1;
for i=1:size(allrawdata,1) 
    
    %Percentile range
    percentiles = prctile(cop_myx1{i,1},[0.06 99.95]); 
 
    outlierIndex =cop_myx1{i,1} < percentiles(1) | cop_myx1{i,1} > percentiles(2);
    
    %Remove outlier values
    cop_myx1{i,1}(outlierIndex) = [];
end
 
% Check if there is any lead reversal?
for i=1:size(allrawdata,1) 
 
% Find the location and take the mean of data>0
r=find(cop_myx1{i,1}>0);
pos_avg=mean(cop_myx1{i,1}(1,r));
 
% Find the location and take the mean of data<0
v=find(cop_myx1{i,1}<0);
neg_avg=mean(cop_myx1{i,1}(1,v));
 
% Flag 0 for correct , flag 1 for reversed
if abs(pos_avg)> abs(neg_avg)
    count1(i,:)=[i 0];    
else
    count1(i,:)=[i 1];   
end
end
 
%% Correct the sign of the signals based on the count1 values
 
final_alldata={};
for i=1:size(allrawdata,1) 
 
    if count1(i,2)==0
    final_alldata(i,1)={myx1{i,1}};
    else
    final_alldata(i,1)={-myx1{i,1}};
    end
end
 
 
%% New labels
% label=[zeros(5050,1); ones(738,1); 2*ones(2456,1); 3*ones(284,1)];


% Remove invalid data

cop_final_alldata=final_alldata;
 
for i=1:size(final_alldata,1)
    percentiles = prctile(cop_final_alldata{i,1},[0.1 99]); 

    outlierIndex =cop_final_alldata{i,1} < percentiles(1) | cop_final_alldata{i,1} > percentiles(2);
    
    %Remove outlier values
    cop_final_alldata{i,1}(outlierIndex) = []; 
end


%% Feature Extraction


% Feature 1: For noise Crest factor
Px1=[];
crest=[];
for i=1:8528
    
    % crest(i,:)=peak2rms(cop_final_alldata{i,1});
    [Px,Fx]=pwelch(cop_final_alldata{i,1},hamming(sfreq),50,sfreq,sfreq);
    Px1(:,i)=20*log10(Px);

end

% Normalize PSD
rat=[];
for i=1:8528
    
    % Px1(:,i)=Px(:,i)/(max(Px(:,i)));
    ind=find(Fx(:,1)>=5 & Fx(:,1)<=40);
    ind2=find(Fx(:,1)>=0 & Fx(:,1)<5);
    rat(i,:)=sum(Px1(ind,i))./sum(Px1(ind2,i));
end


% Feature 2: For noise- (Power in-band 5-40 Hz) / (Out-of-band spectral power)

% Feature 3: Segmented Sample entropy of each series 
 N = sfreq; % 1 sec 
 samp={};
 s_f=[];
 
% Find available no. of local workers
% myCluster = parcluster('local');
mypool=parpool(11);
 
for i=1:size(allrawdata,1)
    
   i
    L=length(cop_final_alldata{i,1});
    parfor j=1:L-(N)
        %     s=msentropy(final_filtered_all{i,1}(1,j:(j+N))',[],[],[],[],[],minM,maxM,1,minR,maxR);
        s=sampenc(cop_final_alldata{i,1}(1,j:(j+N))',3,0.25);
        s=s(1,1);
        s_f=[s_f s];
    end
    samp(i,:)={s_f};
    s_f=[];
    
end
 
delete(mypool)

%  save('Sampen_validation_Fs300.mat','samp')

% Feature 3 : Descriptive features on SE values
f31=[];
for i=1:8528
       samp{i,1}(isinf(samp{i,1})) = 0;
       data= samp{i,1};
       f31(i,1)=mean(data); 
       f31(i,2)=std(data);
       f31(i,3)=kurtosis(data);
       f31(i,4)=skewness(data);
       f31(i,5)=max(data)-min(data);
       f31(i,6)=median(data);
end


% Feature 4 : Descriptive features on ECG values
f41=[];
for i=1:8528
      
       data= cop_final_alldata{i,1};
       f41(i,1)=mean(data); 
       f41(i,2)=std(data);
       f41(i,3)=kurtosis(data);
       f41(i,4)=skewness(data);
       f41(i,5)=max(data)-min(data);
       f41(i,6)=median(data);
       f41(i,7)=mad(data);
       f41(i,8)=mad(data,1);
end


% Feature 5 : Descriptive features on RR interval data
f51=[];
for i=1:8528
       
       [~,qrs_i_raw,~]=pan_tompkin(new{i,3},sfreq,0);
       data= diff(qrs_i_raw)/sfreq;
       f51(i,1)=mean(data); 
       f51(i,2)=std(data);
       f51(i,3)=kurtosis(data);
       f51(i,4)=skewness(data);
       temp=max(data)-min(data);
       if isempty(temp)
           f51(i,5)=0;
       else
           f51(i,5)=temp;
       end
       f51(i,6)=median(data);
       f51(i,7)=mad(data);
       f51(i,8)=mad(data,1);
       
end

rr={};
for i=1:8528
       
       [~,qrs_i_raw,~]=pan_tompkin(new{i,3},sfreq,0);
       rr(i,1)= {diff(qrs_i_raw)/sfreq};
             
end

% Feature 6 : Descriptive features on delta RR interval data
f61=[];
deltarr1={};
for i=1:8528
       
       [~,qrs_i_raw,~]=pan_tompkin(new{i,3},sfreq,0);
       data= diff(qrs_i_raw)/sfreq;
       deltarr=diff(data);
       deltarr1(i,:)={deltarr};
       f61(i,1)=mean(deltarr); 
       f61(i,2)=std(deltarr);
       f61(i,3)=kurtosis(deltarr);
       f61(i,4)=skewness(deltarr);
       temp=max(deltarr)-min(deltarr);
       if isempty(temp)
           f61(i,5)=0;
       else
           f61(i,5)=temp;
       end
       f61(i,6)=median(deltarr);
       
end

% Feature 7: Descriptive features on 2nd order delta RR interval data
f71=[];
for i=1:8528
       
       [~,qrs_i_raw,~]=pan_tompkin(new{i,3},sfreq,0);
       data1= diff(qrs_i_raw)/sfreq;
       data=diff(diff(data1));
       f71(i,1)=mean(data); 
       f71(i,2)=std(data);
       f71(i,3)=kurtosis(data);
       f71(i,4)=skewness(data);
       temp=max(data)-min(data);
       if isempty(temp)
           f71(i,5)=0;
       else
           f71(i,5)=temp;
       end
       f71(i,6)=median(data);
       
end

% Feature 8 : Quantile based mean of less than and above 5% data

f81=[];
for i=1:8528
       
       [~,qrs_i_raw,~]=pan_tompkin(new{i,3},sfreq,0);
       data= diff(qrs_i_raw)/sfreq;
       if isempty(data)
            f81(i,1:4)=0;
       else    
           ran=quantile(data,10);
           ind1=find(data<ran(1,1)); 
           f81(i,1)=mean(data(:,ind1));
           ind2=find(data>ran(1,10)); 
           f81(i,2)=mean(data(:,ind2)); 
           f81(i,3)=min(data);
           f81(i,4)=max(data);
       end
end
save Features_200fs.mat

% Feature 9 : Wavelet based 
    
  f91=[];  
for i=1:8528
          [c,l] = wavedec(new{i,3},5,'sym4');

          [cd1,cd2,cd3,cd4,cd5] = detcoef(c,l,[1 2 3 4 5]);
           cd_2(i,1)={cd2};

          f91(i,[1:3])=[var(cd3) var(cd4) var(cd5)];
          y1=autocorr(cd3);
          y2=autocorr(cd4);
          y3=autocorr(cd5);
          f91(i,[4:6])=[var(y1) var(y2) var(y3)];
%           f91(i,5)=min(cd2)/max(cd2);
%           f91(i,6)=min(cd1)/max(cd1);
%           f91(i,[7:10])=[mad(cd3) mad(cd4) mad(cd5,1) mad(cd2,1)];
end

% Feature 10 : PSPR based on ECG
  f10=[];  
  % %% Load the PAFE reference database
load('Ref_dB_8Hz.mat');
Ref_db=patterndb_ref;


% Load the reference pattern info
load('Ref_patterninfo_8Hz.mat');
matchedpatterninfo_ref=matchedpatterninfo;

%% Resample data to 8 Hz for PSPR
cop_final_alldata1={};

for i=1:8528
    
   cop_final_alldata1(i,1)={resample(cop_final_alldata{i,1},8,200)}; 
    
end


for i=1:8528
i
   [f10(i,1) f10(i,2)]=pspr(cop_final_alldata1{i,1}',Ref_db,matchedpatterninfo_ref);

end


% Feature 11 : PSPR based on RR intervals
%   f11=[];    
%   
% for i=1:8528
%     i
%       f11(i,1:8)=pspr_rr(rr{i,1});
% 
% end

% % Feature 12 : mscoherence
%   f12=[];    
% for i=1:8528
%     i
%    [Cx,f]=mscohere(cop_final_alldata{1,1}(1,1:2700),cop_final_alldata{i,1}(1,1:2700),hamming(300),50,300,300);
%    ind =find(f(:,1)>=5 & f(:,1)<=25);
%    f12(i,1)=mean(Cx(ind,1));
% end

% grp=[zeros(1,5050-25) ones(1,738-25) 2*ones(1,2456-25) 3*ones(1,284-25)]' ;
% T= table([real(f31([26:5050,5076:5788,5814:8244,8270:8528],:)) real(f41([26:5050,5076:5788,5814:8244,8270:8528],:)) real(f51([26:5050,5076:5788,5814:8244,8270:8528],:)) ...
%     real(f61([26:5050,5076:5788,5814:8244,8270:8528],:)) real(f71([26:5050,5076:5788,5814:8244,8270:8528],:)) real(f81([26:5050,5076:5788,5814:8244,8270:8528],:)) ....
%     real(f91([26:5050,5076:5788,5814:8244,8270:8528],[6:10]))  real(f11([26:5050,5076:5788,5814:8244,8270:8528],:)) ...
%     real(rat([26:5050,5076:5788,5814:8244,8270:8528],:)) real(grp)]);

%     real(f10([26:5050,5076:5788,5814:8244,8270:8528],:)) 
    
% T= table([real(f31) real(f41) real(f51)  real(f61) real(f71) real(f81) real(f91) real(f10) nnfeatures real(rat) real(grp)]);


% save analysis.mat
% 
% features= [real(grp) real(f31) real(f41) real(f51) real(f61) real(f71) real(f81) real(f91(:,1:2)) real(rat) ];
% save('features_phase1.mat','features')
                                            % remove  3,9,19,24,25 and 30th features
      
                                            %% Second set of features
% features1=[real(grp) real(f31(:,[1:2,4:end])) real(f41(:,[1:2,4:end])) real(f51) real(f61(:,2:5)) real(f71(:,2:5)) real(f81) real(f91) real(rat) nnfeatures];
% save('features1_phase1.mat','features1')

% features2=[real(grp) real(f31(:,[1:2,4:end])) real(f41(:,[1:2,4:end])) real(f51) real(f61(:,2:5)) real(f71(:,2:5)) real(f81) real(f91) real(rat) nnfeatures];
% 
% 
% try1=[real(f31([26:5050,5076:5788,5814:8244,8270:8528],:)) real(f41([26:5050,5076:5788,5814:8244,8270:8528],:)) real(f51([26:5050,5076:5788,5814:8244,8270:8528],:)) ...
%     real(f61([26:5050,5076:5788,5814:8244,8270:8528],:)) real(f71([26:5050,5076:5788,5814:8244,8270:8528],:)) real(f81([26:5050,5076:5788,5814:8244,8270:8528],:)) ....
%     real(f91([26:5050,5076:5788,5814:8244,8270:8528],:)) real(f10([26:5050,5076:5788,5814:8244,8270:8528],:)) real(f11([26:5050,5076:5788,5814:8244,8270:8528],:)) ...
%     real(rat([26:5050,5076:5788,5814:8244,8270:8528],:)) real(grp)];
% 
% 
% view(trainedClassifier3.ClassificationEnsemble.Trained{3,1},'Mode','graph');
%  3     6     9    23    28    29    34    36    39    40    41    42    43    
%  3     9    12    15    20    23    28    29    34  `  37    42    59
% 3     7     9    12    15    20    23    28    29    34    42    44    59    60
%  3     4     9    10    14    16    19    21    22    23    24    27    28    29    30    33    34    36    37    38    39    40    41    45    59    60

% j=1;
% for i=5051:5060
%        subplot(10,1,j); plot(rr{i,1});
%         j=j+1;
% end
% 
% figure;
% j=1;
% for i=5791:5800
%        subplot(10,1,j); plot(rr{i,1});
%         j=j+1;
% end
% 
%   f12=[];  
% for i=1:8528
% i
%          nor_PSD= Px1(:,i)./max(Px1(:,i));
%          f12(i,1) = wentropy(nor_PSD,'log energy');
% 
% end
% 
% figure;
% cdfplot(rr{19,1})
% hold on;
% cdfplot(rr{20,1});hold on;
% cdfplot(rr{21,1});hold on;
% 
% cdfplot(rr{5056,1});hold on;
% cdfplot(rr{5057,1});hold on;
% cdfplot(rr{5058,1});


%% for standard histogram analysis
% newrr=deltarr1;
% l=cellfun('length',deltarr1);
% in=find(l<30);
% finalrr=[newrr num2cell([zeros(1,5050) ones(1,738) 2*ones(1,2456) 3*ones(1,284)]')];
% finalrr(in,:)=[];
% 
%% Make the reference database

% for 50 NN
% refNN_rr=finalrr(1:25,1);
% figure;
% for i=1:20
% subplot(20,1,i)
%     hist(refNN_rr{i,1}(:,1:30)*1000,20);
% end
% 
% % for 50  AF
% refAF_rr=finalrr(3701:3725,1);
% figure;
% for i=1:20
% subplot(20,1,i)
%     hist(refAF_rr{i,1}(:,1:30)*1000,20);
% end
% 
% % for 50 OT
% refOT_rr=finalrr(4270:4294,1);
% figure;
% for i=1:20
% subplot(20,1,i)
%     hist(refOT_rr{i,1}(:,1:30)*1000,20);
% end
% 
% % for 50 NO
% refNO_rr=finalrr(6055:6079,1);
% figure;
% for i=1:20
% subplot(20,1,i)
%     hist(refNO_rr{i,1}(:,1:30)*1000,20);
% end
% 
% grp1=[zeros(1,3675) ones(1,544) 2*ones(1,1760) 3*ones(1,108)]';
% T= table([real(f12) real(grp1)]);
%% Compare the CDF with reference CDF

% h=[];
% p1=[];p2=[]; p3=[]; p4=[];
% f12=[];
% k=1;
% for i=[26:3700,3726:4269,4295:6054,6080:6187]
%     for j=1:20
%         data= rr{i,1}*1000;
%         if isempty(data)
%             f12(k,1:4)=[0];
%             f12(k,5:8)=[0];
%         else
%             [h(j,1),p1(j,1)] = kstest2(data, refNN_rr{j,1}(:,1:30)*1000);
%             [h(j,1),p2(j,1)] = kstest2(data, refAF_rr{j,1}(:,1:30)*1000);
%             [h(j,1),p3(j,1)] = kstest2(data, refOT_rr{j,1}(:,1:30)*1000);
%             [h(j,1),p4(j,1)] = kstest2(data, refNO_rr{j,1}(:,1:30)*1000);
%         end
%     end
%         f12(k,1:4)=[mean(p1) mean(p2) mean(p3) mean(p4)];
%         f12(k,5:8)=[median(p1) median(p2) median(p3) median(p4)];
%         k=k+1;
%         p1=[];
%         p2=[]; p3=[]; p4=[];
% end
% 
% newtry1=try1;
% newtry1(in,:)=[];

% for 50 NN
refNN_rr=rr(1:50,1);
save('refNN_rr_200Hz.mat','refNN_rr');
% figure;
% for i=1:20
% subplot(20,1,i)
%     hist(refNN_rr{i,1}(:,1:30)*1000,20);
% end

% for 50  AF
refAF_rr=rr(5051:5100,1);
save('refAF_rr_200Hz.mat','refAF_rr');
% figure;
% for i=1:20
% subplot(20,1,i)
%     hist(refAF_rr{i,1}(:,1:30)*1000,20);
% end

% for 50 OT
refOT_rr=rr(5789:5838,1);
save('refOT_rr_200Hz.mat','refOT_rr');
% figure;
% for i=1:20
% subplot(20,1,i)
%     hist(refOT_rr{i,1}(:,1:30)*1000,20);
% end

% for 50 NO
refNO_rr=rr(8245:8294,1);
save('refNO_rr_200Hz.mat','refNO_rr');
% figure;
% for i=1:20
% subplot(20,1,i)
%     hist(refNO_rr{i,1}(:,1:30)*1000,20);
% end
% 
h=[];
p1=[];p2=[]; p3=[]; p4=[];
f12=[];
k=1;
for i=1:8528
    for j=1:50
        data= rr{i,1}*1000;
        if isempty(data)
            f12(k,1:4)=[0];
            f12(k,5:8)=[0];
        else
            [h(j,1),p1(j,1)] = kstest2(data, refNN_rr{j,1}*1000);
            [h(j,1),p2(j,1)] = kstest2(data, refAF_rr{j,1}*1000);
            [h(j,1),p3(j,1)] = kstest2(data, refOT_rr{j,1}*1000);
            [h(j,1),p4(j,1)] = kstest2(data, refNO_rr{j,1}*1000);
        end
    end
        f12(k,1:4)=[mean(p1) mean(p2) mean(p3) mean(p4)];
        f12(k,5:8)=[median(p1) median(p2) median(p3) median(p4)];
        k=k+1;
        p1=[];
        p2=[]; p3=[]; p4=[];
end

% [26:5050,5076:5788,5814:8244,8270:8528]
% 
% grp1=[zeros(1,5050-25) ones(1,738-25) 2*ones(1,2456-25) 3*ones(1,284-25)]';
% 
% try1=[real(f31([26:5050,5076:5788,5814:8244,8270:8528],:)) real(f41([26:5050,5076:5788,5814:8244,8270:8528],1:6)) real(f51([26:5050,5076:5788,5814:8244,8270:8528],1:6)) ...
%     real(f61([26:5050,5076:5788,5814:8244,8270:8528],:)) real(f71([26:5050,5076:5788,5814:8244,8270:8528],:)) real(f81([26:5050,5076:5788,5814:8244,8270:8528],:)) ...
%     real(f91([26:5050,5076:5788,5814:8244,8270:8528],1:6)) real(f10([26:5050,5076:5788,5814:8244,8270:8528],:)) real(f12([26:5050,5076:5788,5814:8244,8270:8528],:)) ...
%      real(rat([26:5050,5076:5788,5814:8244,8270:8528],:))  x([26:5050,5076:5788,5814:8244,8270:8528],1) real(grp1)];
% save latestfeaturesJuly6.mat

% 
% x=(cent(:,1)+cent(:,2))/2;
% 
grp1=[zeros(1,5050) ones(1,738) 2*ones(1,2456) 3*ones(1,284)]';
% 
% try1=[real(f31) real(f41(:,1:6)) real(f51(:,1:6)) ...
%     real(f61) real(f71) real(f81) ...
%     real(f91(:,1:6)) real(f10) real(f12) ...
%      real(rat)  real(R) real(per)  real(perc) real(grp1)];
% 
%  
 
 %% KERNEL DENSITY ESTIMATION
% [f,xi] = ksdensity(rr{5788,1});
% figure
% plot(xi,f);


%% CV model on the 80% data
% try2=try1([1:4040,5051:5640,5789:7752,8245:8471],:);
% try2_test=try1([4041:5050,5641:5788,7753:8244,8472:8528],:);
% 
% grp3=[zeros(1,1010) ones(1,148) 2*ones(1,492) 3*ones(1,57)]';
% 
% yfit = predict(trainedClassifier.ClassificationEnsemble,try2_test(:,1:51)) ;
% yfit1 = predict(trainedClassifier1.ClassificationEnsemble,try2_test(:,1:51)) ;
% yfit2 = predict(trainedClassifier2.ClassificationEnsemble,try2_test(:,1:51)) ;
% yfit3 = predict(trainedClassifier3.ClassificationEnsemble,try2_test(:,1:51)) ;
% yfit4 = predict(trainedClassifier4.ClassificationEnsemble,try2_test(:,1:51)) ;
% 
% acc=confusionmat(grp3,yfit)
% 100*sum(diag(acc))/sum(sum(acc))

%% F13 Bispectrum analysis

% morder = 12;
% figure;
% pyulear(cop_final_alldata{5826, 1},morder,[],100);hold on;
% pyulear(cop_final_alldata{1, 1},morder,[],100);hold on;
% pyulear(cop_final_alldata{100, 1},morder,[],100);



% bispeci(cop_final_alldata{8500,1},21,50,0,'unbiased',100,1);
% 

% a=[];
% for i=1:6000
%     i 
% %     figure;
%     p=arorder(cop_final_alldata{i,1},3);
% %    
%     [~,Bspec,val,f1,f2] = qpctor(cop_final_alldata{i,1},18,p,50);
%     a(i,:)=f1;
% %    ln_bi(i,:)=mean2( abs(Bspec(1:20,1:20)));
% %     bisp(i,:)=val;
% %  [~,~,val] =bicoher (cop_final_alldata{i,1}, 50,100, 100, 50);
% %   bicoh(i,:)=val;
% end
 

%% Poincare plot
% cent=[];
% for i=1:8528
% %     figure
% %     for i=1:length(rr{j,:})-1
% %        plot(rr{j,:}(1,i),rr{j,:}(1,i+1),'b.');
% % %         lsline    % the diagonal line
% %        hold on;
%         if isempty(rr{i,:}) || length(rr{i,:})<=1
%                 cent(i,1:2)=[0];
%         else
%                [idx,C]= kmeans([rr{i,:}(1,1:length(rr{i,:})-1)',rr{i,:}(1,2:length(rr{i,:}))'],1);
%                cent(i,1:2)=[C];
%         end
% %        plot(rr{i,:}(1,1:length(rr{i,:})-1),rr{i,:}(1,2:length(rr{i,:})),'b.');hold on;
% %        plot(C(:,1),C(:,2),'kx');
% %     end
% end

%% Coefficient of determination from Poincare plot
R=[];
for i=1:8528
    
            if isempty(rr{i,:}) || length(rr{i,:})<=2
%                 cod(i,:)=0;
                R(i,:)=0;
            else
                  R12= corrcoef([rr{i,:}(1,1:length(rr{i,:})-1)',rr{i,:}(1,2:length(rr{i,:}))']);
                  r12=R12(1,2);
                  R(i,:)=r12;
%                   cod(i,:)=r12^2;
            end
end


%% PCA
perc=[];
for i=1:8528
     if isempty(rr{i,:}) || length(rr{i,:})<=2
        perc(i,:)=0;
     else
          [COEFF, SCORE, LATENT,TSQUARED, EXPLAINED] = pca([rr{i,:}(1,1:length(rr{i,:})-1)',rr{i,:}(1,2:length(rr{i,:}))']);
 %         d=cumsum(latent)./sum(latent);
%         per(i,:)=d(1,1);
         perc(i,:)=(EXPLAINED(1,1));
%         figure()
% pareto(explained)
% xlabel('Principal Component')
% ylabel('Variance Explained (%)')
% figure
% biplot(pc(:,1:2),'Scores',score(:,1:2),'VarLabels',...
% 		{'X1' 'X2'})
     end
end


%% SLope
% slop=[];
% for i=1:8528
%     
%             if isempty(rr{i,:}) || length(rr{i,:})<=2
%                 slop(i,:)=0;
%             else
%                 coefficients= polyfit(rr{i,:}(1,1:length(rr{i,:})-1),rr{i,:}(1,2:length(rr{i,:})),1);
%                 slop(i,:)=coefficients(1);
%                  
%             end
% end

%% Percentage of rr
per=[];
for i=1:8528
    ind=find(rr{i,:} >=0.6 & rr{i,:} <=1.2 );
    per(i,:)=length(ind)/length(rr{i,:})*100;
end

% try2=try1(:,[2,27,29,31,33,43,51,52,53]);
% y=grp1';
% x=try2';

%% Heart rate range
HR=[];
for i=1:8528
hr=60./(rr{i,:});
HR(i,:)=mad(hr);
end

try1=[real(f31) real(f41(:,1:6)) real(f51(:,1:6)) ...
real(f61) real(f71) real(f81) ...
real(f91(:,1:6)) real(f10) real(f12) ...
real(rat)  real(R) real(per) real(perc) real(HR) real(grp1)];

try2=try1(:,[5 6 7 12 14 16 18 19 21 22 23 24 26 27 28 29 30 33 34 37 41 42 44 45 47 48 49 50 52 54 55 56]);
save Features_200fs.mat

newrr=[];
for i=1:8528
  for j=2:3:length(rr{i,:})-1
      
      newrr(:,j)=1;   % Intialize with category 1

      if (rr{i,:}(:,j) < 0.6 && 1.8*rr{i,:}(:,j)<rr{i,:}(:,j-1))  % if C1
       newrr(:,j)=3;
        pulse=1;
      
            for k=3:3:length(rr{i,:})-1
                 if ((rr{i,:}(:,k-1)< 0.7) && (rr{i,:}(:,k) < 0.7) && (rr{i,:}(:,k+1)< 0.7) || (rr{i,:}(:,k+1)+rr{i,:}(:,k)+rr{i,:}(:,k-1) <1.7))

                    newrr(:,j)=3;
                    pulse=pulse+1;

                 end 
            end
                  if pulse<4
                     newrr(:,j)=1;
                  end
             
       
      elseif ((1.5*rr{i,:}(:,j) < rr{i,:}(:,j-1) && 1.15*rr{i,:}(:,j)<rr{i,:}(:,j+1)) || ((abs(rr{i,:}(:,j-1)-rr{i,:}(:,j))<0.3) && (rr{i,:}(:,j-1)<0.8 || rr{i,:}(:,j)<0.8) && (rr{i,:}(:,j+1)>1.2*(rr{i,:}(:,j) +rr{i,:}(:,j))./2)) ...
            || ((abs(rr{i,:}(:,j)-rr{i,:}(:,j+1))<0.3) && (rr{i,:}(:,j)<0.8 || rr{i,:}(:,j+1)<0.8) && (rr{i,:}(:,j-1)>1.2*(rr{i,:}(:,j) +rr{i,:}(:,j+1))./2))) % if C3,C4,C5
        
                newrr(:,j)=2;
        
      elseif ((2.2< rr{i,:}(:,j-1) <3.0) && (abs(rr{i,:}(:,j-1)-rr{i,:}(:,j))<0.2 || abs(rr{i,:}(:,j)-rr{i,:}(:,j+1))<0.2) )
                newrr(:,j)=4;
        
    end  
  end
end


% qqplot(rr{5700,:}(1,1:length(rr{5700,:})-1),rr{5700,:}(1,2:length(rr{5700,:})));

figure
ellipse_t = fit_ellipse( rr{2000,:}(1,1:length(rr{2000,:})-1),rr{2000,:}(1,2:length(rr{2000,:})),0)

figure
ellipse_t = fit_ellipse( rr{5500,:}(1,1:length(rr{5500,:})-1),rr{5500,:}(1,2:length(rr{5500,:})),0)

figure
ellipse_t = fit_ellipse( rr{7000,:}(1,1:length(rr{7000,:})-1),rr{7000,:}(1,2:length(rr{7000,:})),0)

figure
ellipse_t = fit_ellipse( rr{8500,:}(1,1:length(rr{8500,:})-1),rr{8500,:}(1,2:length(rr{8500,:})),0)

poin=[];
for i=1:8528
    if length(rr{i,:})<=5 || isempty(rr{i,:})
          poin(i,:)=0;
    else
         ellipse_t = fit_ellipse( rr{i,:}(1,1:length(rr{i,:})-1),rr{i,:}(1,2:length(rr{i,:})))  ; 
%             if isempty(ellipse_t.short_axis)
%                 poin(i,:)=0;
%             else
                poin(i,:)=(ellipse_t.long_axis)/(ellipse_t.short_axis);
%             end
    end
end

euc=[];
for i=1:8528
 
         euc(i,:) = dist(rr{i,:}(1,1:length(rr{i,:})-1),rr{i,:}(1,2:length(rr{i,:})))  ; 

end


%% For neural network

% change the output labels
yy=zeros(4,8528);
yy(1,y==0)=1;
yy(2,y==1)=1;
yy(3,y==2)=1;
yy(4,y==3)=1;

[ddd yyhat]=max(network3_outputs);
confusionmat(yy',yyhat')

%% Deep learning 

%Split data into 80 % training -20 % test
try1Train=try1([1:4040,5051:5640,5789:7753,8245:8471],1:56);
try1Test=try1([4041:5050,5641:5788,7754:8244,8471:8528],1:56);

%% Normalize the data

try1Train=normc(try1Train);
try1Test=normc(try1Test);


yyTrain=yy(:,[1:4040,5051:5640,5789:7753,8245:8471]);
yytest=try1([4041:5050,5641:5788,7754:8244,8471:8528],56);

% Develop the model
autoenc1 = trainAutoencoder(try1Train', 30, 'L2WeightRegularization', 0.001, 'SparsityRegularization', 4, 'SparsityProportion',0.05, 'MaxEpochs', 1000);

% autoenc2 = trainAutoencoder(try1(:,1:53)', 30, 'L2WeightRegularization', 0.001, 'SparsityRegularization', 4, 'SparsityProportion',0.05, 'MaxEpochs', 500);

%Extract the features from the network
features1 = encode(autoenc1, try1Train');

%Train a softmax classifier 
softnet = trainSoftmaxLayer(features1,yyTrain,'LossFunction','crossentropy');

%Create a classifier with the encoder and softmax layer 
deepnet = stack(autoenc1,softnet);

% Fine tune the network
deepnet = train(deepnet,try1Train',yyTrain);



% Run the prediction of the network in the test
output = deepnet(try1Test');

% Show the results in a confusion matrix
% plotconfusion(yytest,output)

% Squeeze it further
autoenc2 = trainAutoencoder(features1, 15, 'L2WeightRegularization', 0.001, 'SparsityRegularization', 4, 'SparsityProportion',0.05,'ScaleData',false);

% Extract the features
features2 = encode(autoenc2, features1);

% Train a softmax classifier
softnet = trainSoftmaxLayer(features2,yyTrain,'LossFunction','crossentropy');

% Stack the autoencoders
deepnet = stack(autoenc1,autoenc2,softnet);

% Fine tune the network
deepnet = train(deepnet,try1Train',yyTrain);

% % another layer
% % Squeeze it further
% autoenc3 = trainAutoencoder(features1, 30, 'L2WeightRegularization', 0.001, 'SparsityRegularization', 4, 'SparsityProportion',0.05,'ScaleData',false);
% 
% % Extract the features
% features3 = encode(autoenc3, features1);
% 
% % Train a softmax classifier
% softnet = trainSoftmaxLayer(features3,yyTrain,'LossFunction','crossentropy');
% 
% % Stack the autoencoders
% deepnet = stack(autoenc1,autoenc2,autoenc3,softnet);
% 
% % Fine tune the network
% deepnet = train(deepnet,try1Train',yyTrain);
% 
% % Reclassify 
% output = deepnet(try1Test');
% % Show the results in a confusion matrix
plotconfusion(yytest,output)
view(deepnet) 

%% conduct the t test on the bi coherenac plots obtained above

% %% RR indices
% [trash idx]=sort([allrawdata{:,2}],'ascend');
% new=allrawdata(idx,:);
% 
% 
% rr_ind={};
% for i=1:8528
%        
%        [~,qrs_i_raw,~]=pan_tompkin(new{i,3},300,0);
%        rr_ind(i,:)= {qrs_i_raw};
%             
% end
% 
% 
% cell2csv('rr_ind.csv',rr_ind(1:end,:))

 yfit =predict( trainedClassifier6.ClassificationEnsemble,try1Test(:,1:55)) ;
 confusionmat(yytest, yfit)