样本熵的matlab程序

转自http://www.ilovematlab.cn/thread-89427-1-1.html MATLAB下的动态样本熵计算

1. function SampEnVal = SampEn(data, m, r)
   
2. %SAMPEN  计算时间序列data的样本熵
   
3. %        data为输入数据序列
   
4. %        m为初始分段，每段的数据长度
   
5. %        r为阈值
   
6. % $Author: lskyp
   
7. % $Date:   2010.6.20
   
8. % Orig Version: V1.0--------分开计算长度为m的序列和长度为m+1的序列
   
9. %                           这一版的计算有些问题，需要注意两个序列总数都要为N-m
   
10. % Modi Version: V1.1--------综合计算，计算距离时通过矩阵减法完成，避免重循环
    
11. % V1.1 Modified date: 2010.6.23
    
12. data = data(:)';
    
13. N = length(data);
    
14. Nkx1 = 0;
    
15. Nkx2 = 0;
    
16. % 分段计算距离，x1为长度为m的序列，x2为长度为m+1的序列
    
17. for k = N - m:-1:1
    
18.     x1(k, :) = data(k:k + m - 1);
    
19.     x2(k, :) = data(k:k + m);
    
20. end
    
21. for k = N - m:-1:1
    
22.     % x1序列计算
    
23.     % 统计距离，由于每行都要与其他行做减法，因此可以先将该行复制为N-m的矩阵，然后
    
24.     % 与原始x1矩阵做减法，可以避免两重循环，增加效率
    
25.     x1temprow = x1(k, :);
    
26.     x1temp    = ones(N - m, 1)*x1temprow;
    
27.     % 可以使用repmat函数完成上面的语句，即
    
28.     % x1temp = repmat(x1temprow, N - m, 1);
    
29.     % 但是效率不如上面的矩阵乘法
    
30.     % 计算距离，每一行元素相减的最大值为距离
    
31.     dx1(k, :) = max(abs(x1temp - x1), [], 2)';
    
32.     % 模板匹配数
    
33.     Nkx1 = Nkx1 + (sum(dx1(k, :) < r) - 1)/(N - m - 1);
    
34.     
    
35.     % x2序列计算，和x1同样方法
    
36.     x2temprow = x2(k, :);
    
37.     x2temp    = ones(N - m, 1)*x2temprow;
    
38.     dx2(k, :) = max(abs(x2temp - x2), [], 2)';
    
39.     Nkx2      = Nkx2 + (sum(dx2(k, :) < r) - 1)/(N - m - 1);
    
40. end
    
41. % 平均值
    
42. Bmx1 = Nkx1/(N - m);
    
43. Bmx2 = Nkx2/(N - m);
    
44. % 样本熵
    
45. SampEnVal = -log(Bmx2/Bmx1);