【熵与特征提取】从近似熵，到样本熵，到模糊熵，再到排列熵，究竟实现了什么？（第四篇）——“排列熵”及其MATLAB实现

今天讲排列熵，之前用了三篇文章分别讲述了近似熵、样本熵和模糊熵：

排列熵在原理上与前三种方法差异较大，所以理论部分要一定程度上抛弃惯性思维，接受新的算法理念。下面开始这个熵系列的最后一篇吧（也许）。

一、排列熵

排列熵（Permutation Entropy, PE）是由Bandt和Pompe[1]提出的一种检测时间序列随机性和动力学突变行为的方法，具有计算简单、快速，抗噪能力强，适合在线监测等优点，已经被广泛应用于肌电信号和心率信号分析，气温复杂度以及机械故障检测等。

1.1 排列熵的含义

假设待分析信号序列为x(1),x(2),…,x(N)。

（1）【序列定义】给定模式维数m，可以构造出一组m维矢量： $X(i)=[x(i),x(i+\tau),...x(i+\tau(m-1))]$

需要注意，上述序列有一个参数 $\tau$ ，被称为时间延迟，取值需为正整数。其实这个参数可以理解为对序列的降采样，例如当 $\tau=3$ 时就是每3个数据点进行一次采样，当 $\tau=1$ 时该序列就和近似熵的序列定义相同了。

其实在近似熵、样本熵、模糊熵中也有把这个时间延迟概念引入的情况，但是通常 $\tau$ 设置为1。

（2）【序列排序】对X(i)序列中的各个元素进行排序，例如X(i)中的元素为{20,60,10}时，排序就可以写成（2,3,1），排序中的数字就代表了每个元素排在第几个。（也有另外一种理解方式，会将排序写成（3,1,2），这样写就代表排在该顺序的是原序列的第几个数。这两种理解不影响最终的计算结果）。

（3）【排序出现的概率】对于m维的向量，可能出现的排列方式有 m！种。定义每一种排序出现的概率如下： $P(i)=\frac{Num(X(i))}{N-(m-1)\lambda}$

（4）【排列熵】根据信息熵定义，计算排列熵： $H_{PE}(m)=-\sum_{i=1}^{N-(m-1)\lambda}{P(i)log_{2}P(i)}$

注意，上边求对数我专门标出来了底数是2，因为此处在有些论文中取了以e为底，方法提出者也只是写了log而未标底数，但是笔者核验了原始论文中的案例，确认方法提出者使用的底数为2。不过在下一步归一化后，底数是几其实也无所谓了。

（5）【归一化】为了将熵值的范围调整到 0 到 1 的范围内，进行数据归一化： $H_{NPE}(m) = \frac{H_{PE}(m)}{log_{2}(m!)}$

1.2 排列熵的特点

排列熵具有算法简单，计算速度快、抗噪声能力强等优点，在机械故障诊断、生物信号处理等领域应用较为广泛。

1.3 排列熵的参数选择

排列熵的计算主要涉及m和 $\tau$ 这两个参数。

（1）m的选择。与近似熵、样本熵、模糊熵中m通常取1或2不同，排列熵的m推荐值为3~7。当嵌入维数 m 过小时，重构序列的概率模式较少。当嵌入维数 m过大，会导致计算时间增加，计算效率下降。换句话说，如果能够接受程序运行时间久，可以将m设置得尽量大一些。

（2） $\tau$ 的选择。如果信号数据量不是特别大，没有太大必要对信号降采样，通常设置为1即可。当然也有深入研究相空间重构对取不同值做影响研究的，这种情况该值需视情设置。

2.MATLAB代码实现

排列熵的代码在网上可以找到第三方版本，MATLAB官方库还没有收录。笔者在可靠版本上进行了修改和注释，对照上述理论做了部分修正，并逐行添加了注释，与上述文章中的算法步骤做了对照讲解。我将其命名为kPermutationEntropy，如下：

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

熵特征值共7个——功率谱熵、奇异谱熵、能量熵、近似熵、样本熵、 排列熵、模糊熵

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

比如提取数据“排列熵”就可以像这样写：

fea = genFeatureEn(data,{'PeEn'}) %对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一一对应
%

需要上边这个函数文件以及测试代码的同学，可以在下述链接获取：

特征提取（熵特征指标） | 工具箱文档

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一一对应

需要上边这个函数文件以及测试代码的同学，可以在下述链接获取：

特征提取（时域、频域指标） | 工具箱文档

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

参考

^[1] Bandt C , Pompe B . Permutation Entropy: A Natural Complexity Measure for Time Series[J]. Physical Review Letters, 2002, 88(17):174102.

文章出处登录后可见！

已经登录？立即刷新