基于递归图分析的腭裂语音鼻漏气的识别方法与流程

本发明涉及语音信号处理领域，尤其是一种基于递归图分析的腭裂语音鼻漏气的识别方法。

背景技术：

先天性腭裂畸形最明显的功能障碍之一是语音障碍，由于腭裂丧失了口鼻腔之间正常的骨性和肌性分离，导致患者在发音过程中无法控制气流的大小和方向。腭裂患者由于腭咽闭合不全使发音时腭咽瓣不能完全关闭，气流同时经过了鼻腔和口腔，产生异常共振、鼻漏气、高鼻音等典型的腭裂语音。发音者在发辅音时，有气流从鼻腔漏出，这种气流通过鼻腔的声音有时在发音过程中可以被听见就称为鼻漏气。鼻漏气主要出现于塞音、塞擦音及擦音等压力性辅音，其音质发生变化。

近年来，国内外学者致力于研究更为客观的数字化语音处理与识别方法，但现有技术中较少有用于鼻漏气语音的诊断与评估的研究，并且缺乏对鼻漏气语音的有效自动识别检测。

技术实现要素：

本发明针对现有技术存在的上述技术问题，提供一种基于递归图分析的腭裂语音鼻漏气的识别方法，利用递归图趋势分析得出的改进参数，再针对递归图做区域分块处理得到最小区域矩阵计算相应参数，以及针对递归图做图像分析处理得到图像矩阵，将这些参数作为语音信号的特征参数，采用分类器实现了对腭裂语音鼻漏气的自动识别。

本发明采用以下技术方案：

一种基于递归图分析的腭裂语音鼻漏气的识别方法，包括以下步骤：

(1)语音信号预处理，对输入的语音信号进行降采样、归一化分帧、预加重及幅值归一化；

(2)对预处理后的语音信号求取语音信号的递归图矩阵，将系统中的时间序列提取出来，重现递归信号，选择合适的嵌入维数m和延迟时间τ将一维的非线性时间序列：

{s(i)，i＝1，2，...}(2)

重构出向量:

si＝[s(i)，s(i+τ)，...，s(i+(m-1)τ)](3)，

系统的m维相空间轨道由有时间标记的向量序列构成：

{si，i＝1，2，...，n}(4)

然后用这些相空间上的点作为行和列构成n×n的矩阵递归图；

图中的每个节点ri，j由对应的行、列向量点之间的距离来表示：

ri，j＝θ(ε-||si-sj||)，i，j＝1，2，...，n(5)，

其中ε是根据先前情况设定的阈值常数，用来表示临界距离；符号||*||表示取向量的欧几里德(euclidean)范数；θ(x)是核函数；若ri，j的值为1，则对应递归图中位置点(i,j)就表示为一个黑点，若ri，j的值为0，则对应递归图中位置点(i,j)就表示为一个白点，得到二维递归矩阵rp；

(3)针对递归图做趋势分析，采用定量递归分析的5种针对递归图分析的量化参数，包括确定性det、递归度rrl、最长对角线ll、熵entr和趋势分析rt，将对角线在递归图中分布的概率密度，对应在频域上转化为对角线特征llf函数，并将频域上每一个对角线下区域对应的熵值转化为频域熵值pentr；

(4)针对递归图直接做区域分块处理，对信号递归图做矩阵处理，求取行列值，选取合适的分块系数，根据统计特性来确定分块系数，将递归图分为相应区域后，提取出递归图四个角的最小区域部分矩阵如下：

其中rp为递归矩阵，m为所求的分块系数，分块系数m的计算即将整个递归图分为m*m个部分，根据矩阵所有点数总和，按每个区域1000个点划分，计算公式为：

提取出递归图四个角(左上、左下、右上、右下)的递归矩阵rpa，每个分块后小矩阵的标号用a表示；此时四个角为最小区域部分的矩阵作为特征参数，再针对每一个最小区域部分做处理，计算区域均方值，模拟当一维信号转换到二维图像上分布密集区域的权值，计算公式如下：

