模糊熵、分布熵、近似熵、样本熵理论相关知识与代码实现

本篇为《信号处理》系列博客的第八篇，该系列博客主要记录信号处理相关知识的学习过程和自己的理解，方便以后查阅。

模糊熵

理论基础

模糊熵(Fuzzy Entropy, FuzzyEn, FE) 衡量的也是新模式产生的概率大小，测度值越大，新模式产生的概率越大，即序列复杂度越大。

算法描述如下：

代码实现

%%% 模糊熵计算函数 %%%
function fuzzyen = Fuzzy_Entropy( dim, r, data, tau )
% FUZZYEN Fuzzy Entropy
%   calculates the fuzzy entropy of a given time series data

% Similarity definition based on vectors' shapes, together with the
% exclusion of self-matches, earns FuzzyEn stronger relative consistency
% and less dependence on data length.

%   dim     : embedded dimension 
%   r       : tolerance (typically 0.2 * std)
%   data    : time-series data
%   tau     : delay time for downsampling (user can omit this, in which case
%             the default value is 1)
%

if nargin < 4, tau = 1; end
if tau > 1, data = downsample(data, tau); end

N = length(data);
Phi = zeros(1,2);

for m = dim:dim+1
    
    Ci = zeros(1,N-m+1);
    dataMat = zeros(m,N-m+1);
    
    % setting up data matrix - form vectors
    for j = 1:m
        dataMat(j,:) = data(j:N-m+j);
    end
    
    % baseline
    U0 = mean(dataMat);
    % remove baseline and calculate the absolute values
    Sm = abs(dataMat - repmat(U0,m,1));
    
    % Given vector Si, calculate the similarity degree between its'
    % neighboring vector Sj
    for i = 1:N-m+1
        
        Sm_tmp = Sm;
        Sm_tmp(:,i) = []; % excluded self-matches
        % maximum absolute difference of the corresponding scalar components
        % of Si and Sj (j≠i)
        dij = max(repmat(Sm(:,i),1,N-m) - Sm_tmp);
        % similarity degree
        Aij = exp(-log(2)*(dij/r).^2);
        % averaging all the similarity degree of its neighboring vectors Sj
        Ci(i) = sum(Aij)/(N - m);
        
    end
    
    % summing over the counts
    Phi(m-dim+1) = sum(Ci)/(N-m+1); % φ_m and φ_m+1
    
end

fuzzyen = log(Phi(1))-log(Phi(2));  % fuzzyen = ln(φ_m)-ln(φ_m+1)

end

关于一些参数的选择：

Detection of epileptic seizure based on entropy analysis of short-term EEG

分布熵

理论基础

分布熵(DistEn) 较适应于短数据序列

分布熵解决了样本熵和近似熵的参数依赖性以及非
鲁棒性的缺点， 充分利用分析时间序列的状态空间对应项来量化向量间距离的分布特
性

算法描述如下：

代码实现

%%% 分布熵计算函数 %%%
function DistEn = Distribution_Entropy(dim, data, B, tau)
% DISTEN Distribution Entropy
%   calculates the distribution entropy of a given time series data

%   dim     : embedded dimension 
%   r       : tolerance (typically 0.2 * std)
%   data    : time-series data
%   tau     : delay time for downsampling (user can omit this, in which case
%             the default value is 1)
%   B       : 直方图的参数，分割数据的份数

if nargin < 4, tau = 1; end
if tau > 1, data = downsample(data, tau); end

m = dim;

N = length(data);
dataMat = zeros(m,N-m+1);

% setting up data matrix - form vectors
for j = 1:m
	dataMat(j,:) = data(j:N-m+j);
end

Sm = dataMat;
% Given vector Si, calculate the similarity degree between its'
% neighboring vector Sj
for i = 1:N-m+1
AE
	Sm_tmp = Sm;
	Sm_tmp(:,i) = []; % excluded self-matches
	% maximum absolute difference of the corresponding scalar components
	% of Si and Sj (j≠i)
	dij = max(repmat(Sm(:,i),1,N-m) - Sm_tmp);

    D(i,:)  = dij;% 生成N-m+1*N-m的矩阵
end

Dij=reshape(D,1,(N-(m-1))*(N-(m-1)-1));% 将D重置为1维矩阵，为啥呢？？？

% Dij：一维的向量
% B  ：分割的数量
% n  ：在一个分割区域内，Dij数据出现的次数
% x  ：一个分割区域的中心线位置坐标
[n,x]=hist(Dij,B);
f=n/length(Dij);% 计算概率密度

