样本熵

MATLAB下的动态样本熵计算

样本熵是时间序列复杂度的一种度量，在上个世纪末期由几位非线性动力学研究者提出。样本熵比近似熵更具有相对一致性，在分析生物信号序列的复杂度分析中已经获得成功应用。

文献中给出的样本熵算法步骤如下：

根据文献中提供的算法步骤，我们可以通过C语言或者MATLAB 编程实现。

这里有一个实现方法，这是在ying_327的提问帖（程序样本熵求助）中出现的，现摘抄如下：

复制内容到剪贴板

代码:

function [shang]=ybs(xdate)

% clc,clear,

% xdate = randn(10,1);

m=2;

n=length(xdate);

r=0.2*std(xdate);

cr=[];

gn=1;

gnmax=m;

% while gn<=gnmax

x2m=zeros(n-m+1,m);%存放变换后的向量

d=zeros(n-m+1,n-m);% 存放距离结果的矩阵

cr1=zeros(1,n-m+1);%存放

k=1;

for i=1:n-m+1

for j=1:m

x2m(i,j)=xdate(i+j-1);

end

end

x2m;

for i=1:n-m+1

for j=1:n-m+1

if i~=j

d(i,k)=max(abs(x2m(i,:)-x2m(j,:)));%计算各个元素和响应元素的距离

k=k+1;

end

end

k=1;

end

d;

for i=1:n-m+1

[k,l]=size(find(d(i,:)

cr1(1,i)=l;

end

cr1;

cr1=(1/(n-m))*cr1;

sum1=0;

for i=1:n-m+1

sum1=sum1+cr1(i);

end

end

cr1=1/(n-m+1)*sum1;

cr(1,gn)=cr1;

gn=gn+1;

m=m+1;

end

cr;

shang=-log(cr(1,1)/cr(1,2));

我们先不看程序完成的是否正确，仅从程序结构本身就可以看出，这个程序有着浓厚的C风格，在MATLAB强大的矩阵处理能力下，循环

操作显得异常刺眼，而且效率异常低。我在回复中给出另外一种解法：复制内容到剪贴板

代码:

function SampEnVal = SampEn(data, m, r)

%SAMPEN 计算时间序列data的样本熵

% data为输入数据序列

% m为初始分段，每段的数据长度

% r为阈值

% $Author: lskyp

% $Date: 2010.6.20

data = data(:)';

N = length(data);

% 分段计算距离

for k = N - m:-1:1

xm(k, :) = data(k:k + m - 1);

end

for k = N - m:-1:1

for p = N - m:-1:1

% 统计距离，k=p时距离为0，计算时舍去一个即可

d(k, p) = max(abs(xm(p, :) - xm(k, :)));

end

end

% 模板匹配数以及平均值计算

Nk = 0;

for k = N - m:-1:1

Nk = Nk + (sum(d(k, :) < r) - 1)/(N - m - 1);

end

Bm = Nk/(N - m);

clear xm d Nk

% m值增加1，重复上面的计算

m = m + 1;

% 分段计算距离

for k = N - m:-1:1

xm(k, :) = data(k:k + m - 1);

end

for k = N - m:-1:1

for p = N - m:-1:1

% 统计距离，k=p时距离为0，计算时舍去一个即可

d(k, p) = max(abs(xm(p, :) - xm(k, :)));

end

end

% 模板匹配数以及平均值计算

Nk = 0;

for k = N - m:-1:1

Nk = Nk + (sum(d(k, :) < r) - 1)/(N - m - 1);

end

Bmadd1 = Nk/(N - m);

% 熵值

SampEnVal = -log(Bmadd1/Bm);

刚又分析了这个程序，对比算法来说，这个程序有些问题，注意到两个序列（长度为m的以及长度为m+1）的数量应该都为N-m个，而上面的程序由于中间m进行了加1操作，因此第二次计算的序列数量比第一次少了一个。两个序列应该同时计算。

另外，上面的程序使用了两重循环，经过分析，可以将内层循环矩阵化，矩阵化之后的程序仅进行一次循环，效率增加了很多，但是程序看起来可能不如第一版程序好理解。

复制内容到剪贴板

代码: