【熵与特征提取】从近似熵,到样本熵,到模糊熵,再到排列熵,包络熵,散布熵,究竟实现了什么?(第六篇)——“散布熵”及其MATLAB实现

server/2024/9/20 15:35:24/

今天讲散布熵,之前用了几篇文章分别讲述了功率谱熵、奇异谱熵、能量熵、近似熵、样本熵、模糊熵、排列熵、包络熵这8种类型的熵:

Mr.看海:【熵与特征提取】基于“信息熵”的特征指标及其MATLAB代码实现(功率谱熵、奇异谱熵、能量熵)

Mr.看海:【熵与特征提取】从近似熵,到样本熵,到模糊熵,再到排列熵,究竟实现了什么?(第一篇)——“近似熵”及其MATLAB实现

Mr.看海:【熵与特征提取】从近似熵,到样本熵,到模糊熵,再到排列熵,究竟实现了什么?(第二篇)——“样本熵”及其MATLAB实现

Mr.看海:【熵与特征提取】从近似熵,到样本熵,到模糊熵,再到排列熵,究竟实现了什么?(第三篇)——“模糊熵”及其MATLAB实现

Mr.看海:【熵与特征提取】从近似熵,到样本熵,到模糊熵,再到排列熵,究竟实现了什么?(第四篇)——“排列熵”及其MATLAB实现

Mr.看海:【熵与特征提取】从近似熵,到样本熵,到模糊熵,再到排列熵,包络熵,散布熵,究竟实现了什么?(第五篇)——“包络熵”及其MATLAB实现

散布熵和以上这些类型的熵一样,均是一种衡量时间序列复杂程度的指标,来刻画时间序列的不确定性。

一、散布熵

Rostaghi 和 Azami于 2016 年提出了散布熵[1]。该算法克服近似熵、样本熵与排列熵的部分缺陷,具有计算速度快、受突变信号影响较小等优点,在滚动轴承、齿轮箱等旋转机械特征提取及故障诊断中得到较好应用。

1.1 算法流程

散布熵算法流程是比较复杂的,如果大家对详细的理论讲解有兴趣,可以看一下提出该算法的论文原文[1]。

在这里,我们用相对通俗的表述大致描述一下计算过程,作为参考。

  1. 使用正态累积分布函数(NCDF)将时间序列x中的每个元素映射到0到1之间,得到一个新的时间序列y。这一步的目的是将时间序列x归一化,使其均值为0,标准差为1。
  2. 将时间序列y中的每个元素四舍五入到最近的整数,并将结果乘以类别数c,再加上0.5,得到一个新的时间序列z。这一步将y中的元素映射到1到c之间的整数,其中c是预先设定的类别数。
  3. 初始化一个长度为N-(m-1)*d的数组dp,用于存储每个嵌入向量对应的散布模式。
  4. 对于每个位置i (1 <= i <= N-(m-1)*d),提取一个长度为m的嵌入向量z_emb。嵌入向量z_emb包含z中从位置i开始,每隔d个元素取一个,共取m个元素。
  5. 将嵌入向量z_emb映射为一个散布模式,并将结果存储在dp(i)中。散布模式是一个m位的c进制数,每一位对应嵌入向量中的一个元素。例如,当m=3,c=6时,嵌入向量[1,4,5]对应的散布模式为(5-1)*6^2 + (4-1)*6^1 + (1-1)*6^0 = 870。
  6. 统计dp中每种散布模式出现的次数,得到一个长度为c^m的数组count。数组count的每个元素表示对应散布模式出现的次数。
  7. 将count除以散布模式的总数(即count的元素和),得到散布模式的概率分布p。
  8. 在计算散布熵之前,先移除概率为0的散布模式,以提高数值稳定性。
  9. 根据Shannon熵的定义,计算散布熵DE。散布熵的表达式为:
    DE = -sum(p .* log(p))
    其中,p是散布模式的概率分布,log是自然对数,sum表示对所有非零概率求和。

1.2 参数说明

由上述散布熵定义可知,与散布熵相关的主要参数为嵌入维数 m、类数 c 和时延 d,其相应取值会对散布熵计算结果产生一定影响[2]。

1.对于嵌入维数m而言,如果m太小,就很难检测到信号的动态行为。相反,如果m设置得太大,散布熵方法虽然可以得到更可靠的结果,但是它不能观察到小的变化,而且比较耗时。一般m取2或3居多。

