CN102429662B - 家庭环境中睡眠呼吸暂停综合征的筛查系统 - Google Patents

Abstract

本发明公开一种家庭环境中睡眠呼吸暂停综合征的筛查系统，其中包括鼾声数据获取模块，鼾声处理模块，判别模块。数据获取模块采集鼾声数据；鼾声信号处理模块中首先对鼾声信号进行分帧、预加重、加窗等语音信号的预处理过程，然后采用双门限检测方法完成鼾声的端点检测，端点检测过程中以鼾声样本的前十帧数据作为背景噪声的第一测度，并将鼾声段数据传递给能量比求解模块。最终求解鼾声段的能量比特征，并将得到的能量比特性传递给判别模块中。判别模块中以能量比作为分类依据，能量比高于某一门限值即为OSAS患者，反之为单纯打鼾者。本发明是一种对患者睡眠干扰小、成本低、在家庭环境中亦能达到较高筛查精度的系统。

Description

家庭环境中睡眠呼吸暂停综合征的筛查系统

技术领域

本发明涉及生物医学领域中睡眠呼吸障碍疾病的筛查，特别涉及家庭环境中的阻塞性睡眠呼吸暂停综合征的筛查系统。

背景技术

睡眠呼吸障碍是睡眠过程中出现的呼吸异常，包括睡眠呼吸暂停综合征、低通气综合征、慢性肺部及神经肌肉疾患引起的有关睡眠呼吸障碍等，其中以阻塞性睡眠呼吸暂停综合征(OSAS)为主。调查发现，美国的中年人群中约有9％的女性和24％的男性患有不同程度的OSAS，OSAS患者占澳大利亚总人口的3-5％，印度则有5.4-7.5％人群受OSAS困扰。研究表明，OSAS会造成白天嗜睡，头昏，头疼，记忆力衰退，乏力，反应迟钝，睡眠行为异常等症状。长期患有OSAS可引起高血压、冠心病、心衰、中风等多种疾病。医学界对此疾病的研究十分重视，且已取得了重大成果。其中，多导睡眠仪(PSG)是公认的检测OSAS的金标准，已被广泛地应用到OSAS诊断中，但其设备有限、检测费用昂贵、对患者睡眠干扰较大以及需专业人员监视等不足极大地制约了PSG的普及化使用，超过90％的患者得不到及时诊治。

研究发现：每5个成年人中就有一个存在打鼾现象，而鼾声是OSAS最直接的临床表征，也就是说，“打鼾”这种现象很有可能潜伏着OSAS或者其他病症。此外，鼾声是由呼吸气流激励上气道产生的响应，而OSAS又是由上气道塌陷引起，因此可以通过分析鼾声特性得到上气道的特性，从而实现OSAS筛查。

基于以上原因，近年来，各国专家、学者已经对基于鼾声特性的OSAS筛查进行了探索性研究：例如通过统计鼾声时间间隔，分析基音周期，功率谱特性，小波域中鼾声幅度、密度和能量分布情况，共振峰特性等特性实现OSAS检测。上述已有方法均以实现OSAS准确诊断为最终目的。但是，我们的研究发现，仅利用鼾声这一单一特性几乎不可能达到与PSG同样或相近的检测精度，无法完全替代PSG检测。此外，现有研究对鼾声样本录制条件要求较高，而对录制条件的要求直接限制了方法的普及使用。总的来说，现有方法不适用于家庭环境中的检测。

发明内容

本发明的主要目的是公开一种对患者睡眠干扰小、成本低、在家庭等环境中亦能达到较高筛查精度的方法。

一种家庭环境中睡眠呼吸暂停综合征的筛查系统，

包括鼾声获取模块：采集鼾声并记录鼾声数据，采样率设为8k Hz、精度为16bit鼾声信号；

鼾声处理模块：

鼾声处理模块包括预处理部、端点检测部、能量比求解部，

其中预处理部，提供鼾声信号的分帧、预加重以及加窗过程，完成鼾声信号的预处理，

端点检测部，区分有声段和无声段，完成鼾声信号中鼾声段的检测，

能量比求解部，求解鼾声信号的子带能量比值，得到子带能量比值；

鼾声判决模块：

将上述鼾声处理模块中求得的子带能量比值与门限值比较。

其中鼾声处理模块的预处理部，在处理之前需要对其分帧，帧长取128ms，帧移为64ms，利用一阶FIR滤波器H(z)＝1-0.98z^-1实现预加重，在对鼾声信号进行频域处理之前使用如公式(1)的窗函数来减小频谱泄露，其中n为采样点：

$w\left(n\right)=\left\{\begin{array}{cc}0.54-0.46*\cos (n-1)*2*\mathrm{\&pi;}/512& 0<n\&le;256\\ 1& 256<n\&le;768\\ 0.54-0.46*\cos (n-1)*2*\mathrm{\&pi;}/512& 769<n\&le;1024\end{array}\right.---\left(1\right)$

其中鼾声处理模块的端点检测部，采用短时能量和短时过零率双门限分析法，具体过程如下：

1)计算每帧的短时能量：

