CN103876736B - 一种基于功率谱划分的复杂度谱脑电信号分类识别方法 - Google Patents

Abstract

本发明提出了一种基于功率谱划分的复杂度谱的脑电信号分类识别方法，属于脑电信号的模式识别与分类领域。本发明主要包括：一种基于功率谱划分的脑电信号的复杂度谱定义、分析和提取方法和一种非线性罗辑斯蒂(Logistic)复杂度谱参考模型构建方法。首先对脑电信号定义了一个基于功率谱划分的复杂度谱，并给出其计算方法，然后利用本发明对此映射生成的数据序列进行计算，据此建立一个对脑电信号分析复杂度的谱参考模型，并分析各谱线序列的大小、数量和分布所反映的物理意义，画出该映射基于功率谱划分的复杂度谱参考空间分布模型。本发明可以对脑电信号进行模式识别和分类研究。

Description

一种基于功率谱划分的复杂度谱脑电信号分类识别方法

技术领域

本发明属于脑电信号(EEG)分类和识别方法的领域，涉及一种脑电信号处理中的分析、建模和计算方案，尤其是一种能够对脑电信号进行参考模型化非线性处理以及对脑电信号的复杂度细微结构进行分析的方法。

背景技术

目前，脑科学的研究是科学研究领域中的一个热点。其中，脑电信号的处理共有两类方法：一类是线性分析方法，另一类是非线性方法。线性方法属于传统信息分析方法；非线性方法属于现代信息处理方法。线性方法主要有时域分析、频域分析、时频分析等，非线性方法主要包括非线性动力学方法、信息熵、复杂度等。虽然线性方法在脑电信号的研究领域有着不可替代的作用，也没有证据表明它有被取代的可能，但是在研究中确实发现存在许多不可避免的缺陷。例如，脑电信号在处理中的非线性关联信息量在分析中丢失，使得其复杂性得不到充分的研究等。在此种情况下，出现了脑电信号的非线性分析方法，并且快速发展起来，已经取得许多重要成果。其中，复杂度的方法就是发展起来比较快的一种方法。

20世纪60年代，柯尔莫哥洛夫(Kolmogorov)等人把复杂性定义为：一个系统的复杂程度与该系统的空间结构或时间序列所表示的变化行为的最小描述有关，一般称作柯尔莫哥洛夫复杂性。由于柯尔莫哥洛夫复杂性不可计算，1976年A.Lempel和J.Ziv提出一种度量符号序列复杂性的简单算法，称为Lempel-Ziv复杂度。直到1987年，才由卡斯帕(Kaspar)和舒斯特(Schuster)提出了该算法的计算机实现方法。解幸幸、李舒等也在《Lempel-Ziv复杂度在非线性检测中的应用研究》一文中认为：某事物的算法复杂度等于产生该事物的图形结构或符号序列的最短程序长度与该图形结构或符号序列本身大小之比的极限。

脑电信号是一种非线性信号，复杂度能够客观有效地测度和评价了信号的复杂性和混乱程度。1991年，吴祥宝和徐京华等首次将复杂度分析引入脑电信号的研究中。

1998年陈芳、顾凡及、徐京华等提出C₀复杂度的概念和算法，其主要思想就是：把信号分解成规则成分和不规则成分两部分，C₀复杂度就定义为不规则成分在原信号里所占的比例。随后在2004年，蔡志杰、顾凡及、沈恩华对该算法作了进一步的改进，主要是改进了划分规则部分与不规则部分的方法,采用幅值平方的功率谱代替了幅度谱,然后再计算功率谱的平均值,以此作为分界线来区分规则部分与不规则部分，进而计算出复杂度值。本文即是在此工作的基础上的进一步改进和拓展。