此处矩阵点ri，j∈rpa，na表示该最小区域矩阵的点数总和；

(5)对递归图进行图像处理，将格式转换后的图像矩阵，依次进行二值化处理和两次特定结构模板的滤波处理后得到处理后的矩阵信号特征用于分类识别；

(6)利用分类器对信号进行分类识别，得到自动识别结果。

所述步骤3具体包括以下步骤：

(1.1)语音信号降采样，按照采样频率44100降采样到15000或8000的比例对数据进行处理；

(1.2)归一化：对每一个语音数据求取绝对值最大的幅值max，然后将语音信号中的每一个值分别除以该最大值，进行归一化，相应的计算公式为：

s＝xx/max|xx|(1)

其中xx为输入的语音信号，s为处理后的序列。

所述步骤3中的量化参数具体为：

由于递归图关于主对角线对称，则右下部分等间距平行于主对角线，用rrl表示递归图第l区域中递归度，计算公式如下：

而趋势分析rt用递归度变化速率来表征，计算公式如下：

其中，为{rrl，l＝1，2，...，k}的平均值，因此rt也是{rrl，l＝1，2，...，k}序列关于{l＝1，2，...，k}的线性回归斜率，斜率越小就表示该系统越稳定，而斜率越大也代表了系统的非稳定性或动力学特性的突变；

确定性det的计算公式如下：

其中p(l)表示对角线在信号递归图中分布的概率密度，p(l)计算公式如下：

lmin是对角线长度的统计初值，满足2≤lmin≤nl-1，p表示概率计算，由于信号长度不同，归一化后对最长对角线ll计算公式如下：

其中lmax是对应lmin的值，由于lmin是对角线长度的统计初值，因此lmax是统计长度中对角线的最大值，而ll就是对对角线最大值进行归一化；

熵entr是用来描述信号递归图的复杂程度，其复杂程度越高则熵值越大，假设递归图中所有线长都为同一个值，则熵为0，根据语音熵计算公式：

所述步骤3中的将对角线在递归图中分布的概率密度，对应在频域上转化为对角线特征llf函数，并将频域上每一个对角线下区域对应的熵值转化为频域熵值pentr具体为：

将对角线在递归图中分布的概率密度，对应在频域上转化为对角线特征llf函数，计算公式如下：

其中，β为修正因子，防止超出频域计算范围，根据实际情况选定，

将频域上每一个对角线下区域对应的熵值转化为频域熵值pentr，计算公式如下：

所述步骤5具体为：

(5.1)令rp(x,y)表示了大小为n×n像素的rp递归图矩阵，其中x＝0，1，2，...，n-1，y＝0，1，2，...，n-1；

(5.2)快速傅里叶变换：

令rp(u,v)表示的rp(x,y)的二维离散傅里叶变换，计算公式为：

式中，u＝0，1，2，...，nl-1，v＝0，1，2，...，nl-1；

(5.3)特定模板滤波处理过程：

对于快速傅里叶变换后的递归图rp(u,v)，采用自适应法确定二值化的阈值，对图像二值化；

对图像二值化后进行第一次滤波，选取第一特定结构元进行开操作；

然后进行第二次滤波进行模糊图像，去掉一些细节，最后选取第二特定结构元进行闭操作，去除不关心的非密集成分，得到最终处理后的图像矩阵结果，并放入分类器进行特征区分。

所述步骤6具体为：

(6.1)计算距离

根据给定的测试对象，计算它与训练集中的每个对象的距离，通常计算欧式距离d1(x，y)和曼哈顿距离d2(x，y)，

(6.2)找邻近对象，按照距离的递增关系进行排序。圈定距离最近的k个训练对象，作为测试对象的近邻；

(6.3)确定前k个点所在类别的出现频率；

(6.4)返回前k个点中出现频率最高的类别作为测试数据的预测分类。

综上所述，由于采用了上述技术方案，本发明的有益效果是：

1、本发明提供的基于递归图分析的腭裂语音鼻漏气的识别方法，可以基于获取的腭裂语音鼻漏气语音信号进行全自动检测，利用语音信号处理的方法提取信号特征，相比于现有技术中，本发明的检测方法检测结果客观准确，实现较高程度的自动测量；