for i=1:B
    if f(i)==0
        P(i)=0;% 如果概率密度为零，则求和值置零
    else
        P(i)= (f(i))*log2(f(i));
    end
end
    
DistEn = (-1/log2(B))*sum(P);

end

近似熵

理论基础

近似熵(Approximate Entropy, ApEn) 是一种用于量化时间序列波动的规律性和不可预测性的非线性动力学参数，它用一个非负数来表示一个时间序列的复杂性，反映了时间序列中新信息发生的可能性，越复杂的时间序列对应的近似熵越大

算法描述如下：

代码实现

%%% 近似熵计算函数 %%%
function apen = Approximate_Entropy( dim, r, data, tau )
%ApEn
%   dim  : embedded dimension
%   r    : tolerance (typically 0.2 * std)
%   data : time-series data
%   tau  : delay time for downsampling

%   Changes in version 1
%       Ver 0 had a minor error in the final step of calculating ApEn
%       because it took logarithm after summation of phi's.
%       In Ver 1, I restored the definition according to original paper's
%       definition, to be consistent with most of the work in the
%       literature. Note that this definition won't work for Sample
%       Entropy which doesn't count self-matching case, because the count
%       can be zero and logarithm can fail.
%
%       A new parameter tau is added in the input argument list, so the users
%       can apply ApEn on downsampled data by skipping by tau.
%---------------------------------------------------------------------
% coded by Kijoon Lee,  kjlee@ntu.edu.sg
% Ver 0 : Aug 4th, 2011
% Ver 1 : Mar 21st, 2012
%---------------------------------------------------------------------
if nargin < 4, tau = 1; end
if tau > 1, data = downsample(data, tau); end

N = length(data);
result = zeros(1,2);

for j = 1:2
    m = dim+j-1;
    phi = zeros(1,N-m+1);
    dataMat = zeros(m,N-m+1);
    
    % setting up data matrix
    for i = 1:m
        dataMat(i,:) = data(i:N-m+i);
    end
    
    % counting similar patterns using distance calculation
    for i = 1:N-m+1
        tempMat = abs(dataMat - repmat(dataMat(:,i),1,N-m+1));
        boolMat = any( (tempMat > r),1);
        phi(i) = sum(~boolMat)/(N-m+1);
    end
    
    % summing over the counts
    result(j) = sum(log(phi))/(N-m+1);
end

apen = result(1)-result(2);

end

样本熵

理论基础

样本熵(Sample Entropy, SampEn, SE) 是基于近似熵(ApEn)的一种用于度量时间序列复杂性的改进方法，在评估生理时间序列的复杂性和诊断病理状态等方面均有应用

由于样本熵是近似熵的一种改进方法，因此可以将其与近似熵联系起来理解.

算法描述如下：

代码实现

%%% 样本熵计算函数 %%% 
function SampEn = Sample_Entropy( dim, r, data, tau )
% SAMPEN Sample Entropy
%   calculates the sample entropy of a given time series data

%   SampEn is conceptually similar to approximate entropy (ApEn), but has
%   following differences:
%       1) SampEn does not count self-matching. The possible trouble of
%       having log(0) is avoided by taking logarithm at the latest step.
%       2) SampEn does not depend on the datasize as much as ApEn does. The
%       comparison is shown in the graph that is uploaded.

%   dim     : embedded dimension
%   r       : tolerance (typically 0.2 * std)
%   data    : time-series data
%   tau     : delay time for downsampling (user can omit this, in which case
%             the default value is 1)
%

if nargin < 4, tau = 1; end
if tau > 1, data = downsample(data, tau); end

N = length(data);
result = zeros(1,2);

for m = dim:dim+1% 该循环用于实现算法的第六步
    Bi = zeros(1,N-m+1);
    dataMat = zeros(m,N-m+1);
    
    % setting up data matrix
    for i = 1:m
        dataMat(i,:) = data(i:N-m+i);
    end
    
    % counting similar patterns using distance calculation
    for j = 1:N-m+1
        % calculate Chebyshev distance, excluding self-matching case
        dist = max(abs(dataMat - repmat(dataMat(:,j),1,N-m+1)));
        % calculate Heaviside function of the distance
        % User can change it to any other function
        % for modified sample entropy (mSampEn) calculation
        D = (dist <= r);
        % excluding self-matching case
        Bi(j) = (sum(D)-1)/(N-m);
    end
    
    % summing over the counts
    result(m-dim+1) = sum(Bi)/(N-m+1);
    
end

SampEn = -log(result(2)/result(1));

end