近年来，有关复杂度的研究许多都和混沌时间序列结合在一起，出现许多新的研究成果。主要分为几方面：(1)侧重关注研究混沌时间序列复杂度的计算方法；(2)侧重关注研究混沌时间序列复杂度的强度计算；(3)关注研究混沌时间序列复杂度的稳定性分析；(4)研究关注混沌时间序列复杂度的相关性分析。所有这些研究有一个共同点就是都只关注一个单一指标的复杂度的计算和分析，没有提出对其细微结构进行进一步的分析研究，没有对多复杂度概念指标进行研究，没有提出复杂度谱的定义，分析和计算方法。虽然有人提出结构复杂度的概念，但还是没有对时间序列复杂度的组成、结构成分以及数量分布的进行研究。现有的脑电信号都具有较强的非线性，仅仅分析计算单一复杂度及其性质不能满足实际存在的类型繁多的脑电信号的模式分类、特征和性质。因此，提出一种在统一的非线性背景下从规则周期，经分叉混沌周期，再到随机序列的演化过程的复杂度谱的整个演变过程，并能够计算时间序列复杂度，计算各成分的数量、大小和分布的方法和分析其性质及物理意义，对于分析和研究脑电信号有着很重要的参考意义。

综上所述，现有方法和研究成果虽然比较多，但对于脑电信号的分类与识别研究远远不够，复杂度的研究虽然给这个研究带来较大的发展，但其单一复杂度指标现状难以满足脑电信号强非线性的性质的实际需求，尤其对多种类型的脑电信号并存的状况，给脑电信号的分类和识别带来更大的难度。

发明内容

为解决实际技术中单一复杂度指标不能满足对多种类型脑电信号的分类和识别的问题，本发明提供一种基于功率谱划分的复杂度谱的脑电信号的分类和识别方法，这种方法不仅能够对脑电信号进行分析、建模和计算，还能对脑电信号进行参考模型化非线性处理以及对脑电信号的复杂度细微结构进行分析计算。

本发明采用如下的技术方案。

一种基于功率谱划分的复杂度谱的脑电信号分类和识别方法，其特征在于包括以下步骤：

步骤一：基于功率谱划分的复杂度谱模型的定义及计算：本发明在现有技术和研究的基础上扩展和拓宽C₀复杂度的概念，提出复杂度谱的定义和计算方法，并给出了通过功率谱划分确定复杂度谱长度的方法。具体方法如下：

(1)基于功率谱划分的复杂度谱模型的定义及计算，方法如下：

1)对于给定的脑电信号序列X＝{x(n),n＝0,1,2,...,N-1}，计算其傅立叶(Flourier)变换，得到功率谱序列，记为F；

2)利用公式计算序列F在f(k)点处的功率谱值。得功率谱值序列:

P＝{p(k),k＝0,1,2,...,N-1}；

3)寻找序P的最大值和最小值：P_max＝max(P)，P_min＝min(P)，计算其差值ΔP＝P_max-P_min，将ΔP平均分成m段，其中：2≤m≤N的正整数,N为P的元素个数，并从下到上依次记为ΔP₁,ΔP₂,…,ΔP_m，记新序列为则按如下规则生成新序列：

得到新序列

4)对于序列进行傅立叶(Flourier)逆变换，得新的时间序列为：

$\widetilde{x^j}\left(n\right)=\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}\widetilde{f_N^j}\left(k\right)W_N^{\mathrm{nk}},k=\mathrm{0,1,2},...,N-1,j=\mathrm{1,2},...,m---\left(2\right)$

表示的傅立叶变换项；

5)定义如下求和比：

$C_X\left(j\right)=\frac{\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}{\left|\widetilde{x^j}\left(n\right)\right|}^2}{\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}{\left|x\left(n\right)\right|}^2},j=\mathrm{1,2},...,m---\left(3\right)$

为C_X(j)为序列X的第j复杂度，称序列s＝{C_X(j),j＝1,2,...,m}，为脑电信号时间序列X的基于功率谱划分的复杂度谱；

(2)确定基于功率谱划分的复杂度谱长度m，方法如下：

1)构造关于m的函数：

$E_{CX}\left(m\right)={\left(\underset{i=1}{\overset{m}{\Σ}}{\left(C_X\left(j\right)\right)}^2\right)}^{\frac{1}{2}}---\left(4\right)$

作为目标函数，m是整数，2≤m≤N，N为P的元素个数。

2)作出目标函数曲线，对每一个m求其线性拟合曲线：

y＝a_mm+c (5)

其中：a_m是取m时的斜率，c截距。得到一个斜率序列A＝{a₂,a₃,…,a_m}；

3)计算A中每两个相邻项之间差值的绝对值序列：

Y＝{Δa_i＝|a_i+1-a_i|，i＝1,2,3,4,…,m-1} (6)

Y表示相邻两项的差序列；