$M_m=\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{\left|x_m\left(n\right)\right|}^2---\left(2\right)$

其中，N为帧长1024，m为帧号，n为采样点，M_m为短时能量，x_m(n)为鼾声时域波形。

2)假定前10帧为无声帧，计算各帧的短时能量和短时过零率，统计前10帧噪声的平均过零率Mzc，平均过零率方差Vzc，平均短时能量E。

确定门限值，能量门限值如公式(3)，其中，EL为低门限值，EH为高门限值，平均短时能量E：

$\left\{\begin{array}{c}\mathrm{EL}=2*E\\ \mathrm{EH}=5*\mathrm{EL}\end{array}\right.---\left(3\right)$

过零率门限zc如公式(4)，其中前10帧噪声的平均过零率Mzc，平均过零率方差Vzc：

zc＝Mzc+2*Vzc                            (4)

3)先根据能量门限EL和EH确定出鼾声段的初始起止点[N₁，N₂]，再利用过零率门限zc修正并获得最终起止点[N_B，N_E]，N_B为鼾声起始点，N_E为鼾声终止点，至此，完成鼾声信号的预处理及端点检测过程；

其中鼾声处理模块的能量比求解部，为了用数字信号处理方法对鼾声信号进行处理，首先需要建立鼾声信号产生的数学模型。其中上气道模型的建立和求解是通过线性预测编码(LPC)实现的，LPC分析过程如下所示：

当前值用s(n)表示，表示预测样值，则预测样值表示为：

其中a_i为线性预测系数，P为预测阶数，n为采样点，则预测误差e(n)定义为：

$e\left(n\right)=s\left(n\right)-\hat{s}\left(n\right)=s\left(n\right)-\underset{i=1}{\overset{P}{\mathrm{\&Sigma;}}}{a}_{i}s(n-i),$ $\begin{array}{c}{a}_0=1\\ i=\mathrm{0,1},\&CenterDot;\&CenterDot;\&CenterDot;P\end{array}---\left(5\right)$

按最小均方差准则确定预测系数，使得均方误差E[e²(n)]最小，E[e²(n)]定义为：

$E_n=\underset{n}{\mathrm{\&Sigma;}}{e}^2\left(n\right)=\underset{n}{\mathrm{\&Sigma;}}{[s\left(n\right)-\underset{i=1}{\overset{p}{\mathrm{\&Sigma;}}}{a}_{i}s(n-i)]}^2---\left(6\right)$

要求得最小E[e²(n)]，需要对E[e²(n)]关于a_i求偏导：

$\frac{\&PartialD;E\left(e^2\left(n\right)\right)}{\&PartialD;a_i}=2E\left[e\left(n\right)\frac{\&PartialD;e\left(n\right)}{\&PartialD;a_i}\right]=0,$ i＝0，1，…P                                    (7)

代入 $\frac{\&PartialD;e\left(n\right)}{\&PartialD;a_i}=s(n-i)$ 得：

E[e(n)s(n-i)]＝0                                    (8)

公式(5)代入公式(8)可得：

$E\left[e\left(n\right)s(n-i)\right]=E[s\left(n\right)s(n-i)-\underset{l=0}{\overset{P}{\mathrm{\&Sigma;}}}{a}_{l}s(n-l)s(n-i)]=0,$ $\begin{array}{c}{a}_0=1\\ i=\mathrm{0,2},\&CenterDot;\&CenterDot;\&CenterDot;P\end{array}---\left(9\right)$

因为信号自相关序列为：

R(i-l)＝E[s(n-l)s(n-i)]， $\begin{array}{c}{a}_0=1\\ l,i=\mathrm{0,2},\&CenterDot;\&CenterDot;\&CenterDot;P\end{array}---\left(10\right)$

则公式(9)可写为：

$R\left(i\right)-\underset{l=0}{\overset{P}{\mathrm{\&Sigma;}}}{a}_{l}R(i-l)=0,$ $\begin{array}{c}{a}_0=1\\ i=\mathrm{0,2},\&CenterDot;\&CenterDot;\&CenterDot;P\end{array}---\left(11\right)$

公式(11)称为线性预测的标准方程，写成矩阵形式为：

$\left[\begin{array}{cccc}R\left(0\right)& R\left(1\right)& ...& R(P-1)\\ R\left(1\right)& R\left(0\right)& ...& R(P-2)\\ ...& ...& ...& ...\\ R(P-1)& R(P-2)& ...& R\left(0\right)\end{array}\right]\left[\begin{array}{c}{a}_1\\ a_2\\ .\\ .\\ .\\ a_P\end{array}\right]=\left[\begin{array}{c}R\left(1\right)\\ R\left(2\right)\\ .\\ .\\ .\\ R\left(p\right)\end{array}\right]---\left(12\right)$

公式(11)就是尤利-沃克方程，解方程可得到p阶线性预测系数a_i和预测误差值，本发明利用莱文森-德宾递推算法求解尤利-沃克方程，该算法流程如下：

步骤1：定义初始条件：

E₀＝R(0)，a₁ ⁽⁰⁾＝1

步骤2：定义反射系数

