样本熵(SampEn)的C/C++代码实现与优化

已于 2023-03-07 13:40:41 修改
阅读量352 收藏 4
点赞数
分类专栏： 杂项 文章标签： c语言 c++
于 2023-01-30 00:10:07 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/PeaZomboss/article/details/128795851
版权

正文

本文不介绍什么是样本熵，具体推荐看此文https://blog.csdn.net/Cratial/article/details/79742363，写的很好，里面的示例也被我拿来测试代码写的对不对。

本文所有代码可以在此处找到https://gitcode.net/PeaZomboss/miscellaneous，文件夹231424-sampen。

前言

一开始有人找我帮忙写些C++代码，实现一些算法，这其中呢就有一个是样本熵。

当时粗写了一个，后来似乎性能不太够，当然实际上性能是没问题的，虽然确实不如现在优化的，但也不会很差，原因是代码被改了一部分，当然这个和算法的优化没关系，就不展开了。

不过既然都要优化了，那不如在算法层面也做一些改进，让速度再快一点，于是就想了想办法提升了性能。于是在此之后就整理了一下代码，并严谨测试，然后与大家分享。

开头提到的那篇文章里有用python写的代码，还有一个matlab写的。但是不管是python还是matlab都不是拿来实际应用的，一般追求性能的模块会用C/C++，当然用Java这类的应该也不会太差。

代码实现

首先是一般情况下很容易想到的代码：

static double step(double *X, int N, int m, double r)
{
    double sum = 0;
    for (int i = 0; i <= N - m; i++) {
        int Bi = 0;
        for (int j = 0; j <= N - m; j++) {
            if (i != j) {
                /* 找出最大差值 */
                double D = fabs(X[i] - X[j]);
                for (int k = 1; k < m; k++) {
                    double t = fabs(X[i + k] - X[j + k]);
                    if (D < t)
                        D = t;
                }
                if (D <= r)
                    Bi++;
            }
        }
        sum += 1.0 * Bi / (N - m); // 会有累加误差
    }
    return sum / (N - m + 1);
}

double SampEn(double *X, int N, int m, double r)
{
    double B = step(X, N, m, r);
    if (B == 0) // 尽管大部分时候不会是0
        return 0;
    double A = step(X, N, m + 1, r);
    if (A == 0)
        return 0;
    return -log(A / B);
}

这里的所有步骤基本都是按照原算法描述的内容进行编写的，所以非常好理解啊。而且我已经故意避开了诸如fmax之类的函数，优化找最大差值的逻辑，都是为了提高性能。

不过这个step函数的sum += 1.0 * Bi / (N - m);在N较大的时候累加误差会增加，而实际上这个Bi变量完全可以累加到最后在进行运算，不但减少误差还能降低运算量，所以可以改成：

static double step(double *X, int N, int m, double r)
{
    int Bi = 0;
    for (int i = 0; i <= N - m; i++) {
        for (int j = 0; j <= N - m; j++) {
            if (i != j) {
                double D = fabs(X[i] - X[j]);
                for (int k = 1; k < m; k++) {
                    double t = fabs(X[i + k] - X[j + k]);
                    if (D < t)
                        D = t;
                }
                if (D <= r)
                    Bi++;
            }
        }
    }
    return 1.0 * Bi / (N - m) / (N - m + 1);
}

这样子几乎就是按照原来的算法思路进行了一些简单处理，但提升并不明显。

代码优化

仔细看这个step()函数，那复杂度可是O(n^2)啊，而且还要跑两次，那损耗肯定惊人。但是仔细看这两次计算，只有第一次的m变成了m+1，那能不能把两次运算整合到一起呢？

注意到m+1和m的区别在于循环次数少了一次，寻找最大差值的向量多了一个，其余是没有区别的，那我们完全可以在处理m的时候顺便把m+1的情况也放一起就行了。

请看：

double FastSampEn(double *X, int N, int m, double r)
{
    int Ai = 0, Bi = 0;
    int LoopsSub1 = N - m; // 循环次数减一，因为算法中的表述是从1到N-m+1，而我们是从0开始的
    for (int i = 0; i <= LoopsSub1; i++) {
        for (int j = 0; j <= LoopsSub1; j++) {
            if (i != j) {
                // 这里一样找m个的最大差值
                double D = fabs(X[i] - X[j]);
                for (int k = 1; k < m; k++) {
                    double t = fabs(X[i + k] - X[j + k]);
                    if (D < t)
                        D = t;
                }
                if (D <= r)
                    Bi++;
                // 对于m+1的情况，当到达m维数的边界的时候显然是不行的
                // 所以我们只要限制边界情况就行了
                if (i != LoopsSub1 && j != LoopsSub1) {
                    double t = fabs(X[i + m] - X[j + m]);
                    if (D < t) // 判断最后一个是不是最大值
                        D = t;
                    if (D <= r)
                        Ai++;
                }
            } // i!=j
        } // j
    } // i
    double B = 1.0 * Bi / (N - m) / (N - m + 1);
    double A = 1.0 * Ai / (N - m - 1) / (N - m);
    if (B == 0 || A == 0)
        return 0;
    return -log(A / B);
}

这样修改以后速度可以达到原先的1.5倍，算上循环内部的开销，这样的情况也是挺不错的。

至此，针对样本熵代码的基本优化我能想到的就到这了，不过后面还有进一步的优化。

考虑到实际应用情况，这个m大多数时候都是取值2，所以针对m为2的情况又专门优化了一下：

double FastSampEn_m2(double *X, int N, double r)
{
    int Ai = 0, Bi = 0;
    int LoopsSub1 = N - 2;
    for (int i = 0; i <= LoopsSub1; i++) {
        for (int j = 0; j <= LoopsSub1; j++) {
            if (i != j) {
                double D = fabs(X[i] - X[j]);
                double t = fabs(X[i + 1] - X[j + 1]);
                if (D < t)
                    D = t;
                if (D <= r)
                    Bi++;
                if (i != LoopsSub1 && j != LoopsSub1) {
                    double t = fabs(X[i + 2] - X[j + 2]);
                    if (D < t)
                        D = t;
                    if (D <= r)
                        Ai++;
                }
            }
        }
    }
    double B = 1.0 * Bi / (N - 2) / (N - 1);
    double A = 1.0 * Ai / (N - 3) / (N - 2);
    if (B == 0 || A == 0)
        return 0;
    return -log(A / B);
}

在m取2的情况下可以稍微快一点点。

不过这个代码还是有比较大的提升空间，我找了ChatGPT给我提供优化，但是结果一直出错，后来我分析了一下代码，手动调整了一下，得出了正确的代码，然后再让AI分析。经过一段时间的操作，现在得出了一个效率更高的代码，而且没有任何额外的内存占用。

double FastSampEn_ai(double *X, int N, int m, double r)
{
    int Ai = 0, Bi = 0;
    int LoopTimes = N - m + 1;
    for (int i = 0; i < LoopTimes; i++) {
        for (int j = i + 1; j < LoopTimes; j++) {
            double D = fabs(X[i] - X[j]);
            for (int k = 1; k < m; k++) {
                double t = fabs(X[i + k] - X[j + k]);
                if (t > D)
                    D = t;
            }
            if (D <= r)
                Bi++;
            if (j < N - m) {
                double t = fabs(X[i + m] - X[j + m]);
                if (t > D)
                    D = t;
                if (D <= r)
                    Ai++;
            }
        }
    }
    double B = 1.0 * Bi / ((N - m + 1) * (N - m) / 2);
    double A = 1.0 * Ai / ((N - m) * (N - m - 1) / 2);
    if (A == 0 || B == 0)
        return 0;
    return -log(A / B);
}

简单分析下思路，因为两层循环里面，实际上有一半的计算是重复的，所以直接使用for (int j=i+1;...)可以减少一半的计算量，最后只要再让结果乘个2就好了。其实AI提出了好几种办法，比如使用空间换时间的方法等，但是存在不少bug，最后我采纳了这个比较简单明了且效果好的代码，并进行了微调测试，最终使速度为原来的两倍。

然后再写一个手动替换m为2的情况：

double FastSampEn_ai_m2(double* X, int N, double r)
{
    int Ai = 0, Bi = 0;
    int LoopTimes = N - 1;
    for (int i = 0; i < LoopTimes; i++) {
        for (int j = i + 1; j < LoopTimes; j++) {
            double D = fabs(X[i] - X[j]);
            double t = fabs(X[i+1]-X[j+1]);
            if (t > D)
                D = t;
            if (D <= r)
                Bi++;
            if (j < N - 2) {
                double t = fabs(X[i + 2] - X[j + 2]);
                if (t > D)
                    D = t;
                if (D <= r)
                    Ai++;
            }
        }
    }
    double B = 1.0 * Bi / ((N - 1) * (N - 2) / 2);
    double A = 1.0 * Ai / ((N - 2) * (N - 3) / 2);
    if (A == 0 || B == 0)
        return 0;
    return -log(A / B);
}

这应该是开销最小的速度基本达到极限的单线程方法了。

测试代码

接下来为了测试上述代码是否正确，就用开头提到的那篇文章的数据进行测试看看：

std::vector<double> x;
for (int i = 0; i < 17; i++) {
    x.push_back(85);
    x.push_back(80);
    x.push_back(89);
}
std::cout << SampEn(x.data(), x.size(), 2, 3) << '\n';
std::cout << FastSampEn(x.data(), x.size(), 2, 3) << '\n';
std::cout << FastSampEn_m2(x.data(), x.size(), 3) << '\n';
std::cout << FastSampEn_ai(x.data(), x.size(), 2, 3) << '\n';
std::cout << FastSampEn_ai_m2(x.data(), x.size(), 3) << '\n';

原文给出的结果是0.0008507018803128114，不过由于std::cout的舍入问题，实际输出的是0.000850702，是符合结果的。

然后再测试一下性能：

for (int i = 0; i < 10000; i++) {
    x.push_back(85);
    x.push_back(80);
    x.push_back(89);
}

double se;
clock_t t;

t = clock();
se = SampEn(x.data(), x.size(), 2, 3);
t = clock() - t;
std::cout << se << ", time = " << t << " ms\n";

t = clock();
se = FastSampEn(x.data(), x.size(), 2, 3);
t = clock() - t;
std::cout << se << ", time = " << t << " ms\n";

t = clock();
se = FastSampEn_m2(x.data(), x.size(), 3);
t = clock() - t;
std::cout << se << ", time = " << t << " ms\n";

t = clock();
se = FastSampEn_ai(x.data(), x.size(), 2, 3);
t = clock() - t;
std::cout << se << ", time = " << t << " ms\n";

t = clock();
se = FastSampEn_ai_m2(x.data(), x.size(), 3);
t = clock() - t;
std::cout << se << ", time = " << t << " ms\n";

先是塞了30000个数据，然后用clock()函数计时，依次测试每一个函数，然后看看结果是否一致，时间差距有多少。

如果有兴趣的话，可以去用之前提到那篇文章的python代码跑跑看，30000个数据花了我半个多小时才跑完，而我的这个代码即使没经过优化也不超过一分钟。而matlab的代码我没试过，电脑上没这个软件，也不是学这个的。

完整代码在开头给出的仓库里，打开test.cpp就可以看到。

运行测试

编译器和版本：g++ (x86_64-win32-seh-rev0, Built by MinGW-W64 project) 8.5.0

测试机器CPU：AMD Ryzen 5 4600H

测试操作系统：Windows 11 21H2

  • 编译命令g++ -D ALL_IN_ONE test.cpp -o test_x64_O0 -O0
0.000850702
0.000850702
0.000850702
0.000850702
0.000850702
2.21505e-09, time = 11192 ms
2.21505e-09, time = 6663 ms
2.21505e-09, time = 5560 ms
2.21505e-09, time = 3150 ms
2.21505e-09, time = 2566 ms
  • 编译命令g++ -D ALL_IN_ONE test.cpp -o test_x64_O1 -O1
