function [d,C,stat_RR,pow,AFEv,OriginCount,IrrEv,PACEv,p_ratio,HOS,P,P_absence,f_presence,Noise,wevelet_power]=featureExtract(recordName)
% function [stat_RR]=featureExtract(recordName)
% 
% 
% Written by: Na Liu April 7 2017
%             nkuliuna@163.com
% 
%% 去噪，提取各层小波能量特征
[tm,ecg,fs,siginfo]=rdmat(recordName);
[ecg,wevelet_power] = fliter(ecg,'db5');
%% 时域统计特征
pow = sum(ecg.^2)/(length(ecg)/fs);
ecg2=ecg; ecg2(abs(ecg2)>2)=0;
QRS=qrs_detect2(ecg2,0.25,0.6,fs); 
mean_R=mean(ecg(QRS));
if length(QRS)<7|(QRS(end)-QRS(1))<length(ecg)/2|abs(mean_R)>2
   QRS=findQRS(ecg,1);
end
%     close all
%     plot(ecg);hold on;plot(QRS,ecg(QRS),'ro')

mean_R=mean(ecg(QRS));
std_R=std(ecg(QRS));
if mean_R>0,else ecg=-ecg;end

RR=diff(QRS')/fs;
if isempty(RR),min_RR=0;else min_RR=min(RR);end
mean_RR=mean(RR);std_RR=std(RR);median_RR=median(RR);
CV=std_RR/mean_RR; delta_RR_max = max(RR)-min(RR);
alpha=0.75; beta=1.4; gamma=0.3;
RR2=RR;RR2(RR2<alpha*mean_RR&RR2>beta*mean_RR)=[];
NADev=mean(sum(abs(RR2-mean_RR))/mean_RR);
diff_RR=abs(diff(RR)); 
rms_diff_RR=rms(diff_RR); std_diff_RR=std(diff_RR); pNN50=sum(diff_RR>0.05)/length(diff_RR);
diff_RR(diff_RR>gamma)=[];
NADiff=mean(sum(diff_RR)/mean_RR);
stat_RR=[mean_R,std_R,mean_RR,std_RR,min_RR,median_RR,CV,delta_RR_max,NADev,NADiff,rms_diff_RR,std_diff_RR,pNN50];

%% 谱估计，频域特征：F(f-presence)，Noise
[pxx,f2] = pwelch(ecg,1*fs,[],[],fs);
f=0:0.5:150;
px = interp1(f2,pxx,f,'pchip');
%     plot(f2,pxx,'o-',f,px,'ro-');
px1 = px(f>=3&f<=10); px2 = px(f>=10&f<=125);
f_presence = trapz(f(f>=3&f<=10),px1)/trapz(f,px);
Noise = rms(ecg)*trapz(f(f>=3&f<=10),px1.*log2(px1))/trapz(f(f>=10&f<=125),px2.*log2(px2));

%% 样本熵，AFEv，频域特征P_absence
if length(RR)<6 
    d = 3;
    C = -10;
    AFEv = -99; OriginCount = 0; IrrEv = 0; PACEv = 0;
    p_ratio = 0;
    P_absence = 100;
else
    r0 = std_RR * 0.2;
    [d,r] = sampleEntropy_RR(RR, 1, r0, 1);
    C = d+log(2*r)-log(mean_RR);
%         M = ceil(length(RR)/12);
%         r0=0.03;
%         for i=1:M-1
% %             r0=std(RR(i*12-11:i*12))*0.2;
%             [d(i),r] = sampleEntropy_RR(RR(i*12-11:i*12), 1, r0, 1);
%             C(i) = d(i)+log(2*r)-log(mean(RR(i*12-11:i*12)));
%         end 
% %         r0=std(RR(max(1,end-11):end))*0.2;
%         [d(M),r] = sampleEntropy_RR(RR(max(1,end-11):end), 1, r0, 1);
%          C(M) = d(M)+log(2*r)-log(mean(RR(max(1,end-11):end)));
%          d=mean(d);C=mean(C);
%         C = d+log(2*r)-log(mean_RR)+log(length(RR));
    [AFEv,OriginCount,IrrEv,PACEv] = comput_AFEv(RR);
% 计算P-wave absence
    left1=QRS-0.2*fs;
    num1 = sum(left1>0);
    e = zeros(num1,num1);
    for i = 1:num1-1
       if left1(i)<=0
       else 
           for j = i+1:num1
               e(i,j)=sum((ecg(left1(i):QRS(i))-ecg(left1(j):QRS(j))).^2);
           end
           p(i)=sum(e(i,:))/(num1-i);
       end
    end
    P_absence=sum(p);
end
%%  P波存在率p_ratio，高阶统计量：方差、峰度、偏度，P波6段均值
plocs_1=[];
left = floor(QRS(2:end)-0.205*mean(RR)*fs);
right= ceil(QRS(2:end)-0.06*mean(RR)*fs);
num=sum(left>0);
v=zeros(num,1);%方差
ske=zeros(num,1);%偏度
kur=zeros(num,1);%峰度
p1=zeros(num,1);p2=zeros(num,1);p3=zeros(num,1);p4=zeros(num,1);p5=zeros(num,1);p6=zeros(num,1);

for i=1:length(left) 
    if left(i)<=0
    else
        [pks,locs]=findpeaks(ecg(left(i):right(i)),'MinPeakHeight',0.25*ecg(QRS(i+1)));
        [~,ind]=max(ecg(locs-1+left(i)));
        ploc=locs(ind)-1+left(i);plocs_1=[plocs_1,ploc];
%         9个特征：方差，峰度，偏度和6段均值
        L=length(ecg(left(i):right(i)));
        if L==36,Pwave=ecg(left(i):right(i));
        else Pwave=resample(ecg(left(i):right(i)),36,L);
        end
        Pwave=Pwave-mean(Pwave);%减去均值
        v(i)=var(Pwave);ske(i)=skewness(Pwave);kur(i)=kurtosis(Pwave);
        p1(i)=mean(Pwave(1:6));p2(i)=mean(Pwave(7:12));p3(i)=mean(Pwave(13:18));
        p4(i)=mean(Pwave(19:24));p5(i)=mean(Pwave(25:30));p6(i)=mean(Pwave(31:36));
    end  
end
p_ratio=length(plocs_1)/num;
V=mean(v);SKE=mean(ske);KUR=mean(kur);
HOS=[V,SKE,KUR];
P1=mean(p1);P2=mean(p2);P3=mean(p3);P4=mean(p4);P5=mean(p5);P6=mean(p6);
P=[P1,P2,P3,P4,P5,P6];
end