fMRI功能脑影像数据中的信号复杂度分析

fMRI功能脑影像数据中的信号复杂度分析人脑是目前所知的宇宙中最复杂和最神奇的结构，脑的高级功能是自然界中最复杂的运动形式。脑的基本构成单位是神经细胞（神经元）和神经胶质细胞。人脑的神经元数约达1011个。大脑分为左右两个半球，结构上具有一定的对称性，而功能上各有一定的机能分工。大脑两半球主要由灰质表层、白质和基底神经节组成。大脑皮质（亦称皮层）约由140亿个神经元组成。神经元之间通过突触（包括化学突触和电突触）来传递信息。每一个神经元通常拥有几百个以至几千个突触联结，人脑的全部突触数约达10大脑是一个非常复杂的非线性系统，早期通过解剖和研究脑局部损伤患者来获取脑的结构和功能知识，但这样获得的知识很有限。自从有了医学影像技术后就可以无损伤地测量人脑内部结构。如计算机断层扫描技术（ComputerTomography，CT）、磁共振成像技术（MagneticResonanceImaging，MRI）、弥散张量成像（DiffusionTensorImaging，DTI）都能给出的人脑的结构图像。脑功能成像不同于脑结构成像，它能显示脑内发生功能变化的区域及其时空特征。脑电图（Electroencephalograph，EEG）和脑磁图（Magnetoencephalograph，MEG）记录了脑活动时头皮表面的电信号和磁信号，能够给出脑功能快速活动时的信息。在磁共振成像技术基础上发展起来的功能磁共振成像技术（functional-MRI，fMRI）基于血氧水平依赖性（BloodOxygenLevelDependent，BOLD）效应，通过测量脑活动时脑内各处血流中含氧量的变化，来反映相应的神经细胞活动的变化。通过记录空间上分离的不同脑区的脑电、脑磁等信号，并统计它们在时间上的相关性或因果性等关系，人们可以有效地描绘不同脑区之间的协同工作模式。由于大脑的复杂性，与大脑有关的生物医学信号也表现出高度的非线性特征。随着非线性科学和计算机技术的发展，越来越多的非线性方法被广泛应用于复杂信号的分析之中。在大脑活动的研究中，近几年兴起的是对脑成像数据进行复杂度分析。在非线性动力学发展的过程中，研究者们提出了许多具有不同侧重点和应用范围的复杂度指数，比如，近似熵、排列组合熵、Hurst指数（HE）、分形维数……以上概念中，既有时域熵，又有频域熵，反映出复杂度分析本身的复杂性和多角度，应用它们可以得到关于信号的许多性质。本文的工作基于大脑fMRI数据，通过计算熵对其进行单通道时间序列复杂性分析的研究，借此来探索大脑临界行为的一些特点。功能磁共振（fMRI）长久以来，科学家们就注意到这样的事实：即人脑的功能反映在大脑皮是按空间分区的，在脑次级结构也是按空间分隔的。研究脑功能映射有许多成功的模式。如脑电图和脑磁图可检测大脑对诱发刺激响应的电或磁信号，但很难对活动区作准确的空间定位。在众多的模式中，用于脑功能定位的功能磁共振成像是一种非常有效的研究脑功能的非介入技术，已经成为最广泛使用的脑功能研究手段。它虽然是一种非介入技术，但却能对特定的大脑活动的皮层区进行准确、可靠的定位，空间分辨率达到1mm，并且能以各种方式对物体反复进行扫描。fMRI技术基于血液中氧合血红蛋白和脱氧血红蛋白磁化率不同导致的血液局部磁敏感性变化，和血氧水平依赖效应。大脑的活动是快速的神经元生理和生化变化，是消耗大量能量的过程。脑活动的增加将伴随脑局部血流、血流容量、血氧消耗的增加，而血流量增加明显高于耗氧量的增加。这将使局部氧合血红蛋白增多，脱氧血红蛋白相对减少，脱氧血红蛋白具有较高的顺磁性，主要被限制在细胞膜内，导致局部磁场环境高度不均匀。fMRI信号采集使用高场强结合高梯度场及快速梯度切换率的MR设备，选择对磁化率变化最敏感的扫描序列，常用GRE-EPI序列。可在几分钟内完成一个fMRI试验，并获得较高的空间分辨率。fMRI信号是复杂的多通道、多维数的非线性信号，因此对其复杂度的描述就不能用单一的指标。在本文所述的工作中我们分别计算了样本熵、近似熵、模糊熵、排列组合熵和赫斯特指数。fMRI数据本文计算所采用的是在脑区水平的健康被试的fMRI数据，样本容量301。在246分区模式下，每个脑区产生一个时间序列，每个序列1200个数据点。为减小数据方差和均值对计算结果的影响，我们对每个时间序列进行标准化处理，即去除均值和方差缩放，将数据按期属性（按列进行）减去其均值，并除以其方差。图1.fMRI数据图2.以节点标注的人脑246个分区复杂度指标近似熵近似熵(ApproximateEntropy，ApEn)是由Pincus于1991年提出的一种度量序列复杂性和统计量化的非线性动力学参数。它用一个非负数来表示一个时间序列的复杂性，反映时间序列中新信息发生率，越复杂的时间序列对应的近似熵越大。设一维时间序列为X(i)，i=1，…，N，通过延迟重构得到m和m+1维向量，当m维相比较时，模式序列相接近，那么当m+1维相比较时，仍然保持接近的可能性。保持接近的可能性越大，则ApEn的值越小。重复序列的ApEn值很小，而明显随机行为的ApEn值就很大。近似熵的定义与算法如下。设原始数据为x1按序号连续顺序组成一组m维矢量：Xi=1~N−m+1定义Xi与Xj之间的距离d此时Xi，X(j)中其他对应元素间差值自然都小于d。对每一个i计算Xi与其余矢量Xj(j=1~N−m+1)给定阈值r，对每一个i值统计d[X距离总数N-m的比值，记作CiCi=先将Cimrϕ再把维数加1，变成m+1，重复步骤（1）~（4），得到Cim理论上此序列的近似熵为：ApEn(m,r)=一般言之，此极限值以概率1存在。当N为有限值时按上述步骤得出的是序列长度为N时ApEn的估计值，记作ApEn(m,r)=实践上，一般取m=2,r=0.1~0.2SDxSDx是原始数据x(i),i=1~N的标准差。近似熵只要比较短的数据就能得到比较稳健的估计值，所需数据点大致1000样本熵样本熵(SampleEnropy，SampEn)是由Richman和Moornan提出的一种时间序列复杂性测度方法。样本熵在算法上相对于近似熵算法做了改进：相对于近似熵而言，样本熵计算的则是和的对数。样本熵旨在降低近似熵的误差，与已知的随机部分有更加紧密的一致性。样本熵与近似熵相比优点在于：第一样本熵不包含自身数据段的比较，它是条件概率的负平均自然对数的精确值，因此样本熵的计算不依赖数据长度；第二样本熵具有更好的一致性；第三样本熵对于丢失数据不敏感，即使数据丢失多达到三分之一，对计算值影响依然很小。样本熵是CP严格的自然对数，可用SampEn(m，r，N)来表示。其中，N为长度，r为相似容限，维数为m及m+1。样本熵旨在降低近似熵的误差，与已知的随机部分有更加紧密的一致性。样本熵的定义与算法如下。一般地，对于由N个数据组成的时间序列x（1）按序号连续顺序组成一组m维矢量：Xi=1~N−m+1（2）定义Xi与Xj之间的距离d此时Xi，X(j)中其他对应元素间差值自然都小于d。对每一个i计算Xi与其余矢量Xj(j=1~N−m+1)（3）给定阈值r，对每一个i值统计d[Xi,Xj]于1≤i定义Bi=（4）定义BB（5）再把维数加1，变成m+1，重复步骤（1）~（4），得到Ai,（6）这样，Bm(r)是两个序列在相似容限下匹配m个点的概率，而两个序列匹配m+1个点的概率，此序列的样本熵为：Sam当N为有限值时，可以用下式估计：Sam模糊熵模糊熵（Fuzzyentropy）在计算向量之间相似度时采用了模糊隶属度函数取代了近似熵算法中的二值函数，使得熵值连续平稳。同时，该算法在继承了近似熵优点的基础上，还具有对于不同的参数结果稳定，抗噪性强等特点，比近似熵更适合于度量时间序列的复杂度分析。模糊熵的定义与算法如下。对于一个M点采样序列定义：u按照序号的连续顺序重构生成一组n维矢量，X其中uj,uj+1u定义两个n维矢量Xin和Xj的一个，即（4）用模糊函数μ(dijn,m,r)定义两个矢量XD上式子中函数μdijn（5）定义函数O重复步骤（2）~（5），按照序号顺序重构一组（n+1）维矢量，函数定义如下O模糊熵定义为FuzzyEn(当M值为有限值时，按照上述七步得出的序号长度为M时模糊熵的估计FuzzyEn排列组合熵2002年Bandt等人提出排列熵（Permutationentropy）方法来检测随机信号的突变。作为一种基于复杂性量度的非线性动力学参数，其特点是基于相邻值的比较而无需考虑值的大小，因此计算简单而快速，抗噪能力强，只需要较短的序列长度就能估计出较稳定的系统特征值。可以较好地检测出复杂系统的动力学突变。根据熵值计算方法的不同，排列组合熵有Shannonpermutationentropy（SpE）、Tsallispermutationentropy（TpE）和Renyipermutationentropy（RpE）等，本文采用Shannonpermutationentropy计算排列组合熵，定义与算法如下。（1）设有长度为N的时间序列u1，u2，u3（2）通过将原序列进行重构，将每个子序列Xi表示，其中X（3）然后对每个Xi内部进行递增排序，也就是ui+j1−1L≤ui（4）通过上面的步骤，就将连续的m维子空间用一个符号序列表示了，其中这些符号的个数有m!.将所有符号的概率分布用P1,（5）计算ShannonEntropy，则时间序列u1，u2，u3H当Pj=1/m!，也就是每种符号都有且他们的概率都相等，此时时间序列的复杂度最高，排列熵最大，是lnm!0赫斯特指数1.赫斯特指数（Hurstexponet），以英国水文学家哈罗德·赫斯特命名，起初被用来分析水库与河流之间的进出流量，后来被广泛用于分形分析。被用于度量一段时间序列的长期记忆性和分形规模。赫斯特指数有三种形式：（1）如果H=0.5，表明时间序列可以用随机游走来描述；（2）如果0.5<H<1，表明时间序列存在长期记忆性；（3）如果0≤H<0.5，表明粉红噪声(反持续性)即均值回复过程。也就是说，只要H≠0.5，就可以用有偏的布朗运动(分形布朗运动)来描述该时间序列数据。2.计算方法包括聚合方差法(AggregatedVariancemethod)，R/S分析法(R/Smethod)，周期图法(Periodogrammethod)，绝对值法(AbsoluteValuemethod)，残差方差法(Varianceofresiduals)，小波分析法(Abry-Veitchmethod)，Whittle法(Whittleestimator)。我们使用R/S分析法计算。（1）将时间序列分割成不同的片段。比如取沪深三百日回报率(dailyreturn)一百天内的数据作为时间序列，将其按照以下六种规格分割：a.单个片段大小是整个序列，分成1组；

b.单个片段大小是序列的1/2，分成2组；

c.单个片段大小是序列的1/4，分成4组；

d.单个片段大小是序列的1/8，分成8组；

e.单个片段大小是序列的1/16，分成16组；

f.单个片段大小是序列的1/32，分成32组；（2）计算每个片段的均值(meanvalue)，按照上一步的六种分法总共要计算1+2+4+8+16+32=63个mean值：Mean其中Xi是片段中的元素（3）针对每个片段计算离差序列(deviation)（63个离差序列）:Yi=X其中，Y=每个片段计算出的新的离差序列，X=片段中的元素，M=第2步中计算的每个片段的均值。（4）计算每个离差序列的最大差距(widestdifference)，得到63个最大差距RR其中，R=片段中的元素；Y=离差序列。（5）计算每个片段的标准差(standarddeviation):σ=i=1nX（6）计算每个片段的R/S值：rescaledrange=RSt=R其中，R=第4步中计算的每个片段的最大距离，σ=每个片段的标准差。（7）对每种分割方法，将其各个片段的R/S值求平均得到AverageR/S(ARS)；ARS=1t2−t（8）计算Hurst指数a.将每种分段方法的片段大小（size）和ARS对10取对数

b.这样我们就有了6组对数序列。将lg(ARS)作为被变量Y，lg(Size)作为变量X，线性回归估计斜率H，H就是Hurst指数。参数选取4.1根据文献调研，我们选择了在计算fMRI数据复杂度研究中常用的推荐参数值来启动计算。模糊熵中，取嵌入维数m=2，相似容限r=0.3（由于数据进行了标准化处理，所以标准差为1，下同）；Shannon排列组合熵中，取嵌入维数m=2，延迟时间t=2.4.2近似熵和样本熵的参数调整对于近似熵和样本熵，从通行的定义中我们发现，两者本质上是相同的，样本熵是近似熵的改进。但在实际计算中发现样本熵表现出了更好的稳定性，近似熵则对相似容限r有较高的敏感性，在某些r取值下，全脑近似熵和样本熵甚至会出现完全相反的趋势。图3.Apen和Saen的相反趋势试验发现，在r=0.4附近，Apen会发生这种趋势上的反转，同时，作出全脑平均熵值随r的变化关系。图4.不同r下Apen和Saen的全脑平均值SEQ图4._不同r下Apen和Saen的全脑平均值\*ARABIC1由上图可以看出，随着r的增大，Apen和Saen在全脑的平均值趋于一致。为使二者在趋势上保持一致，同时保留部分数值上的差异，我们选择r=0.6、m=2进行计算。在某个被试下的结果如图5所示。图5.r=0.6，m=2.同一被试246脑区的Apen与Saen大脑同步特性以及临界行为在大脑信号研究中，人们越来越多地发现同步现象（synchronizationphenomena）是不同脑区信息交流的关键特征（Gray等，1989）。大脑信息整合和处理的主要潜在机制可能是由各神经网络的同步振荡引起的。如脑电中的同步现象，尤其是脑电中高频段的同步振荡，是建立起不同脑区间信息交流的一个关键部分。脑的同步活动是脑系统功能区域整合或绑定的潜在机制。从解剖学角度来看，大脑皮层可被划分为若干区域。大量动物和人类的实验表明，小到神经元，大到脑的某个区域内部和不同区域之间，都存在同步现象。认知活动需要整合广泛分布在大脑各区域的大量神经元细胞来实现，这种大范围整合是通过调节振荡神经元在有限时间段内进入特定频段上的精确锁定完成的。人在感知和识别物体时，相关脑区会自动发生同步化神经活动。脑电同步分析方法可以分为线性和非线性分析方法。线性分析方法包括时域分析、时域或者频域的参数化分析以及时域或者频域的非参数化分析；非线性方法包括互信息熵、似然同步、广义同步以及相位同步实际上生物系统或信号间的同步常常没有幅度同步关系，而仅有相位关系，特别是大尺度的跨脑区的同步。在实际中应用最多的是基于希尔伯特变换的相位同步，通过相位同步指数（Synchronization）来表征两个时间序列同步性强弱。在同步指数的基础上，我们还可以计算同步指数的香农熵，称为同步熵（Synchronousentropy，SynE）。相关研究发现，在中高程度同步下，大脑具有最大的同步熵，这表明正常大脑工作在一种临界状态下。图6展示了这方面研究的一种结果。某些同步指数的对应被试在整个样本中具有最大的全脑平均同步熵，也就是具有最大的全脑平均复杂度。提取同步熵极值附近的被试进一步研究他们的大脑行为，发现这些被试大脑的某些行为特征满足临界现象的幂律特性。在这个样本中同步熵随同步指数的分布形态表明在普遍情况下，大脑应该在某种临界状态下工作。本文的工作计算了多种非线性信号在全脑水平和脑区水平下的熵值，比较了在不同尺度上与同步指数之间的变化关系，希望从更多的复杂度指标中寻找大脑临界特性的证据。进一步探究大脑的临界和同步行为有助于我们加深对大脑工作机制的理解，对神经网络的构建、模拟，人工智能的发展，脑疾病的诊断、治疗和预防都有重要的指导意义。图6.Syn-SynE在syn=0.5附近出现synE的极大值Hurst指数我们首先研究Hurst指数随同步指数的变化行为。如前所述，Hurst指数可以用来度量时间序列有无长程记忆，其值介于0到1之间。人脑处于临界态时，脑神经活动具有长程时间自相关性，所以近年来许多fMRI研究利用赫斯特指数来衡量信号的自相关程度并进一步作为临界状态的判据。赫斯特指数越大，说明时间序列的长程自相关程度越高，系统更接近临界点。我们首先计算了每个脑区在301个被试下赫斯特指数的平均值图7.不同脑区赫斯特指数在所有被试下的平均值以及误差同时，我们借助同步熵来表征大脑功能的同步特性。将所有被试的Hurst指数与同步熵做线性相关分析，发现两者之间存在显著的正相关。由于神经活动在临界时趋于有更高程度的长程自相关，即时间序列有更大的赫斯特指数，这种正相关就表示神经系统处于临界态时趋向于拥有更高程度的同步特性。图8.syn-Hurst随着同步指数的增加观察到Hurst指数的上升趋势其他复杂度指标的结果7.1全脑平均水平在赫斯特指数之外，我们又计算了Apen、Saen、Fuzzyen、SpE，以此为依据来判断大脑的临界状态。图9.不同脑区Apen、Saen、fuzzyen和Peen在所有被试下的平均值和误差横轴标注出了脑区编号对应的大脑解剖学部位，分别是额叶（Frontallobe）、颞叶（Temporallobe）、顶叶（Parietallobe）、脑岛（Insula）、扣带回（Cingulategyrus）、枕叶（Occipitallobe）、和皮层下中枢（Subcorticalcenter）。这种分区将在下面的讨论中继续使用，具体内容也将在下文提到。从计算结果可以得知，虽然复杂度指标不同，但是熵值在脑区上的分布呈现出相似的形态，展现出分区聚集的趋势。这令我们很自然地同大脑的功能分区联系起来。下面我们会看到，熵值的这种分布恰好对应大脑的某种功能分区或者生理构造。相关研究发现，Hurst指数在超过0.5时，与同一过程的排列组合熵呈现负相关。图10.Hurst和Peen的相关行为上图引自A.B.Pessa1,∗andHaroldoV.Ribeiro.Characterizingstochastictimeserieswithordinalnetworks.physics.data-an3Oct2019通过作SpE和Apen、Saen、Fuzzyen的相关分析可知，它们之间存在显著的正相关，那么，在Syn和Hurst指数呈现正相关的情形下，Syn应和其他复杂度指标呈现负相关，即在中高程度同步下，其他复杂度指标出现数值的极小，或者随syn呈现下降趋势。图11.不同被试的Apen、Saen、fuzzyen、Peen在所有脑区下的平均值与syn的关系可以看出各熵随syn呈现下降趋势7.2分脑区水平为了更好地显示脑区水平的熵值与同步指数的关系，研究大脑不同部位在临界状态的特点，我们将大脑246个分区划分为七个模块，即额叶、颞叶、顶叶。脑岛、扣带回、枕叶和皮层下中枢。其中，额叶（Frontallobe）（脑区编号1-68）包括额上回（1-14）、额中回（15-28）、额下回（29-40）、眶回（41-50）、中央前回（53-64）和中央旁小叶（65-68）；颞叶（Temporallobe）（69-124）包括颞上回（69-80）、颞中回（81-88）、颞下回（89-102）、梭状回（103-108）、海马旁回（109-120）和颞后上沟（121-124）；顶叶（Parietallobe）（125-162）包括顶叶上小叶（125-134）、顶叶下小叶（135-146）、楔前叶（147-154）和中央后回（155-162）；脑岛（Insula）分区编号为163-174；扣带回（Cingulategyrus）分区编号为175-188；枕叶（Occipitallobe）（189-210）包括枕部中央皮质（189-198）和枕外侧皮质（199-210）；皮层下中枢（Subcorticalcenter）（211-246）包括杏仁核（211-214）、海马体（215-218）、基底神经节（219-230）和丘脑（231-246）.计算每个被试在不同脑区上的熵的平均值，同样关注其与同步指数之间的关系，可以看出，即使是在脑区水平下，也是在中高程度同步下具有熵的极值，也就是说大脑的不同区域的工作机制具有一定的共性，但不同脑区间也存在一定的差异。图12.不同被试的Apen、Saen、fuzzyen、Peen在不同脑区下的平均值与syn的关系除了Hurst指数之外，熵在不同脑区下的平均值同样随syn呈现下降趋势，不同脑区的熵随syn的变化趋势存在一定程度差异7.3熵值的脑地形图为了更直观地展现不同的熵在大脑不同区域的分布情况，我们借助BrainNetViewer软件、使用颜色梯度在整个大脑上对熵作了标注。图13.Apen图14.Saen图15.fuzzyen图16.Peen图17.Hurst可以直观地看处，在前额叶、颞叶、扣带回和脑岛，近似熵、样本熵、模糊熵和排列组合熵数值较大，在枕叶数值较小；而赫斯特指数则正好相反。参考文献杨福生.近似熵_一种适用于短数据的复杂性度量陈慧.近似熵与样本熵在静息态fMRI数据中的对比研究张艳艳.近似熵理论及应用南姣芬.基于近似熵静息态fMRI的大脑功能性别差异研究吕杰.基于样本熵的生物医学信号复杂度分析方法及其优化算法刘慧.生物时间序列的近似熵和样本熵方法比较相洁.基于模糊熵的时间序列非线性检测方法孙克辉.模糊熵算法在混沌序列复杂度分析中的应用武薇.基于排列熵理论的非线性系统特征提取研究ChristophBandtandBerndPompe，PermutationentropyAnaturalcomplexitymeasurefortimeseries.physicsreviewletters88（17）29April2002A.B.Pessa1,∗andHaroldoV.Ribeiro.Characterizingstochastictimeserieswithordinalnetworks.physics.data-an3Oct2019赵晓军.时间序列的相关性及复杂性研究李凌.利用相位同步聚类方法研究大脑功能连通性王娜.脑电信号相位同步特性研究屈景怡.大脑皮层神经网络模型同步问题研究陈清坚.静息态下大脑临界行为及动力学研究

附录计算近似熵clc;Apen2=[];files=dir('E:\HCP_ROI_SIGNAL\ROI_signal');mn=size(files,1);fork=1:mn-2cd('E:\HCP_ROI_SIGNAL\ROI_signal');Q=k+2;str=files(Q).name;cd(str);load('ROI_246_RS.mat');data=ROI_246_RS;[H,L]=size(data);fora=1:Ldata=ROI_246_RS(:,a);%每一列算一个熵值sigma=std(data);data=(data-mean(data))./sigma;%数据标准化N=length(data);result=zeros(1,2);dim=2;%近似熵的嵌入维数一般取为2r=0.6;%r=0.2sigma~0.35sigma,不要小于0.2sigma，数据标准化后标准差为1forj=1:2m=dim+j+1;phi=zeros(1,N-m+1);dataMat=zeros(m,N-m+1);fori=1:mdataMat(i,:)=data(i:N-m+i);endfori=1:N-m+1tempMat=abs(dataMat-repmat(dataMat(:,i),1,N-m+1));boolMat=any((tempMat>r),1);phi(i)=sum(~boolMat)/(N-m);endresult(j)=sum(log(phi))/(N-m+1);%算上了和自身的距离，应是sum(log(phi))-1%endApen2(k,a)=result(1)-result(2);endend计算样本熵clc;Saen2=[];files=dir('E:\HCP_ROI_SIGNAL\ROI_signal');mn=size(files,1);fork=1:mn-2cd('E:\HCP_ROI_SIGNAL\ROI_signal');Q=k+2;str=files(Q).name;cd(str);load('ROI_246_RS.mat');data=ROI_246_RS;[H,L]=size(data);fora=1:Ldata=ROI_246_RS(:,a);%每一列算一个熵值sigma=std(data);data=(data-mean(data))./sigma;%数据标准化N=length(data);result=zeros(1,2);dim=2;%样本熵的嵌入维数一般取为2r=0.6;%r=0.2sigma~0.35sigma,不要小于0.2sigma，可取0.2.*sigma~0.35.*sigma，数据标准化后标准差为1form=dim:dim+1;Bi=zeros(1,N-m+1);dataMat=zeros(m,N-m+1);fori=1:mdataMat(i,:)=data(i:N-m+i);endfori=1:N-m+1dist=max(abs(dataMat-repmat(dataMat(:,i),1,N-m+1)));D=(dist<r);Bi(i)=(sum(D)-1)/(N-m);%算上了和自身的距离，应是sum(D)-1endresult(m-dim+1)=sum(Bi)/(N-m+1);endSaen2(k,a)=-log(result(2)/result(1));endend计算模糊熵%data：时间序列数据%dim:嵌入维数，可取为2或3;%r：相似容限0.3.*std~0.35.*std%tau：下采样延迟时间，一般取为1clc;fuzzyen=[];files=dir('E:\HCP_ROI_SIGNAL\ROI_signal');mn=size(files,1);fork=1:mn-2cd('E:\HCP_ROI_SIGNAL\ROI_signal');Q=k+2;str=files(Q).name;cd(str);load('ROI_246_RS.mat');data=ROI_246_RS;[H,L]=size(data);fora=1:Ldata=ROI_246_RS(:,a);%每一列算一个熵值sigma=std(data);data=(data-mean(data))./sigma;%数据标准化N=length(data);Phi=zeros(1,2);dim=2;%近似熵的嵌入维数一般取为2r=0.3;%r=0.2sigma~0.35sigma,不要小于0.2sigma，数据标准化后标准差为1tau=1;form=dim:dim+1Ci=zeros(1,N-m+1);dataMat=zeros(m,N-m+1);forj=1:mdataMat(j,:)=data(j:N-m+j);endU0=mean(dataMat);Sm=abs(dataMat-repmat(U0,m,1));fori=1:N-m+1Sm_tmp=Sm;Sm_tmp(:,i)=[];dij=max(repmat(Sm(:,i),1,N-m)-Sm_tmp);Aij=exp(-log(2).*(dij/r).^2);Ci(i)=sum(Aij)/(N-m);endPhi(m-dim+1)=sum(Ci)/(N-m+1);endfuzzyen(k,a)=log(Phi(1))-log(Phi(2));%fuzzyen=ln(φ_m)-ln(φ_m+1);endend计算排列组合熵clc;Peen=[];files=dir('E:\HCP_ROI_SIGNAL\ROI_signal');mn=size(files,1);fork=1:mn-2cd('E:\HCP_ROI_SIGNAL\ROI_signal');Q=k+2;str=files(Q).name;cd(str);load('ROI_246_RS.mat');data=ROI_246_RS;[H,L]=size(data);forb=1:Ldata=ROI_246_RS(:,b);%每一列算一个熵值sigma=std(data);data=(data-mean(data))./sigma;%数据标准化m=5;%m为嵌入维数，应如何判断其合适取值？t=2;%t为延迟时间，应如何判断其合适取值？y=data;ly=length(y);permlist=perms(1:m);c(1:length(permlist))=0;forj=1:ly-t*(m-1)[a,iv]=sort(y(j:t:j+t*(m-1)));%即y(j:t:j+t*(m-1)，1),j到j+t*(m-1)，步长为tforjj=1:length(permlist)if(abs(permlist(jj,:)-iv'))