申玉静，王 寻，唐 闽，吴 伟

（1.中国医学科学院阜外医院心律失常中心，北京100037；2.上海电机学院智能制造（中德）学院，上海201306；3.中国科学院声学研究所语言声学与内容理解重点实验室，北京100190）

0 引言

近年来，随着我国城市化进程的加快以及人口老龄化问题的日趋严峻，心血管疾病患者日益增多，因心脏疾病而死亡的病例也逐年增加。尽早发现心脏异常有利于降低心血管疾病的死亡率[1]。目前，存在磁共振成像[2]、计算机断层显像[3]和超声成像[4]等多种心脏检查方法，但这些方法所需设备费用高，操作复杂，且需要专业人员才能操作，难以在基层及社区推广。心电检查所用设备虽然相对简单[5]，对人体损害很小，但难以诊断瓣膜疾病。心音是心脏跳动过程中由瓣膜开闭和血流冲击所产生的，能够反映心脏部分组织的健康状况。当心脏瓣膜发生病变后，心音中往往会有杂音出现。同时，采集心音无需昂贵的设备，这使得心音听诊长期以来被用于诊断心脏瓣膜疾病。一个心动周期通常可以分为4个部分：第一心音（S1）、第二心音（S2）、S1～S2间期和S2～next S1间期[6]，如图1所示。其中，S1和S1～S2属于收缩期，S2和S2～next S1属于舒张期[7]。无论是医生的听诊教学，还是依据心音进行自动诊断系统的模型训练，都需要大量的心音数据。然而一些研究条件欠缺的机构难以获取大量不同类型的心音，因此不利于基于心音的心脏疾病诊断方法的推广。

图1 一个完整的心动周期

为了解决上述问题，很多研究人员提出了合成心音的方法。如高佼[8]使用模拟电路模型类比心血管系统中的腔室、动脉、静脉以及心室收缩和外周阻力的生理控制机制，计算得出各瓣膜两侧血压差和瓣膜关闭时刻。然后将瓣膜振动类比为弹簧振子的阻尼受迫振动，将此前由电路模型计算出的瓣膜两侧血压差用Ramp函数拟合，求解瓣膜振动方程，进而得到瓣膜振动波形，并结合瓣膜关闭时刻合成S1与S2。姜斌[9]设计了基于心脏生理参数的复合心音合成模型，其中包括心脏仿真和心音合成2个子模块，前者为后者提供运行参数。心音合成子模块又分为2个部分，分别产生左侧心音和右侧心音。通过将心脏活动与电路类比，建立心脏活动方程，并使用四阶Runge-Kunta算法求解，得到左、右侧心音轨迹图和波形图，并在时域上将两者叠加，从而得到完整的合成心音时域波形。上述模型均基于解析形式的数学公式，虽然可以模拟出S1与S2，但由于生物信号较为复杂，这些模型所合成的心音与真实心音仍然存在一定的差距。此外，这些模型难以模拟各类心脏杂音，通用性有限。

目前，很多研究机构致力于发展基于大数据的医疗软件平台[10]。本文通过构建心音数据库，并以此为基础设计心音合成软件，为心音听诊的教学和自动诊断中的数据生成提供便利。

1 心音数据采集

心音合成软件需要建立在大量心音数据的基础上，因此心音采集是非常重要的环节。医学上常用的听诊区有5个，包括主动脉瓣听诊区、肺动脉瓣听诊区、主动脉瓣第二听诊区、三尖瓣听诊区和二尖瓣听诊区[11]。常见的瓣膜疾病有主动脉瓣狭窄、主动脉瓣返流、二尖瓣狭窄、二尖瓣返流和三尖瓣返流5种。本研究使用长虹九心格电子听诊器（NC-STH-002）采集了184位志愿者5个听诊区的心音，每个听诊区采集时长1 min左右。数据分布详见表1。

表1 采集的心音数据分布

所有采集的心音均存储为4 kHz采样率、16 bit的.wav格式音频文件。采集时由电子听诊器通过蓝牙将心音发送至搭载Android系统的平板计算机内保存。采集完成后，再通过蓝牙将全部心音数据发送至台式计算机存储。

2 数据处理

2.1 心音降噪

本研究中，电子听诊器使用压电陶瓷传感器将心音振动转换为电信号，这些电信号中往往包含环境噪声和电噪声[12]。在将采集的心音数据收录进数据库之前，需要先对心音数据进行降噪处理。目前，常用的心音降噪方法有谱减法[13]和小波阈值降噪[14]等。本研究使用最优改进对数谱幅度（optimally modified log-spectral amplitude，OMLSA）估计方法对心音进行降噪。OMLSA是2002年Cohen等研究人员提出的降噪方法[15-18]，目标是最小化纯净音频与估计的降噪后音频之间的差异。降噪后音频频谱可表示为

