CN1303247A - 电岛记录的频率估计 - Google Patents

Abstract

直接确定EIG事件的基本频率,而不用首先检测事件本身。利用包含多于一个事件的分析窗口(3-5)。利用与在语音中用来检测和估计音调的那些技术有些类似的技术(1-8)。本发明基于两个重要的生物学发现。第一个是EIG由功能步产生器产生。第二个是EIG信号在大多数时间是准周期的。音调检测算法的使用基于EIG信号的基本准周期特性。准周期是指,在连续事件之间的间隔不是准确地相同,而是可以稍微改变;连续事件的振幅和形状也可以呈现某些变化。

Description

电岛记录的频率估计

本发明一般涉及一种用来处理电生理信号的方法和系统，更具体地说，涉及与血糖浓度有关的电活动的检测和分析。

血糖水平监视对于糖尿病患者很重要。葡萄糖水平的连续监视能大大地减少由代谢失衡引起的医学并发症。

胰岛(胰岛)位于胰腺内，并且负责人体中胰岛素的制造。岛是一个多细胞串。在岛内的β细胞响应电活动冲击中的葡萄糖。在血糖监视器中由供给动物的胰腺导出的这种岛的使用，公开在美国专利5,101,814中。电岛记录(EIG)是胰岛电活动的测量。本发明利用EIG提供一种连续的血糖水平传感器。

研究证明在EIG信号的基本频率、与在围绕岛的介质中的葡萄糖水平之间的清楚相关。因此，电岛记录的基本频率的估计是有意义的实际值。

基本频率定义为EIG的“事件”的频率。认为“事件”代表在岛中细胞的同步电活动。通过模拟ECG，在EIG中的“事件”可与心脏周期(PQRST复合波)相比较。

尽管相信EIG处理不能由任何实物进行，但本专利申请的所有者试图直接检测事件，并且然后计算基本频率。与该手段有关的问题是，EIG“事件”的形状和尺寸变化很大，并因此难以实现可靠和健全的事件检测。

按照本发明，直接确定EIG事件的基本频率，而不首先检测各个事件本身。所有这些算法必须使用包含多于一个事件的分析窗口。本发明利用与用来检测和估计语音处理中音调的那些类似的技术。

本发明基于两个重要的生物学发现。第一个是EIG由功能步产生器产生。第二个是EIG信号在大多数时间是准周期的。因为EIG信号的基本准周期特性，使用音调检测算法。用准周期我们指，(1)在连续事件之间的间隔不是准确相同，而是稍微改变；(2)连续事件的振幅和形状也可以呈现某些变化。

试验了几种音调检测算法。他们中的三种达到了良好结果：自相关、分段自相关及谐波波峰分析。本发明的最佳实施例集中在基于自相关方法的算法上。

算法的最佳形式包括如下步骤：

●非EIG信号的检测。该信号可以包括非EIG段，如人为段和无声段。扫描该信号，并且标记和忽略非EIG段。

●把信号划分成重叠分析窗口，每个窗口四秒长，并且与相邻窗口具有75％的重叠。分析窗口包含多于一个事件。

●把用于在语音处理中的音调检测的类型的自相关变换的改进形式应用于单个分析窗口。对于每个分析窗口重复该步骤。

●由分析窗口的自相关值导出基本频率。通常基本频率由最大自相关值指示。对于每个分析窗口重复该步骤。

●后处理器用来“滤除”所有各个分析窗口的结果。在该阶段添加以前标记的非EIG段(人为段和无声段)。

为了产生修改的自相关转换，基于算法的现存自相关按如下方式适合于EIG：

●改进算法设计成用来由自相关值确定音调。算法通常选择最高的自相关波峰(值)。在EIG中我们发现，有时真正的音调不由波峰代表，而由在几个相邻波峰之间的波谷代表。一种算法设计成确定这些情形，并且正确地估计音调。我们称这种现象为“火山”形自相关曲线，因为发现“山”的中心在较低的地面上。

