基于MATLAB的心音信号处理.doc

基于MATLAB的心音信号处理

目录

1.概述 2
1.1 心音信号及MATLAB简介 2
1.2 设计目的 5
1.3 设计任务 5
2.设计流程 5
2.1 设计方案选择 5
2.1.1预处理方法的选择 5
2.1.2时域分析 7
2.1.3频域分析 7
2.2 MATLAB仿真结果及分析 8
2.2.1预处理的结果及分析 9
2.2.2时域分析的结果及分析 12
2.2.3频域分析的结果及分析 13
3.总结 16
4.设计心得与体会 17

1.概述
1.1 心音信号及MATLAB简介
心音信号是指在心动周期中，心肌收缩、心脏瓣膜启闭，心室壁、大动脉瓣等被血流冲撞，引起机械振动发出的声音。产生的声音信号通过周围组织传导到胸壁，可以通过听诊器来听诊，也可以利用传感器将心音信号转换为电信号，进行心音信号的显示和存储。心音作为人体的一种重要的心脏、心血管生理信号，体现着和心脏等生理情况相关的信息，能够反映心脏等的生理和病理信息。心音对于有关心脏疾病和心血管疾病的诊断具有重要的诊断价值，是评估一个人心脏功能情况的重要依据[6]。
心音的频率一般在5-600HZ左右，杂音频率可达1500HZ左右，但是人的听觉系统能够接受的声音频率为16HZ-20KHZ,而且仅对频率在1000~5000HZ频率的声音最为敏感。心脏听诊的不足可以用心音图弥补。心音图可以直观的显示正常和异常心音的各段持续的时间，进一步提高心脏疾病诊断的准确性。

图1.正常的心音信号的时域图
正常心音按出现时间先后可以分为四个部分。第一心音（S1)、第二心音（S2)是最常见的，此外，还有第三心音（S3)、第四心音（S4)。临床上一般听到的是第一心音和第二心音，一些儿童、青少年以听到第三心音，有些老年人也可以听到第四心音。第一心音（S1)产生的原因是房室瓣关闭、血流急速冲击房室瓣，主动脉壁和肺动脉壁被心室喷射出的血液撞击，引起振动。一般在心尖搏动处最强，持续时间一般为0.1-0.12s。在音调方面，比较低沉。第二心音（S2)在T波的末端出现。主要产生原因是主动脉瓣关闭和肺动脉瓣关闭。持续时间一般为0.08s，相对较短[4]。第三心音（S3)主要是血流流向心室时，血流速度过急，在心室内又突然减速引起室壁的振动而产生。一般在第二心音后的0.1-0.2s后产生。在儿童、青少年的心脏上可以听到生理性的第三心音。第四心音是由于心房收缩后，血流流进心室的速度很快，心室壁振动而引起。心脏听诊是诊断心脏和心血管疾病不可或缺的方法，该方法简便、有效。
MATLAB是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。Matlab是一个高级的矩阵/阵列语言，它包含控制语句、函数、数据结构、输入和输出和面向对象编程特点[5]。用户可以在命令窗口中将输入语句与执行命令同步，也可以先编写好一个较大的复杂的应用程序（M文件）后再一起运行。新版本的MATLAB语言是基于最为流行的C++语言基础上的，因此语法特征与C++语言极为相似，而且更加简单，更加符合科技人员对数学表达式的书写格式。使之更利于非计算机专业的科技人员使用。而且这种语言可移植性好、可拓展性极强，这也是MATLAB能够深入到科学研究及工程计算各个领域的重要原因。MATLAB由一系列工具组成。这些工具方便用户使用MATLAB的函数和文件，其中许多工具采用的是图形用户界面。包括MATLAB桌面和命令窗口、历史命令窗口、编辑器和调试器、路径搜索和用于用户浏览帮助、工作空间、文件的浏览器。随着MATLAB的商业化以及软件本身的不断升级，MATLAB的用户界面也越来越精致，更加接近Windows的标准界面，人机交互性更强，操作更简单。而且新版本的MATLAB提供了完整的联机查询、帮助系统，极大的方便了用户的使用。
简单的编程环境提供了比较完备的调试系统，程序不必经过u编译就可以直接运行，而且能够及时地报告出现的错误及进行出错原因分析。在MATLAB环境进行信号处理可以方便的调用MATLAB提供的函数，运用简单的语句就可以实现极为复杂的运算，加快了信号处理的进程，MATLAB还提供GUI的功能，便于用户设计友好的交付界面。
1.2 设计目的
1.掌握医学信号处理的整个过程
2.学会使用MATLAB进行信号处理
3.学会分析信号及信号处理后得到的结果
4.完成对心音信号的处理及分析
1.3 设计任务
根据心音信号的噪声特性及时频特性，选用合适信号处理方法，运用MATLAB完成对心音信号的预处理获得干净的心音信号，并进行时域分析、频域分析，获得其相应的病理信息，并利用MATLAB提供GUI，设计显示界面。
2.设计流程
2.1 设计方案选择
2.1.1预处理方法的选择
心音信号是微弱低频生理信号，容易混入多种多样的噪声，许多因素（呼吸状态、心脏的血液流动动力学状态以及环境等）都会影响心音信号。由于心音的频率一般在5—600HZ左右。使用巴特沃斯滤波器进行滤波，其具有通频带内的频率响应曲线最大限度平坦，没有起伏，而在阻频带则逐渐下降为零的特点。所以本设计采用巴特沃斯滤波器进行带通滤波所以设计的带通滤波器的通带频率为5-600hz ，通带最大衰减3 dB和阻带最小衰减18dB。
利用巴特沃斯滤波器滤波的信号，还有一定噪声在其中，本设计设计了小波阈值滤波的方法再次进行去噪。
小波变换是用数学变换分析信号的一种方法[5]。在分析低频长时信号时，时间分辨率很低，频域分辨率很高。小波变换的局部化特性是比较好的，这种局部特性反映在频率域和空间域两方面。由于心音的频率一般在5—600HZ左右，杂音频率可达1200HZ左右。根据相关文献db3 小波对心音信号处理的效果较好，所以对心音信号釆用db3小波进行5层分解，然后采用小波软阈值去噪方式进行去噪，最后对各层小波分解系数进行重构。软阈值处理是把小波系数大于阈值的变为该点与阈值的差值。
小波软阈值去噪[4]的过程如下：
（1)选择小波和小波分解的层数j,计算含噪声 信号的小波分解系数;
(2)对每层系数选择一个阈值,并且对高频系 数用阈值处理;
(3)根据第j层的低频系数和从第一层到第j 层的高频系数,计算信号的小波重构。软阈值处理的数学表示为
djk-λ,djk≥λ
djk'= 0, |djk|<λ
djk+λ,djk≤-λ