式中，Y（k，l）表示l时刻k频率处的频谱幅度值；GH1和Gmin分别表示有效信号存在和不存在时的增益；p（k，l）为有效信号存在的后验概率。p（k，l）的计算公式如下：

式中，q（k，l）表示第l帧频率为k的频带上有效信号不存在的概率；ξ（k，l）为音频先验信噪比。ξ（k，l）的表达式为

式中，α为权重因子常数，用于控制降噪程度和有效信号保留之间的平衡，本文设定α=0.9；γ（k，l）为音频后验信噪比；GH1（k，l-1）为上一帧的对数谱增益函数。

γ（k，l）的表达式为式中，λd（k，l）为音频信号的噪声谱估计，是由改进的最小值控制递归平均（improved minima controlled recursive averaging，IMCRA）算法[19]计算的。v（k，l）的表达式为

同时，GH1满足[20]：

本课题组此前的研究表明，在OMLSA方法降噪的基础上，进一步使用小波阈值方法降噪，总体上可以获得更好的降噪效果[21]。因此，本研究结合OMLSA方法与小波阈值方法对采集到的心音降噪，并将降噪后心音存放在数据库中。小波阈值降噪时选用coif5小波，分解层数为6层，使用软阈值方式，阈值选择准则为rigrsure。降噪前后的心音波形如图2所示。

2.2 数据库建立

图2 主动脉瓣狭窄患者降噪前后的心音波形

医学常用的数据库包括Qracle[22]、SQL Server[23]和Access等。由于目前本研究所用心音均存储在本地，因此选用Office 2013中的Access数据库存储心音信息。新建心音数据库并在数据库中建立4个表格，依次命名为S1、S12、S2和S21，分别用于存储心音4个部分的信息。图3显示了S1表格中存储的部分心音信息。其中，“听诊区”属性下1～5分别表示该编号的心音采集于主动脉瓣听诊区、肺动脉瓣听诊区、主动脉瓣第二听诊区、主动脉瓣第三听诊区和二尖瓣听诊区。“健康状况”属性下序号1～6与表1中序号1～6逐一对应。如“健康状况”属性为2，则表示该心音采集于主动脉瓣狭窄的患者。

图3 心音数据库中S1表格的部分心音信息

3 软件设计

MATLAB GUI（graphical user interface）是一种常用的用户可视化交互工具，已经被应用于语音信号处理[24]、地球科学[25]和统计[26]等领域的软件设计。它具有控件丰富、自由度高等优点，可供研究人员快速创建应用程序[27]，因此本研究选用MATLAB 2014b的GUI设计心音合成软件。

如图4所示，软件分为心音数据选择、心音合并以及查看文件和软件信息3个部分。数据选择主要有4个部分，分别用于选择S1、S1～S2、S2和S2～next S1的心音。以S1为例说明心音数据选择的操作方法，其数据选择与波形预览界面如图5所示。首先需要在时长、听诊区和健康状况的Edit控件中输入想要查找的条件，然后点击“查找”按钮查询，查询到的心音显示在表格中。点击表格任意一行，会在下方绘图区域显示被选中心音的波形。其他心音成分选择方法与此类似。

图4 心音合成软件界面

当逐个选择4个部分的心音并点击“添加”按钮后，需要分别输入S1～S2、S2和S2～next S1心音的开始时间以及4个部分心音的放大倍数，合成的波形和时频图如图6所示。由于软件中默认将S1设置为从零时刻开始，所以不需要手动设置S1的开始时间。输入完成后点击“合并”按钮，即可将所选4个部分的心音合并为一个完整心动周期的心音。如果对心音波形较为满意，可以点击“播放”试听。在“持续周期数”后的文本框中输入想要生成的心音波形的周期数，然后点击“生成”按钮，即可生成多个周期的心音.wav文件。当生成心音文件后，点击图4中“查看生成文件”按钮即可打开生成的心音文件所在文件夹。此文件夹中保存了所有生成的心音文件，用户可以随时查看和试听。

图5 S1心音数据选择与波形预览界面

图6 心音波形生成界面

4软件实现

4.1 心音数据选择

心音数据选择部分主要功能为选择4个部分的心音。本节以S1的选择和显示为例，叙述程序设计。先在时长、听诊区和健康状况的Edit控件中输入文本作为查找信息，然后点击“查找”Button控件在Access数据库中查询。点击“查找”时，回调函数S1_search_Callback会首先获取上述3个Edit控件中的文本内容：

然后使用MATLAB内建函数database连接数据库

再使用Select语句查询数据库中信息：

其中，S为当前查找的表格名（S1、S12、S2或S21），order by…asc语句是为了让查询结果按照数值从小到大排序。当前排序的优先级为听诊区＞健康状况＞时长＞编号。

