互近似熵文档格式.docx
《互近似熵文档格式.docx》由会员分享，可在线阅读，更多相关《互近似熵文档格式.docx（17页珍藏版）》请在冰豆网上搜索。
MIX（P）=（1-Z）*X+Z*Y
其中，X为正弦信号序列，Y为均布随机数。
Z为0和1序列，Z取1的概率为P，P越大，则该序列理论上越复杂。
构造2000点MIX（0.1）和MIX（0.9）过程进行试验。
%%MIXstochasticprocesses
%usingtotestthelackofconsistencyofapen
N
=2000;
fs=50;
t
=0:
1/fs:
（N-1）/fs;
x
=sqrt
（2）*sin（2*pi*t/12）;
y
=unifrnd（-sqrt（3）,sqrt（3）,[1,N]）;
p1
=unifrnd（0,1,[1,N]）;
p9
=p1;
p1（p1>
=0.9）=1;
p1（p1~=1）
=0;
p9（p9>
=0.1）=1;
p9（p9~=1）
MIX1=（1-p1）.*x+p1.*y;
MIX9=（1-p9）.*x+p9.*y;
m=2;
r=logspace（log10（0.01）,0,20）;
ApEn1=zeros（size（r））;
ApEn9=zeros（size（r））;
fork=1:
length（r）
ApEn1（k）=apen（MIX1,m,r（k））;
ApEn9（k）=apen（MIX9,m,r（k））;
end
semilogx（r,ApEn1,'
o'
）
holdon
semilogx（r,ApEn9,'
p'
ylim（[04]）
set（gca,'
XTickLabel'
[0.010.11]）;
结果如下图所示：
aa.jpg
上图中“o”序列为MIX（0.1）过程，另外为MIX（0.9）过程。
可以明显看出，这两个序列的近似熵值随阈值r的选择出现交叉现象：
当r取值较小时，获得的结果是过程MIX（0.1）复杂度较高，随着r的取值增大，过程MIX（0.9）的复杂度高于MIX（0.1）。
这种不一致性影响了我们的判断。
原创——MATLAB近似熵计算函数
熵是系统无序的度量，近似熵可以仅仅依靠少量的数据，便可以对一个序列的无序性进行有效评估（PincusSM1991）。
近似熵的计算方法如图所示（徐安2005）：
snap1.jpg
snap2.jpg
snap4.jpg
根据上述步骤，在MATLAB中可以很方便的编程实现，但是问题在于，如果没有一定的编程技巧，上述步骤编写出的程序运行速度将会非常之慢，以下为一组数据：
（1）直接按照上述算法编写，利用嵌套循环形式，计算长度为2000点的序列ApEn，其计算时间为：
259.920109秒；
（2）对循环使用一定的技巧进行优化，尽量减少嵌套循环，并对循环中出现的变量预开辟空间，仍然计算上述长度为2000的序列的近似熵，其计算时间缩短到3.891304秒！
（3）使用matlab中提供的一些特殊函数，再次优化循环，仍然计算上述数据序列的近似熵，其计算时间为：
1.578624秒！
时间的节省效果是非常显著的，如果试验中计算100个数据序列，每个数据序列长度为2000点，那么，使用
（1）的方法，所需时间大约为7小时；
如果采用优化循环后的方法
（2），需要的时间大概为7分钟；
使用再次优化后的（3），使用时间大概不超过2分钟。
因此，在当前计算机硬件迅猛发展的前提下，即使不使用快速算法，这些程序的计算时间仍然可以接受的，只是需要在编程时增加一些优化，多一些技巧。
最后，给出一个试验，计算均匀独立同分布时间序列的近似熵，每一个序列长度为5000（使用最优化版程序计算速度大概为8秒一个序列），绘制近似熵结果与阈值r的选取之间的关系，r值取从0.01到1直接按照对数间隔取20个点，因此共需要计算20次，总共时间大约3分钟左右。
试验程序如下：
以下内容为程序代码:
rnddat=unifrnd（0,1,5000,1）;
rnddat=rnddat./std（rnddat）;
%标准差归一化
ApEn=zeros（size（r））;
ApEn（k）=apen（rnddat,m,r（k））;
semilogx（r,ApEn,'
ylim（[06]）
实验结果：
snap5.jpg
注意横轴绘制的对数坐标。
关于程序优化的技巧，请参考：
《MATLABN个实用技巧》
关于本程序下载以及问题讨论，请在21世纪电子论坛matlab讨论组发帖进行：
%求近似熵函数
function
resultapen=APEN（input_signal,input_m,input_r）;
%input_signal为要求近似熵的序列，input_m为给定的模式的维数，input_r为相似容限
x=input_signal;
m=input_m;
r=input_r;
[input_row,input_col]=size（x）;
%输入的信号序列一般为1行n列
n=input_col;
%信号的点数
rownum=n-m+1;
%构造的m维矢量的行数
signalmatrix=zeros（rownum,m）;
%构造信号矩阵
fori=1:
1:
rownum
forj=1:
m
signalmatrix（i,j）=x（i+j-1）;
end
sumC=0;
forii=1:
nnum=0;
forjj=1:
U=zeros（m,1）;
%存放
d=0;
forkk=1:
U（kk,1）=abs（signalmatrix（ii,kk）-signalmatrix（jj,kk））;
d=max（U）;
%X（i）与X（j）的距离的最大值
ifd<
r
nnum=nnum+1;
%如果距离小于阈值，则个数+1；
C=（nnum-1）/（rownum-1）;
ifC>
sumC=sumC+log（C）;
resultapen=sumC/rownum;
return;
样本熵近似熵代码应用于matlab
matlab,样本,代码,应用
Ñ
ù
±
¾
ì
Ø
½
ü
Ë
Æ
´
ú
Â
ë
Ó
¦
Ó
Ã
Ú
matlab
matlab,Ñ
´
Ó
Ö
Ê
Ç
´
ó
¼
Ò
¿
ª
µ
×
Å
Ï
Û
Î
Ä
º
ò
Á
~~
Õ
â
¸
ö
Ð
È
Ý
è
·
ß
Ô
§
Ü
¿
É
²
function[shang]=jss（xdate）
m=2;
n=length（xdate）;
r=0.2*std（xdate）;
cr=[];
gn=1;
gnmax=m;
whilegn<
=gnmax
x2m=zeros（n-m+1,m）;
%´
æ
ä
»
d=zeros（n-m+1,n-m）;
%´
à
À
½
á
¹
û
cr1=zeros（1,n-m+1）;
k=1;
n-m+1
forj=1:
x2m（i,j）=xdate（i+j-1）;
x2m;
ifi~=j
d（i,k）=max（abs（x2m（i,-x2m（j,））;
%¼
ã
÷
Í
k=k+1;
d;
[k,l]=size（find（d（i,<
r））;
%½
«
RÐ
¡
ý
ø
L
cr1（1,i）=l;
cr1;
cr1=（1/（n-m））*cr1;
sum1=0;
ifcr1（i）~=0
sum1=sum1+log（cr1（i））;
cr1=1/（n-m+1）*sum1;
cr（1,gn）=cr1;
gn=gn+1;
m=m+1;
cr;
shang=cr（1,1）-cr（1,2）;
function[shang]=ybs（xdate）
sum1=sum1+cr1（i）;
shang=-log（cr（1,1）/cr（1,2））;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%存放变换后的向量
%存放距离结果的矩阵
d（i,k）=max（abs（x2m（i,
-x2m（j,
））;
%计算各个元素和响应元素的距离
[k,l]=size（find（d（i,
<
%将比R小的个数传送给L
基于脑电信号的认知动力学系统研究——线性/非线性方法及动态时—频—空分析
升级会员 