(4)式中,λ为阈值,djk为小波系数,djk为处理后的小波系数。
2.1.2时域分析
本部分主要包括对心音信号提取其包络（运用MATLAB提供的函数实现）、心率的计算、第一、二心音的间隔时间。
正常的心率是45-100次/分，S1与S2的间隔是0.05-0.5s。
心率的计算主要通过在获取心音信号的峰值信息以后，根据两个S1波之间的时间间隔即为心动周期的基本原理，通过计算得到。S1与S2的间隔时间，是在求出多个心音信号间隔以后，在求平均值得到的结果。
2.1.3频域分析
快速傅氏变换，是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶虚、实等特性，对离散傅立叶变换的算法进行改进获得的。FFT算法可分为按时间抽取算法和按频率抽取算法[5]。
小波分解在心音的身份识别、心音的诊断检测、心音的特征提取方面应用方面尤为突出。不同的心脏疾病的心音信号在各分解频带上具有不同的能量分布，可将小波包系数的能量作为心音信号的特征向量。本处理采用处理心音信号常用的db3小波对心音信号进行分层取，获得相关频带的信号。

图2. 小波分解示意图
2.2 MATLAB仿真结果及分析
本设计采用的信号是来自师兄给的一个心音音频信号，该心音信号的采样频率为11025HZ,采样位数是8位，采用点数为 262144。其使用MATLAB显示后得到下图：
图3.心音信号波形
有图3可看出信号的噪声较大，分析其频域信息，可得到其频域信息，发现其信号（包括噪声在内）主要集中在低频成份。