2、本发明提供的基于递归图分析的腭裂语音鼻漏气的识别方法，提出了基于通过递归图提取相应的递归趋势分析特征，该递归特征能更好反映语音信号由于发声机理的改变而导致的更加细微的数据差异，又将该差异转换到频域做进一步处理，利用特定的模板滤波得到相应的图像矩阵特征。该特征的提取基于人类听觉特性，在突出鼻漏气语音及正常语音细微差异的同时提高了运算速度。

3、本发明提供的基于递归图分析的腭裂语音鼻漏气的识别方法，基于在递归图的频域中提取递归矩阵特征的同时，还根据信号的特点选取合适的特定模板做进一步的处理，得到的处理后的特征能更好地体现语音之间的细微差别。

4、本发明提供的基于递归图分析的腭裂语音鼻漏气的识别方法，能够更好符合语音信号特性，因为递归图分析是基于非线性动力学的分析方法，而语音信号是非线性信号，递归图能对其有很好的表征。而在非线性动力学的表现上，人声和空气声有明显差异，故而能有效表征两者之间的特征差异。

附图说明

图1是本发明实施例提供的腭裂鼻漏气自动检测方法框图。

图2是正常语音和有鼻漏气语音的信号递归图。

图3是转换后的频域熵值pentr示意图。

图4是对转换后的频域熵值pentr示意图各个部分的截取图。

图5是对正常语音和有鼻漏气语音信号递归图的观察示意图。

图6是递归图的分块矩阵示意图。

图7是对递归图进行图像分析处理的处理流程图。

图8是正常语音和有鼻漏气语音的递归图变换后的示意图。

图9是对变换后的图像做滤波处理的流程图。

图10是二值化后的正常语音和有鼻漏气语音示意图。

图11是第一次滤波后的正常语音和有鼻漏气语音示意图。

图12是第二次滤波后的正常语音和有鼻漏气语音示意图。

具体实施方式

为了使本领域的人员更好地理解本发明的技术方案，下面结合本发明的附图，对本发明的技术方案进行清楚、完整的描述，基于本申请中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的其它类同实施例，都应当属于本申请保护的范围。

需要说明的是，在本实施例中采用的语音数据是按照普通话构音测量表进行录制得到，语音数据前期由专业的语音师进行判听，由人工对鼻漏气数据进行详尽的标注，根据鼻漏气发音特点，选取含辅音/c//ch//d//f//k//g//j//s//sh//t//x//p//q//z//zh/等的字样本。

进一步需要说明的是，本发明以及相应实施例提出的信号特征处理方法、特征提取方法、信号识别与分类方法都仅仅是对信号的处理和识别方法本身进行研究和改进，虽然针对的为鼻漏气信号，实现的自动分类识别结果可以作为评估参考，但在临床或医疗领域其评估结果也仅仅是一个辅助性的评估，对于具体的治疗方法仍需要并主要依赖于医生的临床经验和医生提供的治疗方法。

如图1所示，本实施例提供的基于递归图分析的腭裂语音鼻漏气的识别方法如下：

(1)语音信号预处理，对输入的语音信号进行降采样及幅值归一化：

(1.1)语音信号降采样，按照采样频率44100降采样到15000或8000的比例对数据进行处理；

(1.2)归一化：语音信号采集的时候，由于设备差异会存在信号采集条件不统一，也为了避免声音音量大小不一，噪声等问题，对每一个语音数据求取绝对值最大的幅值max，然后将语音信号中的每一个值分别除以该最大值，进行归一化，相应的计算公式为：

s＝xx/max|xx|(1)

其中xx为输入的语音信号，s为处理后得到的序列。

(2)求取语音信号的递归图矩阵，对于混沌系统将一种非稳态轨道形成的递归现象作为其基本特征。由于混沌吸引子在它的相空间中体积是有限的，因此构成吸引子的非稳态轨道在有限的吸引子空间中不断近似逼近又分叉远离，递归图法就是运用系统动力学的原理，将系统中的时间序列提取出来，重现递归信号，选择合适的嵌入维数m和延迟时间τ就可以将一维的非线性时间序列：

