多尺度熵---UnderstandingMultiscaleEntropy
⼀、引⾔
与其他熵测量⽅法⼀样，多尺度熵的⽬标是评估时间序列的复杂性。使⽤多尺度熵的主要原因之⼀是不知道时间序列中相关的时间尺度。例如，在分析语⾳信号时，在单词时间尺度下统计信号的复杂度会⽐统计整个语⾳⽚段的复杂度更加有效。但如果你不知道⾳频信号代表语⾳，甚⾄对语⾳概念没有任何了解，你就不知道应该运⽤什么时间尺度以从原始信号中获得更多有⽤的信息。因此，通过多个时间尺度来分析问题将会得到更多信息。在脑电图中，潜在的脑电模式是未知的，因此相关的时间尺度也是未知的。所以，需要通过多尺度样本熵来分析哪个尺度对特定场合下脑电信号的分析更有⽤.
⼆、计算多尺度熵
多尺度熵的基本原理包括对时间序列进⾏粗粒化或下采样 - 主要是在越来越粗略的时间分辨率下分析时间序列。具体操作如下：
假设有以采样频率为1kHz采样得到的时间序列(原始时间尺度T为1 ms)
$$x_1，x_2，x_3，...... x_N$$
缘分有take2
粗粒化（ Coar Graining）数据意味着对不同数量的连续点取平均值，以创建不同尺度（或分辨率）的信号.
1. 当scale为1时，粗粒化数据的结果就是原始时间序列.
2. 当scale为2时，粗粒化后的时间序列是通过计算两个连续时间点的平均值形成的，如下图(A)所⽰。即定义y1 = (x1 + x2)/2; y2 = (x3 + x4)/2，以此类推.
3. 当scale为3时，粗粒化时间序列为三个连续时间点的平均值构成，如下图(B)所⽰。即定义y1 = (x1 + x2 + x3)/3; y2 = (x4 + x5 +
x6)/3，依此类推.
图⽰来⾃⽂献[1].
粗粒化过程分为两种形式：⼀是⾮重叠式，每次跳跃τ个数据，取τ个数据做平均以产⽣新的数据；⼆是重叠式，每次跳跃1-τ个数据，取τ个数据做平均。
上述粗粒化过程为⾮重叠式，其数学定义为
其中，τ表⽰时间尺度.
然后，计算与每个尺度或分辨率对应的样本熵，并绘制样本熵－尺度曲线图。曲线下的⾯积是所选尺度范围内样本熵值之和，是多尺度熵的度量。⼀个波动较⼤的时间序列会产⽣较⼤的熵值，因此可以认为是具有较⾼复杂度的信号。同样，具有⾼度规律性的信号，其熵值也较低。
三、代码实现
备注：⽰例程序来⾃Björn Herrmann的个⼈⽹站上有关多尺度熵的介绍，对四种不同噪声信号的多尺度熵进⾏了详细分析，说明了采⽤传统的样本熵进⾏分析可能存在的问题，. 在此对Doctor Björn所做贡献表⽰感谢！
考研政治教材（为节省篇幅，此处仅给出了主程序，其中涉及的⼦函数请参见⽂末附录，或直接下载源代码.）
% matlab 主程序
awgclear all;
rand('ed',10)
%% generate signals
Sf = 1000; % Sampling frequency
dur = 30; % Time duration
y = colored_noi(Sf,dur,0); % white noi
%% Plot time cour
t = (0:length(y)-1)/Sf; % time vector
figure, t(gcf,'Color',[1 1 1]), hold on
t(gca,'FontSize',12)
plot(t,y,'-','Color',[0,0,0])
ylim([-3 3])
xlabel('Time (s)')
ylabel('Amp. (a.u.)')
title(['White' ' noi (1/f^{' num2str(0) '})'])
%% Plot spectra
[F xf] = spectra(y, Sf, Sf*10, @hann);
figure, t(gcf,'Color',[1 1 1]), hold on
t(gca,'FontSize',12)
plot(xf,abs(F),'-','Color',[0,0,0])
% xlim([0 400])
ylim([-0.01 0.06])
xlabel('Frequency (Hz)')
ylabel('Amplitude (a.u.)')
title(['White ' 'noi (1/f^{' num2str(0) '})'])
interestrate
%% Compute MultiScale Entropy
% [m sf] = MSE_Costa2005(x,nSf,m,r)
%
% x - input signal vector (e.g., EEG signal or sound signal)
% nSf - number of scale factors
% m - template length (epoch length); Costa ud m = 2 throughout
% r - matching threshold; typically chon to be between 10% and 20% of % the sample deviation of the time ries; when x is z-transformed:
% defined the tolerance as r times the standard deviation
%
% m - multi-scale entropy
% sf - scale factor corresponding to m
[m,sf] = MSE_Costa2005(y,20,2,std(y)*0.15);
%% Plot m
figure, t(gcf,'Color',[1 1 1]), hold on, pp = []; labfull = {};主要看气质是什么意思
t(gca,'FontSize',12)
pp = plot(sf,m,'-','Color',[0,0,0],'LineWidth',2);
labfull = ['White' ' noi (1/f^{' num2str(0) '})'];
ylim([0 3])
xlabel('Scale')
ylabel('SampEn')
title('Multi-scale entropy')
legend(pp,labfull)
legend('boxoff')
四、多尺度熵在脑电分析中的应⽤
由于MSE计算不同时间尺度下信号的熵，因此对于理解像EEG等⽣物信号在不同时间尺度下的复杂度变化⼗分有⽤. 熵与时间尺度的曲线可能会产⽣⼀个峰值，表明在该时间尺度下存在最⼤熵，该时间尺度可能具有更⼤的相关性。
事实上，Escudero等⼈[2]的研究表明，即使在10个电极的⼤时间尺度上，MSE也能发现阿尔茨海默症患者和对照组之间的显著差异，⽽且与对照组相⽐，阿尔茨海默症患者的脑电图活动没有那么复杂.
在Partk等⼈[3]的另⼀项研究中，他们对正常受试者（Normal）、阿尔茨海默病患者（Severe AD）和轻度认知障碍(MCI)患者的脑电信号进⾏了MSE分析。结果再次显⽰，阿尔茨海默病(AD)患者脑电信号的复杂性降低，其熵值较低，如图所⽰。正常和MCI被试在尺度为6、7时具有局部最⼤熵值，之后熵值下降。他们的分析还表明，在时间尺度因⼦为1时，轻度认知障碍、阿尔茨海默病和正常受试者的样本熵在统计学上是⽆法区分的，因此，在这种情况下，样本熵或近似熵等熵指标⽆法区分健康受试者和病理受试者.
在另⼀项研究中，Catarino等⼈[4]对健康受试者和⾃闭症谱系障碍(ASD)患者在执⾏社交和⾮社交任务(包括⾯部和椅⼦构成的视觉刺激)时的脑电信号进⾏了MSE分析。他们的结果显⽰，在枕叶和顶叶区域，与对照组相⽐，⾃闭症组的脑电图复杂性显著下降(p值统计结果).
图⽚来⾃⽂献[4]
因此，尽管以⾮特异性的⽅式，各种疾病状态都体现为MSE测量值下降。
MSE中的参数
MSE只是简单地将样本熵度量扩展到不同的时间尺度，因此，基础参数m(⽐较的线段长度)和r(两个线段之间的距离度量)是相同的。在前⾯的中讨论过，这些参数对于MSE的性能⾮常重要。关于如何选择这些参数，并没有标准的指导准则。对于数据长度m, 研究表明在MSE应⽤中，需要保证每个时间尺度下有⾜够的数据量。Costa等⼈[5]的研究表明，当⽩噪声和 1/f 噪声的数据点减少时，样本熵的均值(超过30个尺度)会发散。特别是在 1/f 噪声情况下，由于其⾮平稳性，与⽩噪声相⽐，发散速度更快。要查看参数m，r和数据长度对样本熵的影响，请参阅我们之前的以及Costa等⼈关于这些问题的另⼀篇优秀⽂章[5].
关于r，需要记住的是，为了避免噪声对样本熵估计的显著贡献，r必须⼤于⼤部分的信号噪声。选择r的另⼀个标准是基于信号的动态特性（signal dynamics）。
然⽽，最重要的⽅⾯是计算r的⽅式，以及所选择的距离测量（通常是欧⼏⾥德距离）是否真正与时间序列相关。例如，对于语⾳信号，欧⼏⾥德距离可能不是两个单词之间距离的最精确度量。脑电图完全有可能也是这种情况。
参考⽂献
[1] Busa, Michael A., and Richard EA van Emmerik. “Multiscale entropy: A tool for understanding the complexity of postural control.” Journal of Sport and Health Science 5.1 (2016): 44-51.
[2] Escudero, J., et al. “Analysis of electroencephalograms in Alzheimer’s dia patients with multiscale
entropy.” Physiological measurement 27.11 (2006): 1091.
[3] Park, Jeong-Hyeon, et al. “Multiscale entropy analysis of EEG from patients under different pathological
conditions.” Fractals15.04 (2007): 399-404.
[4] Catarino, Ana, et al. “Atypical EEG complexity in autism spectrum conditions: a multiscale entropy analysis.” Clinical neurophysiology 122.12 (2011): 2375-2383.
[5] Costa, Madalena, Ary L. Goldberger, and C-K. Peng. “Multiscale entropy analysis of biological signals.” Physical review
E 71.2 (2005): 021906.
附录
colored_noi.m
function [y] = colored_noi(Sf,dur,b)
% [y] = colored_noi(Sf,dur,b)
%
% Sf - sampling frequency
% dur - duration
% b - exponent of the 1/f^b power law function
% b = 0 (white), 1 (pink), 2 (brown), -1 (blue), -2 (violet)
instruction是什么意思%
% y - noi time domain signal
% (y is normalized such that mean(y)==0 and std(y)==1)
%
% Description: The script calculates a colored noi.
%
%
% ---------
%
% Copyright (C) 2015
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public Licen as published by
% the Free Software Foundation, either version 3 of the Licen, or
% (at your option) any later version.
%
同事用英语怎么说
% This program is distributed in the hope that it will be uful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU General Public Licen for more details.
%
% You should have received a copy of the GNU General Public Licen
% along with this program. If not, e </licens/>.
%
% -----------------------------------------------------------------------
% B. Herrmann, Email: herrmann., 2015-08-27
% get number of samples
N = round(Sf * dur);
if mod(N,2)
M = N + 1;
el
M = N;
end
% get fourier spectrum and the corresponding frequency axis
新会计准则F = fft(rand([M 1]));
xf = (0:M/2-1)*(Sf/M);
% power law function
p = 1./xf(2:end).^b;
% get scaled magnitudes and pha
R = sqrt((abs(F(2:length(xf))).^2).*p');
P = angle(F(2:length(xf)));
% get symmetrical spectrum
F(2:length(xf)) = R .* (cos(P) + 1i.*sin(P));
F(length(xf)+2:end) = real(F(length(xf):-1:2)) - 1i*imag(F(length(xf):-1:2));
% IFFT
y = ifft(F);
% get correct number of samplesshabby什么意思
y = real(y(1:N));
% ensure unity standard deviation and zero mean value
demeter
y = y - mean(y);
yrms = sqrt(mean(y.^2));
y = y / yrms;
spectra.m
function [F xf] = spectra(X, Sf, nfft, win)
% [F xf] = spectra(X, Sf, nfft, win)
%
% Obligatory inputs:
% X - signal matrix (fft is done along the first dimension)
% Sf - sampling frequency