0.000850702
0.000850702
0.000850702
0.000850702
0.000850702
2.21505e-09, time = 2483 ms
2.21505e-09, time = 1619 ms
2.21505e-09, time = 1143 ms
2.21505e-09, time = 637 ms
2.21505e-09, time = 515 ms
  • 编译命令g++ -D ALL_IN_ONE test.cpp -o test_x64_O2 -O2
0.000850702
0.000850702
0.000850702
0.000850702
0.000850702
2.21505e-09, time = 2703 ms
2.21505e-09, time = 1622 ms
2.21505e-09, time = 1159 ms
2.21505e-09, time = 692 ms
2.21505e-09, time = 523 ms
  • 编译命令g++ -D ALL_IN_ONE test.cpp -o test_x64_O3 -O3
0.000850702
0.000850702
0.000850702
0.000850702
0.000850702
2.21505e-09, time = 2641 ms
2.21505e-09, time = 1614 ms
2.21505e-09, time = 1160 ms
2.21505e-09, time = 727 ms
2.21505e-09, time = 359 ms
  • 编译命令g++ -D ALL_IN_ONE test.cpp -o test_x64_Os -Os
0.000850702
0.000850702
0.000850702
0.000850702
0.000850702
2.21505e-09, time = 2969 ms
2.21505e-09, time = 1639 ms
2.21505e-09, time = 1156 ms
2.21505e-09, time = 803 ms
2.21505e-09, time = 591 ms

32位g++未作测试，但是理论上速度肯定是会更快的。

  • 使用VC++2022在x64的Release模式下编译
0.000850702
0.000850702
0.000850702
0.000850702
0.000850702
2.21505e-09, time = 2950 ms
2.21505e-09, time = 1854 ms
2.21505e-09, time = 1071 ms
2.21505e-09, time = 852 ms
2.21505e-09, time = 506 ms
  • 使用VC++2022在x86的Release模式下编译
0.000850702
0.000850702
0.000850702
0.000850702
0.000850702
2.21505e-09, time = 3882 ms
2.21505e-09, time = 2343 ms
2.21505e-09, time = 1177 ms
2.21505e-09, time = 1026 ms
2.21505e-09, time = 505 ms

可以看到提升还是显著的，尤其是AI优化的部分，比我手动优化的还快一倍。

更新记录

  • 2023-02-04：修正错别字，优化部分表述。
  • 2023-03-07：新增AI的优化方案以及测试结果，删去不必要的内容。
Peacoor Zomboss
关注
  • 0
    点赞
  • 踩
  • 4
    收藏
    觉得还不错? 一键收藏
  • 2
    评论