2.对于类数c的影响,如果c太小,两个非常不同的振幅值很可能被分为同一类。如果c太大,则散布熵具有更大的计算负担,并且对噪声更敏感。c可以取6附近的整数值。

3.对于延时d的选择,当d大于1时,可能会出现混叠。另外,和排列熵一样,时延d对散布熵的估计影响很小。因此,根据上述研究结论,通常选择d为1,这样既能保证较高的计算效率,又能提供可靠的分析。

2.MATLAB代码实现

3.1 散布熵编程实现

散布熵目前网上还未找到特别靠谱的MATLAB代码,不过好在算法流程还是很清楚的,且算法提出者的论文中列举了案例,我们可以按照流程重新编写代码,并对照论文案例即可验证正确性。

对此,我编写了名为kDispersionEn的函数:

function DE = kDispersionEn(x, m, c, d)
% 散布熵(Dispersion Entropy),笔者编写,算法介绍见:https://zhuanlan.zhihu.com/p/694155930/
% 输入:
%   x: 单变量时间序列,一个长度为N的行向量
%   m: 嵌入维数
%   c: 类别数(通常建议取6)
%   d: 时间延迟(通常建议取1)
% 输出:
%   DE: 散布熵

然后拿论文中的数据验证一下:

导入数据上图中的数据x,并设置d=1,m=2,c=3。

% 参数设置
m = 2; % 嵌入维数
c = 3; % 类别数
d = 1; % 时间延迟x = [9,8,1,12,5,-3,1.5,8.01,2.99,4,-1,10];
% 计算散布熵
DE1 = dispersionEn(x, m, c, d);

在求散布熵过程中有一个关键步骤,即计算散布概率p。

论文中计算得到该概率值分别为:

在程序中,计算得到的每种散布模式的出现次数如下图:

这些数值即论文中每个p的分子,他们的和就是p的分母11。

至此基本验证成功,代码正确。将p值带入信息熵公式,得到散布熵为:1.8462

这里还发现了论文的一个小错误,论文中将该熵算成了1.8642,不过无伤大雅也就是了。

2.2 封装函数

为了特征提取代码的易用性,笔者对一系列熵特征提取进行了封装,包括上边添加注释的代码都集中到一起。由于搞科研写论文时,对特征提取的需要往往是集中性的、多种类的、需求各异的,所以我把之前介绍过的熵特征值和后边将会降到的集中熵特征进行了打包(上边的kDispersionEn函数也打包在其中了):

熵特征值共9个——功率谱熵、奇异谱熵、能量熵、近似熵、样本熵、排列熵、模糊熵、包络熵、散布熵

以上9种全都集中到一个封装函数里,实现一行代码完成特征提取。

比如提取数据“包络熵”就可以像这样写:

fea = genFeatureEn(data,{'enveEn'}) %对data求包络熵

如果提取数据“功率谱熵、奇异谱熵、能量熵、近似熵、样本熵、排列熵、模糊熵、包络熵、散布熵”这全部9种特征,就可以这样写:

fea =genFeatureEn(data,{'psdE','svdpE','eeE','ApEn', 'SpEn','FuzzyEn','PeEn','enveEn','DE'});  
%调用genFeature函数,完成特征提取,算出的特征值会保存在fea变量里

也就是说需要提取哪个特征,在函数中直接指定就可以了。输出的fea变量里就会得到相应的这些特征值,顺序也是与输入的排序保持一致的。

此外,针对同学们写论文要水图(哦不,科研)的需求,程序运行完后还可以画图这样的柱状图:

如果输入一维数据,得到的是柱状图

函数也可以输入二维数据,如果输入二维数据,则就是逐行求取熵特征。比如下图,导入的是一个有三行数据的二维数组的绘图结果:

如果输入的是二维数据,得到的是折线图,横坐标是数据组数

这个函数的介绍如下:

fea = genFeatureEn(Data,featureNamesCell,option);  %调用genFeature函数,完成特征提取,算出的特征值会保存在fea变量里,%fea变量的长度和featureNamesCell中指定的特征量一致,且顺序一一对应%程序运行完成后,在MATLAB的工作区,双击fea变量,可以查看求得的具体数值% function fea = genFeatureEn(data,featureNamesCell,options)
% 熵相关算法的信号特征提取函数
% 输入:
% data:待特征提取的时域信号,可以是二维数据,维度为m*n,其中m为数据组数,n为每组数据的长度。即每行数据为一组。行列方向不可出错
% options:其他设置,使用结构体的方式导入。目前可设置变量包括:
%   -svdpEn:即奇异值的窗口长度。
%   -Apdim:近似熵参数,Apdim为近似熵的模式维度
%   -Apr:近似熵参数,Apr为近似熵的阈值
%   -Spdim:样本熵参数,Spdim为样本熵的模式维度
%   -Spr:Spr为样本熵的阈值
%   -Fuzdim:模糊熵参数,Fuzdim为模糊熵模式维度
%   -Fuzr:模糊熵参数,Fuzr为模糊熵的阈值
%   -Fuzn:模糊熵参数,Fuzn为模糊熵权重
%   -Pedim:排列熵参数,Pedim为排列熵模式维度
%   -Pet:排列熵参数,Pet为排列熵的时间延迟
%   -DEm: 散布熵参数,DEm为散布熵模式维度
%   -DEc: 散布熵参数,DEc为散布熵类别数(通常建议取6)
%   -DEd: 散布熵参数,DEd为散布熵时间延迟(通常建议取1)
%   -fs:采样频率,采样频率即每秒钟采集的数据点数,按照实际情况设置,该参数目前在包络熵特征采集中用到
% featureNamesCell:拟进行特征提取的特征名称,该变量为cell类型,其中包含的特征名称为字符串,特征名称需要在下边列表中:
% 目前支持的特征(2024.4.23,共9种):
%      psdE:功率谱熵
%      svdpE:奇异谱熵
%      eeE:能量熵
%      ApEn:近似熵
%      SpEn:样本熵
%      FuzzyEn:模糊熵
%      PeEn:排列熵
%      enveEn:包络熵
%      DE:散布熵
% 
% 输出:
% fea:数据data的特征值数组,其特征值顺序与featureNamesCell一一对应

需要上边这个函数文件以及测试代码的同学,可以在公众号 khscience(看海的城堡)中回复“特征提取”获取。

3.其他:时域、频域特征提取的MATLAB代码实现

除了上述熵特征的提取,笔者还对之前文章中讲到过的时域和频域特征进行了代码实现,具体包括:

有量纲特征值8个——最大值、最小值、峰峰值、均值、方差、标准差、均方值、均方根值(RMS) 无量纲特征值6个——峭度、偏度、波形因子、峰值因子、脉冲因子、裕度因子 频域特征值5个——重心频率、均方频率、均方根频率、频率方差、频率标准差 谱峭度特征4个——谱峭度的均值、谱峭度的标准差、谱峭度的偏度、谱峭度的峭度

以上23种全都集中到一个封装函数里,实现一行代码完成特征提取。

比如提取数据“重心频率”就可以像这样写:

fea = genFeatureTF(data,{'FC'}) %对data数据求重心频率

如果提取数据“最大值、最小值、峰峰值、均值、方差、标准差、均方值...”这全部22种特征,就可以这样写:

fea =genFeatureTF(data,{'max','min','mean','peak','arv','var','std','kurtosis',...'skewness','rms','waveformF','peakF','impulseF','clearanceF',...'FC','MSF','RMSF','VF','RVF',...'SKMean','SKStd','SKSkewness','SKKurtosis'});  %调用genFeature函数,完成特征提取,算出的特征值会保存在fea变量里

也就是说需要提取哪个特征,在函数中直接指定就可以了。输出的fea变量里就会得到相应的这些特征值,顺序也是与输入的排序保持一致的。

这个函数的介绍如下:

function fea = genFeatureTF(data,fs,featureNamesCell)                                           
% 时域、频域相关算法的信号特征提取函数
% 输入:
% data:待特征提取的时域信号,可以是二维数据,维度为m*n,其中m为数据组数,n为每组数据的长度。即每行数据为一组。行列方向不可出错
% fs:采样频率,如果不提取频域特征,fs值可以设置为1
% featureNamesCell:拟进行特征提取的特征名称,该变量为cell类型,其中包含的特征名称为字符串,特征名称需要在下边列表中:
% 目前支持的特征(2022.5.23,共23种):
% max :最大值
% min :最小值
% mean :平均值
% peak :峰峰值
% arv  :整流平均值
% var  :方差
% std  :标准差
% kurtosis  :峭度
% skewness  :偏度
% rms       :均方根
% waveformF :波形因子
% peakF     :峰值因子
% impulseF  :脉冲因子
% clearanceF:裕度因子
% FC:重心频率
% MSF:均方频率
% RMSF:均方根频率
% VF:频率方差
% RVF:频率标准差
% SKMean:谱峭度的均值
% SKStd:谱峭度的标准差
% SKSkewness:谱峭度的偏度
% SKKurtosis:谱峭度的峭度
% 
% 输出:
% fea:数据data的特征值数组,其特征值顺序与featureNamesCell一一对应