●来自语音处理的一种发声/不发声决定机构是适用的。把“不发声”EIG段定义为非信号(未决定段)。一个后处理器用来决定这些未决定段的音调。尽管不发声语音段的确存在，但“不发声”EIG段是一种虚拟的非信号，而不真正存在。

●设计一种专用预处理算法。信号经受卷积以便增大事件的宽度。这不象目的在于增大高振幅部分和/或滤除共振峰的语音预处理。这种预处理技术称作“变肥”。

●使用4秒非常长的分析窗口，以便找到在从0.25 Hz至5 Hz范围的频率。在语音中，习惯使用约30毫秒的窗口，以便找到在从80 Hz至300 Hz范围的频率。

本发明也考虑了分段自相关算法的使用。这认为是用于EIG的最好方法，但其他的方法也是适用的。分段自相关在“使用分段自相关的语音信号的音调检测(Pitch detection of speech signals usingSegmented Autocorrelation)”/I.A.Atkinson,A.M.Kondoz,B.G.Evans/Electronics Letters Vol.31 No.7 pp.533-535/1995年3月中描述。

参照附图由当前最佳实施例的如下详细描述，将更完整地理解以上简要描述以及本发明的进一步目的、特征、及优点，在附图中：

图1是功能方块图，表示最佳实施例的过程的概要；

图2是流程图，描述预检测器、无声检测及人为段检测器；

图3是流程图，描述一个单独分析窗口的基本频率估计；

图4是图3的块25之一的流程图扩展，并且描述基本频率估计器的决定算法；及

图5是流程图，描述后处理相位算法，这些算法与来自相邻分析的信息相结合以便纠正局部估计误差；

参照附图现在将详细地描述本发明的最佳实施例。首先用于基本频率估计的过程和系统的最佳实施例的示意概要呈现在图1中。算法包括如下步骤：

a．非EIG信号的检测。该信号可以包括非EIG段，如人为段和无声段。扫描该信号，并且标记和忽略非EIG段。使用一个无声检测器1和一个人为段检测器2。“人为段”是严重污染EIG信号的主要干扰噪声信号。

b．把信号划分成重叠分析窗口。每个窗口四秒长，并且与相邻窗口具有75％的重叠。分析窗口包含多于一个事件。考虑分析窗口能多达在连续事件之间间隔的四十倍。

c．把自相关变换应用于单个分析窗口。对于每个分析窗口重复该步骤。预处理3是选择性的，并且在最佳实施例中不使用。

d．由分析窗口5的自相关值导出基本频率。通常基本频率由最大自相关值指示。对于分析窗口的估计基本频率计算确定性大小6。对于每个分析窗口重复该步骤。

e．一个后处理器7用来“滤除”所有各个分析窗口的结果。在该阶段添加以前标记的非EIG段(人为段和无声段)。

f．输出是信号的估计基本频率8。

本发明的最佳实施例包括最好在使用Windows 95操作系统的奔腾PC上运行的软件。算法在“The Mathworks Inc.”的“Matlab”软件上实施。算法以Matlab语言写成，并且在Matlab程序外壳。对于“The Mathworks Inc.”的“Matlab”也需要“数字信号处理工具箱”。

图2是描述预检测器的流程图。无声检测器(11、12、13)包括两个步骤。在第一步骤11，在每秒期间计算振幅相关度量，及在第二步骤12，把振幅相关度量与一个固定阈值相比较。“Sig”是原始数字输入信号。它是从岛记录的、在模数转换之后的原始信号。在块11中，检查每个一秒窗口。“Max(Sig)”是在检查的一秒窗口内的最大样本值。“Min(Sig)”是在检查的一秒窗口内的最小样本值。由于“Sig”总是实数，所以计算实数最大值和最小值。如果度量在阈值以下大于5秒，那么把该段分类为无声段13。该阈值取决于测量设备和环境。它也取决于数字化尺度。在最佳实施例中，EIG信号的最大振幅接近4000振幅单位，从而使用100振幅单位的阈值。“振幅单位”是在模数转换器的两个相邻量化级之间的振幅差。在最佳实施例中，一个“振幅单位”对应于原始电信号的0.25微伏特。“原始”是指在放大之前。