专栏目录
c语言中计算熵的函数,如何用c++编程实现各种信源熵的计算
weixin_31231693的博客
05-20 1306
算术编码图象数据压缩标准(如jpeg,jbig)扮演重角色.算术编码,消息01之间实数进行编码.算术编码两基本参数:符号概率它编码间隔.信源符号概率决定压缩编码效率,也决定编码过程信源符号间隔,而这些间隔包含01之间.编码过程间隔决定符号压缩后输出.算术编码需输入符号,各符号概率还需编码符号序列,根据概率可算出初始编码间隔,先设几变量后面可:High——当前编码限,Low——当前编码限,high—...
各种熵的代码（包括模糊熵、样本熵、相似熵等等）
05-04
本代码真实有效，可以运行，包括各种熵的matlab代码，样本熵、近似熵、模糊熵等等。
脑电样本熵
05-11
利用matlab计算样本熵 数据选自矩阵形式
matlab直方图隐藏信息代码-hctsa-1:高度比较的时间序列分析
05-27
matlab直方图隐藏信息代码hctsa ，一个高度比较的时间序列分析代码存储库 hctsa是一个用于运行高度比较时间序列分析的软件包，使用（完全支持R2014b或更高版本；用于python cf.）。 该软件提供了一个代码框架，可从时间序列（或时间序列数据集）中提取成千上万的时间序列分析功能，以及用于对数据进行规范化和聚类，生成数据的低维表示的工具，识别不同时间序列类别之间的区分特征，使用大量时间序列特征集学习多元分类模型，找到与感兴趣时间序列最接近的匹配项，以及一系列其他可视化和分析功能。 所有这些类型的分析都在我们的随附资料中进行了描述。 提供了有关hctsa的全面文档，可以在线阅读或以pdf，epub或mobi格式下载。 任何反馈是非常有帮助的（），并且具体地，涉及该代码的任何改进将是非常在或形式理解。 下载资料库 对于不熟悉git的用户，只需单击“下载.zip”按钮即可下载存储库的当前版本。 建议将存储库与git一起使用。 为此，请使用以下代码将其克隆到本地计算机，然后将其设置为与主存储库保持同步： git remote add upstream git://github.c
SampleEntropy.rar_matlab_sampleEntropy_信号 熵_样本熵_样本熵的参数
07-15
matlab计算信号样本熵，参数.m文件中有说明
matlab信息熵代码-ML-in-Diagnose-:ML诊断
05-19
matlab信息熵代码 校创-利用机器学习实现“烟雾病”诊断 说明 该仓库用于 兰州大学-物理科学与技术学院-ML校创项目 使用 成员：钮博恒、魏弘量、吴若岑、王悦晨、张峰瑞 本次校创旨在 从具体的项目中，相互学习。大家应该都是第一次参与校创/使用GitHub合作，如有任何建议或问题，均可以在仓库专门区域，或者是qq群中提出:beaming_face_with_smiling_eyes: 4月2日更新 我跟俞老师联系了一下，他表示确实需要计算多种熵，所以目前主要在于搜集熵计算方法，不仅要了解算法，还要给出代码实现。我们最后拿到的是3+1维（空间+时间）的时间序列，熵算法分为时域和频域，所以还涉及到数值傅里叶变换，这个后面再说。我大概查了下，表格中给出的算法基本上在谷歌上可以查到具体的代码，但是读懂原理，理解代码可能还是比较费劲。 现分配如下： 钮博恒和王悦晨：Time domain entropy全部 吴若岑：Renyi entropy 魏弘量：Spectral entropy 张峰瑞：Tsallis entropy 文章我已经查了一些了，在article文件夹里，可以直接参考。大家查的时候依据自己搜索习惯来查吧，知网啥的我没去看。 最后以md文
matlab toolbox 向量,RVM_matlabToolBox 相关向量机(RVM)的 源程序，包含快速算法，内含代码使用说明。 RVM采取 支持 238万源代码下载- www.pudn.c...
weixin_36048947的博客
03-17 272
文件名称: RVM_matlabToolBox下载 收藏√ [5 4 3 2 1]开发工具: matlab文件大小: 157 KB上传时间: 2013-11-21下载次数: 8提 供 者: 小芝麻详细说明：相关向量机(RVM)的matlab源程序，包含快速算法，内含代码使用说明。 RVM采取是与支持向量机相同的函数形式稀疏概率模型，对未知函数进行预测或分类。 优点： (1) 不仅仅...
样本熵(Sample Entropy)
BIG橙橙滴
10-03 9193
样本熵 Sample Entropy计算理解概念计算方法代码实现C语言版功能快捷键如何改变文本的样式插入链接与图片如何插入一段漂亮的代码片生成一个适合你的列表创建一个表格设定内容居中、居左、居右SmartyPants创建一个自定义列表如何创建一个注脚注释也是必不可少的KaTeX数学公式新的甘特图功能，丰富你的文章 计算理解 图片来源：Multiscale entropy analysis of ...
【Python】近似熵，样本熵，模糊熵计算高效版
Explore_Zhou
10-26 7336
最近在学习机器学习，发现对于与生物医学信号相关的机器学习任务，在选定特征时，各种针对时间序列的熵是绕不开的重要特征，诸如近似熵，样本熵，模糊熵等。因为它们所包含的信息要远比均值方差等特征要多得多，通过写python程序实现的过程中收获了不少，这里简单总结一下。计算各种熵的关键还是在于对计算方式的理解，如果博客说法不一，那就去查找文献，如果文献说法不一，那就去找提出这个熵的论文。
模糊熵、分布熵、近似熵、样本熵理论相关知识与代码实现
梁小憨憨的博客
04-03 7232
本篇为《信号处理》系列博客的第八篇，该系列博客主要记录信号处理相关知识的学习过程和自己的理解，方便以后查阅。 模糊熵、分布熵、近似熵、样本熵理论相关知识与代码实现模糊熵理论基础代码实现分布熵理论基础代码实现近似熵理论基础代码实现样本熵理论基础代码实现 模糊熵 理论基础 模糊熵(Fuzzy Entropy, FuzzyEn, FE) 衡量的也是新模式产生的概率大小，测度值越大，新模式产生的概率越大，即序列复杂度越大。 算法描述如下： 代码实现 %%% 模糊熵计算函数 %%% function fuzzye
信号处理算法（2）：样本熵（SampEn）
热门推荐
雨寒sgg的博客
07-10 6万+
  本文主要介绍样本熵的物理意义、计算过程以及Matlab里代码实现。 1.物理意义   样本熵（Sample Entropy，SampEn）与近似熵的物理意义相似，都是通过度量信号中产生新模式的概率大小来衡量时间序列复杂性，新模式产生的概率越大，序列的复杂性就越大。   样本熵的值越低，序列自我相似性就越高；样本熵的值越大，样本序列就越复杂。 2.计算过程   样本熵的计算过程如下： ...