需要上边这个函数文件以及测试代码的同学,可以在公众号 khscience(看海的城堡)中同样回复“特征提取”获取。

上述2个函数(熵特征提取函数“genFeatureEn”和时频特征提取函数“genFeatureTF”)会持续更新,有哪些想要加进去的特征指标,同学们可以在评论区留言,笔者会考虑纳入到这个“特征提取指标全家桶”中。


http://www.ppmy.cn/server/13791.html

相关文章

使用mac自带服务器(一行命令即可启动)

打开终端&#xff0c;开启Apache: 开启apache: sudo apachectl start 重启apache: sudo apachectl restart 关闭apache: sudo apachectl stop 启动后地址&#xff1a; http://127.0.0.1/ mac下Apache服务器的文件路径&#xff1a; 点击Finder 然后按住快捷键CommandShiftG 输入…

XiaodiSec day031 Learn Note 小迪安全学习笔记

XiaodiSec day031 Learn Note 小迪安全学习笔记 记录得比较凌乱&#xff0c;不尽详细 day31 上传漏洞 前置 基础内容在 ctfshow 中演示 中间件 cms 中的文件上传 开始 文件上传一般配合抓包 前台验证, 在前台改就可上传成功 php 后缀的文件有 php 后门&#xff0c;可连…

如何链接多个modbus_tcp设备,并将设备数据写入同一个modbusSlave,以便外部客户端获取所有链接设备的数据。

在modbus通信中&#xff0c;一个modbus服务器一次只能链接一个客户机&#xff0c;那么&#xff0c;外部客户端要获取多个设备的modbus数据&#xff0c;就需要使用链接一个专用的mosbus服务器&#xff0c;一下就是详细解决方法。 第一步&#xff1a;创建modbus客户端&#xff0…

LeetCode 1052. 爱生气的书店老板

题目链接 https://leetcode.cn/problems/grumpy-bookstore-owner/description/?envTypedaily-question&envId2024-04-23 先把最初的满意人数累加算出来&#xff0c;然后使用滑动窗口来模拟连续 minutes分钟不生气&#xff0c;计算不生气minutes分钟最大的满意数 class S…

LESS vs. SCSS:选择何种CSS预处理器?

CSS 预处理器如 LESS 和 SCSS 已成为现代 Web 开发中不可或缺的工具&#xff0c;它们为原生 CSS 提供了强大的变量、嵌套、混合&#xff08;Mixins&#xff09;、运算等高级特性&#xff0c;显著提升了样式表的可维护性和开发效率。尽管功能相似&#xff0c;LESS 和 SCSS 在语法…

计算机网络原原理学习资料分享---第一章/第一节(为有梦想的自己加油!)

计算机网络原理 课程知识框架 计算机网络原理 课程框架 第一章 计算机网络概述 重点 第二章 网络应用 重点 第三章 传输层 重点 难点 第四章 网络层…

云从科技AI智能体云月亮相中国铁建GSF项目展示中心

近日&#xff0c;中国铁建大湾区科学论坛永久会址项目综合展示体验中心&#xff08;以下简称“中国铁建GSF项目展示中心”&#xff09;迎来了一位特别的客服——云月数智人。云月是云从从容多模态大模型的融合承载体——AI智能体&#xff08;AI-Agent&#xff09;&#xff0c;她…

Mysql全局优化总结

Mysql全局优化总结 从上图可以看出SQL及索引的优化效果是最好的&#xff0c;而且成本最低&#xff0c;所以工作中我们要在这块花更多时间 服务端系统参数 官方文档&#xff1a;https://dev.mysql.com/doc/refman/8.0/en/server-system-variables.html#sysvar_max_connections…

生成对抗网络(GAN)

一、介绍 生成对抗网络&#xff08;Generative Adversarial Network&#xff0c;GAN&#xff09;是一种机器学习模型&#xff0c;由生成器和判别器两部分组成&#xff0c;用于生成逼真的图像、文本或音频等内容。GAN的作用主要包括以下几个方面&#xff1a; 图像生成&#xff…