4)在Y中寻找这样的项Δa_g，使得Y中所有i≥g的项，都有Δa_i≤Δa_g，则第g项Δa_g就是序列Y开始收敛或稳定变化的起始点，从g项后所有序列项就是收敛和稳定变化的区域，此时，Y序列相对Δa_g的稳定或收敛的变化精度为：

$\mathrm{\δ}=\frac{{\mathrm{\Δa}}_g}{max\left\{Y\right\}}---\left(7\right)$

其中：0＜δ≤0.05，g即是相对于δ开始收敛点，也是相对于精度δ的最小划分数目。

步骤二：建立非线性罗辑斯蒂(Logistic)映射的复杂度谱参考演化模型：本发明利用步骤一定义的方法，对经典的非线性罗辑斯蒂映射序列的数据进行计算建立复杂度参考演化模型，主要包括计算复杂度谱和标定模型，以及复杂度计算窗口的最小长度和稳定复杂度值的确定方法，具体方法如下：

(1)建立非线性罗辑斯蒂映射的复杂度谱参考演化模型,方法如下：

1)根据罗辑斯蒂(Logistic)映射公式，生成控制参数u∈[0.75,2]的演化区间作为建立脑电信号复杂度谱从周期经混沌向随机演化的标准数据，按照步骤一的方法计算其每一时间序列的复杂度谱：

$S^i=\{C_X^i\left(j\right),j=1,2,\mathrm{...},m\},i\∈\[0.75,2\]---\left(8\right)$

Sⁱ表示i对应复杂度谱。

2)按照给定的数据步长，逐步加大计算窗口的长度并重复计算S,直到各复杂度谱中的每个复杂度值收敛或变化稳定为止，此时计算和统计稳定区间的如下过程数据：

i)计算最小数据窗口的长度：

N_min＝min{N₁,N₂,N₃,…N_i} (9)

N_i表示在收敛或变化稳定是各计算窗口长度；

ii)计算S在稳定区域的各谱线的近似平均值：

$\stackrel{\‾}{S}=\{\stackrel{\‾}{C_X^i\left(j\right)},j=1,2,\mathrm{...},m,i\∈\[0.75,2\]\}---\left(10\right)$

表示S的近似平均值序列；

iii)计算S中各复杂度值的变化范围：

表示脑电信号X在控制参数u＝i时第j复杂度在不同长度的窗口计算得到的复杂度值序列。

3)根据复杂度谱序列最后计算值，并对其校验后在三维坐标图中作出三维空间模型图，并标出变化范围的边界。作出各序列每个复杂度谱的复杂度值的窗口长度三维模型。

(2)确定最小计算窗口长度和稳定区平均复杂度值，算法如下：

1)脑电信号X在u＝i时第j复杂度在以L/30的步长(L为数据长度)来增加计算窗口长度,计算得到的复杂度值序列为：

v表示起始计算时的数据窗口长度，为1≤v≤L的正整数；d表示每次增加的计算步长，为1≤d≤L的正整数。

2)计算每相邻两项的中位值：

$t_{v+id}=\frac{1}{2}(Z_{v+id}+Z_{v+(i+1)d})---\left(13\right)$

并构成新的序列：

T＝{t_v+id，i＝0,1,2,…} (14)

T表示新的平均序列；

3)以新序列T代替重复迭代步骤二中(2)的第2)步10次，取其最后一项作为近似平均复杂度值t_end；

4)将序列的每一项Z_v+id与t_end作差值：

ΔZ_v+id＝Z_v+id-t_end (15)

取序列中ΔZ_v+id≥0项元素的下标组成新的序列H_big，取序列中ΔZ_v+id≤0项元素的下标组成新的序列H_small；

5)计算序列H_big中的最小值：

b_big＝min{H_big} (16)

计算H_small中的最小值：

b_small＝min{H_small} (17)

然后求h_big和h_small两者中的最大值,记为N₁，其值为：

N₁＝max{b_big,b_small} (18)

从而，得到生成新序列的起始点；

6)在步骤二中(2)的第1)步的复杂度谱序列中，用H表示其下标，对该序列下标H≥N₁的所有项N₁,N₁+d,N₁+2d,…,将其组成一个新序列G：

G＝{Z_N1,Z_N1+d,Z_N1+2d,...} (19)

设G有r项，计算序列G的平均值：