$k_i=[R\left(i\right)-\underset{j=1}{\overset{i-1}{\mathrm{\&Sigma;}}}{a}_j^{(i-1)}R(i-j)]/E_{i-1}=a_i^{\left(i\right)};$

a_j ⁽ⁱ⁾＝a_j ^(i-1)-k_ia_i-j ^(i-1)；

E_i＝(1-k_i ²)E_i-1

该步从低阶开始递推，直到i＝p时截止，循环结束后给出各阶次的所有参数。

步骤3：线性预测的系数：

a_j＝a_j ^(p)1≤j≤p

经过以上过程，求得线性预测器的预测系数a_i，并建立p阶全极点模型AR(p)构建时变数字滤波器H(z)模拟上气道结构，其中U(z)为激励信号，S(z)为输出信号：

$H\left(z\right)=\frac{S\left(z\right)}{U\left(z\right)}=\frac{G}{1-\underset{i=1}{\overset{p}{\mathrm{\&Sigma;}}}{a}_i{z}^{-i}}---\left(13\right)$

通过上述模型求出上气道的LPC对数谱，鼾声受上气道的共振作用，呼吸气流的不同频率能量重新分配，其中最高能量对应频率即为第一共振峰。

选取的两通频带范围为：第一子带200-350Hz，第二子带500-650Hz，第一子带对应正常呼吸时鼾声第一共振峰集中频段，第二子带对应发生低通气时鼾声第一共振峰集中频段。

方案A直接求和子带能量比方法，确定低通气状态和正常状态下鼾声第一共振峰对应的频段后，利用AR模型构建预处理后鼾声信号的能量谱，分别计算所选两子带频段内频谱能量，即频域中实现带内能量谱累加求和，第一子带内能量为E₁，第二子带内能量为E₂，然后取E₂与E₁的比值，直接求和子带能量比如公式(14)，其中f表示频率，N为帧长，sub-band1表示第一子带频率范围，sub-band2表示第二子带频率范围，S_N(f)表示信号的短时LPC谱包络。

$\mathrm{Er}=E_2/E_1=\underset{f\&Element;\mathrm{sub}-\mathrm{band}2}{\mathrm{\&Sigma;}}{\left|S_N\left(f\right)\right|}^2/\underset{f\&Element;\mathrm{sub}-\mathrm{band}1}{\mathrm{\&Sigma;}}{\left|S_N\left(f\right)\right|}^2---\left(14\right)$

方案B改进的求和子带能量比方法，对直接求和能量比方法进行改进，通过加权系数对带内鼾声能量谱加权求和，其中第一子带加权能量和为E₁′，第二子带加权能量和为E₂′，加权求和子带能量比如公式(15)，其中B₁(f)表示第一子带加权系数，B₂(f)表示第二子带加权系数，其它表示同公式(14)。

${\mathrm{Er}}^{\&prime;}=E_2^{\&prime;}/E_1^{\&prime;}=\underset{f\&Element;\mathrm{sub}-\mathrm{band}2}{\mathrm{\&Sigma;}}{\left|S_N\left(f\right)\right|}^2*B_2\left(f\right)/\underset{f\&Element;\mathrm{sub}-\mathrm{band}1}{\mathrm{\&Sigma;}}{\left|S_N\left(f\right)\right|}^2*B_1\left(f\right)---\left(15\right)$

方案B借鉴梅尔滤波器的设计思想，以汉明窗形式给出加权系数，加权系数曲线如图1所示。

本发明以两个频段的能量比作为特征对OSAS进行筛查，这样，我们研究的不是某个单频点，而是整个频率段的信息。即使第一共振峰频率由于噪声的影响而出现偏差，只要偏差不是很大，仍可落在设定的频段内，算法仍然可以做出正确的判决。具有更强的鲁棒性和针对性，更适用于家庭等环境中OSAS筛查。

                                        (15)

附图说明

图1(a)是低频段加权系数曲线

图1(b)是高频段加权系数曲线

图2(a)是正常鼾声和低通气鼾声的时域图

图2(b)是正常鼾声和低通气鼾声的频域图

图3(a)是本发明提出的方案A中各步骤的流程图

图3(b)是本发明提出的方案B中各步骤的流程图

图4为本发明提出两种方法的箱线图对比

图5为本发明提出两种方法性能评价的ROC曲线

具体实施方式

如附图所示本发明的家庭环境中睡眠呼吸暂停综合征的筛查系统，

包括鼾声获取模块：所用鼾声信号利用非接触单向麦克风(20-20000Hz，SM81)置于患者口鼻上方18-20cm处录制6-8小时鼾声信号，录制时采样率设为8k Hz，精度为16bit，并利用电脑(IBM Think Pad 400)记录并保存数据。

鼾声处理模块：

鼾声处理模块主要包括预处理部、端点检测部、能量比求解部，

其中预处理部，提供鼾声信号的分帧、预加重以及加窗过程，完成鼾声信号的预处理。

研究发现：99％的鼾声持续时间分布在0.5s到1.8s之间。98.5％的鼾声间隔时间分布在1.4s到4.0s。在一个鼾声持续时间内，鼾声的时域参数和频率参数是基本保持不变的。

