特征提取（复合多尺度熵）

最后更新于：2024-06-28 23:32:15

一、代码运行环境

MATLAB2018a及更新版本。

二、程序介绍

注：图标代表该m文件为脚本文件，可以直接运行；图标代表函数文件，在没有输入变量的情况下无法直接运行。更详细的解释可以看这里。

1.testGenFeaCMSE.m文件

特征提取函数genFeatureMSE的测试文件，可以直接运行。程序运行完成后，在MATLAB的工作区，双击mseFea变量，可以查看求得的具体数值。

可以提取的复合多尺度熵类型包括：复合多尺度样本熵、复合多尺度模糊熵、复合多尺度排列熵、复合多尺度包络熵、复合多尺度功率谱熵、复合多尺度能量熵、复合多尺度奇异谱熵、复合多尺度散布熵（完整版20240628更新）。可以根据需求灵活选取计算所有或者部分复合多尺度熵。

代码运行后，会绘制不同尺度粗粒度信号图（完整版独有）、不同尺度粗粒度信号图-取平均后（完整版独有）和各类型熵的多尺度计算结果，如下图：

2.genFeatureCMSE.m文件

复合多尺度熵特征提取函数。9种全都集中到一个封装函数里，只需要输入待提取数据、采样频率和特征类型，实现一行代码完成特征提取。

需注意该函数文件不能直接运行，需参照testGenFeaCMSE的演示案例进行调用。此外需注意，因复合多尺度熵本身特征限制，该函数仅对一维数据求取复合多尺度熵。

function cmseFea = genFeatureCMSE(data, scales, featureNamesCell, options)
    % 复合多尺度熵特征提取函数
    % 输入：
    % data：待特征提取的时域信号，必须是一维数据
    % scales：要分析的尺度数组
    % featureNamesCell：拟进行特征提取的特征名称，该变量为cell类型，其中包含的特征名称为字符串，特征名称需要在下边列表中：
    % 目前支持的特征（2024.6.28，共9种）：
    %      psdE：功率谱熵
    %      svdpE：奇异谱熵
    %      eeE：能量熵
    %      ApEn：近似熵
    %      SpEn：样本熵
    %      FuzzyEn：模糊熵
    %      PeEn：排列熵
    %      enveEn：包络熵
    %      DE：散布熵（需要完整版）
    % 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：采样频率，采样频率即每秒钟采集的数据点数，按照实际情况设置，该参数目前在包络熵特征采集中用到
    %   -figflag：是否画图，'on'为画图，'off'为不画图
    
    % 
    % 输出：
    % cmseFea：不同尺度下的特征值，是二维数据，每行代表一种熵的计算结果
    % cmseFea的行数和featureNamesCell中指定的特征量一致，且顺序一一对应

该文件中还包含了功率谱熵、奇异谱熵、能量熵、近似熵、样本熵、排列熵、模糊熵、包络熵、散布熵（完整版20240628更新）的子函数，函数介绍如下：

function ie = kInformationEntopy(sig,SegmentNum)

% 计算信号的信息熵
% 参考《矿用带式输送机托辊远程故障诊断系统》
% 输入：
% sig：输入信号
% SegmentNum：拟分组数，如果不输入，则自动使用斯特格斯（Sturges）经验公式计算。
% 输出：
% ie：信息熵求解结果

function svdpE = kSingularSpectrumEntropy(data,n)

% 求信号的奇异谱熵
% 参考《矿用带式输送机托辊远程故障诊断系统》
% 输入：
% data：待分析信号
% n：窗口长度，注意：2<=n<=lenght(data)-1
% 输出：
% svdpE：奇异谱熵值

function psdE = kPowerSpectrumEntropy(data)

% 求信号的功率谱熵
% 功率谱使用周期图法计算得到
% 参考《矿用带式输送机托辊远程故障诊断系统》
% 输入：
% data：待分析信号
% 输出：
% psdE：功率谱熵值

function eE = kEnergyEntropy(data)

% 求信号基于emd分解算法的能量熵
% 参考《面向高铁走行部故障诊断算法的研究与实现》
% 输入：
% data：待分析信号
% 输出：
% eE：能量熵值

function ApEn = kApproximateEntropy(data, dim, r)


%  计算近似熵(ApEn)，在网上流传的ApproximateEntropy函数基础上修改
%  输入：
%  data  - 待分析数据，需要是一维数据
%  dim - 模式维度
%  r  -  阈值大小，一般选择r=0.1~0.25，（网上流传的ApproximateEntropy需要再乘以data数据的标准差，本程序不需要，只用输入0.1-0.25间的数字即可）
%  输出：
%  ApEn：近似熵值大小
%  (网上流传的ApproximateEntropy还有一个降采样参数tau，此函数中不进行降采样，删去该参数)
%  论文：Pincus S M . Approximate entropy as a measure of system complexity[J]. Proceedings of the National Academy of Sciences ,1991,88(6):2297—2301.
% 理论讲解：https://zhuanlan.zhihu.com/p/494761890/

function kSampleEnValue = kSampleEn(data, dim, r)