${\stackrel{\‾}{Z}}_1=\frac{1}{r}\underset{i=0}{\overset{r-1}{\Σ}}{Z}_{N_1+id}---\left(20\right)$

表示G的平均值；

7)用替代t_end，重复步骤二中(2)的第3)～5)步，直到计算得到的N_i值与上一次计算的N_i-1值相同为止。记最后一次计算的N_i为N_min，为

这样，最后计算得到的N_min和即为所求的最小计算窗口长度和稳定区平均复杂度值。

步骤三：计算已获取多种脑电信号的复杂度谱序列特征，在参考模型中标定，方法如下：

(1)对于已收集到的多种脑电信号，按照步骤一和步骤二的计算方法，计算N_min、和指标；

(2)并将这些指标标定在复杂度谱三维空间模型图和最小窗口三维模型图中，以此确定这几类脑电信号指标在这两个模型中的空间位置和分布范围的大小。

步骤四：对于新采集的脑电未知信号计算其特征指标，并将其与参考模型的特征空间做相似度对比，计算其和各类空间之间的距离dis(dis表示在特征空间的欧氏距离，为dis＞0的正实数),最后依据dis的大小进行排序，取最小值的dis的类作为分类识别输出。

本发明的有益效果是，本方法不仅能够对脑电信号进行分析、建模和计算，还能对脑电信号进行参考模型化非线性处理以及对脑电信号的复杂度细微结构进行分析计算。为研究者为提供一种新的脑电信号分析和计算方法。

附图说明

图1为基于功率谱划分的复杂度谱的脑电信号分类识别方法流程图；

图2为计算基于功率谱划分的脑电信号的复杂度谱流程图；

图3为确定功率谱划分m的流程图；

图4为某一脑电信号的基于功率谱划分的复杂度谱；

图5为随控制参数u的演化功率谱划分数m变化曲线；

图6为精度变化和划分数m关系曲线；

图7为计算演化参考模型流程图；

图8为计算最小稳定窗口长度和谱线稳定值流程图；

图9为罗辑斯蒂(Logistic)映射复杂度谱线稳定值演化参考模型图(u∈[0.75,2])；

图10为罗辑斯蒂(Logistic)映射复杂度谱线变化幅度值分布图(u∈[0.75,2])；

图11为罗辑斯蒂(Logistic)映射复杂度谱线最小稳定计算窗口长度分布图(u∈[0.75,2])。

具体实施方式

下面结合附图和实施示例对本发明作进一步说明。本实施示例针对已获取的某一脑电信号进行分类和识别的示范过程。

图1是本发明所述的方法流程图，具体包括以下步骤：

步骤一：计算基于功率谱划分的脑电信号复杂度谱，方法如图2所示，包括以下内容：

(1)根据公式(1)～(3),对脑电信号计算基于功率谱划分的复杂度谱序列。图4是某一脑电信号功率谱划分为36时的脑电复杂度谱，其中m值是根据步骤一的(2)按照精度0.02确定的。

(2)根据图3所示，按照公式(4)～(7)，计算不同精度下的功率谱划分需要的最小复杂度的长度。本发明已经计算得出了罗辑斯蒂(Logistic)映射在0.75～2之间序列演化时，精度为0.5～0.01之间需要的最小划分数m，如图5所示。图6为随计算精度变化所需的最小划分数m变化曲线。

步骤二：建立非线性罗辑斯蒂映射的复杂度谱演化参考模型，方法如下：

(1)根据公式(8)，按照图7的流程计算不同窗口的复杂度值，直到稳定为止。构成复杂度谱序列。

(2)根据公式(9)～(20)，按照图7的流程计算最小稳定窗口长度和稳定复杂度谱，以及谱线的变化幅度。

(3)根据计算的稳定复杂度谱值绘制复杂度参考模型图，如图9所示；根据谱线值变化幅度绘制参考模型变化分布图，如图10所示；根据最小稳定窗口长度绘制参考模型最小窗口变化分布图，如图11所示。

步骤三：计算各类脑电信号的复杂度谱序列并提出特征在参考模型中标定，方法如下：

(1)利用步骤一计算多种类型脑电信号的复杂度谱，根据复杂图谱距离参考模型的谱线的欧氏距离排序，将其标定在参考模型中距其最小谱线所在的区域。