人为段检测器(14、15、16、17)包括两个步骤。在第一步骤14，把信号与一个适合阈值相比较。使用300秒长的记录通过如下计算，计算适合阈值：

把一个高基准定义为300秒长记录的振幅值的直方图的百分之99.833。

把一个低基准定义为300秒长记录的振幅值的直方图的百分之0.167。

把振幅基准定义为高基准减低基准。

把一个高阈值定义为高基准加振幅基准。

把一个低阈值定义为高基准减振幅基准。

通过简单地分类300秒长记录的所有样本值(“振幅”)产生在块14中使用的“直方图”。从最小值至最大值进行分类。根据这些分类值导出两个基准。“高基准”是分类值的百分之99.833。“低基准”是分类值的百分之0.167。术语“直方图”是指值的分类，并且能由术语“分类振幅值”代替。

第二步骤15包括检查是否有其中信号的振幅高于高阈值或低于低阈值的样本。如果发现这样的样本，那么把他们分类为包含人为段16。剩余样本17不分类为人为段或无声段。

图3包含描述在每个分析窗口上进行的步骤的流程图。每个分析窗口是4秒长。使用重叠分析窗口。在两个连续窗口之间有3秒重叠。不分析包含无声段或人为段的分析窗口。分析窗口在块21中描述，并且由如下注释数学表示：

x(n)是信号的原始数字样本。我们较早称该原始输入信号为“Sig”。这是相同的信号，只是这里使用更严格的数学注释。分析窗口包含4秒的信号。在最佳实施例中的抽样频率(“FS”)是100Hz，因此分析窗口包含400个样本。这用“N(窗口长度)=400”标记。在分析窗口中的第一样本是x(l)，其中“l”是样本指数。因此在分析窗口中的最后样本是x(l+399)=x(l+N-1)。因此分析窗口由样本描述，它包含从x(l)至x(l+N-1)。

可以使用预处理22，但在最佳实施例中没有实施。

矩形窗口23应用于样本，并且计算归一化偏置自相关转换24。它用如下公式数学表示： $\mathrm{\φ}\left(m\right)=\frac{1}{{\mathrm{\Σ}}_{k=0}^{N-1}{\left[x(k+1)w\left(k\right)\right]}^2}{\mathrm{\Σ}}_{n=o}^{N-\left|m\right|-1}\left[x(n+1)w\left(n\right)\right]\left[x(n+1+m)w(n+m)\right]---\left(1\right)$

其中“w(n)”是开窗口函数。由于“w(n)”是矩形窗口，如在块23中描述的那样，所以它等效于分析窗口内的“1”：

w(n)=1  0≤n≤N-1

0  其他

自相关转换在L.R.Rabiner和R.W.Schafer,1978年的书“语音信号的数字处理(Digital Processing of Speech Signals)”第141-164页中进一步描述。基本频率决定算法在图4中更详细地描述。

我们现在参照图4，并且以后返回到图3。

图4是流程图，描述基本频率决定算法。实际上它在基本周期上决定，由此能导出基本频率。为了由自相关值确定基本频率设计一种改进的算法。在典型的算法中，通常选择最高自相关波峰(值)。在EIG中，发现有时真正的单调不由波峰代表，而是由在几个相邻波峰之间的波谷代表。一种算法设计成确定这些情况，并且正确地估计音调。我们称这种现象为“火山”形自相关曲线，因为发现“山”的中心在较低的地面上。