因为鼾声信号具有短时平稳性，在处理之前需要对其分帧，帧长取128ms，帧移为64ms。利用一阶FIR滤波器H(z)＝1-0.98z^-1实现预加重，其目的是为了去除口唇辐射的影响，增加鼾声信号的高频分辨率。在对鼾声信号进行频域处理之前使用如公式(1)的窗函数来减小频谱泄露，其中n为采样点：

$w\left(n\right)=\left\{\begin{array}{cc}0.54-0.46*\cos (n-1)*2*\mathrm{\&pi;}/512& 0<n\&le;256\\ 1& 256<n\&le;768\\ 0.54-0.46*\cos (n-1)*2*\mathrm{\&pi;}/512& 769<n\&le;1024\end{array}\right.---\left(1\right)$

端点检测部，区分有声段和无声段，完成鼾声信号中鼾声段的检测。

采用短时能量和短时过零率双门限分析法，具体过程如下：

1)计算每帧的短时能量：

$M_m=\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{\left|x_m\left(n\right)\right|}^2---\left(2\right)$

其中，公式编辑N为帧长1024，m为帧号，n为采样点，M_m为短时能量，x_m(n)为鼾声时域波形。

2)假定前10帧为无声帧，计算各帧的短时能量和短时过零率，统计前10帧噪声的平均过零率Mzc，平均过零率方差Vzc，平均短时能量E。

确定门限值，能量门限值如公式(3)，其中，EL为低门限值，EH为高门限值，平均短时能量E：

$\left\{\begin{array}{c}\mathrm{EL}=2*E\\ \mathrm{EH}=5*\mathrm{EL}\end{array}\right.---\left(3\right)$

过零率门限zc如公式(4)，其中前10帧噪声的平均过零率Mzc，平均过零率方差Vzc：

zc＝Mzc+2*Vzc                    (4)

3)先根据能量门限EL和EH确定出鼾声段的初始起止点[N₁，N₂]，再利用过零率门限zc修正并获得最终起止点[N_B，N_E]，N_B为鼾声起始点，N_E为鼾声终止点，至此，完成鼾声信号的预处理及端点检测过程；

能量比求解部，求解患者的子带能量比值，以子带能量比值作为特征实现OSAS检测。

为了用数字信号处理方法对鼾声信号进行处理，首先需要建立鼾声信号产生的数学模型。因为鼾声也是一种声音，其模型与语音信号类似。本发明是针对OSAS疾病的检测，而鼾声是呼吸气流激励上气道产生的响应，因此，我们只关心患者的上气道模型。上气道模型的建立和求解是通过线性预测编码(LPC)实现的，LPC分析过程如下所示：

当前值用s(n)表示，

表示预测样值，则预测样值表示为：

其中a_i为线性预测系数，P为预测阶数，n为采样点，则预测误差e(n)定义为：

$e\left(n\right)=s\left(n\right)-\hat{s}\left(n\right)=s\left(n\right)-\underset{i=1}{\overset{P}{\mathrm{\&Sigma;}}}{a}_{i}s(n-i),$ $\begin{array}{c}{a}_0=1\\ i=\mathrm{0,1},\&CenterDot;\&CenterDot;\&CenterDot;P\end{array}---\left(5\right)$

按最小均方差准则确定预测系数，使得均方误差E[e²(n)]最小，E[e²(n)]定义为：

$E_n=\underset{n}{\mathrm{\&Sigma;}}{e}^2\left(n\right)=\underset{n}{\mathrm{\&Sigma;}}{[s\left(n\right)-\underset{i=1}{\overset{p}{\mathrm{\&Sigma;}}}{a}_{i}s(n-i)]}^2---\left(6\right)$

要求得最小E[e²(n)]，需要对E[e²(n)]关于a_i求偏导：

$\frac{\&PartialD;E\left(e^2\left(n\right)\right)}{\&PartialD;a_i}=2E\left[e\left(n\right)\frac{\&PartialD;e\left(n\right)}{\&PartialD;a_i}\right]=0,$ i＝0，1，…P                            (7)

代入 $\frac{\&PartialD;e\left(n\right)}{\&PartialD;a_i}=s(n-i)$ 得：

E[e(n)s(n-i)]＝0                                    (8)

公式(5)代入公式(8)可得：

$E\left[e\left(n\right)s(n-i)\right]=E[s\left(n\right)s(n-i)-\underset{l=0}{\overset{P}{\mathrm{\&Sigma;}}}{a}_{l}s(n-l)s(n-i)]=0,$ $\begin{array}{c}{a}_0=1\\ i=\mathrm{0,2},\&CenterDot;\&CenterDot;\&CenterDot;P\end{array}---\left(9\right)$

因为信号自相关序列为：

R(i-l)＝E[s(n-l)s(n-i)]， $\begin{array}{c}{a}_0=1\\ l,i=\mathrm{0,2},\&CenterDot;\&CenterDot;\&CenterDot;P\end{array}---\left(10\right)$

则公式(9)可写为：