(2)逐步积累大量的样本点，最后在参考模型中形成不同类型的脑电信号模型区域。通过这些区域分析其空间分布关系，进而分析其内在联系、演化关系和物理意义。

步骤四：对新的脑电未知信号，计算其基于功率谱划分的复杂度谱，并将其与参考模型空间的各类脑电信号区域做相似度对比，计算其和各类空间之间的距离dis,最后依据dis的大小进行排序，取最小值的dis的类作为分类和识别输出。例如：图4的复杂度谱序列s₁在演化模型空间中，经计算得到与Φ特征空间、与Θ特征空间以及与Π特征空间的综合欧氏距离分别为0.15、0.05和0.3。对这三个数进行排序，最小者为0.05，它对应着Θ特征空间。据此，输出的分类和识别结果为：s₁为Θ类型脑电信号的复杂度谱。

Claims (3)

1.一种基于功率谱划分的复杂度谱的脑电信号分类识别方法，其特征在于包括以下步骤：

步骤一：基于功率谱划分的复杂度谱模型的定义及计算，方法如下：

(1)对于给定的脑电信号序列X＝{x(n),n＝0,1,2,...,N-1}，计算其傅立叶(Flourier)变换，得到功率谱序列，记为F；

(2)利用公式计算序列F在f(k)点处的功率谱值，得到功率谱值序列:P＝{p(k),k＝0,1,2,...,N-1}；

(3)寻找P的最大值和最小值：P_max＝max(P)，P_min＝min(P)，计算其差值ΔP＝P_max-P_min，将ΔP平均分成m段，其中：2≤m≤N的正整数,N为P的元素个数，并从下到上依次记为ΔP₁,ΔP₂,…,ΔP_m，记新序列为则按如下规则生成新序列：

得到新序列

(4)对于序列进行傅立叶逆变换，得新的时间序列为：

$\widetilde{x^j}\left(n\right)=\underset{k=0}{\overset{N-1}{\Σ}}\widetilde{f_N^j}\left(k\right)W_N^{nk},k=0,1,2,...,N-1,j=1,2,...,m,$

表示的傅立叶变换项；

(5)定义如下求和比：

$C_X\left(j\right)=\frac{\underset{n=0}{\overset{N-1}{\Σ}}|\widetilde{x^j}\left(n\right){|}^2}{\underset{n=0}{\overset{N-1}{\Σ}}|x\left(n\right){|}^2},j=1,2,...,m,$

C_X(j)为序列X的第j复杂度值，称序列S＝{C_X(j),j＝1,2,...,m}，为脑电信号时间序列X的基于功率谱划分的复杂度谱；

步骤二：建立非线性罗辑斯蒂(Logistic)映射的复杂度谱参考演化模型，方法如下：

(1)根据罗辑斯蒂映射公式，生成控制参数u∈[0.75,2]的演化区间作为建立脑电信号复杂度谱从周期经混沌向随机演化的标准数据，按照步骤一的方法计算脑电信号每一时间序列的复杂度谱：

$S^i=\{C_X^i\left(j\right),j=1,2,...,m\},i\∈\[0.75,2\],$

Sⁱ表示i对应复杂度谱；

(2)按照给定的数据步长，逐步加大计算窗口的长度并重复计算S,直到各复杂度谱中的每个复杂度值收敛或变化稳定为止，此时计算和统计稳定区间的如下过程数据：

1)计算最小数据窗口的长度：

N_min＝min{N₁,N₂,N₃,…N_i}，

N_i表示在收敛或变化稳定时各计算窗口长度；

2)计算S在稳定区域各谱线的近似平均值：

$\stackrel{\‾}{S}=\{\stackrel{\‾}{C_X^i\left(j\right)},j=1,2,...,m,i\∈\[0.75,2\]\},$

表示S的近似平均值序列；

3)计算S中各复杂度值的变化范围：

表示在时间序列X的控制参数u＝i时第j复杂度在不同长度的窗口计算得到的复杂度值序列；

(3)根据复杂度谱序列最后计算值，并对其校验后在三维坐标图中作出三维空间模型图，并标出变化范围的边界，作出各序列每个复杂度谱的窗口长度三维模型；

步骤三：计算已获取多种脑电信号的复杂度谱序列特征，并在参考模型中标定，方法如下：

