排列熵matlab实现

排列熵算法

熵可以表征信号的复杂性以及度量信息的不确定性，适用于处理非线性问题。排列熵（Permutation entropy，PE）通过比较相邻时间序列的值，检测时间序列的动态变化。

排列熵算法基本原理

对一长度为N的时间序列 { X ( i ) , i = 1 , 2 , ⋅ ⋅ ⋅ , N } \left \{ X\left ( i \right ),i= 1,2,\cdot \cdot \cdot ,N\right \} {X(i),i=1,2,⋅⋅⋅,N}，进行相空间重构
$$\left[\begin{array}{c}x\left(1\right)x\left(1+\tau \right)\cdot \cdot \cdot x\left[1+\left(m-1\right)\tau \right]\\ x\left(2\right)x\left(2+\tau \right)\cdot \cdot \cdot x\left[2+\left(m-1\right)\tau \right]\\ x\left(j\right)x\left(j+\tau \right)\cdot \cdot \cdot x\left[j+\left(m-1\right)\tau \right]\\ ⋮\\ x\left(K\right)x\left(K+\tau \right)\cdot \cdot \cdot x\left[K+\left(m-1\right)\tau \right]\end{array}\right],j=1,2,\cdot \cdot \cdot ,K$$
m为嵌入维数，$\tau$为延迟时间，K=N-(m-1)$\tau$。
将$X\left(i\right)$重构矩阵中的第 j 个重构分量 $x\left(j\right),x\left(j+\tau \right),\cdot \cdot \cdot ,x\left[j+\left(m-1\right)\tau \right]$按照数值大小进行升序排列，$j_1,j_2,\cdot \cdot \cdot j_m$表示重构分量中各元素所在列的索引，即$$x\left[i+\left(j_1-1\right)\tau \right]⩽x\left[i+\left(j_2-1\right)\tau \right]⩽\cdot \cdot \cdot ⩽x\left[i+\left(j_m-1\right)\tau \right]$$
若重构分量中存在相等的值，则按照$j_1,j_2$的大小排序，即当$j_1<j_2$时，有$x\left[i+\left(j_1-1\right)\tau \right]⩽x\left[i+\left(j_2-1\right)\tau \right]$。任意时间序列$X\left(i\right)$经相空间重构所得的重构矩阵的每一行都能得到一组符号序列：$S\left(l\right)=\left(j_1,j_{2,\cdot \cdot \cdot j_m}\right)$，$l=1,2,\cdot \cdot \cdot ,k$，且$k⩽m!$。
$m$维相空间映射有$m!$种不同的符号序列，$S\left(l\right)$只是其中一种符号序列。若将每一种符号序列出现的概率记为 $P_1,P_2,\cdot \cdot \cdot P_k$，按照信息熵的定义，时间序列$X\left(i\right)$的$k$种不同符号序列的排列熵定义为
$$H_{PE}\left(m\right)=-\sum _{j=1}^k{P}_{j}ln{P}_j$$

matlab代码：

function pe= pec(y,m,t)
%  Calculate the permutation entropy
%  Input:   y: time series;
%           m: order of permuation entropy
%           t: delay time of permuation entropy, 
% Output: 
%           pe:    permuation entropy
ly = length(y);
permlist = perms(1:m);
c(1:length(permlist))=0;
    
 for j=1:ly-t*(m-1)
     [a,iv]=sort(y(j:t:j+t*(m-1)));
     for jj=1:length(permlist)
         if (abs(permlist(jj,:)-iv))==0
             c(jj) = c(jj) + 1 ;
         end
     end
 end
c=c(find(c~=0));
p = c/sum(c);
pe = -sum(p .* log(p));

取材于matlab官方交换文档，可根据自己需求对程序进行相应的修改。
链接: https://ww2.mathworks.cn/matlabcentral/fileexchange/37289-permutation-entropy.

嵌入维数 m 和延迟时间 $\tau$ 的选取

使用平均互信息（AMI）估计相空间重建的延迟时间，使用虚假最近邻点（FNN）算法估计相空间重构的嵌入维数。

[~,eLag,eDim] = phaseSpaceReconstruction(Y);

eDim为估计的嵌入维数，eLag为估计的延迟时，phaseSpaceReconstruction函数为matlabR2019a的Predictive Maintenance Toolbox自带的相空间重构函数。
phaseSpaceReconstruction: https://ww2.mathworks.cn/help/predmaint/ref/phasespacereconstruction.html.