样本熵代码
09-02
样本熵代码
样本熵程序代码
04-28
样本熵的Matlab程序代码，基于睡眠脑电信号的非线性特征的提取，可用于睡眠的分期
香农模式的内容和意义
09-22
香农模式的内容和意义，对学习有一定的帮助，可以供大家借鉴一下。
基于信息熵的属性约简算法c/c++代码
08-25
基于信息熵的属性约简算法 将数据存保存在excel文件中（新建文件默认有三个表单，删除后两个，只保留第一个。） 然后文件另存为txt格式，保存类型为 文本文件(制表符分隔)(*.txt) 元组最大个数为200，属性个数为50.（可以修改源代码） 测试的例子一个是论文里的实例，一个是网上的实例。 网页保存为pdf在文件夹里 源代码打包在压缩包
SampEn_熵_sampleentropy_样本sampen_样本熵_
10-01
样本熵是通过度量信号中产生新模式的概率大小来衡量时间序列复杂性，新模式产生的概率越大，序列的复杂性就越大。
SampEn.rar_SampEn_matlab sampen_样本熵_样本熵快速_熵 矩阵
07-15
通过matlab矩阵处理功能实现样本熵的快速计算。
SampEn样本熵计算
03-05
计算样本熵的Matlab程序，可用于机械故障诊断，生物医学信号处理（例如脑电、肌电、心电等）
样本熵：样本熵 (SampEn) 的可承受的矢量化实现。-matlab开发
05-29
此函数根据 Richman, JS 和 Moorman, JR (2000) 的建议计算样本熵 (SampEn) 算法。 该脚本是可承受的、压缩的和矢量化的。 因此，计算成本是最小的。 此外，还考虑了未定义 SampEn 的特殊情况： - 如果 B = 0，则...
python 样本熵sampen
最新发布
04-01
Python中可以用pyeeg库计算样本熵（sample entropy）。 安装pyeeg库： ```python !pip install pyeeg ``` 计算样本熵： ```python import pyeeg data = [1, 2, 3, 4, 5, 6, 5, 4, 5, 6, 5, 4, 3, 2, 1] m = 2 # 子序列长度 r = 0.5 * pyeeg.std(data) # 阈值 sampen = pyeeg.sampen(data, m, r) print(sampen) ``` 其中，`data`是待计算样本熵的数据，`m`是子序列长度，`r`是阈值。 输出结果： ``` 0.6931471805599453 ``` 参考文献： [1] https://pyeeg.readthedocs.io/en/latest/entropy.html#sample-entropy-sampen

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 2
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值