(1)对于已收集到的多种脑电信号，按照步骤一和步骤二的计算方法，计算N_min、和指标；

(2)并将这些指标标定在窗口长度三维模型中，以此确定这几类脑电信号指标在这两个模型中的空间位置和分布范围的大小；

步骤四：对于新采集的脑电未知信号计算其特征指标，并将其与参考模型的特征空间做相似度对比，计算其和各类空间之间的距离dis,最后依据dis的大小进行排序，取最小值的dis的类作为分类识别输出。

2.根据权利要求1所述的一种基于功率谱划分的复杂度谱的脑电信号分类识别方法，其特征在于，步骤一中确定功率谱划分数目m的方法如下：

(1)构造关于m的函数

$E_{CX}\left(m\right)={\left(\underset{i=1}{\overset{m}{\Σ}}{\left(C_X\left(j\right)\right)}^2\right)}^{\frac{1}{2}},$

作为目标函数，m是整数，2≤m≤N，N为P的元素个数；

(2)作出目标函数曲线，对每一个m求其线性拟合曲线：

y＝a_mm+c，

其中：a_m是取m时的斜率，c截距，得到一个斜率序列A＝{a₂,a₃,…,a_m}；

(3)计算A中每两个相邻项之间的差值的绝对值序列：

Y＝{Δa_i＝|a_i+1-a_i|，i＝1,2,3,4,…,m-1}，

Y表示相邻两项的差序列；

(4)在Y中寻找这样的项Δa_g，使得Y中所有i≥g的项，都有Δa_i≤Δa_g，则第g项Δa_g就是序列Y开始收敛或稳定变化的起始点，从g项后所有序列项就是收敛和稳定变化的区域，此时，Y序列相对Δa_g的稳定或收敛的变化精度为：

$\mathrm{\δ}=\frac{{\mathrm{\Δa}}_g}{max\left\{Y\right\}},$

其中：0<δ≤0.05，g即是相对于δ开始收敛点，也是相对于精度δ的最小划分数目。

3.根据权利要求1所述的一种基于功率谱划分的复杂度谱的脑电信号分类识别方法，其特征在于，在权利要求1步骤二中所提到的计算最小数据窗口的长度和计算S在稳定区域各谱线的近似平均值，计算方法如下：

(1)脑电信号序列X在u＝i时第j复杂度在以一定步长增加的窗口下,计算得到的复杂度谱序列为：

v表示起始计算时的数据窗口长度，d表示每次增加的计算步长；

(2)计算每相邻两项的中位值：

$t_{v+id}=\frac{1}{2}(Z_{v+id}+Z_{v+(i+1)d}),$

并构成新的序列：

T＝{t_v+id，i＝0,1,2,…}，

T表示新的平均序列；

(3)以新序列T代替重复迭代权利要求3中的第(2)步10次，取其最后一项的值作为近似平均复杂度值t_end；

(4)将序列的每一项Z_v+id与t_end作差值：

ΔZ_v+id＝Z_v+id-t_end，

取序列中ΔZ_v+id≥0项的元素下标组成新的序列H_big，取序列中ΔZ_v+id≤0项的元素下标组成新的序列H_small；

(5)计算序列H_big中的最小值：

b_big＝min{H_big}，

计算H_small中的最小值：

b_small＝min{H_small}，

然后求h_big和h_small两者中的最大值,记为N₁，其值为：

N₁＝max{b_big,b_small}，

从而，得到生成新序列的起始点；

(6)在权利要求3中第(1)步的复杂度谱序列中，用H表示其下标，对该序列下标H≥N₁的所有项N₁,N₁+d,N₁+2d,…,将其组成一个新序列G：

G＝{Z_N1,Z_N1+d,Z_N1+2d,...}，

设G有r项，计算序列G的平均值：

${\stackrel{\&OverBar;}{Z}}_1=\frac{1}{r}\underset{i=0}{\overset{r-1}{\&Sigma;}}{Z}_{N_1+id},$

表示G的平均值；

(7)用替代t_end，重复权利要求3中的第(3)～(5)步，直到计算得到的N_i值与上一次计算的N_i-1值相同为止，记最后一次计算的N_i为N_min，为

这样，最后计算得到的N_min和即为所求的最小数据窗口长度和S在稳定区域各谱线的近似平均值。