$R\left(i\right)-\underset{l=0}{\overset{P}{\mathrm{\&Sigma;}}}{a}_{l}R(i-l)=0,$ $\begin{array}{c}{a}_0=1\\ i=\mathrm{0,2},\&CenterDot;\&CenterDot;\&CenterDot;P\end{array}---\left(11\right)$

公式(11)称为线性预测的标准方程，写成矩阵形式为：

$\left[\begin{array}{cccc}R\left(0\right)& R\left(1\right)& ...& R(P-1)\\ R\left(1\right)& R\left(0\right)& ...& R(P-2)\\ ...& ...& ...& ...\\ R(P-1)& R(P-2)& ...& R\left(0\right)\end{array}\right]\left[\begin{array}{c}{a}_1\\ a_2\\ .\\ .\\ .\\ a_P\end{array}\right]=\left[\begin{array}{c}R\left(1\right)\\ R\left(2\right)\\ .\\ .\\ .\\ R\left(p\right)\end{array}\right]---\left(12\right)$

公式(11)就是尤利-沃克方程，解方程可得到p阶线性预测系数a_i和预测误差值。本发明利用莱文森-德宾递推算法求解尤利-沃克方程，该算法流程如下：

步骤1：定义初始条件：

E₀＝R(0)，a₁ ⁽⁰⁾＝1

步骤2：定义反射系数

$k_i=[R\left(i\right)-\underset{j=1}{\overset{i-1}{\mathrm{\&Sigma;}}}{a}_j^{(i-1)}R(i-j)]/E_{i-1}=a_i^{\left(i\right)};$

a_j ⁽ⁱ⁾＝a_j ^(i-1)-k_ia_i-j ^(i-1)；

E_i＝(1-k_i ²)E_i-1

该步从低阶开始递推，直到i＝p时截止，循环结束后给出各阶次的所有参数。

步骤3：线性预测的系数：

a_j＝a_j ^(p)  1≤j≤p

经过以上过程，求得线性预测器的预测系数a_i，并建立p阶全极点模型AR(p)构建时变数字滤波器H(z)模拟上气道结构，其中U(z)为激励信号，S(z)为输出信号：

$H\left(z\right)=\frac{S\left(z\right)}{U\left(z\right)}=\frac{G}{1-\underset{i=1}{\overset{p}{\mathrm{\&Sigma;}}}{a}_i{z}^{-i}}---\left(13\right)$

通过上述模型求出上气道的LPC对数谱，鼾声受上气道的共振作用，呼吸气流的不同频率能量重新分配，其中最高能量对应频率即为第一共振峰。当患者上气道塌陷或者阻塞而发生低通气时，鼾声最高能量出现在高频段，反之出现在低频段。

基于上述原理，本发明借鉴梅尔倒谱中子带能量的思想，提出一种新的、能减小噪声影响的基于鼾声子带能量比的OSAS筛查方法，并针对鼾声特性，简化了梅尔滤波器组的设计，仅使用两个滤波器的能量特性实现OSAS筛查。根据文献提供数据，并对本数据库内鼾声样本统计分析，本发明选取的两通频带范围为：第一子带200-350Hz，第二子带500-650Hz。前者对应正常呼吸时鼾声第一共振峰集中频段，后者对应发生低通气时鼾声第一共振峰集中频段。因为无病打鼾者只存在正常打鼾状态，其鼾声第一共振峰基本落在低频段，因此，第一子带能量高于第二子带能量，后者与前者的比值通常较小；反之二者的比值通常较大，根据此特性，可以达到区分OSAS患者和无病打鼾者的目的。

本发明基于上述原理提出两种子带能量比值的求解方案，直接求和子带能量比方法和改进求和子带能量比方法，如下所述：

方案A直接求和子带能量比方法

确定低通气状态和正常状态下鼾声第一共振峰对应的频段后，本发明利用AR模型构建预处理后鼾声信号的能量谱，分别计算本发明所选两子带频段内频谱能量，即频域中实现带内能量谱累加求和，第一子带内能量为E₁，第二子带内能量为E₂，然后取E₂与E₁的比值，最后根据子带能量比的特性实现OSAS筛查。直接求和子带能量比如公式(14)，其中f表示频率，N为帧长，sub-band1表示第一子带频率范围，sub-band2表示第二子带频率范围，S_N(f)表示信号的短时LPC谱包络。

$\mathrm{Er}=E_2/E_1=\underset{f\&Element;\mathrm{sub}-\mathrm{band}2}{\mathrm{\&Sigma;}}{\left|S_N\left(f\right)\right|}^2/\underset{f\&Element;\mathrm{sub}-\mathrm{band}1}{\mathrm{\&Sigma;}}{\left|S_N\left(f\right)\right|}^2---\left(14\right)$

方案B改进的求和子带能量比方法

利用梅尔倒谱的思想，本发明对直接求和能量比方法进行改进，通过加权系数对带内鼾声能量谱加权求和，其中第一子带加权能量和为E₁′，第二子带加权能量和为E₂′，加权求和子带能量比如公式(15)，其中B₁(f)表示第一子带加权系数，B₂(f)表示第二子带加权系数，其它表示同公式(14)。