至决定算法的输入是一个包含自相关系数的向量31。步骤32和33求出在允许频率范围内的最高极点。“火山”效应检测器34首先检查在可疑“火山”区域中是否发现最高波峰之一。如果发现这样一种波峰，则其指数(i)通过以便检查与最高波峰的振幅相比较的波峰振幅35。如果结果是肯定的，那么检测“火山”36，并且估计基本周期是与最高波峰对应的自相关滞后的一半。自相关由以上公式(1)表示。m是“自相关滞后”。它是在对其计算相关的两个段之间的滞后。第一段是x(n+l)至x(n+l+N-｜m｜-1)，而第二段是x(n+l+m)至(x+l+m+N-｜m｜-1)。同样，见L.R.Rabiner和R.W.Schafer,1978年“语音信号的数字处理(Digital Processing of Speech Signals)”第141-164页。

一种音调平分检测算法37首先检查自相关系数滞后的比值37，并且然后检查其振幅38。如果最高波峰之一都满足这些条件，那么把它选择为估计基本周期39。如果“火山”检测器和音调平分检测器都产生否定的结果，那么假定估计基本周期是最高自相关波峰的滞后40。图4的基本周期决定算法的输出是估计的基本周期。该输出返回到图3中的块26。

“发声”决定机构26用来决定在当前分析窗口中是否应该估计基本频率。如果在允许滞后(20至333个样本的滞后)范围内的最大自相关系数低于0.3，那么把分析窗口分类为“不发声”，并且不估计基本频率27。在该情况下，把零确定性等级送给分析窗口。

发声/不发声决定机构适用于来自语音处理。把“不发声”EGI段定义为非信号(未决定段)。后处理器在后阶段中用来决定这些未决定段的基本频率。注意不发声语音段的确存在，而“不发声”EIG段是一种虚拟的非信号，并且不真正存在。

如果当前分析窗口没有分类为“不发声”，那么估计过程继续，并且对于当前分析窗口计算确定性等级28。较高的确定性等级指示有较高的概率正确地估计基本频率。对于每个分析窗口返回几个输出29：

1．估计基本周期-如果把分析窗口分类为“不发声”，则不返回该值；

2．“发声”/“不发声决定”-规定在分析窗口中是否估计基本周期；及

3．确定性等级-如果把分析窗口分类为“不发声”，则不返回该值。

图5描述估计过程的最后阶段：后处理器。后处理器把所有获得数据41组合成原始信号的基本频率估计。后处理器在原始信号的300秒长记录上操作。第一步骤42是把所有估计基本周期值转换成估计基本频率(“音调”)值。该阶段的输出是以赫兹对每个“发声”分析窗口的基本频率估计。下个步骤43把以前检测的无声段1插入到基本频率估计结果中。在仅包含无声段的所有窗口中分配基本频率一个0(Hz)值。包含非相象基本频率估计的所有分析窗口在下个步骤44标记为“不发声”。

下个步骤45把一个基本频率(“音调”)值分配给无声段。后处理器在所有分析窗口上循环一次。它把音调值分配给所有“不发声”窗口(包含人为段的窗口也处理成“不发声”)。从相邻窗口的音调值的平均值导出音调值。能添加一些小数据相关改进，以便产生一种较灵活的平均方法，但这些是与数据非常相关的小变化。然后应用一个音调平分/加倍误差校正机构46。该算法在所有分析窗口上循环一次。如果当前窗口的音调是相邻窗口平均的一半或两倍，那么当前音调取作相邻的灵活平均。术语“灵活平均”是指全体数值在除去极值之后的平均。例如：我们按如下方式定义8个数值的“灵活平均”：

1．除去最大和最小数值及

2．计算剩余6个数值的平均值。

使用这种“灵活平均”，以便得到更稳定和可靠的结果。

最后频率47寻找具有与其相邻显著不同的音调的分析窗口。该算法在所有分析窗口上循环一次。如果当前窗口的音调大于或小于相邻窗口平均值的35％，那么当前音调取作相邻的灵活平均。

整个基本频率估计过程的最终输出48是对于每个分析窗口的基本频率值。

