一种基于matlab的ECG频谱分析.docx
《一种基于matlab的ECG频谱分析.docx》由会员分享，可在线阅读，更多相关《一种基于matlab的ECG频谱分析.docx（17页珍藏版）》请在冰豆网上搜索。
一种基于matlab的ECG频谱分析
1引言
1.1ECG处理的意义
生物医学信号属于强噪声背景下的低频微弱信号，它是由复杂的生命体发出的不稳定的自然信号。
作为一种对判断人体生命状况极其重要的生理信号，处理心电信号就显得很有必要，尤其在临床诊断上的应用。
根据心电信号，可以绘制出心电图。
心电图是心脏兴奋的发生、传播及过程的客观指标，是反映心脏兴奋的电活动过程，它对心脏基本功能及其病理研究方面，具有重要的参考价值。
心电图可以分析与鉴别各种心律失常；也可以反映心肌受损的程度和发展过程；心房、心室的功能结构情况。
在指导心脏手术进行及指示必要的药物处理上有参考价值。
心电图的检查必须结合多种指标和临床资料，进行全面综合分析，才能对心脏的功能结构做出正确的判断。
在临床上，分析心电信号，可以确诊心肌梗塞及急性冠状动脉供血不足，协助诊断慢性冠脉供血不足、心肌炎、心肌病及心包炎，判定有无心房、心室肥大，从而协助某些心脏病的病因学诊断，例如风湿性、肺源性、高血压性和先天性心脏病等，观察某些药物对心肌的影响，包括治疗心血管疾病的药物（如洋地黄、抗心律失常药物）及可能对心肌有损害的药物。
此外，对某些电解质紊乱（如血钾、血钙的过高或过低），心电信号不仅有助于诊断，还可以对指导治疗有重要参考价值。
1.2ECG处理的发展与现状
心电信号处理是目前信号处理领域中的研究热点之一,其真正实现将有力地促进医疗事业的发展和人们健康水平的提高,也将是现代信号处理理论与技术在医疗领域中应用的重大突破。
心电信号处理系统的研究内容广泛,涉及的基础理论和关键技术繁多,是一个多学科交叉的庞大课题。
到目前为止,现有的心电信号的分析方法中还存在着诸多的缺陷和不足,在理论研究和实际应用方面仍有许多改进和创新的空间。
由于计算机和微电子技术的问世，20世纪50年代末、60年代初就已经采用了计算机进行心电信号处理，随后随着计算机以及信号检测、分析与处理技术的不断完善，20世纪后半叶心电图的计算机辅助分析应用日益广泛，在心电监护、心电Holter、心电图自动分析仪、心电体表标测以及高频心电、心室晚电位及希氏束、胎儿心电及心率变异性等许多方面都取得了进展。
如今，随着数字信号处理技术的发展和各种数字信号处理器性能的提高，用数字信号处理器来进行心电图辅助分析也成为了可能。
本文主要介绍心电信号的预处理去工频干扰、平滑滤波和QRS复波检测的方法及演示结果。
2Matlab处理ECG的优势
2.1Matlab简介
Matlab是由美国MathWorks公司20世纪80年代中期推出的数学软件。
MATLAB是“MatricLaboratory”的缩写，意为“矩阵实验室”，优秀的数值计算能力和卓越的数据可视化能力使其很快在数学软件中脱颖而出。
Matlab已经发展成为多学科、多种工作平台的功能强大的大型软件。
在欧美的高校和研究机构中，MATLAB是一种非常流行的计算机语言，许多重要的学术刊物上发表的论文均是用MATLAB来分析计算以及绘制出各种图形。
Matlab是一完整的并可扩展的计算机环境，是一种进行科学和工程计算的交互式程序语言。
它的基本数据单元是不需要指定维数的矩阵，它可直接用于表达数学的算式和技术概念，而普通的高级语言只能对一个个具体的数据单元进行操作。
因此，解决同样的数值计算问题，使用MATLAB要比使用Basic、Fortran和C语言等提高效率许多倍。
许多人赞誉它为万能的数学“演算纸”。
MATLAB采用开放式的环境，你可以读到它的算法，并能改变当前的函数或增添你自己编写的函数。
2.2Matlab的优势
友好的工作平台和编程环境。
MATLAB由一系列工具组成。
这些工具方便用户使用MATLAB的函数和文件，其中许多工具采用的是图形用户界面。
包括MATLAB桌面和命令窗口、历史命令窗口、编辑器和调试器、路径搜索和用于用户浏览帮助、工作空间、文件的浏览器。
随着MATLAB的商业化以及软件本身的不断升级，MATLAB的用户界面也越来越精致，更加接近Windows的标准界面，人机交互性更强，操作更简单。
而且新版本的MATLAB提供了完整的联机查询、帮助系统，极大的方便了用户的使用。
简单的编程环境提供了比较完备的调试系统，程序不必经过编译就可以直接运行，而且能够及时地报告出现的错误及进行出错原因分析。
简单易用的程序语言。
Matlab一个高级的矩阵/阵列语言，它包含控制语句、函数、数据结构、输入和输出和面向对象编程特点。
用户可以在命令窗口中将输入语句与执行命令同步，也可以先编写好一个较大的复杂的应用程序（M文件）后再一起运行。
新版本的MATLAB语言是基于最为流行的C＋＋语言基础上的，因此语法特征与C＋＋语言极为相似，而且更加简单，更加符合科技人员对数学表达式的书写格式。
使之更利于非计算机专业的科技人员使用。
而且这种语言可移植性好、可拓展性极强，这也是MATLAB能够深入到科学研究及工程计算各个领域的重要原因。
出色的图形处理功能。
Matlab自产生之日起就具有方便的数据可视化功能，以将向量和矩阵用图形表现出来，并且可以对图形进行标注和打印。
高层次的作图包括二维和三维的可视化、图象处理、动画和表达式作图。
可用于科学计算和工程绘图。
新版本的MATLAB对整个图形处理功能作了很大的改进和完善，使它不仅在一般数据可视化软件都具有的功能（例如二维曲线和三维曲面的绘制和处理等）方面更加完善，而且对于一些其他软件所没有的功能（例如图形的光照处理、色度处理以及四维数据的表现等），Matlab同样表现了出色的处理能力。
同时对一些特殊的可视化要求，例如图形对话等，MATLAB也有相应的功能函数，保证了用户不同层次的要求。
另外新版本的MATLAB还着重在图形用户界面（GUI）的制作上作了很大的改善，对这方面有特殊要求的用户也可以得到满足。
应用广泛的模块集合工具箱。
MATLAB对许多专门的领域都开发了功能强大的模块集和工具箱。
一般来说，它们都是由特定领域的专家开发的，用户可以直接使用工具箱学习、应用和评估不同的方法而不需要自己编写代码。
目前，MATLAB已经把工具箱延伸到了科学研究和工程应用的诸多领域，诸如数据采集、数据库接口、概率统计、样条拟合、优化算法、偏微分方程求解、神经网络、小波分析、信号处理、图像处理、系统辨识、控制系统设计、LMI控制、鲁棒控制、模型预测、模糊逻辑、金融分析、地图工具、非线性控制设计、实时快速原型及半物理仿真、嵌入式系统开发、定点仿真、DSP与通讯、电力系统仿真等，都在工具箱家族中有了自己的一席之地。
3ECG特征参数及分析
3.1心电信号的特点
心电信号属于强噪声背景下的低频微弱信号，它是由复杂的生命体发出的不稳定的自然信号，由于受到人体诸多因素的影响，因而有着一般信号所没有的特点。
信号弱。
例如从母体腹部取到的胎儿心电信号仅为10μv,成人的心电信号范围也仅为5mv.
噪声强。
由于人体自身信号弱，加之人体又是一个复杂的整体，因此信号易受噪声的干扰。
如胎儿心电混有很强噪声，它一方面来自肌电、工频等干扰，另一方面，在胎儿心电中不可避免地含有母亲心电，母亲心电相对我们要提取的胎儿心电则变成了噪声。
随机性强。
心电信号信号不但是随机的，而且是非平稳的。
正是因为生物医学信号的这些特点，使得心电信号处理成为当代信号处理技术最可发挥其威力的一个重要领域。
3.1.1心电信号的特征参数
图3.1
如图3.1为完整的心电信号波形图，分别由P波、P—R段、P—R间期、QRS复合波、S—T段、T波和U波组成。
P波。
心脏的兴奋发源于窦房结，最先传至心房，故心电图各波中最先出现的是代表左右两心房兴奋过程的P波。
兴奋在向两心房传播过程中，其心电去极化的综合向量先指向左下肢，然后逐渐转向左上肢。
如将各瞬间心房去极的综合向量连结起来，便形成一个代表心房去极的空间向量环，简称P环。
P环在各导联轴上的投影即得出各导联上不同的P波。
P波形小而圆钝，随各导联而稍有不同。
P波的宽度一般不超过0.11秒，电压（高度）不超过0.25毫伏。
P—R段。
是从P波终点到QRS波起点之间的曲线，通常与基线同一水平。
P－R段由电活动经房室交界传向心室所产生的电位变化极弱，在体表难于记录出。
P—R间期。
是从P波起点到QRS波群起点的时间距离，代表心房开始兴奋到心室开始兴奋所需的时间，一般成人约为0.12～0.20秒，小儿稍短。
超过0.21秒为房室传导时间延长。
QRS复波。
代表两个心室兴奋传播过程的电位变化。
由窦房结发生的兴奋波经传导系统首先到达室间隔的左侧面，以后按一定路线和方向，并由内层向外层依次传播。
随着心室各部位先后去极化形成多个瞬间综合心电向量，在额面的导联轴上的投影，便是心电图肢体导联的QRS复合波。
典型的QRS复合波包括三个相连的波动。
第一个向下的波为Q波，继Q波后一个狭高向上的波为R波，与R波相连接的又一个向下的波为S波。
由于这三个波紧密相连且总时间不超过0.10秒，故合称QRS复合波。
QRS复合波所占时间代表心室肌兴奋传播所需时间，正常人在0.06～0.10秒之间。
S—T段。
由QRS波群结束到T波开始的平线，反映心室各部均在兴奋而各部处于去极化状态，故无电位差。
正常时接近于等电位线，向下偏移不应超过0.05毫伏，向上偏移在肢体导联不超过0.1毫伏，在单极心前导程中V1，V2，V3中可达0.2～0.3毫伏；V4，V5导联中很少高于0.1毫伏。
任何正常心前导联中，ST段下降不应低于0.05毫伏。
T波。
是继QRS波群后的一个波幅较低而波宽较长的电波，反映心室兴奋后再极化过程。
心室再极化的顺序与去极化过程相反，它缓慢地从外层向内层进行，在外层已去极化部分的负电位首先恢复到静息时的正电位，使外层为正，内层为负，因此与去极化时向量的方向基本相同。
连接心室复极各瞬间向量所形成的轨迹，就是心室再极化心电向量环，简称T环。
T环的投影即为T波。
再极化过程同心肌代谢有关，因而较去极化过程缓慢，占时较长。
T波与S－T段同样具有重要的诊断意义。
U波。
在T波后0.02～0.04秒出现宽而低的波，波高多在0.05毫伏以下，波宽约0.20秒。
一般认为可能由心舒张时各部产生的负后电位形成，也有人认为是浦肯野氏纤维再极化的结果。
3.2心电信号噪声的来源及特点
心电信号在经过采集、数模转换过程中，不可避免的受到各种类型的噪声干扰，这些干扰使得得到的心电信号的信噪比较低,甚至淹没了心电信号。
通常心电信号中主要包括以下3种噪声：
①工频干扰
主要包括50HZ电源线干扰及高次谐波干扰。
由于人体分布电容的存在使入体具有天线效应以及较长的导联线暴露在外，50HZ的工频干扰在心电信号中是常见的，依情况不同，其干扰幅度达心电信号峰一峰值的0—50%。
②肌电干扰
由于病人的紧张或寒冷刺激，以及因某些疾病如甲状腺机能亢进等，都会产生高频肌电噪声，其产生是众多肌纤维分时随机收缩时引起的，频率范围很广（DC-1000V）,谱特性接近白噪声，其频率一般在5HZ—2KHZ之间。
③基线漂移
这种噪声是因呼吸、肢体活动或运动心电图测试所引起的。
稍微剧烈的肢体运动将引起心电信号波形发生改变，严重地破坏了心电信号分析的准确性。
上下波动和扭曲的心电图也令医师眼花缭乱，影响诊断，其频率一般在0.05—2HZ之间。
4ECG频谱分析的ＧＵＩ界面设计
4.1ＧＵＩ简介
图形用户界面GUI（GraphicalUserInterfaces）是由窗口光标、按键、菜单、文字说明等对象构成的一个用户界面。
用户通过一定的方法（如鼠标或键盘）选择、激活这些图形对象，使计算机产生某种动作或变化，比如实现计算、绘图等。
在Matlab中可以很方便地使用用户界面开发环境创建GUI应用程序。
设计的GUI布局代码会存储在FIG文件中，同时Matlab会产生一个M文件用于存储调用函数。
通过点击控件，Matlab会运行该控件下的Callback函数。
因此，编写GUI应用程序，可简单地认为是布局控件，编写callback函数以及解决函数间数据传递问题。
GUI可很好地进行技术，方法的演示，并且它还可以被编译器编译成应用程序。
4.2ＧＵＩ实现的功能
为了能实现人机交换的功能，使得对心电信号的处理的操作更加简洁方便，需要制作一个GUI界面及GUI程序，以实现以下的功能：
1以菜单栏的方式载入心电数据。
如点击菜单，出现open的选项，然后弹出一个对话框，可以选择所需要的格式的文件并载入。
2显示处理后的数据波形图。
设计axes控件，可以显示被处理的心电信号波形图以及处理后的心电数据的波形图。
3调整参数，以达到处理心电信号的最佳效果。
设置EditText控件，可以任意输入想要的数据，调控数据的显示范围以及滤波器的参数。
4数据存储。
设计了一个可以保存处理后的心电数据的pushbutton控件，实现一键保存数据的功能。
4.3GUI设计
4.3.1GUI布局
图4.3
如图4.3所示,显示部分有4个axes控件显示载入的原始心电信号波形、原始心电信号的频谱图、处理后的心电信号图、处理后的心电信号频谱图以及用QRS复波检测方法设计的滤波器滤出来的波形。
Input为原始信号输入并进行可视化的部分。
该部分GUI界面上设置有6个可编辑的文本框，可以在里面输入适当的参数，并设有两个按钮，分别可以显示原始心电信号的波形以及经傅里叶变换后的频谱图。
Filtering是滤波控制的部分。
该部分设有三个滤波器的选择按钮和相应的编辑文框。
另外还有综合了三个滤波器功能的滤波控件，可同时调用三个滤波器进行滤波。
QRS是QRS复波检测部分，该部分可以控制QRS复波检测。
4.3.2编写控件代码
Callback为控件被触发时调用的函数。
取值是字符串，可以是某个M文件或一小段Matlab语句。
当用户激活控件对象（如用鼠标点击控件对象图标或移动滚动条的滑块）时，应用程序就运行该属性定义的子程序。
具体程序代码见附录。
Backgroundcolor定义控件对象区域的背景色，缺省颜色是系统定义的浅灰色。
具体程序代码参见附录。
5ECG信号处理
现代医学表明，心电信号（ECG）含有临床诊断心血管疾病的大量信息，ECG的检测与分析在临床诊断中具有重要价值，是了解心脏的功能与状况、辅助诊断心血管疾病、评估各种治疗方法有效性的重要手段。
但由于实际检测工况的非理想，ECG信号中往往含有工频噪声及电极极化等引起的各种随机噪声。
噪声的存在降低了诊断的准确性，其中影响最大的是50Hz工频干扰和基线漂移噪声。
因此，在ECG信号检测过程中，如何抑制工频干扰和基线漂移是必须解决的问题。
心电信号是微弱低频人体生理电信号，通常频率在0.05—100Hz，幅值不超过4mv，通过安装在皮肤表面的电极来获取。
由于人体是一个复杂的生命系统，存在50H工频干扰及基线漂移等其他生理电信号的干扰。
噪声可能会影响到医生的临床诊断，因此，需对心电信号进行滤波，即必须做好前端数据采集的软硬件设计以保证心电数据的可靠和准确。
5.1平滑处理
要对心电信号进行平滑处理，就必须要平滑滤波。
平滑滤波是低频增强的空间域滤波技术。
它一般的目的有两类：
一类是模糊；另一类是消除噪音。
最常碰到的信号处理任务之一是平滑数据以抑制高频噪声。
高频噪声源包括60Hz电网噪声，运动伪迹和量化误差。
求几个数据点的平均值的方法是减弱高频噪声的一种简单方法。
这种滤波器被称之为移动平均滤波器。
本次对心电信号进行平滑处理的任务是保留信号中有用的频率的成分,滤掉不想要的成分,即为滤掉常说的噪声,留下想要的成分。
在本次毕业设计的制作过程中，对心电数据平滑处理使用的移动平均滤波器主要有海宁滤波器以及巴特沃斯滤波器。
5.1.1海宁滤波器
海宁移动平均滤波器是最简单的平滑滤波器之一。
正如其差分方程所示，海宁滤波器计算的是一组加权移动平均值，其中间数据点的权系数好似其余两个数据点的2倍。
y（nT）=1/4[x（nT）+2x（nT-T）+x（nT-2T）]
一旦有了差分方程（实现数字滤波器的数值算法），即可确定传递函数。
利用离散值和z域值间的关系类推，传递函数能够完全代表滤波器的性能。
x（nT）和y（nT）是当前采样时刻相联系的输入和输出序列中的点，他们分别对应于无延时z域数值X（z）和Y（z）。
同样，x（nT-T）是前一采样点的输入数值，
对应于滞后一个采样点的z域输入值，或X（z）z-1。
若将输出的Y（z）视为输入X（z）的函数，则其z变换的方程式为：
图5.1.1（a）海宁滤波器流程图
图5.1（a）所示的是用功能块直接实现方程中的各项的流程图。
z的-2次幂决定了设计中需要两个延时器。
为了与因子2和1/4相乘，必须有两个乘法器，还需要两个加法器来将各项相加。
此方程的传递函数为：
5.1.1（b）海宁滤波器零-极点分布图
海宁滤波器有两个零点，均位于z=-1；有两个极点，均位于z=0.因为它们对幅值响应的所有频率影响相同，故图中未画出。
海宁滤波器的实现是通过一段计算机的程序来完成的。
程序代码如下：
m=length（x）;
y2=zeros（1,m）;
fori=1:
m
x0=x（i）;
ifi==1
x1=0;
else
x1=x（i-1）;
end
ifi<=2
x2=0;
else
x2=x（i-2）;
end
y2（i）=1/4（x0+2*x1+x2）;
end
在for循环内，x（nT）的值从x（i）中获得，x（nT-T）的值从x（i-1）中获得，x（nT-2T）的值从x（i-2）中获得，i的值在1~m中取。
然后输入数据在延时器中移动。
下次输入之前，前一采样点的数据x（nT-T）移至输入x（nT-2T），即x（i-2）中。
当前输入x（nT）在时间点上倒退一个点，代替x（nT-T）。
在下一个循环中，重新获得x（nT）的新值，重复此过程直至m的数值全部被处理为止。
在Matlab程序坏境中运行改程序得到的效果如图5.1.1（c）和5.1.1（d）所示：
原始数据波形及频谱：
图5.1.1（c）
滤波后的波形及频谱：
5.1.1（d）
如上图所示，经海宁滤波器滤波之后，原始的心电信号波形变得平滑了许多，减少了许多毛边。
从频谱图也能发现这一变化。
5.1.2巴特沃斯滤波器
本次任务主要是滤去心电信号中的高频噪声，因此需要设计一个低通滤波器。
所谓低通滤波器，就是只能让低频分量通过。
在低通滤波器的设计中，有许多不同的逼近方法：
巴特沃斯逼近，切比雪夫逼近，椭圆逼近方法等。
其中，最简单的逼近为巴特沃斯逼近。
巴特沃斯滤波器的特点是通频带内的频率响应曲线最大限度平坦，没有起伏，而在阻频带则逐渐下降为零。
因此对于处理心电信号是比较合适的。
巴特沃斯低通滤波器的幅度平方函数
用下式表示：
式中，N称为滤波器的阶数。
当
=0时，
；
=
时，
，
是3dB截止频率。
在
=
附近，随着
加大，幅度迅速下降。
幅度特性与
和N的关系如图5.1.2（a）所示。
幅度下降的速度与阶数N有关，N越大，通带越平坦，过渡带越窄，过渡带与阻带幅度下降的速度越快，总的频响特性与理想低通滤波器的误差越小。
图5.1.2（a）巴特沃斯低通滤波器幅度特性与
和N的关系
以s代替j
将幅度平方函数
写成s的函数：
复变量s=
此式表明幅度平方函数有2N个极点，极点
用下式表示：
2N个极点等间隔分布在半径为
的圆上（该圆为巴特沃斯圆），间隔是
rad。
例如N=3，极点间隔为
rad，如图5.1.2（b）所示。
图5.1.2（b）三阶巴特沃斯滤波器极点分布图
为形成因果稳定的滤波器，2N个极点中只取s平面坐平面的N个极点构成
，而右半平面的N个极点构成
。
的表达式为：
得出
后即能得到数字化的巴特沃斯滤波器。
用Matlab信号处理工具箱函数设计巴特沃斯滤波器，其中有五种用格式，在本次设计任务中，用到的是格式之一的[N,wc]=buttord（wp,ws,Rp,As）和[B,A]=butter（N,wc,’ftype’）。
其中，[N,wc]=buttord（wp,ws,Rp,As）用于计算巴特沃斯数字滤波器的阶数N和3dB截止频率wc。
调用参数wp，ws分别为数字滤波器的通带、阻带截止频率的归一化值，要求：
0≤wp≤1，0≤ws≤1。
1表示数字频率pi。
Rp，As分别为通带最大衰减和组带最小衰减（dB）。
当ws≤wp时，为高通滤波器；当wp和ws为二元矢量时，为带通或带阻滤波器，这时wc也是二元向量。
N，wc作为butter函数的调用参数。
具体程序代码见附录。
根据计算要求的设计指标为：
通带最大衰减Rp＝0.2dB，阻带最小衰减As＝15dB，通带边界频率wp＝100Hz，采样频率ws＝500Hz。
得出的效果如图5.1.2（c）和5.1.2（d）所示：
原始数据波形及频谱：
图5.1.2（c）
滤波后的波形及频谱：
图5.1.2（d）
如上图所示，经巴特沃斯滤波器滤波之后，原始的心电信号波形变得平滑而干净，滤去了许多由于高噪声干扰而产生的描边。
而且，从频谱图也能发现这一变化。
5.2去工频干扰
升级会员 