${\mathrm{Er}}^{\&prime;}=E_2^{\&prime;}/E_1^{\&prime;}=\underset{f\&Element;\mathrm{sub}-\mathrm{band}2}{\mathrm{\&Sigma;}}{\left|S_N\left(f\right)\right|}^2*B_2\left(f\right)/\underset{f\&Element;\mathrm{sub}-\mathrm{band}1}{\mathrm{\&Sigma;}}{\left|S_N\left(f\right)\right|}^2*B_1\left(f\right)---\left(15\right)$

方案B借鉴梅尔滤波器的设计思想，以汉明窗形式给出加权系数，加权系数曲线如图1所示。

本发明以两个频段的能量比作为特征对OSAS进行筛查，这样，我们研究的不是某个单频点，而是整个频率段的信息。即使第一共振峰频率由于噪声的影响而出现偏差，只要偏差不是很大，仍可落在设定的频段内，算法仍然可以做出正确的判决。具有更强的鲁棒性和针对性，更适用于家庭等环境中OSAS筛查。

鼾声判决模块：

判决模块中，本发明采用线性分类器将上一模块中求得的子带能量比分类，高于门限值的即为OSAS患者，低于门限值的为单纯打鼾者。最终，利用ROC曲线和箱线图对本发明中提出两种方案的分类特性进行简要评估。AUC是评价诊断试验整体性能的定量指标，本发明提出两种方案的曲线下面积分别为：方案A的AUC＝0.829，方案B的AUC＝0.938。方案B优于方案A是因为，方案B考虑到第一共振峰在各频率点出现概率不同，在出现概率较高的中心频率附近赋予较大的权重系数，在出现概率较低的边缘频率点赋予较小的权重系数，使得改进方法可以灵活的修正鼾声第一共振峰的偏差，有效减小谐波频率处于子带边缘的噪声干扰。

经过以上处理本发明得到的筛查灵敏度为92.86％，特异性为92.74％，已基本满足OSAS筛查的实际应用要求，有望成为一种类似于常规体检的筛查方法。

图2(a)是正常鼾声和低通气鼾声的时域图，上图为单纯打鼾者鼾声时域曲线，下图为OSAS患者鼾声时域曲线。明显发现单纯打鼾者鼾声片断时域曲线表现为多个振幅、间隔大致相仿的复合波；OSAS患者鼾声的时域曲线表现为多个振幅、间隔不规则的复合波。图2(b)是正常鼾声和低通气鼾声的频域图，上图为单纯打鼾者鼾声LPC频谱包络图，下图为OSAS患者鼾声频谱包络图。频谱包络中最大幅值所对应的频率值即为鼾声的第一共振峰频率。图中可明显发现正常鼾声的能量集中在低频段，低通气鼾声能量集中在高频段。本发明通过提取鼾声的能量特性实现OSAS的筛查。

描述本发明工作流程的示例如图3所示：

1.鼾声的采集模块中，我们利用非接触式麦克风(距患者口鼻18-20cm)录制鼾声数据，并记录在计算机硬盘中。录制时设置信号的采样率为8k Hz，采样精度为16bit。

2.实验所用样本是患者整晚6-8小时鼾声数据，与语音相似，鼾声也是非平稳过程，采用语音信号处理的方法对其进行预处理。与语音不同的是鼾声音素单一、有声段无声段变化明显，预处理过程比正常语音简单。包括，取帧长1024点，帧移512点；加入汉明窗减小频谱泄露；端点检测采用短时能量检测方法。经过以上几步完成鼾声信号的预处理。

3.本发明利用AR(p)模型构建预处理后鼾声信号的能量谱如图2(b)所示，上图为正常鼾声的LPC谱包络，下图为低通气鼾声LPC谱包络。

4.能量比求解模块中，方案A中直接计算本发明中所选两频段内的子带内的能量值累加值E₁、E₂，然后计算高频子带与低频子带能量的比值E_r。同理，方案B中计算本发明中所选两频段内的子带能量谱的加权累加值E₁′、E₂′，其中加权系数曲线如图1所示，然后计算高频子带与低频子带能量的比值E_r′，最后将得到的能量比传递到判决部分中。

5.判决部分中，以高低频子带能量比值作为分类依据，分析本发明数据库中所有鼾声样本，找到最优的门限值作为分类值完成OSAS筛查，当患者能量比高于某门限值时，判为OSAS患者，反之为单纯打鼾者。

6.筛查效果评估部分，我们利用ROC曲线和箱线图对本发明中提出两种方案的分类特性进行简要评估。所述箱线图以能量比Er为特征，如图4所示，图4(a)为方案A(直接求和子带能量比法)，图4(b)为方案B(改进求和能量比法)。图中发现，方案B明显优于方案A。所述ROC曲线以1-特异度为横轴，灵敏度为纵轴，画出本发明提出的两种方案的ROC曲线如图5所示，并利用曲线下面积(AUC)对两种方案的筛查结果进行评定。

Claims (5)

