今天讲散布熵，之前用了几篇文章分别讲述了功率谱熵、奇异谱熵、能量熵、近似熵、样本熵、模糊熵、排列熵、包络熵这8种类型的熵：
Mr.看海：【熵与特征提取】基于“信息熵”的特征指标及其MATLAB代码实现（功率谱熵、奇异谱熵、能量熵）
Mr.看海：【熵与特征提取】从近似熵，到样本熵，到模糊熵，再到排列熵，究竟实现了什么？（第一篇）——“近似熵”及其MATLAB实现
Mr.看海：【熵与特征提取】从近似熵，到样本熵，到模糊熵，再到排列熵，究竟实现了什么？（第二篇）——“样本熵”及其MATLAB实现
Mr.看海：【熵与特征提取】从近似熵，到样本熵，到模糊熵，再到排列熵，究竟实现了什么？（第三篇）——“模糊熵”及其MATLAB实现
Mr.看海：【熵与特征提取】从近似熵，到样本熵，到模糊熵，再到排列熵，究竟实现了什么？（第四篇）——“排列熵”及其MATLAB实现
Mr.看海：【熵与特征提取】从近似熵，到样本熵，到模糊熵，再到排列熵，包络熵，散布熵，究竟实现了什么？（第五篇）——“包络熵”及其MATLAB实现
散布熵和以上这些类型的熵一样，均是一种衡量时间序列复杂程度的指标,来刻画时间序列的不确定性。
一、散布熵
Rostaghi 和 Azami于 2016 年提出了散布熵[1]。该算法克服近似熵、样本熵与排列熵的部分缺陷，具有计算速度快、受突变信号影响较小等优点，在滚动轴承、齿轮箱等旋转机械特征提取及故障诊断中得到较好应用。
1.1 算法流程
散布熵的算法流程是比较复杂的，如果大家对详细的理论讲解有兴趣，可以看一下提出该算法的论文原文[1]。
在这里，我们用相对通俗的表述大致描述一下计算过程，作为参考。
- 使用正态累积分布函数(NCDF)将时间序列x中的每个元素映射到0到1之间,得到一个新的时间序列y。这一步的目的是将时间序列x归一化,使其均值为0,标准差为1。
- 将时间序列y中的每个元素四舍五入到最近的整数,并将结果乘以类别数c,再加上0.5,得到一个新的时间序列z。这一步将y中的元素映射到1到c之间的整数,其中c是预先设定的类别数。
- 初始化一个长度为N-(m-1)*d的数组dp,用于存储每个嵌入向量对应的散布模式。
- 对于每个位置i (1 <= i <= N-(m-1)*d),提取一个长度为m的嵌入向量z_emb。嵌入向量z_emb包含z中从位置i开始,每隔d个元素取一个,共取m个元素。
- 将嵌入向量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。
- 统计dp中每种散布模式出现的次数,得到一个长度为c^m的数组count。数组count的每个元素表示对应散布模式出现的次数。
- 将count除以散布模式的总数(即count的元素和),得到散布模式的概率分布p。
- 在计算散布熵之前,先移除概率为0的散布模式,以提高数值稳定性。
- 根据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”）会持续更新，有哪些想要加进去的特征指标，同学们可以在评论区留言，笔者会考虑纳入到这个“特征提取指标全家桶”中。
