排列熵、模糊熵、近似熵、样本熵的原理及MATLAB实现

一、排列熵

1.概念

一种检测动力学突变和时间序列随机性的方法，能够定量评估信号序列中含有的随机噪声。

2.基本原理

(1) 对一组长度为N的时间序列X进行相空间重构，得到矩阵Y为

其中，m为嵌入维数，t为延迟时间，K=N-(m-1)t。
矩阵Y中的每一行都是一个重构分量，共有K个重构分量。

(2) 将每一个重构分量按照升序重新排列，得到向量中各元素位置的列索引构成一组符号序列。
m维相空间映射不同的符号序列总共有m!种。

(3) 计算每一种符号序列出现的次数除以m!种不同的符号序列出现的总次数作为该符号序列出现的概率。
(4) 时间序列X的排列熵的计算公式为：
(5) 排列熵的最大值为ln(m!)，将排列熵值进行归一化处理，即：
排列熵值的大小表示时间序列X的随机程度：熵值越小，说明时间序列越简单、规则；反之，熵值越大，则时间序列越复杂、随机。

说明： 我之前看排列熵的理论一直看不懂，搞不清排列熵的计算过程，直到我认真阅读排列熵算法的代码，并且结合理论，才终于搞清楚。建议读者阅读排列熵的代码时，边调试边思考，有助于理解。

补充说明

关于计算每一种符号序列出现的次数除以m!种不同的符号序列出现的总次数作为该符号序列出现的概率的理解
假设有时间序列
X={x(1),x(2),x(3),x(4),x(5),x(6),x(7),x(8),x(9),x(10),x(11),x(12)}={2,5,3,4,6,7,9,5,3,7,9,2};
计算信号的排列熵时，设延迟时间t为1，嵌入维数m为3，则时间序列重构后的矩阵Y为
{x(1),x(2),x(3);
x(2),x(3),x(4);
x(3),x(4),x(5);
x(4),x(5),x(6);
x(5),x(6),x(7);
x(6),x(7),x(8);
x(7),x(8),x(9);
x(8),x(9),x(10);
x(9),x(10),x(11);
x(10),x(11),x(12)}
矩阵Y中的每一行都是一个重构分量,所以重构分量的个数K=10，则每一个重构分量按照升序重新排列后各元素位置的列索引构成一组符号序列，所以矩阵Y各个元素的索引构成的矩阵为
{1,3,2;
2,3,1;
1,2,3;
1,2,3;
1,2,3;
3,1,2;
3,2,1;
2,1,3;
1,2,3;
3,1,2;}
由于嵌入维数m=3,符号序列共有3!=6种:（(1,2,3)，(1,3,2)，(2,1,3)，(2,3,1)，(3,1,2)，(3,2,1)）。
上述10个重构分量的符号序列中，不同的符号序列出现的次数为：
(1,2,3)为4次，(1,3,2)为2次，(2,1,3)为1次，(2,3,1)为0次，(3,1,2)为2次，(3,2,1)为1次，所以6种符号序列出现的总次数为10次(实际只有5种符号序列，即此时k=5，小于3!=6种，这句话是对第(2)步中k的解释)。

3.MATLAB代码

%% 主函数调用排列熵函数求时间序列的排列熵值
load X;
[m,n]=size(X);  % X为时间序列,一行为一个时间序列。
% 求排列熵:pe为排列熵
M = 3;  % 嵌入维数
T = 1;  % 延迟时间
pe=zeros(1,m);
for i=1:m
    [pe(i),~] = PermutationEntropy(X(i,:),M,T);
end

%% 排列熵算法
function [pe ,hist] = PermutationEntropy(y,m,t)

%  Calculate the permutation entropy(PE)
%  排列熵算法的提出者：Bandt C，Pompe B. Permutation entropy:a natural complexity measure for time series[J]. Physical Review Letters,2002,88(17):174102.

%  Input:   y: time series;
%           m: order of permuation entropy 嵌入维数
%           t: delay time of permuation entropy,延迟时间

% Output: 
%           pe:    permuation entropy
%           hist:  the histogram for the order distribution
ly = length(y);
permlist = perms(1:m);
[h,~]=size(permlist);
c(1:length(permlist))=0;

 for j=1:ly-t*(m-1)
     [~,iv]=sort(y(j:t:j+t*(m-1)));
     for jj=1:h
         if (abs(permlist(jj,:)-iv))==0
             c(jj) = c(jj) + 1 ;
         end
     end
 end
hist = c;
c=c(c~=0);
p = c/sum(c);
pe = -sum(p .* log(p));
% 归一化
pe=pe/log(factorial(m));
end

参考文献

[1]赵昕海,张术臣,李志深,等.基于VMD的故障特征信号提取方法[J].振动、测试与诊断,2018,38(1):11-13.
[2]排列熵matlab实现