Android Native Hook: 原理、方案对比与具体实现

文章目录 一、原理二、方案对比三、具体实现3.1 Inline Hook3.2 PLT/GOT Hook 四、实践案例&#xff1a;在Android应用中Hook open 函数4.1 Inline Hook实现4.2 PLT/GOT Hook实现 五、实践技巧和优化建议六、总结 在Android开发中&#xff0c;Hook技术是一种常用的技巧&#xf…

掌握Intellij-platform plugin插件开发必备知识(二)

插件开发基础 Intellij platform 插件可以使用java和Kotlin两种语言开发&#xff0c;开发插件时需遵守官方发布的特定的代码\开发规范&#xff0c;又由于Intellij平台插件体系非常庞大加之网上资料的匮乏&#xff0c;在某种程度上来讲初次接触还是有些难度的&#xff0c;对初学…

IDEA中添加servlet模板

官方代码链接 #if (${PACKAGE_NAME} && ${PACKAGE_NAME} ! "")package ${PACKAGE_NAME};#end #parse("File Header.java")import javax.servlet.*; import javax.servlet.http.*; import javax.servlet.annotation.*; import java.io.IOException…

redis的持久化

引入&#xff1a; ①MySQL的事务&#xff0c;有四个比较核心的特性&#xff1a; 1.原子性 2.一致性 3.持久性&#xff08;和持久化说的是一回事&#xff08;把数据存在硬盘上>持久 把数据存在内存中>不持久&#xff09;重启进程/重启主机之后&#xff0c;数据是否存…

【Java--数据结构】提升数据处理速度!深入理解Java中的顺序表机制

欢迎关注个人主页&#xff1a;逸狼 创造不易&#xff0c;可以点点赞吗~ 如有错误&#xff0c;欢迎指出~ 目录 两种创建顺序表的方法及区别 认识ArrayList的构造方法 不带参数的构造方法 带参数的构造方法 利用Collection 构造方法 举例 ArrayList 常用方法演示 add addAll remo…

idea连接Docker数据库

我们在docker下创建了数据库&#xff0c;想要更方便的查看和操作该数据库&#xff0c;idea和DataGrip或者其他软件都可以。在数据库连接时需要填写数据库名字&#xff0c;主机&#xff0c;端口&#xff0c;数据库用户名和密码。 输入之后先不要点击OK和按Enter键&#xff0c;我…

前端点击按钮触发复制文本

1. 效果展示&#xff1a; 点击复制小图标进行内容的复制 在这里我们先不考虑适用插件的情况&#xff0c;因为如果只是简单的复制&#xff0c;则不需要插件 2. 绑定事件 这里我们以vue为例子&#xff0c; 原生和react我后面补上 <i slot"prefix" class"i…

Linux系统IO

Linux系统中的IO函数主要包括两大类&#xff1a;标准C库中的函数和Linux系统调用。这些函数可以用于文件操作、网络通信、设备控制等多种IO任务。以下是Linux系统中常用的IO函数和系统调用的概述&#xff1a; 标准C库IO函数 这些函数是高级的、封装好的&#xff0c;并且与操作…

dedebiz文章内页自动给正文图片加了style样式怎么去掉

dede文章内页自动给正文图片加了style样式怎么去掉 打开&#xff1a;/system/archive/archives.class.php 查找&#xff1a;box-shadow&#xff0c;找到如下&#xff1a; margin:20px 0;box-shadow:0 1px 2px rgba(0,0,0,.1)改成下面这样&#xff1a; box-shadow:0 0px 0px rgb…

Dockers数据卷Volume

数据卷 Docker数据卷是什么 在 Docker 中&#xff0c;容器的文件系统是临时的&#xff0c;当容器被删除时&#xff0c;其文件系统也会随之销毁。但有时我们希望容器内的数据能够持久化&#xff0c;即使容器被删除后&#xff0c;数据也不会丢失&#xff0c;这时就可以使用数据…

idea中打印日志不会乱码,但是部署到外部tomcat中乱码了。

问题&#xff1a;如图Tomcat乱码&#xff0c;而且启动时的系统日志不会乱码&#xff0c;webapp中的打印日志才乱码。 idea中的情况如下&#xff1a;正常中文展示。 问题分析&#xff1a;网上分析的原因是Tomcat配置的字符集和web应用的字符集不匹配&#xff0c;网上集中的解决…