基于Hilbert_Huang变换的第一心音信号时频分析_cropped

基于Hilbert_Huang变换的第一心音信号时频分析_cropped Ξ 基于 Hilbert - Huang 变换的第一心音信号时频分析 杜寿昌 ,梁虹 ()云南大学 信息学院 ,云南 昆明 650091 ( ) 摘 要 根据瓣膜原理 ,第一心音 S1是在心脏收缩期由二尖瓣和三尖瓣关闭时引起的振动产生的 ,含有多个频率分量 . 对第一心音信号的分析研究 ,在临床上对心脏疾病的诊断有重要意义 . 本文用一种全新的时频分析方法 : Hilbert - Huang 变换 ( ) HHT,对 30 例心音数据进行心音分析实验 . 实验结果表明 : HHT 方法可以有效的分析心音信号 ; S1 含有二尖瓣 M1 及三尖瓣 T1 两个主要成份 ;异常 S1 的 M1 和 T1 的频率比正常 S1 有升高 . 关键词 经验模态分解 ; Hilbert - Huang 变换 ;模态函数 ;第一心音 ;时频分析 () 【文章编号】1005 —7188 200402 - 0095 - 04 【中图分类号】TN911173 【文献标识码】A 床医学研究所提供 . 采样频率为 5 KHz ,与心音同步 0 引言 记录的还有心电图. 所选的 15 例是二尖瓣或主动脉 瓣经修复后心脏功能正常的心音 ,另 15 例是异常心 正常人的心音随心脏搏动呈周期性变化 ,临床 ( ) ( ) 上分 为 第 一 心 音 S1 , 第 二 心 音 S2 , 第 三 心 音 音 . ( ) ( ) S3,第四心音 S4和心脏杂音 . 多数情况下人耳只 2 原理 能听到 S1 和 S2 . 心音听诊是一种古老的诊断心脏 2 . 1 EMD 方法 疾病的方法 ,通过听取心音 ,听诊者可以获得用以判 6 ,7 断心脏疾病的许多信息 . 然而 ,有资料表明人耳对中 Huang提 出 一 种 新 的 信 号 处 理 方 法 ———经 ( ) 高频段和较大强度的心音可较好的分辨 ,而对那些 验模态分解方法 EMD. EMD 方法的目的是把信号 () 有重要诊断意义的低频心音及强度小的心音常不能 分解为多个模态函数 intrinsic mode function , IMF分1 分辨. 人耳听力敏感局限性和听诊者的临床经验 量 . 接着 ,对分解得到的分量进行希尔伯特变换 ,从 在很大程度上影响了听诊的准确度 ,且不能给出定 而得出时频平面上的能量分布谱图 ———希尔伯特谱 () ( ) 量分析的结果 . 而心音图 phonocardiogram PCG除了 Hilbert Spectra. 3 Huang 提出了满足下面两个条件的数据称为模 ( ) 记录听力敏感区 10Hz 以上以 内 的 信 息 外 , 还 记 态函数 :一是在整个数据序列中 ,极值点的数量与过 3 2 录了大量 10Hz 以下的信息. 对心音图的各部分 零点的数量必须相等 ,或最多相差为 1 ; 二是在任何 进行进一步的定量分析将有助于临床上对心脏疾病 时间点上 ,由极大值和极小值定义的包络的平均值 的诊断提供科学的依据 . 为零. 很 显 然 , 大 多 数 信 号 不 满 足 上 述 条 件 , 不 是 心音的组成成分复杂 ,且心音信号是非平稳信IMF 分量 ,因此 ,必须把信号分解为模态分量 . Huang 3提出了把信号分解 为 模 态 分 量 的 算 法 ———EMD 算 ,各种不同数字信号处理方法被用于对心音进号 3 ,5( ) 法 ,也被称为筛选 Sifting过程. ,时频分析是其中的重要方法和手段 . 本行分析 主要步骤 :文采用一种全新的时频分析方法 - 经验模态分解方 () 1找出原序列中的所有局部极值点 ,其中所有 ( ) 法 EMD Empirical Mode Decomposition及 Hilbert 谱分 ( ) 的局部最大值被一个三次样条 spline连接成为上包 ( )分析第一心音 ,析 ,统称 Hilbert - Huang 变换 HHT () ( ) 络 Upper envelope,包络序列值 xt. 同理 ,局部最 取得较好效果 .max () 小值产生下包络 Lower envelope,该包络序列值 xmin 1 心音数据 ( ) t,上下包络应将所有的数据都包含在它们之间. 本论文分析的心音数据是由加拿大蒙特利尔临 Ξ 收稿日期 :2003 - 02 - 18 () 基金项目 :云南省自然科学基金资助项目 2001C0004M. () 作者简介 :杜寿昌 1970,,男 ,云南永胜人 ,讲师 ,研究生 ,主要从事电子技术教学及实际应用. () ( )2. 2 Hilbert 变换和 Hilbert 谱 2求取每一时刻的瞬时平均值 m t ( ) ( ) ( ) m t= xt+ xt / 2 ( )通过 EMD , 就 可 以 对 每 一 个 IMF 分 量 进 行1 max min () ( ) 3求原始信号与 m t的差值 Hilbert 变换. + ? ( ) ( ) ( )x t- m t= ht ( )1 2 (τ) x 1 ( )( ) ^x t = H[ x t] = τπ?t - ( )- ?9 ( ) 理想情况下 , ht 应是一个 IMF 分量. 但对于 1 τ ( ) ( )x t 3 h t d= 非线性数据 ,包络均值可能不同于真实的局部均值 , 1 一些非对称波仍 可 能 存 在. 在 第 二 次 处 理 中 , 分 量 ( )h t = πt ( ) ht被当作待处理数据 ,重复以上三步 ,于是 :1 ( ) ( ) 以 实 信 号 x t 作 为 实 部 , 实 信 号 x t 的 ( ) ( ) ( )()ht- mt= ht3 1 11 11 ( ) Hilbert 变 换 x^ t 作 为 虚 部 , 构 成 了 一 个 复 信 号( ) ht 是一个可以把处理过程重复 k 次 , 直到 1 k ( ) z t, 称为解析信号或称信号预包络 , 记为 : 模态分量 ,于是 :θ( ) j t( ) ( ) α( ) ( )() 10 = x t + jx^ t = t ez t ( ( )()) ()- mt = ht 4 5h ( ) 1 k - 11 k 1 k 2 2α( ) ( ) ( ) t= x t+ ^x t ( )11 ( ) ( ) ( ) 定义 : ht= ct, 则 ct 就是从原始数据 1 k 1 1 ( )x^ t θ( ) ()12 t= arctan ( ) 中处理得到的第一个模态分量 IMF.( )x t 为了保证模态分量保存足够的反映物理实际的 其瞬时频率按下式计算 : 幅度与频率调制 ,必须确定一个筛选过程停止的准 θ( )1 dt ( )() 13 f t = ) ( π 则 . 可用前后两个 h t 的标准差 sd 的大小作为模 2dt () 对式 8中每一个 IMF 分量进行希尔伯特变换态函数的判据 : T 2 ( ) ( ) 之后 , 可以把数据表示成下面的形式 :t - h t ] [ h k - 1 k ( )5 s d = 2 n n ? t = 0 θ( )( )j( ) h tt j w tdt i i k - 1 ? ( ) α ( ) z t = Re t eα= Re e i i?? i = 1 i = 1 实践表明 ,当 sd 的值介于 0 . 2,0 . 3 时 ,既能保 ( )14 证模态函数的线性和稳定性 ,又能使所得的模态函 式中 , 省略了残余函数 r ,Re 表示取实部 . 数具有相应的物理意义 . ( ) ( ) αθ即得到幅值函数 t 和相位函数 t , 然后 按( ) ( ) 然后用原序列 x t 减 去ct 得 到 剩 余 分 量 1 ( ) 式 13求得瞬时频率 . ( ) rt: 1 () 式 14称为 Hilbert 幅值谱 ,简称 Hilbert 谱. 记作 :( ) ( ) ( )()x t - ct = rt 6 1 1 n ( ) j w tdti ? ( ) α( ) 15 H w , t= Re e 至此 ,提取第一个模态函数过程全部完成 .i? i = 1 ( ) 由于剩余部分 rt仍然包含较长周期分量的 1 进一步可以定义边际谱为 : + ?( ) 信息 , 所以 rt仍被当作新的原数据按以上相同 1 ( ) ( )( ) h w= H w , tdt 16 的处理过程来处理. 该处理过程可对所有的接下来? - ? ( )的剩 余 分 量 rt进 行 处 理 , 依 次 提 取 第 2 , 第 3 , j 以上的 EMD 和相应的 Hilbert 谱分析方法统称 ( ), 直到第 n 个模态函数 ct n 为 Hilbert - Huang 变换. ( ) ( ) ( )( )( ) ( )rt- ct= rtrt- ct= rt ( ) 1 2 r n - 1n n 3 方法 ()7 之后 , 当分量 c或剩余分量 r变成小到比预定 EMD 方法及 HHT 用 Matlab6 . 1 实现 ,在 Microsoftn n 值小时 ; 或当剩余分量 r变成单调函数 , 而从中不 Windows 98 环境中运行 . n ()能再筛选出基本模式分量时 , 筛选停止 . 将等式 6 () 1提取第一心音 : ( ) 与等式 7相加 , 最终得到原序列 : 为了减少计算量 ,消除 1250 Hz 以上的杂音 , 从 n PCG 中 提 取 S1 之 前 , 首 先 将 原 采 样 频 率 降 为( ) ( ) ( )( )x t= ct+ rt8 i n ? i = 1 2500 Hz. 用于临床诊断的心音通常在 1000 Hz 以下 , 到此已经把原始数据分解成 n 个模态分量 , 及() 仍低于新采样频率 2500 Hz的一半. 一个剩余分量 r. n 然后确定 S1 的起止点 . 确定起止点方法在有关 8 文献中已有报道 ,有的以奇异点来定位,有的以心 9 电图或颈动脉脉冲为参考. 本文利用心音图与心 电图的关系 ,从心音图中选择一个较好的第一心音 并将 它 作 为 模 板 , 找 到 心 电 R 波 波 峰 , 截 取 其 后 10 120ms 的范围定位出 S1 . 然后 , 在 PCG 中的其它 各心动周期中挑选与模板心音互相关系数在 0 . 990 以上的八个第一心音求均值 . 求均值后的第一心音 9 即为我们进行分析的第一心音信号. 图 1 为正常 第一心音时域波形 ;图 2 为异常第一心音时域波形 . () 2将提取出来的第一心音进行 EMD 分解得各 IMF 波形如图 3 及图 4 所示 : 从 IMF 波形可看出 , 第一个模态函数 C1 是从 心音中分解出的振幅最大 ,频率最高 ,波长短 ,依次 下去 ,振幅逐渐变小 ,频率逐渐变低 ,波长越来越长 , 各种模态函数的这种分布情况 ,是由其本性决定的 , 它总是把最主要信号先提取出来 ,即 EMD 分解的头 要的频率成份比正常心音的频率升高 .几个内在模态函数集中了原序列中最显著的信息 . () 3Hilbert 变换和 Hilbert 谱分析4 实验结果及分析 图 5 为正常第一心音边际谱波形 ; 图 6 为异常 第一心音边际谱波形 . 边际谱表达了每个频率在全 表 1 所示为实验时读取的 30 例第一心音的数据. () 局上的幅度 或能量贡献 ,代表在统计意义上的全 心音成份分析 : 部数据累加幅度 ,通过边际谱可明显看出第一心音 图 5 ,图 6 及表 1 中显示出 S1 主要由两部分组 的两个主要特征频率 . 根据瓣膜理论 ,它们分别对应 成 ,第一部分是 S1 的二尖瓣成分 M1 . 这部分为二尖 M1 和 T1 两个成份 . 瓣的关闭及左室壁和肌腱索的突然牵张所产生的 , 表 1 第一心音特征频率实验数据峰值可认为是二尖瓣关闭时刻. 较为陡峭的波峰 ,频 ) () (( ) ( ) ( ) ( ) ID 正常F1 HzF2 HzID 异常F1 HzF2 Hz 率低 ,是 S1 的三尖瓣成分 T1 ,该部分为三尖瓣关闭 10006 68 33 10004 100 65 及半月瓣开放 ,血流减速所致 ,峰值可认为是三尖瓣 10010 71 24 10009 110 56 关闭时刻 . 10011 60 32 10013 87 37 异常第一心音 M1 和 T1 均升高的原因是由于 10012 80 32 10015 112 61 生物瓣膜变异 、钙化 、硬化等造成瓣膜硬度增加 ,使 10016 50 34 10027 100 68 瓣膜关闭时产生振动的频率升高. 二尖瓣狭窄导致 10022 67 32 10034 84 47 10033 70 41 10041 80 36 瓣叶瓣强有力的关闭 ,使第一心音响亮清脆 .10064 68 39 10057 110 78 5 结论及展望 10066 73 24 10062 112 50 10069 77 40 10093 95 50 结论 : EMD 方法及 HHT 可以有效的分析 PCG 这 10072 79 42 10097 110 60 样的非平稳信号 ,通过对第一心音 S1 的分析研究 ,可 10074 57 34 20007 90 67 看出 S1 包含二尖瓣 M1 和三尖瓣 T1 两个主要的成 10084 75 32 20016 82 50 分 ;异常 S1 的 M1 和 T1 的频率比正常 S1 有升高. 10095 68 40 20046 100 58 由于心音信号的随机性 、时变性 、复杂性 ,心音 10096 61 40 20064 121 53 研究应建立一个公认的分析指标体系及完善的心音 Mean 68. 3 34. 6 Mean 99. 5 55. 7 数据库 ,还应结合生理 、病理特性进行综合研究 ,同 Standard 8 . 4 5. 7 Standard 12. 8 11 . 4 时需进一步研究心音发声机理.ID 为病例号 . F1 为第一特征频率 . F2 为第二特征频率 . EMD 方法是一种分析非线性 、非平稳信号的较 Mean 为均值 . Standard 为标准偏差 . 按照同样的方法 ,对 30 例心音数据进行心音分 好的方法 ,已在很多领域得到成功应用. 当然 ,它的研 析实验 ,发现利用 EMD 方法及 HHT 可以有效的分 究还刚刚起步 ,许多问题还有待研究. 如信号分解及 析心音信号 ,通过对 M1 和 T1 频率的检测 , 分别读 插值方法的影响 ,模式混淆问题 ,趋势或均值问题等. ( ) 取两个主要谱峰的峰值 结果见表 1,可看出正常 第一心音 M1 频率主要集中在 68 . 3 ?8 . 4 Hz , T1 主 6 致谢 要集中在 34 . 6 ?5 . 7 Hz , 异常心音 M1 的频率主要集 中在 99 . 5 ?12 . 8 Hz , T1 主要集中在 55 . 7 ?11 . 4 Hz. 衷心 感 谢 加 拿 大 蒙 特 利 尔 临 床 医 学 研 究 所 统计意义上 ,异常第一心音 S1 含有 M1 和 T1 两个主 Durand L G. 为我们提供心音数据 .参考文献 : 1 罗健仲 ,罗琳 . 心脏听诊 M . 第 5 版 . 人民卫生出版社 ,2000. 27 - 35. 2 DURAND L G , PIBAROT P. Digital Signal Processing of the Phonocardiogram : Review of the Most Recent AdvancementsJ . Critical () Reviews in Biomedical Engineering , 1995 ,233/ 4:163 - 167. () 王衍文 ,王海滨 ,程敬之. 心音信号的识别与分类J . 生物医学工程学杂志 , Vol . 16 3, 1999 :382 - 387.3 () 吴延军 ,徐泾平 ,赵艳 ,等 . 心音信号的时频分析 J . 西安交通大学学报 ,1997 ,31 8:25 - 28.4 () 李桥 ,邵庆余 . 心音信号分析方法述评 J . 山东生物医学工程 , 1994 ,13 2:52 - 56.5 ()下转第 115 页 参考文献 : () SHAPIRO L G. Matching three - dimensional objects using a relational paradigmJ . Pattern Recognition , 1984 ,17 4:385 - 405. 1 SHAPIRO L G , HARAL ICK R M. Structural description and inexact matchingJ . IEEE Trans. Patt . Anal . Mach. Intel . ,1981 ,PAMI - 2 () 1 2: ,173 - 184. FENG Jian2sheng. Structured three - dimensional object recognitionD . England , The University of London ,1987. 3 FENG Jian2sheng. A part - distance measure for inexact matching A . Proc . AMSE Int . conf . on Signals &Systems C . Dalian , 4 1989. 139 - 145. FENG Jian2sheng , He Peng2ling. The partially consistent la - beling problemA . Proc . CIE & IEEE Int . Conf . On Signal Processing 5 C . Beijing ,1990. 827 - 830 Improving of strulture description and matching ZHANG Xue2qin , L I J uan ( )The Second Northwest Institute for Nationalities , Yinchuan 750021 ,China Abstract :Based on the description of relationship data structure and the equdl attitude of a part of elements , the au2 thors put forward a kind of top - down method improving relationship data structure of the object and matching result . Key words : data structure ; top - down ; matching () 责任编辑 王 琳()上接第 98 页 6 HUANG N E. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non - stationary time series analysisC . Proc . Roy. Soc . London , 1998 , Ser . A 454 ,903 - 995 . HUANG N E , SHEN Z , LONG S R. A new view of nonlinear water waves : The Hilbert spectrumR . Annu. 7 Rev. Fluid Mech. , 1999 . 417 - 457 . () 吴效明 ,岑人经 ,吴云鹏 ,等. 心音信号检测的一种新方法J . 生物医学工程学杂志 ,1996 ,13 4:349 - 352.8 WANG W L , GUO Z , YANG J ,et al . Analysis of the First Heart Sound Using the Matching Pursuit Method J . 9 () Med. & Bio . Eng. & Comput ,2001 ,39 6:644 . 10 刘毅 ,邵庆余 ,李桥 ,等 . 老年人第三第四心音信号小波变换时频分析研究 J . 北京生物医学工程 , () 2000 ,19 1:42 - 44 . Time - frequency analysis of the first heart sound based on the hilbert - huang transform DU Shou2chang , L IANG Hong ( )College of Information , Yunnan University , Kunming 650091 ,China Abstract :According to the valvular theory , the vibrations affected by the mitral and tricuspid valves closure in early ( ) systole produce the first heart sound S1. S1 usually includes many frequency components. Analysis of the first heart sound will be helpful to diagnose the clinic heart diseases. In this paper , a new method using the Hilbert - Huang Trans2 ( ) form HHTwas used to decompose S1 into a series of time - frequency atoms. By examining 30 cases of heart sound da2 ta , It has being proved the HHT is an efficiency method to study the heart sound , S1 consists of two major high - frequency components , M1 and T1 ,the abnormal bioprosthetic mitral valve generates higher frequency components than the normal one. Key words : EMD ; HHT ; IMF ; the first heart sound ; time - frequency analysis ()高 飞责任编辑