查询到的已排序数据存储于search_result矩阵，使用set(handles.S1_uitable,'Data',search_result)函数将查询结果显示在S1_uitable表格控件中。S1_uitable的回调函数S1_uitable_CellSelectionCallback会在点击其中任意一个单元格的时候触发。为了在点击某一行时显示该行数据对应的心音波形，回调函数中编写的主要代码如下：

其中，eventdata.Indices为当前选中的单元索引，将其行索引赋值给hangIndex_S1，newData_S1为通过Get函数获取到的当前显示表格中的数据，heart_id_S1为当前所点击的心音编号，heart_wav_S1为读入的心音数据。Axes函数用于获取当前画图的Axes控件，然后使用Plot函数画图，波形如图5中蓝色曲线所示。

4.2 心音生成

依次选择4个部分的心音，点击“添加”按钮，然后输入各个部分心音的开始时间和放大倍数，点击“生成”按钮，即可将所选4个部分的心音合并为一个完整心动周期的心音。生成心音时，先定义全零数组composed_wav，数组长度需足以包含4段心音：composed_wav=zeros(1,(floor(fs*S21_begin_t)+length(S21_readed)));

其中，S21_begin_t为S2～next S1段的开始时刻，fs为采样率4 kHz。然后根据每段心音的开始时刻给composed_wav赋值。例如S2段的赋值程序为composed_wav((floor(fs*S2_begin_t))+k2)=composed_wav((floor(fs*S2_begin_t))+k2)+S2_readed(k2).*S2_amp;即将composed_wav中从floor(fs*S2_begin_t)开始的一段心音波形数据加上S2乘以放大倍数S2_amp后的值，生成的心音波形和时频图已在图6中给出。时频图是对时域信号进行短时傅里叶变换（short-time Fourier transform，STFT）得到的[28]，计算公式如下：

式中，X（τ，ω）为变换后得到的时频信号；x（t）为时域心音信号；h（t-τ）为窗函数，本文选用Hanning窗，帧长为256，帧移为32。点击“播放”按钮时，软件调用sound函数。在“持续周期数”后的文本框中输入待生成心音的心动周期数，点击“生成”按钮时会调用函数audiowrite(['created_',S1_id,'_',S12_id,'_',S2_id,'_',S21_id,'_',datestr(now,'yyyy-mm-dd_HH-MM-SS'),'.wav'],composed_wav_all,4000)保存心音。其中composed_wav_all为多周期合并后的心音数据，datestr()函数用于生成字符串格式的当前时刻，4 000为当前心音的采样率（单位为Hz）。

4.3 查看文件和软件信息

当生成心音后，点击图7（a）中“查看生成文件”按钮即会触发回调函数打开预设的文件夹。在此文件夹中可查看此前生成的.wav格式的心音。点击“使用说明”按钮会调用回调函数open（'instruction.fig'）弹出软件介绍，如图7（b）所示。软件介绍预先在Text控件中写好，内容为软件的使用方法。

5 教学效果分析

图7 心音文件与软件信息查看

为了研究心音合成软件在教学方面的效果，由专业医务人员在中国医学科学院阜外医院进行实验。实验前将患者及家属分为A、B 2组，每组40人。经过调查，确认所有测试者均非医学专业。采用不同的方式给上述2组人员授课，教授心音与心脏瓣膜疾病间的关系。对A组人员授课时，使用华盛顿大学医学系的临床高级物理诊断学习与教学（advanced physical diagnosis learning and teaching at the bedside，APDLTB）数据库中的心音，对B组人员授课时采用心音合成软件。授课结束后，A组每位人员可以任意选择数据库中的心音反复听10 min，B组人员自行操作心音合成软件10 min，按照个人意愿合成各类心音并试听。教学完成后，让测试者填写问卷以判断测试者对知识的掌握程度。问卷中包含10个关于心音和心脏瓣膜疾病关系的选择题，2组测试者问卷填写结果详见表2。

表2 2组测试者调查问卷填写结果统计单位：人

由表2可以看出，B组获得优秀和良好的人数明显高于A组；B组获得中等和较差的人数少于A组。可见使用心音合成软件的培训效果优于使用APDLTB数据库。这可能是由于本软件自由度较高，能够激发被测试者主动探索的兴趣，使用户在操作过程中对心音与心脏瓣膜疾病的关系有了较为深刻的理解。

6 结语

本文以MATLAB GUI为平台设计了一套心音合成软件。首先采集大量健康人和心脏疾病患者的心音，然后对心音进行降噪和分割处理，构建Access数据库并实现心音合成软件。本研究中没有采用传统的基于解析公式的心音合成的方案，而是以大量数据为基础，这样可以生成更加真实的心音。此外，本软件界面直观、页面跳转少、使用流畅，为心音听诊教学提供了便利，具有较高的实用价值和进一步开发的潜力。且与使用APDLTB数据库中的心音进行教学相比，本软件能够达到更好的教学效果。未来的研究计划包括扩大心音的采集范围以包含更多的疾病、进一步完善软件功能、推动软件的实用化等。