1.一种家庭环境中睡眠呼吸暂停综合征的筛查系统，其特征在于：

包括鼾声获取模块：采集鼾声并记录鼾声数据，采样率设为8kHz、精度为16bit鼾声信号；

鼾声处理模块：

鼾声处理模块包括预处理部、端点检测部、能量比求解部，

其中预处理部，提供鼾声信号的分帧、预加重以及加窗过程，完成鼾声信号的预处理，

端点检测部，区分有声段和无声段，完成鼾声信号中鼾声段的检测，

能量比求解部，求解鼾声信号的子带能量比值，得到子带能量比值；

鼾声判决模块：

将上述鼾声处理模块中求得的子带能量比值与门限值比较；

其中鼾声处理模块的能量比求解部；

为了用数字信号处理方法对鼾声信号进行处理，首先需要建立鼾声信号产生的数学模型；其中上气道模型的建立和求解是通过线性预测编码LPC实现的，LPC分析过程如下所示：

当前值用s(n)表示，表示预测样值，则预测样值表示为：

其中a_i为线性预测系数，P为预测阶数，n为采样点，则预测误差e(n)定义为：

$e\left(n\right)=s\left(n\right)-\hat{s}\left(n\right)=s\left(n\right)-\underset{i=1}{\overset{P}{\mathrm{\&Sigma;}}}{a}_{i}s(n-i),\begin{array}{c}{a}_0=1\\ i=\mathrm{0,1},...P\end{array}---\left(5\right)$

按最小均方差准则确定预测系数，使得均方误差E[e²(n)]最小，E[e²(n)]定义为：

$E_n=\underset{n}{\mathrm{\&Sigma;}}{e}^2\left(n\right)=\underset{n}{\mathrm{\&Sigma;}}{[s\left(n\right)-\underset{i=1}{\overset{p}{\mathrm{\&Sigma;}}}{a}_{i}s(n-i)]}^2---\left(6\right)$

要求得最小E[e²(n)]，需要对E[e²(n)]关于a_i求偏导：

$\frac{\&PartialD;E\left(e^2\left(n\right)\right)}{\&PartialD;a_i}=2E\left[e\left(n\right)\frac{\&PartialD;e\left(n\right)}{\&PartialD;a_i}\right]=0,i=\mathrm{0,1},...P---\left(7\right)$

代入 $\frac{\&PartialD;e\left(n\right)}{\&PartialD;a_i}=s(n-i)$ 得：

E[e(n)s(n-i)]＝0      （8）

公式（5）代入公式（8）可得：

$E\left[e\left(n\right)s(n-i)\right]=E[s\left(n\right)s(n-i)-\underset{l=0}{\overset{P}{\mathrm{\&Sigma;}}}{a}_{l}s(n-l)s(n-i)]=0,\begin{array}{c}{a}_0=1\\ i=\mathrm{1,2},...P\end{array}---\left(9\right)$

因为信号自相关序列为：

$R(i-l)=E\left[s(n-l)s(n-i)\right],\begin{array}{c}{a}_0=1\\ l,i=\mathrm{1,2},...P\end{array}---\left(10\right)$

则公式（9）可写为：

$R\left(i\right)-\underset{l=0}{\overset{P}{\mathrm{\&Sigma;}}}{a}_{l}R(i-l)=0,\begin{array}{c}{a}_0=1\\ i=\mathrm{1,2},...P\end{array}---\left(11\right)$

公式（11）称为线性预测的标准方程，写成矩阵形式为：

$\left[\begin{array}{cccc}R\left(0\right)& R\left(1\right)& ...& R(P-1)\\ R\left(1\right)& R\left(0\right)& ...& R(P-2)\\ ...& ...& ...& ...\\ R(P-1)& R(P-2)& ...& R\left(0\right)\end{array}\right]\left[\begin{array}{c}{a}_1\\ a_2\\ .\\ .\\ .\\ a_P\end{array}\right]\left[\begin{array}{c}R\left(1\right)\\ R\left(2\right)\\ .\\ .\\ .\\ R\left(p\right)\end{array}\right]---\left(12\right)$

公式（11）就是尤利-沃克方程，解方程可得到p阶线性预测系数a_i和预测误差值，利用莱文森-德宾递推算法求解尤利-沃克方程，该算法流程如下：

步骤1：定义初始条件：

E₀＝R(0),a₁ ⁽⁰⁾＝1

步骤2：定义反射系数

$k_i=[R\left(i\right)-\underset{j=1}{\overset{i-1}{\mathrm{\&Sigma;}}}{a_j}^{(i-1)}R(i-j)]/E_{i-1}={a_i}^{\left(i\right)};$

a_j ⁽ⁱ⁾＝a_j ^(i-1)-k_ia_i-j ^(i-1);

E_i＝(1-k_i ²)E_i-1

本步骤从低阶开始递推，直到i＝p时截止，循环结束后给出各阶次的所有参数；

步骤3:线性预测的系数：

a_j＝a_j ^(p)  1≤j≤p

经过以上过程，求得线性预测器的预测系数a_i，并建立p阶全极点模型AR(p)构建时变数字滤波器H(z)模拟上气道结构，其中U(z)为激励信号，S(z)为输出信号：

