求一个近似熵快速算法的matlab程序

快速近似熵算法：

function [ApEn_value,Cmr,Cmr_1] = fast_ApEn(signal,r_factor)

% Estimate the Aproximate Entropy (ApEn) of a signal, using a fast % algorithm, for the ApEn parameter "m" equal to 2

% The pattern length "m" for which this routine was implemented is 2. For another values of "m", the instructions marked with a (*) in the end must be changed.% m=1 or m=2

% r between 0.1*STD and 0.25*STD, where STD is the signal standard deviation % N (signal length) between 75 and 5000

% [ApEn_value] = fast_ApEn(signal,r_factor)% Input variables:% signal - signal% r_factor - factor of the criterion of similarity r_factor*std(signal)% Output variables:

% ApEn_value - ApEn calculated from the signal

% Optional output variables:% C_m% C_m_1

% if length(signal)<75 | length(signal)>5000

% slength=input('Signal length inappropriate. Continue anyway? (y/n)','s')% if strcmpi(slength,'y')==0% return% end% end

% if r_factor<0.1 | r_factor>0.25

% r_factor_flag=input('Value for r parameter is inappropriate. Continue anyway? (y/n)','s')

% if strcmpi(r_factor_flag,'y')==0% return

% end% end% Initial variables definition.

signal=signal(:)

N=length(signal)

Cmr_ij=[]Cmr_i=[]Cmr=[]

Cmr_ij_1=[]Cmr_i_1=[]Cmr_1=[]

% D and S matrixes computation.

D=abs(signal*ones(1,N)-ones(N,1)*signal')

S=zeros(N,N)S(find(D<=r_factor*std(signal)))=1

% C's computation for "m" and "m+1" patterns.

m=2S(N+1,(m+1):N)=0% necessary for the loop to be possible ("artificial" definitions)

for k=1:N-(m-1)% m pattern.

Cmr_ij=S(k,1:N-1).*S(k+1,2:N)% (*)

Nm_i=sum(Cmr_ij)

Cmr_i=Nm_i/(N-(m-1))

Cmr=[CmrCmr_i]% m+1 pattern.

Cmr_ij(end)=[]

Cmr_ij_1=Cmr_ij.*S(k+2,3:N)

Nm_i_1=sum(Cmr_ij_1)

Cmr_i_1=Nm_i_1/(N-m)

Cmr_1=[Cmr_1Cmr_i_1]

end

Cmr_1(end)=[]% the last C value for the "m+1" pattern is artificial.% Phi’s computation.

phi_m=mean(log(Cmr))

phi_m_1=mean(log(Cmr_1))% ApEn final calculation.

ApEn_value=[phi_m-phi_m_1]% phi_m% phi_m_1

%求近似熵函数function resultapen=APEN(input_signal,input_m,input_r)%input_signal为要求近似熵的序列，input_m为给定的模式的维数，input_r为相似容限x=input_signalm=input_mr=input_r[input_row,input_col]=size(x)%输入的信号序列一般为1行n列n=input_col%信号的点数

rownum=n-m+1%构造的m维矢量的行数signalmatrix=zeros(rownum,m)%构造信号矩阵for i=1:1:rownumfor j=1:1:msignalmatrix(i,j)=x(i+j-1) endendsumC=0

for ii=1:1:rownumnnum=0 for jj=1:1:rownumU=zeros(m,1)%存放d=0 for kk=1:1:m U(kk,1)=abs(signalmatrix(ii,kk)-signalmatrix(jj,kk)) endd=max(U)%X(i)与X(j)的距离的最大值if d<r nnum=nnum+1%如果距离小于阈值，则个数+1； endendC=(nnum-1)/(rownum-1) if C>0sumC=sumC+log(C) endendresultapen=sumC/rownumreturn

function shang(x)

[n,m]=size(x)

k=1/log(n)

X=zeros(n,m)

for j=1:m

for i=1:n

c=sort(x(:,j))

big=x(n,j)

small=x(1,j)

X(i,j)=(x(i,j)-small)/(big-small)+1

end

---