{s(i)，i＝1，2，...}(2)

重构出向量:

si＝[s(i)，s(i+τ)，...，s(i+(m-1)τ)](3)，

系统的m维相空间轨道由有时间标记的向量序列构成：

{si，i＝1，2，...，n}(4)

然后用这些相空间上的点作为行和列构成了n×n的矩阵递归图。图中的每个节点ri，j由对应的行、列向量点之间的距离来表示：

ri，j＝θ(ε-||si-sj||)，i，j＝1，2，...，n(5)

上式中ε是根据先前情况设定的阈值常数，用来表示临界距离。符号||*||表示取向量的欧几里德(euclidean)范数。θ(x)是核函数。若ri，j的值为1，则对应递归图中位置点(i,j)就表示为一个黑点，同理，若ri，j的值为0，则对应递归图中位置点(i,j)就表示为一个白点。所以，递归图表示将一个m维相空间的轨道关系映射到一个二维图上。用rp表示最后得到的二维矩阵，也就是递归矩阵。因为语音信号为一维信号，则m＝1。如图2所示图中分别为正常语音/sha/和有鼻漏气语音/sha/的递归图。

(3)针对递归图做趋势分析，由于递归图反映了动力学系统相空间轨迹运行关系的二维映射图，为了量化递归图中表现出来的系统递归现象，采用定量递归分析(recurrencequantificationanalysis，roa)的5种针对递归图分析的量化参数，包括确定性det、递归度rrl、对角线特征llf、频域熵entr和趋势分析rt。不同的roa参数描述了系统不同的动力学行为。例如，因为确定性描述了动力学轨迹周期递归的程度，那么确定性就与随机性相对，值越大就表明确定性强，相反则表示随机性强。

因为递归图关于主对角线对称，则右下部分等间距平行于主对角线，用rrl表示递归图第l区域中递归度，计算公式如下：

其中，nl＝n×n；

由于递归图中的点分布反应了语音信号中语音段和非语音段的内在结构。因为语音信号的动力学存在差异，语音段在递归图中的点密集比例会比非语音段小。因此将递归度rrl作为主要区分参数。鼻漏气在主观听觉感知上表现为，气流声音代替了正常的发声，故鼻漏气语音会在辅音阶段存在一段气流声，以此作为和语音段的区分。

而趋势分析rt用递归度变化速率来表征，计算公式如下：

其中，为{rrl，l＝1，2，...，k}的平均值，因此rt也表示{rrl，l＝1，2，...，k}序列关于{l＝1，2，...，k}的线性回归斜率，斜率越小就表示该系统越稳定，而斜率越大也代表了系统的非稳定性或动力学特性的突变。

确定性det的计算公式如下：

其中p(l)表示对角线在信号递归图中分布的概率密度，p(l)计算公式如下：

其中lmin是对角线长度的统计初值，满足2≤lmin≤nl-1。p表示概率计算。由于本实施例所分析信号长度不同，归一化后对最长对角线ll计算公式如下：

其中lmax是对应lmin的值，由于lmin是对角线长度的统计初值，因此lmax是统计长度中对角线的最大值，而ll就是对对角线最大值进行归一化。

本实施例改进p(l)计算方法的原因在于一般采用的算法是根据实验信号特性设计的，故在此处要根据语音信号的特点来对其算法改进，而语音信号递归图呈现一定对称性，所以用对角线为主要分割进行概率密度的计算。

另一方面，由于目前针对于递归图的趋势分析仅在时域上，本实施例通过研究发现在频域上的递归分析也有一定作用，并做进一步的具体分析如下。

将对角线在递归图中分布的概率密度，对应在频域上可以转化为对角线特征llf函数，计算公式如下：

其中，β为修正因子，防止超出频域计算范围，根据实际情况选定，

熵entr是用来描述信号递归图的复杂程度，其复杂程度越高则熵值越大，假设递归图中所有线长都为同一个值，则熵为0，根据语音熵计算公式：