% 计算data序列的样本熵
% 输入：
% data：一维数据
% dim：重构维数，一般选择1或2，优先选择2，一般不取m>2
% r：阈值大小，一般选择r=0.1~0.25，（网上流传的SampEn需要再乘以data数据的标准差，本程序不需要，只用输入0.1-0.25间的数字即可）
% 输出：
% kSampleEnValue：样本熵值大小
% 论文：Richman J S, Moorman J R. Physiological time-series analysis using approximate entropy and sample entropy[J]. American Journal of Physiology-Heart and Circulatory Physiology, 2000.
% 理论讲解：https://zhuanlan.zhihu.com/p/519042619/

function FuzEn = kFuzzyEntropy(data,dim,r,n)


% 计算模糊熵，在网传FuzzyEntropy.m基础上修改并注释
% 输入：
% data：一维数据
% dim：重构维数，一般选择1或2，优先选择2，一般不取m>2
% r：阈值大小，一般选择r=0.1~0.25，（网上流传的SampEn需要再乘以data数据的标准差，本程序不需要，只用输入0.1-0.25间的数字即可）
% n：权重，取较小的整数，如2或者3。推荐设置为2
% 输出：
% FuzEn：模糊熵值大小

% 论文：Chen W , Wang Z , Xie H , et al. Characterization of Surface EMG Signal Based on Fuzzy Entropy[J]. IEEE Transactions on Neural Systems and Rehabilitation Engineering, 2007, 15(2):266-272.
% 理论讲解：https://zhuanlan.zhihu.com/p/519633544/

function pe = kPermutationEntropy(data,m,t)


% 计算排列熵，在网传PermutationEntropy.m基础上修改并注释
% 输入：
% data：一维数据
% dim：重构维数，推荐值为3~7，如果能够接受程序运行时间久，可以将m设置得尽量大一些。
% t：时间延迟，取整数，设置为1的情况较多。
% 输出：
% pe：排列值大小

% 论文：Bandt C，Pompe B. Permutation entropy:a natural complexity measure for time series[J]. Physical Review Letters,2002,88(17):174102.
% 理论讲解：https://zhuanlan.zhihu.com/p/519809254/

function enveEn = kEnvelopeEntropy(data,fs)


% 计算包络熵
% 输入：
% data：一维数据
% fs：采样频率，采样频率即每秒钟采集的数据点数，按照实际情况设置
% 输出：
% enveEn：包络熵值
% 论文《自适应最大相关峭度解卷积方法及其在轴承早期故障诊断中的应用》

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

三、快速开始

1.运行测试脚本

先在MATLAB里打开下载好的文件夹，然后运行testGenFeaCMSE.m程序，程序运行完毕后如果没报错，且工作区成功得到mseFea变量，图像正常绘制，则说明运行环境正常，程序正确。

2.导入数据

复制一个testGenFeaCMSE.m的文件副本，在这个副本里做如下修改：

根据你的文件类型的不同（excel,txt,csv等），将数据导入MATLAB的方法有所不同。同学们可以看博主针对常用文件的导入方法的这个教程，教程上未包含的数据类型，大家可以再参考这个文档。

如果你已经实现了数据导入，这时候应该拥有了一个一维数据变量，这时候就可以调用函数进行特征提取了。

3.实现特征提取

参照testGenFeaCMSE.m文件中第二步，genFeatureCMSE函数的调用方法，根据实际采样频率和特征提取类型进行修改设置，运行程序即可。

四、关于完整版与公开版代码

五、获取公开版程序（需使用电脑浏览器打开）

注：公开版代码需使用MATLAB2022a及以上版本。

六、获取完整版程序（使用电脑浏览器或者手机浏览器打开）

获取通道一（淘宝）：点击此处获取完整版程序

获取通道二（本页面）：点击下面“立即支付”按钮，付款后获取完整版代码下载链接和售后联系方式~本通道处于测试阶段，使用该通道可以额外优惠（仅需53元）。付款完成后刷新一下本页面即可看到下载链接。

（注意支付跳转失败的话，请使用浏览器打开本页面）

七、完整版代码重要更新

20240628 更新了复合多尺度散布熵（需完整版代码）

20240217 完成初版代码

八、常见问题

1.程序运行了很长时间，一直跑不出结果。

原因：这通常是由于分析数据的数据量比较大，并且对奇异谱熵（svdpE）特征值进行提取导致的。由于程序运行奇异值分解时非常吃电脑资源，所以如果提取这个特征，运行时间可能会非常长。

建议：1.如果不是必须求奇异谱熵这个特征，那么建议在程序运行时间过长时不再提取这个特征，这样做可以立竿见影，也不会影响其他特征的提取。2.如果一定要计算奇异谱熵，那么建议对数据进行重采样或者数据截段，减小待分析数据的数据量。

2.求出来的能量熵（eE）的值为0

原因：可能是分析对象数据是经emd分解后的imf分量。需要注意，请不要对数据进行emd分解，然后把emd分解得到的某一imf分量求能量熵，因为在能量熵计算函数里会对数据进行emd分解，imf分量再次进行emd分解将只能得到它本身这一个分量，而对于一个常数求信息熵得到的结果就是0.

建议：请不要对经emd分解后的单一imf分量求能量熵。