【熵与特征提取】从近似熵，到样本熵，到模糊熵，再到排列熵，究竟实现了什么？（第二篇）——“样本熵”及其MATLAB实现
上一篇文章中介绍了近似熵的概念和代码实现，同时也指出了近似熵的缺陷。
样本熵是对近似熵算法的改进，也是当前运用比较广泛的一种熵特征值计算方法。 样本熵也是时间序列复杂度的 一种度量，样本熵越小，时间序列复杂度越小，自相似性越高。通常用于故障诊断和生物时间序列分析当中。
一、样本熵
样本熵（Sample Entropy）是Richman等人提出的一种与近似熵不同的不计数自身匹配的统计量[1]。
样本熵是在近似熵基础上做的一种改进方法，所以在理论上存在继承关系。这里推荐同学们先读一下上一篇文章：
在理解近似熵的理论后再看下述介绍。
1.1 样本熵的含义
（1）【序列定义】样本熵的序列定义和近似熵是相同的，这里不再赘述。
（2）【近似的含义】与近似熵不同，遍历所有的X(i)和X(j)的组合时，i的范围是从1到N-m，且不包括i与j相等的情况，所以总共会有N-m-1种组合[1]。
（3）【近似比例】将近似数量与总数量的比值记为 B_{i}^{m}(r)，注意此时的总数量是N-m-1
（4）此时再求得：B^{m}(r)=\frac{1}{N-m}\sum_{i=1}^{N-m}{B_{i}^{m}}(r) ，注意，相对于近似熵的步骤，这里没有求对数，这样就避免了出现ln0的情况。
（5）将维数加1变为m+1，重复以上步骤，得到 A^{m}(r)
（6）【样本熵】求得信号序列的样本熵为：SampEn(m,r)=lnB^{m}(r)-lnA^{m}(r)
1.2 样本熵的特点
相对于近似熵，样本熵体现出以下更好的性质[2]：
（1）样本熵计算近似时，不包含与自身数据段的比较，计算误差较小，且不依赖数据长度。
（2）样本熵具有更好的一致性。如果一时间序列比另一时间序列有较高的SampEn值，那么对于其他的r和m值，也具有较高的SampEn值。
（3）样本熵对于丢失数据不敏感，即使数据丢失多达1/3，对SampEn计算值影响很小。
1.3 样本熵的参数选取
样本熵的计算主要涉及m和r这两个参数。
（1）m的选择。m可以被称为模式维数，它是用来计算样本熵的时序列的窗口长度。通常m取1或者2，且取m=2更为多见。
（2）r的选择。r即相似容限阈值。当r的值较大时，会丢失较多的信息，而当r的值较小时，又会不能理想的估计出系统的统计特性。与近似熵相同，通常r的取值在0.1和0.25SD(x)(SD(x)是序列的标准差)之间。
2.MATLAB代码实现
目前MATLAB官方还没有样本熵函数，网上流传的样本熵代码版本比较多，笔者对比了一下运算出来的结果虽然相差不多但是多少都有一些差别。所以笔者自己写了一个版本，是比较严格地按照上述理论流程写的，并逐行添加了注释，并与上述文章中的算法步骤做了对照讲解。我将其命名为kSampleEn，如下：
为了特征提取代码的易用性，笔者对一系列熵特征提取进行了封装，包括上边添加注释的代码都集中到一起。由于搞科研写论文时，对特征提取的需要往往是集中性的、多种类的、需求各异的，所以我把之前介绍过的熵特征值和后边将会讲到的几种熵特征进行了打包：
熵特征值共7个——功率谱熵、奇异谱熵、能量熵、近似熵、样本熵、排列熵、模糊熵
以上7种全都集中到一个封装函数里，实现一行代码完成特征提取。
比如提取数据“样本熵”就可以像这样写：
fea = genFeatureEn(data,{'SpEn'}) %对data求样本熵如果提取数据“功率谱熵、奇异谱熵、能量熵、近似熵、样本熵、排列熵、模糊熵”这全部7种特征，就可以这样写：
fea =genFeatureEn(data,{'psdE','svdpE','eE', 'ApEn', 'SpEn','PeEn','FuzzyEn'});
%调用genFeature函数，完成特征提取，算出的特征值会保存在fea变量里也就是说需要提取哪个特征，在函数中直接指定就可以了。输出的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为排列熵的时间延迟
% featureNamesCell：拟进行特征提取的特征名称，该变量为cell类型，其中包含的特征名称为字符串，特征名称需要在下边列表中：
% 目前支持的特征（2022.5.23，共7种）：
% psdE：功率谱熵
% svdpE：奇异谱熵
% eE：能量熵
% ApEn：近似熵
% SampleEn：样本熵
% FuzzyEn：模糊熵
% PerEn：排列熵
%
% 输出：
% 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”）会持续更新，有哪些想要加进去的特征指标，同学们可以在评论区留言，笔者会考虑纳入到这个“特征提取指标全家桶”中。