再将其映射到频域上，将频域上每一个对角线下区域对应的熵值转化为频域熵值pentr，计算公式如下：

如图3所示为频域熵值pentr的整体图，在图3中表征为更加细微的数据点，进一步将图3的各个部分放大截图得到如图4所示，因此频域熵值pentr表征了在频域下的点密集程度。最后将改进后的确定性det、递归度rrl、对角线特征llf、频域熵entr和趋势分析rt五个参数作为鼻漏气语音的信号特征。

(4)针对递归图直接做区域分块处理，通过对有无鼻漏气语音递归图的观察，发现正常语音递归图主要分为四个区域，而稳定性良好且熵值较小，如图5的左图所示。有鼻漏气语音递归图主要分为八个区域，熵值较正常语音大，在递归图四个角区域存在一些密集的点分布，如图5的右图所示。

对信号递归图做矩阵处理，求取行列值，选取合适的分块系数。根据统计特性来确定分块系数，将递归图分为相应区域后，提取出递归图四个角的最小区域部分矩阵。分块系数m的计算是通过将整个递归图m*m个部分，根据矩阵所有点数总和，假设其均匀分布，按每个区域1000个点划分，计算公式为：

故用公式来表示分块后的区域矩阵和原递归矩阵的关系如下：

其中rp为递归矩阵。

提取出递归图四个角(左上、左下、右上、右下)的递归矩阵rpa，每个分块后小矩阵的标号用a表示，如图6所示。此时，四个角为最小区域部分的矩阵作为特征参数。有鼻漏气语音的最小区域部分矩阵比正常语音点分布密集，可以作为有效的信号特征来判断语音中鼻漏气是否存在。

再针对每一个最小区域部分做处理，计算区域均方值，模拟了当一维信号转换到二维图像上分布密集区域的权值，计算公式如下：

此处矩阵点ri，j∈rpa，na表示该最小区域矩阵的点数总和。

通过上式的计算，我们可以得到四个最小区域部分矩阵对应的权值ea值，将其作为鼻漏气语音的特征参数。

(5)对递归图进行图像分析处理，处理流程如图7所示，

(5.1)将rp递归图矩阵以unit8格式存储，读入matlab时默认jpg格式，则此时令rp(x,y)表示了一副大小为n×n像素的数字图像，其中x＝0，1，2，...，n-1，y＝0，1，2，...，n-1。

(5.2)快速傅里叶变换：

令rp(u,v)表示的rp(x,y)的二维离散傅里叶变换，计算公式为：

式中，u＝0，1，2，...，nl-1，v＝0，1，2，...，nl-1。使用确定频率变量u和v可以得到x和y，也可以将指数项展开为正弦函数和余弦函数。

变换后的图像如图8所示，通过观察发现，有鼻漏气语音中央区域点更密集，且存在十字空白区，表现了有鼻漏气语音在空气动力学上存在语音信号与非语音信号的差异，也就是鼻腔空气流动声音代替了发音。因此从频域上证明了时域递归分析的结果，而且频域上更能明显反应差异之处。

(5.3)特定模板滤波处理过程：

如图9所示，将做完快速傅里叶变换的递归图rp(u,v)进行格式转换，保存为图片格式读入matlab中，由于递归图点分布有些地方较密集，有些地方较稀疏，故采用自适应法确定二值化的阈值，对图像二值化。结果如图10所示。

对图像二值化后进行第一次滤波，该第一次滤波为中值滤波，选取特定结构元进行开操作，

综合递归图的点特性，选取b为结构元可以有效的将密集点分布聚集，合成一块区域，更易于区分特征，结果如图11所示。

然后进行第二次滤波，该第二次滤波为高斯滤波，模糊图像，去掉一些细节，最后选取特定结构元(半径15的圆盘结构)进行闭操作，去除不关心的非密集成分，得到最终结果如图12所示。

通过处理结果明显发现，有鼻漏气语音在图像中央部分，存在异于正常语音的点密集分布区域，在经过特定滤波变换后，更加突出明显。