尽管为了说明目的已经公开了本发明的最佳实施例，但熟悉本专业的技术人员将认识到多种添加、改进及替代是可能的，而不脱离本发明的范围和精神。例如，本发明的原理和应用不限于描述的具体生物学现象。他们同样好地能用来测量，例如来自肌肉活细胞的光学信号或神经元的电信号、或者由类似或其他生物学微观结构产生的各种信号。特别是，认为本发明能在其中细胞或其他生物微观结构在其自然环境外运动的任何用途中发现实用性和用作生物传感器。

Claims (20)

1．一种用来监视产生以事件为特征的可检测信号的生物微观结构活动的方法，包括由可检测信号估计事件发生基本频率的步骤，而不用检测各个事件的发生。

2．根据权利要求1所述的方法，其中在横跨多于一个事件的分析窗口期间分析事件。

3．根据权利要求1或2任一项所述的方法，其中事件是由从其原始环境变位的生物微观结构产生的信号。

4．根据权利要求1或2任一项所述的方法，其中事件是由活细胞产生的信号。

5．根据权利要求1或2任一项所述的方法，其中事件是由胰岛产生的信号。

6．根据权利要求1或2任一项所述的方法，其中估计步骤包括自相关运算。

7．根据权利要求6所述的方法，进一步包括如下步骤之一：

根据位于在几个相邻波峰中的较低自相关值，估计基本频率；

把可检测信号的“不发声”段处理成关于音调未决定的段，并且经以后处理估计这些段的音调；

寻求估计在.25至5赫兹范围内的基本频率；

利用在几秒范围内的分析窗口持续时间；

进行具有增大事件有效持续时间的效果的预处理运算；及

利用进行分段自相关的自相关过程。

8．根据权利要求1所述的方法，进一步包括把在预定时间内的“不发声”窗口数量用作病人血液成分水平的度量的步骤。

9．根据权利要求1所述的方法，进一步包括把确定性等级用作病人血液或组织成分水平的度量的步骤。

10．根据权利要求1所述的方法，进一步包括把基本频率作为病人血液或组织成分水平的度量的步骤。

11．根据权利要求1所述的方法，其中血液成分水平是在生物微观结构附近的血糖水平。

12．根据权利要求1所述的方法，其中事件由用作病人体内的测头的胰岛中的活细胞产生，并且基本频率估计用作病人血糖水平的度量。

13．根据权利要求2任一项所述的方法，其中分析窗口横跨高达连续事件之间间隔40倍的持续时间。

14．在一种用来监视产生以事件为特征的可检测信号的生物微观结构活动的系统中，一个传感器能够接收可检测信号，并且一个处理器包括一个用来由可检测信号估计事件发生基本频率的模块，而不用首先检测各个事件的发生。

15．根据权利要求14所述的系统，其中处理器进一步包括一个分析窗口，在该分析窗口期间分析事件，该分析窗口横跨多于一个事件。

16．根据权利要求14或15任一项所述的系统，其中用于估计的模块包括进行自相关运算的多个元件。

17．根据权利要求16所述的系统，用于估计的模块进一步包括如下子模块之一：

一个子模块，根据位于在几个相邻波峰中的较低自相关值，估计基本频率；

一个子模块，把可检测信号的“不发声”段标识成关于音调未决定的段，并且这些段的音调由一个以后处理子模块估计；

一个子模块，把基本频率估计控制在.25至5赫兹的范围内；

一个子模块，把分析窗口控制成具有在几秒范围内的持续时间；及

一个子模块，进行具有增大事件有效持续时间的效果的预处理运算。

18．根据权利要求14或15任一项所述的系统，其中处理器建造成完成一个分段自相关过程。

19．根据权利要求14或15任一项所述的系统，其中传感器是一个能够检测由胰岛中的活细胞发送的信号的测头，频率估计是带有这些细胞的病人的血糖水平指示。

20．根据权利要求15任一项所述的系统，其中分析窗口横跨高达连续事件之间间隔40倍的持续时间。

Images