图4.信号的频谱特性
2.2.1预处理的结果及分析
1.带通滤波
采用巴特沃斯滤波器以后的得到信号如下图所示：
图5.滤波过后的信号
经过带通滤波过后的信号不仅没有变得干净，反而混入了许多的高频成分进去。分析带通滤波器的响应函数，得到下图：
图6.滤波器响应函数
分析获得信号的频谱特性，经过FFT变换得到下图，有图可看出经过滤波后的信号混入了许多毛刺，产生更多的噪声。具体原因难以得知，估计跟信号的采样频率有关，鉴于接下来还须用到小波的阈值去噪，带通滤波的作用不大，取消带通滤波器的使用。
图7.滤波得到的信号的频谱特性
2.50hz陷波器
为滤除工频干扰设计了，50hz的陷波器，其响应的函数图像如下图：
图8.50hz陷波器响应函数
获得的信号经过FFT变换得到下图的频谱特性，可以看出50hz的工频干扰已被滤除，最终得到的信号如图10：
图9.FFT变换过后的频谱特性

图10.经过陷波器过后心音信号
3.小波软阈值滤波
采用db3小波将信号进行8层分解，由于信号的采样频率为11025hz，将1、2、3层信号舍弃，将4到8层的信号的高频成分进行软阈值滤波后进行重构，各层频率范围为
分解层数
近似分量
细节分量
第1层
0~2756
2756~5512
第2层
0~1378
1378~2756
第3层
0~689
689~1378
第4层
0~348
348~689
第5层
0~172
172~348
第6层
0~86
86~172
第7层
0~43
43~86
第8层
0~21.5
21.5~43
表1 db3小波8层分解后各层分量对应的频段范围 单位：HZ

得到的经过小波软阈值滤波以后的信号如下图，经过小波去噪以后，可以看出信号的基线变窄了，信号更为清晰。

图11.小波软阈值去噪过后的心音信号
2.2.2时域分析的结果及分析
利用MATLAB提供的提取包络的函数，对处理过的信号提取包络的到下图，为更为清晰的看到心音信号，只取了信号的10000到100000个点的信号。并计算出的心跳次数为87次/分， S1与S2间隔为0.28s，属于正常范围。

图12.信号的包络图
2.2.3频域分析的结果及分析
1. 频域分析
利用MATLAB提供的FFT函数，对信号进行频域分析，得到信号频谱图如下图所示。可以看出信号的频率主要集中在0到100hz：
图13.心音信号频谱特性
2. 小波分解
小波分解在心音的身份识别、 心音的诊断检测、心音的特征提取面应用尤其突出。不同心脏疾病的心音信号在各分解频带上具有不同的能量分布， 可将小波包系数的能量作为心音信号的特征向量。将希
图14.心音信号小波分解的其中三层细节信息
l
尔伯特变换和小波分析结合起来， 可提高提取信号包络信息的精度[2]。对心音分裂信号进行分解， 然后再特定曾进行重构， 可以检测出分裂不明显的心音分裂信号。用离散小波变换对心音信号进行分解然后从样本中提取特征值再将选取的特征值输入到前馈型神网络进行训练和识别，可以对心音进行诊断。本设计采用常用于心音信号分析的db3小波进行信号的分解，分解层数为8层。上图所示为其中三层的信息。
3. 短时傅里叶变换
对于心音这种非平稳信号 , 由于不同时刻有不同的频率成分 , 单纯在时域或频域都不足以表示信号特征。为获得心脏的动力学特征 , 需要研究心音信号频率成分随时间的变化。本文采用的短时傅氏变换 ( STFT) 所做出的时频分析 , 是通过对加窗后的信号进行傅氏变换实现的。
图15.短时傅里叶变换得到的三维时频图像
其中横轴表示频率，纵轴表示窗层数。
2.2.4 部分MATLAB函数
1.50hz陷波器的设计
f0=50;
Ts=0.001;
NLen=512;
n=0:NLen-1;
apha=-2*cos(2*pi*f0*Ts);
beta=0.96;
b=[1 apha 1];
a=[1 apha*beta beta^2];
y=dlsim(b,a,xy);
layer=8;
2.fft变换
L=length(x);
NFFT=2^nextpow2(L);
Y=fft(x,NFFT)/L;
f=fz/2*linspace(0,1,NFFT/2+1);
axes(handles.axes1);
plot(f,abs(Y(1:(NFFT/2)+1)));