(6)构建鼻漏气语音信号的特征参数矩阵：

将步骤(3)(4)的一行参数进行拼接，即将步骤(3)计算后的五个递归趋势分析参数和步骤(4)中的四个最小区域矩阵的ea值拼接为单行的参数矩阵。对于步骤(5)中的处理，在处理前将unit格式的矩阵转换成jpg格式图像，在处理后再将其转换成矩阵格式，选取统一的截取标准，即对所有语音数据进行分析后，取其中点数的最小值就作为截取标准。假设为m，则我们可以得到一个m*m的矩阵，对于该矩阵第一行数据m个，前9个为趋势分析参数和ea，其余m-9值用0补齐，最终得到m+1行×m列的特征矩阵，将其放入分类器进行识别。

(7)利用分类器对信号进行分类识别

knn分类器是依据knn算法又称为k近邻分类(k-nearestneighborclassification)算法实现的。knn是通过测量不同特征值之间的距离进行分类的。如果一个样本在特征空间中k个最相似(即特征空间中最邻近)的样本中的大多数属于某一个类别，则该样本也属于这个类别。k通常是不大于20的整数。knn算法中，所选择的邻近都是已正确分类的对象。该方法在定类决策上只依据最邻近的一个或者几个样本的类别来决定待分样本所属的类别。具体实现步骤如下：

(7.1)计算距离

根据给定的测试对象，计算它与训练集中的每个对象的距离，通常计算欧式距离d1(x，y)和曼哈顿距离d2(x，y)。

(7.2)找邻近对象：按照距离的递增关系进行排序。圈定距离最近的k个训练对象，作为测试对象的近邻。

(7.3)确定前k个点所在类别的出现频率；

(7.4)返回前k个点中出现频率最高的类别作为测试数据的预测分类。

(8)通过以上处理步骤可完成本实施例的鼻漏气语音的自动分类与识别，进一步验证该方法得到相应的实验结果如下：

根据变量统一原则，和前文所对应，采样点数n，嵌入维度m，延迟时间τ，临界距离ε。

其中1类为鼻漏气语音，2类为无鼻漏气语音。

腭裂语音鼻漏气自动识别正确率最高可达84.63％。分别从降采样点、延迟时间、临界距离三个因素对识别正确率的影响做分析，由于以上三个因素对腭裂语音信号递归图有决定性作用，降采样点反映了递归图的分辨率，延迟时间和临界距离反映了递归图内部规律，并依据这两个参数生成符合某种特定运动规律的腭裂语音信号递归图。递归图是本文整个算法分析的根本，通过对三个主要因素的讨论，逐步确定最佳取值应用在下一个因素中。

(1)降采样点的影响

设降采样点n(个)，延迟时间τ(ms)，临界距离ε(单位)。设置延迟时间τ＝3，临界距离ε＝5。

表1降采样点数对腭裂语音鼻漏气识别的影响

(2)延迟时间的影响

设置延迟时间分别在n＝8000和n＝15000两种情况下，临界距离ε＝5。

表2延迟时间对腭裂语音鼻漏气识别的影响

(3)临界距离的影响

设置延迟时间分别为n＝8000和n＝15000两种情况，延迟时间τ＝3。

表3临界距离对腭裂语音鼻漏气识别的影响

(4)语音单元的影响

通过对上述三个因素的讨论，只做了一个语音样本对应一个递归图矩阵，产生五个递归参数、四个最小区域矩阵的特征参数。由于一个递归图矩阵表征信息有限，则假设语音单元也会对实验结果产生影响。对语音信号进行分帧，本文采用hamming窗，帧长为200ms，帧移为40ms。选取作为特征参数的帧数为fn。设置降采样点数n＝30000，因为分帧会降低数据计算量，故由表1可知，数据点数越多，鼻漏气识别正确率越高。根据表2、表3，选取延迟时间τ＝1，3，临界距离ε＝5，7，分别进行测试，得到如下结果：

表4语音单元帧数对腭裂语音鼻漏气识别的影响