function classifyResult = challenge(recordName)
%
% Sample entry for the 2017 PhysioNet/CinC Challenge.
%
% INPUTS:
% recordName: string specifying the record name to process
%
% OUTPUTS:
% classifyResult: integer value where
%                     N = normal rhythm
%                     A = AF
%                     O = other rhythm
%                     ~ = noisy recording (poor signal quality)
%
% To run your entry on the entire training set in a format that is
% compatible with PhysioNet's scoring enviroment, run the script
% generateValidationSet.m
%
% The challenge function requires that you have downloaded the challenge
% data 'training_set' in a subdirectory of the current directory.
%    http://physionet.org/physiobank/database/challenge/2017/
%
% This dataset is used by the generateValidationSet.m script to create
% the annotations on your training set that will be used to verify that
% your entry works properly in the PhysioNet testing environment.
%
%
% Version 1.0
%
%
% Written by: Chengyu Liu and Qiao Li January 20 2017
%             chengyu.liu@emory.edu  qiao.li@emory.edu
%
% Last modified by: Ruhi Mahajan Aug1 2017
%                   rmahajan@uthsc.edu 
%

%classifyResult = 'N'; % default output normal rhythm

%% AF determination

%[tm,ecg,fs,siginfo]=rdmat(recordName);
%[QRS,sign,en_thres] = qrs_detect2(ecg',0.25,0.6,fs);
%if length(QRS)<6
%else
%    RR=diff(QRS')/fs;
%    AFEv = comput_AFEv(RR);
%    if AFEv>1
%        classifyResult = 'A';
%    end
%end

%% You can add the determination rules for other rhythm and noisy recordings here

%UTHSC-ORNL CBMI TEAM
%% Prepare New ECG for Classification
[~,ecg,~,~]=rdmat(recordName);

%% This code is for the Physionet 2017 challenge to distinguish Normal, Atrial fibrillation, other and noisy ECG rhtyhms.
%   It loads following functions:
        % sampenc- to calculate sample entropy
        %
%   Authors Dr. Oguz Akbilgic and Dr.Ruhi Mahajan
%   Last Modified on April 7,2017
%   UTHSC-ORNL Center for Biomedical Informatics, Memphis,TN


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load('Model31compact.mat');
load('refNN_rr_200Hz.mat');
load('refAF_rr_200Hz.mat');
load('refOT_rr_200Hz.mat');
load('refNO_rr_200Hz.mat');

% %% 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;
%% Data preprocessing

allrawdata=1000*ecg';
%% Resampe and filter the data for outlier detection

sfreq=200;
df=designfilt('lowpassfir','FilterOrder',35,'PassbandFrequency',1,'StopbandFrequency',35,'SampleRate',sfreq);

% % Resample data to 200 Hz
myx1=resample(allrawdata,sfreq,300);

% Filter data
myx= filtfilt(df,myx1);% filter

%% Check if the ECG signals are reversed

% Check for the outlier data and remove the outliers
 copy_myx=myx;
 
%Percentile range
percentiles = prctile(copy_myx,[0.06 99.95]); 
outlierIndex =copy_myx < percentiles(1) |copy_myx > percentiles(2);

%Remove outlier values
copy_myx(outlierIndex) = [];

% Check if there is any lead reversal?

% Find the location and take the mean of data>0
r=find(copy_myx>0);
pos_avg=mean(copy_myx(1,r));

% Find the location and take the mean of data<0
v=find(copy_myx<0);
neg_avg=mean(copy_myx(1,v));

% Correct the sign of the signals based on the count1 values
if abs(pos_avg)<= abs(neg_avg)
    final_alldata=-myx;  
else
    final_alldata=myx; 
end
   
    cop_final_alldata=final_alldata;
 

    percentiles = prctile(cop_final_alldata,[0.1 99]); 

    outlierIndex =cop_final_alldata < percentiles(1) | cop_final_alldata > percentiles(2);
    
    %Remove outlier values
    cop_final_alldata(outlierIndex) = []; 
    
%% Feature extraction

% % Feature 1: For noise- (Power in-band 5-40 Hz) / (Out-of-band spectral power)
%     
%     [Px,Fx]=pwelch(cop_final_alldata,hamming(sfreq),50,sfreq,sfreq);
%      Px1=20*log10(Px);
%     ind=find(Fx(:,1)>=5 & Fx(:,1)<=40);
%     ind2=find(Fx(:,1)>=0 & Fx(:,1)<5);
%     rat=sum(Px1(ind,1))./sum(Px1(ind2,1));
% 

% Feature 2: Segmented Sample entropy of each series 
 N = sfreq; % 1 sec 
 samp=[];
 s_f=[];

    L=length(cop_final_alldata);
   for j=1:L-(N)
        s=sampenc(cop_final_alldata(1,j:(j+N))',3,0.25);
        s=s(1,1);
        s_f=[s_f s];
    end
    samp=s_f;


% Feature 3 : Descriptive features on SE values
f31=[];

       samp(isinf(samp)) = 0;
       data= samp;
%        f31(1,1)=mean(data); 
%        f31(1,2)=std(data);
%        f31(1,3)=kurtosis(data);
%        f31(1,4)=skewness(data);
       f31(1,1)=max(data)-min(data);
       f31(1,2)=median(data);



% Feature 4 : Descriptive features on ECG values
f41=[];
     
       data= cop_final_alldata;
       f41(1,1)=mean(data); 
%        f41(1,2)=std(data);
%        f41(1,3)=kurtosis(data);
%        f41(1,4)=skewness(data);
%        f41(1,5)=max(data)-min(data);
       f41(1,2)=median(data);



% Feature 5 : Descriptive features on RR interval data
f51=[];

       [~,qrs_i_raw,~]=pan_tompkin(allrawdata,sfreq,0);
       datarr1= diff(qrs_i_raw)/sfreq;
%        f51(1,1)=mean(datarr1); 
       f51(1,1)=std(datarr1);
%        f51(1,3)=kurtosis(datarr1);
       f51(1,2)=skewness(datarr1);
%        temp=max(datarr1)-min(datarr1);
%        if isempty(temp)
%            f51(1,5)=0;
%        else
%            f51(1,5)=temp;
%        end
       f51(1,3)=median(datarr1);
       



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


% Feature 7: Descriptive features on 2nd order delta RR interval data
f71=[];

       
       [~,qrs_i_raw,~]=pan_tompkin(allrawdata,sfreq,0);
       data2= diff(qrs_i_raw)/sfreq;
       datarr=diff(diff(data2));
%        f71(1,1)=mean(datarr); 
       f71(1,1)=std(datarr);
       f71(1,2)=kurtosis(datarr);
       f71(1,3)=skewness(datarr);
       temp=max(datarr)-min(datarr);
       if isempty(temp)
           f71(1,4)=0;
       else
           f71(1,4)=temp;
       end
       f71(1,5)=median(datarr);
       


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

f81=[];
       
       [~,qrs_i_raw,~]=pan_tompkin(allrawdata,sfreq,0);
       data= diff(qrs_i_raw)/sfreq;
       if isempty(data)
            f81(1,1:2)=0;
       else    
%            ran=quantile(data,10);
%            ind1=find(data<ran(1,1)); 
%            f81(1,1)=mean(data(:,ind1));
%            ind2=find(data>ran(1,10)); 
%            f81(1,2)=mean(data(:,ind2)); 
           f81(1,1)=min(data);
           f81(1,2)=max(data);
       end



% Feature 9 : Wavelet based 
    
          f91=[];  

          [c,l] = wavedec(allrawdata,5,'sym4');
          [cd3,cd4,cd5] = detcoef(c,l,[3 4 5]);
                  
          f91(1,1)=[var(cd5)];
     
          
% Feature 10 : PSPR based on ECG
           f10=[];  
           % Resample for PSPR features
           cop_final_alldata1=resample(cop_final_alldata,8,200); 
           [f10(1,1) f10(1,2)]=pspr(cop_final_alldata1',Ref_db,matchedpatterninfo_ref);

           
% Feature 12: Based on the KS test  
h=[];
f12=[];

    for j=1:50
        data= datarr1*1000;
        if isempty(data)
            f12(1,1:2)=[0];
            f12(1,3:6)=[0];
            p1=0;
            p2=0;
            p3=0;
            p4=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(1,1:2)=[ mean(p2) mean(p3) ];
        f12(1,3:6)=[median(p1) median(p2) median(p3) median(p4)];
    
%% Feature 13: Coefficient of variation in the Poincare plot
     
             R=[];
            if isempty(datarr1) || length(datarr1)<=2
                 R=0;
            else
                  R12= corrcoef([datarr1(1,1:length(datarr1)-1)',data(1,2:length(datarr1))']);
                  r12=R12(1,2);
                  R=r12;
            end

%% Feature 14: PCA on the Poincare plot  
             
             perc=[];
             if isempty(datarr1) || length(datarr1)<=2
                       perc=0;
             else
                      [COEFF, SCORE, LATENT,TSQUARED, EXPLAINED] = pca([datarr1(1,1:length(datarr1)-1)',datarr1(1,2:length(datarr1))']);
                       perc=(EXPLAINED(1,1));
             end

% %% Feature 15: Percentage of rr interval
%             per=[];
%                 ind=find(datarr1 >=0.6 & datarr1 <=1.2 );
%                 per=length(ind)/length(datarr1)*100;
           
%% Feature 16: Heart rate range
            HR=[];
            hr=60./(datarr1);
            HR=mad(hr);
            

testdata=[real(f31) real(f41) real(f51) real(f61) real(f71) real(f81) real(f91) real(f10) real(f12) real(R) real(perc) real(HR)];
result= predict(trainedClassifier31.ClassificationEnsemble,testdata);
 
if  result==0
    classifyResult= 'N';
elseif result==1
    classifyResult= 'A';
elseif result==2
    classifyResult= 'O';
else result==3
    classifyResult= '~';
    
   
end 