$H\left(z\right)=\frac{S\left(z\right)}{U\left(z\right)}=\frac{G}{1-\underset{i=1}{\overset{p}{\mathrm{\&Sigma;}}}{a}_i{z}^{-i}}---\left(13\right)$

通过上述模型求出上气道的LPC对数谱，鼾声受上气道的共振作用，呼吸气流的不同频率能量重新分配，其中最高能量对应频率即为第一共振峰；

选取的两通频带范围为：第一子带200-350Hz，第二子带500-650Hz，第一子带对应正常呼吸时鼾声第一共振峰集中频段，第二子带对应发生低通气时鼾声第一共振峰集中频段。

2.根据权利要求1所述的家庭环境中睡眠呼吸暂停综合征的筛查系统，其特征在于：其中鼾声处理模块的预处理部

在处理之前需要对鼾声信号分帧，帧长取128ms，帧移为64ms，利用一阶FIR滤波器H(z)＝1-0.98z^-1实现预加重，在对鼾声信号进行频域处理之前使用如公式（1）的窗函数来减小频谱泄露，其中n为采样点：

$w\left(n\right)=\left\{\begin{array}{cc}0.54-0.46*\cos (n-1)*2*\mathrm{\&pi;}/512& 0<n\&le;256\\ 1& 256<n\&le;768\\ 0.54-0.46*\cos (n-1)*2*\mathrm{\&pi;}/512& 769<n\&le;1024\end{array}\right.---\left(1\right).$

3.根据权利要求1所述的家庭环境中睡眠呼吸暂停综合征的筛查系统，其特征在于：其中鼾声处理模块的端点检测部

采用短时能量和短时过零率双门限分析法，具体过程如下：

1）计算每帧的短时能量：

$M_m=\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{\left|x_m\left(n\right)\right|}^2---\left(2\right)$

其中，N为帧长1024，m为帧号，n为采样点，M_m为短时能量，x_m(n)为鼾声时域波形；

2）假定前10帧为无声帧，计算各帧的短时能量和短时过零率，统计前10帧噪声的平均过零率Mzc，平均过零率方差Vzc，平均短时能量E；

确定门限值，能量门限值如公式（3），其中，EL为低门限值，EH为高门限值，平均短时能量E：

$\left\{\begin{array}{c}\mathrm{EL}=2*E\\ \mathrm{EH}=5*\mathrm{EL}\end{array}\right.---\left(3\right)$

过零率门限zc如公式（4），其中前10帧噪声的平均过零率Mzc，平均过零率方差Vzc：

zc＝Mzc+2*Vzc      （4）

3）先根据能量门限EL和EH确定出鼾声段的初始起止点[N₁,N₂]，再利用过零率门限zc修正并获得最终起止点[N_B,N_E]，N_B为鼾声起始点，N_E为鼾声终止点，至此，完成鼾声信号的预处理及端点检测过程。

4.根据权利要求1所述的家庭环境中睡眠呼吸暂停综合征的筛查系统，特征在于：所述能量比求解部采用直接求和子带能量比方法，该方法具体如下：

确定低通气状态和正常状态下鼾声第一共振峰对应的频段后，利用AR模型构建预处理后鼾声信号的能量谱，分别计算所选两子带频段内频谱能量，即频域中实现带内能量谱累加求和，第一子带内能量为E₁，第二子带内能量为E₂，然后取E₂与E₁的比值，直接求和子带能量比如公式（14），其中f表示频率，N为帧长，sub-band1表示第一子带频率范围，sub-band2表示第二子带频率范围，S_N(f)表示信号的短时LPC谱包络。

$\mathrm{Er}=E_2/E_1=\underset{f\&Element;\mathrm{sub}-\mathrm{band}2}{\mathrm{\&Sigma;}}{\left|S_N\left(f\right)\right|}^2/\underset{f\&Element;\mathrm{sub}-\mathrm{band}1}{\mathrm{\&Sigma;}}{\left|S_N\left(f\right)\right|}^2---\left(14\right)$

5.根据权利要求4所述的家庭环境中睡眠呼吸暂停综合征的筛查系统，特征在于：改进的求和子带能量比方法如下：

对直接求和能量比方法进行改进，通过加权系数对带内鼾声能量谱加权求和，其中第一子带加权能量和为E₁′，第二子带加权能量和为E₂′，加权求和子带能量比如公式（15），其中B₁(f)表示第一子带加权系数，B₂(f)表示第二子带加权系数，其它表示同公式（14）。

${\mathrm{Er}}^{\&prime;}={E_2}^{\&prime;}/{E_1}^{\&prime;}=\underset{f\&Element;\mathrm{sub}-\mathrm{band}2}{\mathrm{\&Sigma;}}{\left|S_N\left(f\right)\right|}^2*B_2\left(f\right)/\underset{f\&Element;\mathrm{sub}-\mathrm{band}1}{\mathrm{\&Sigma;}}{\left|S_N\left(f\right)\right|}^2*B_1\left(f\right)---\left(15\right)$

Images