3.sfft变换
N=length(xy);
Nw=3000;
L=Nw/2;
Ts=round((N-Nw)/L)+1;
nfft=128;
TF=zeros(Ts,nfft);
for i=1:Ts
xw=xy((i-1)*L+1:i*L+L);
temp=fft(xw,nfft)/L;
TF(i,:)=temp;
end
axes(handles.axes6);
mesh(abs(TF));

4.小波分解及软阈值滤波
layer=8;
[C,L]=wavedec(y,layer,'db3');
cD1=detcoef(C,L,1);
cD2=detcoef(C,L,2);
cD3=detcoef(C,L,3);
cD4=detcoef(C,L,4);
cD5=detcoef(C,L,5);
cD6=detcoef(C,L,6);
cD7=detcoef(C,L,7);
cD8=detcoef(C,L,8);
axes(handles.axes2);
plot(cD4,'b');
title('小波分解的第四层的高频部分');
axes(handles.axes3);
plot(cD5,'b');
title('小波分解的第五层的高频部分');
axes(handles.axes4);
plot(cD6,'b');
title('小波分解的第六层的高频部分');
Thr4=thselect(cD4,'minimaxi');
Thr5=thselect(cD5,'minimaxi');
thr6=thselect(cD6, 'minimaxi');
thr7=thselect(cD7,'minimaxi');
thr8=thselect(cD8,'minimaxi');
THR=[thr4,thr5,thr6,thr7,thr8];
[s3,CXC,LXC,PERF0,PERFL2]=wdencmp('lvd',y,'db3',5,THR,'s');%wedencmp
3.总结
本设计完成了对心音信号的提取显示，并对信号进行小波软阈值去噪、工频陷波，然后对信号进行了时域分析，包括包络的提取、心率及第一二心音间隔的计算，最后对信号进行了频域分析，包括小波分解、FFT变换、SFFT变换。
本设计还有许多不足之处，还需在后续的学习中继续完善。对于去噪，还不能灵活的使用小波进行去噪，未能达到比较好的处理效果；对于时域分析，提取的特征值太少；对于频域分析，由于未能理解所示图像中包含是生理信息，难以判断信号处理的效果是否良好；平台方面选择了相对比较低级的MATLAB。面对以上不足，还需继续学习信号处理的相关知识，尤其是医学信号的处理。
4.设计心得与体会
通过本次设计，我主要的收获如下：
1.明白并理解了信号处理的整个过程，知道信号处理的基本步骤及基本方法。
2.提高了阅读文献的能力，也提高了自学的能力。
3.进一步熟练掌握了对MATLAB进行信号处理的能力，学会使用了MATLAB的GUI设计功能。
由于知识储备不足，在进行设计过程中，我会不断面对新的问题，需要不断通过各种途径去加以解决，当问题得以解决，便会有小小的成就感。在进行程序的编写调试的过程中，有时一个程序会反复的出错，调试上好几遍都过不了，我这时就会变得十分的烦躁，改来改去反而错误会更多，让同学帮忙看才发现，只是一些很小的错误，只要稍加细心就会发现。可惜只因为粗心大意，浪费了大量的时间，而且由于不够有耐心，造成的时间更加多。在以后的学习中这中粗心大意、不耐烦的习惯一定要努力去改变。
通过本次课程设计不仅仅是完成一项作业，在这个过程我学会的也不仅仅是相关专业知识，更多的是知识以外的东西。最后，对任课老师侯老师的孜孜不倦的教授表示衷心的感谢。

参考文献

[1] 于云之，聂邦.心音的临床意义及研究现状[J].现代医学仪器与 应用，1997,9（3）：9-12
[2] 陶冶薇，张会香.基于LabVIEW的心音多功能分析仪[J].南京邮电大学学报，2012,32（4）：75-79
[3] 张孝桂，何为.基于嵌入式系统的便携式心音分析仪的研究[J].仪器仪表学报，2007,28（2）：302-306
[4] 倪林.小波变换与图像处理[M].合肥，中国科学技术大学出版社2010.6
[5] 薛年喜.MATLAB在数字信号处理中的应用[M].北京清华大学出版社，2008,01
[6] 黄政钦.心音心电采集系统设计与信号处理[D].云南大学，2014.05

16