1、基于MATLAB的心电信号的分析与处理设计河南科技大学课 程 设 计 说 明 书课程名称 医学信号处理 题 目 基于MATLAB的心电信号的分析与处理设计（2） 院 系 医学技术与工程学院班 级 医疗器械工程111班学生姓名 指导教师 侯海燕 宋卫东日 期 2014年9月11号课程设计任务书（指导教师填写）课程设计名称医学信号处理 学生姓名 专业班级医疗器械工程111班设计题目 基于MATLAB的心电信号的分析与处理设计（2）一、课程设计目的1熟练掌握使用MATLAB程序设计方法2掌握数字信号处理的基本概念、理论、方法3掌握序列离散傅里叶变换的MATLAB实现，并进行频谱分析4熟练掌握使用MA
2、TLAB设计IIR或FIR数字滤波器5学会用MATLAB对信号进行分析和处理二、设计内容、技术条件和要求一）设计内容与技术条件1根据给定的一段MITBIH心电信号(101号），画出心电信号的时域波形和频谱图（幅频和相频）；2根据心电信号频率范围及其噪声的频率范围设计2个滤波器（一个IIR，一个FIR）实现对心电信号滤波。滤波器的种类（高通，低通，带通,带阻）,滤波器性能指标（通阻带截止频率，衰减系数）,滤波器的设计方法（IIR有冲击响应不变法和双线性变换法，FIR有窗函数法及频率抽样法)等自行设计.要求输出所设计的滤波器的系统函数，画出滤波器的频率响应(幅频响应和相频响应)曲线；3用该滤波器对
3、心电信号进行滤波，画出滤波以后心电信号的时域波形和频谱（幅频）；分析信号滤波前后心电信号的时域和频域的变化；4两个滤波器滤波效果异同分析;5运用GUI设计一个心电信号处理系统界面。（选作）二）设计要求1根据滤波器的性能指标要求，设计数字滤波器;2程序中按照IIR滤波器的步骤一步步完成设计； 尽可能的少调用MATLAB自带的函数文件；3设计程序要通用性好,整齐易懂，并要求主要语句有注释；4设计结果中的图示要美观，整齐,有标题，有纵横坐标标示；5课程设计报告要有理论依据、设计过程，结果分析。报告要求实事求是、文理通顺、字迹端正.三、时间进度安排9月1日：查阅资料，复习医学信号处理课程的相关内容，熟
4、悉MATLAB的工作环境，掌握运算编程方法。9月2日5日：按照课程设计要求，完成课程设计的程序设计部分。9月8日10日：调试并完善程序，撰写课程设计说明书.9月1112日：答辩。四、主要参考文献参考资料（1）数字信号处理-理论、算法与实现胡广书清华大学出版社(2）数字信号处理及MATLAB实现余成波清华大学出版社（3)MATLAB7。0在数字信号处理中的应用罗军辉机械工业出版社（4）数字信号处理原理及其MATLAB实现从玉良电子工业出版社（5）MATLAB帮助文件相关网站（1)MATLAB中国论坛（2)MATLAB学习网（3）MATLAB中文论坛http：/指导教师签字:侯海燕 2014年8月
5、25日1、设计目的意义12、设计内容12。1设计原理12.2要求12.3内容13、计过程及分析13。1原始心电信号分析13。2 设计滤波器3 4、结果分析75、参考文献 7附录81、设计目的意义本次课程主要达到以下几个目的： (1） 熟练掌握使用MATLAB程序设计方法(2) 掌握序列离散傅里叶变换的MATLAB实现，并进行频谱分析(3) 熟练掌握使用MATLAB设计IIR或FIR数字滤波器 2、设计内容 2.1设计原理(1)频谱分析（傅里叶变换)（2)滤波器设计IIR（冲激响应不变法，双线性变换法） 详解见附录FIR（窗函数法，频率抽样法） 详解见附录（3)信号滤波（时域卷积，频域相乘）2。
6、2要求要求设计出心电数据处理的处理与分析程序。(1) 处理对象：心电数据;（2） 内容:心电数据仿真，心电数据处理（仿真数据，真实数据）；（3） 结果：得到处理结果.2.3内容（1）原始心电信号的时域波形和频谱图(幅频和相频).（2）IIR滤波器的频率响应(幅频响应和相频响应）曲线；原始心电信号经IIR滤波器滤波后的心电信号的时域波形和频谱图(幅频和相频）；(3)FIR滤波器的频率响应（幅频响应和相频响应）曲线；原始心电信号经FIR滤波器滤波后的心电信号的时域波形和频谱图（幅频和相频）；3、设计过程及分析3.1原始心电信号分析用load函数将原心电信号导入b = load（C:UsersAdm
7、inistratorDesktop课程设计101ecg.txt)，并画出心电信号的时域波形和频谱图(幅频和相频）,如图1所示:图1 原始心电信号的时域波形图及频谱图心电信号由于受到人体诸多因素的影响，因而有着一般信号所没有的特点:（1)信号弱。心电信号是体表的电生理信号，一般比较微弱，幅度在10pV5mV，频率为005100Hz。例如从母体腹部收取到的胎儿心电信号仅10zV50IV.（2）噪声强。由于人体自身信号弱,加之人体又是一个复杂的系统，因此信号容易受到噪声干扰.(3）随机性强。心电信号不仅是随机的，而且是非平稳的。同时，在心电图检测过程中极易受到各种噪声源的干扰，从而使图像质量变差,使
8、均匀和连续变化的心电数值产生突变，在心电图上形成一些毛刺.使原本很微弱的信号很难和噪声进行分解。可能出现的噪声有如下的种类：1)工频干扰工频干扰是由电力系统和人体的分布电容引起的，其频率包括50Hz（MIT-BIH数据库数据工频因为是美国标准，所以是60Hz）的基波及其各次谐波，其幅值成分在ECG峰一峰值的0-50范围内变化.2）引起基线漂移的干扰心电信号有时候会出现信号基线起伏不平的现象,造成这样的现象有很多原因，主要的有：呼吸运动人体呼吸时胸腔内器官和组织会发生一定程度的变化，会对在体表记录到的心电图波形的幅度和形态有所影响，表现为基线随呼吸产生周期性或非周期性漂移，从而导致心电波形的幅度
9、随呼气和吸气而分别上抬和下移。呼吸运动是引起心电基线漂移的主要原因。运动伪迹运动伪迹是由于人体轻微运动造成电极与入体的接触电阻发生变化而引入的一种干扰,它的产生原因仅仅是接触电阻的变化，而不是接触的断续.这种干扰同样导致信号基线的变化，但不是基线的跃变。信号记录和处理中电子设备引起的干扰这种干扰对信号影响很大，严重时可完全淹没心电信号或使得基线剧烈漂移，其中导联开路和放大器的热移是主要因素。这种干扰往往无法通过心电分析算法来校正。由于心电波形已经完全畸变，此时对这些数据分析已无太大意义。所以一般跳过此段数据。3)高频噪声心电信号中的高频噪声主要是肌电噪声。肌肉收缩会产生mV级的肌电干扰，表现为
10、心电图上不规则的细小波纹，使心电图模糊不清或产生失真。肌电噪声的特点是频率范围较广，频谱分布非常复杂。3。2 设计滤波器（1）IIR滤波器的设计1）IIR滤波器的设计过程：按照技术要求设计一个模拟滤波器,得到模拟低通滤波器的传输函数H（s），再按一定的转换关系将H（s）转换成数字低通滤波器的系数函数H（z)。这样设计的关键问题就是找到这样的转换关系，将s平面上的H(s）转换成z平面上的H(z）。2）巴特沃斯滤波器分母多项式的因式表示，如表1所示：表1 巴特沃斯滤波器分母多项式的因式表示3）巴特沃斯低通滤波器的阶数公式 N=log10（10（As/10）1)/（10(Rp/10)-1）)/（2*
11、log10（ws/wp）)）4)巴特沃斯低通滤波器函数由巴特沃斯低通滤波器的阶数公式和巴特沃斯滤波器分母多项式的因式表示求出归一化巴特沃斯低通滤波器Has (s ）N=7则Has(s )=1/（s+1)*（s2+0。4450s+1)*（s2+1。247s+1)(s2+1。8022s+1）5）巴特沃斯低通滤波器的频域特性，如图2所示:图2 巴特沃斯低通滤波器的相频和幅频特性6）经过巴特沃斯低通滤波器器后心电信号的时域波形和频谱图,如图3所示:图3 经过巴特沃斯低通滤波器器后心电信号的时域波形图和频谱图对比原始信号的时域波形图和频谱图可得通过低通滤波器后的心电信号波形图可以明显看出波形变得平滑,由
12、工频干扰产生的毛刺被低通滤波器成功滤除。（2)FIR滤波器的设计1）FIR滤波器的设计过程：1给定理想的频率响应函数 Hd(ejw）及技术指标, w;2求出理想的单位抽样响应hd（n);3根据阻带衰减选择窗函数w(n）;4根据过渡带宽度确定N 值N=A/w；5求所设计的FIR滤波器的单位脉冲响应h（n）=hd(n）*w（n)；6计算频率响应Hd(ejw)，验算指标是否满足要求。2）布拉克曼窗低通滤波器的频域特性，如图4所示： 图4 布拉克曼窗低通滤波器的相频和幅频特性3)经过布拉克曼窗低通滤波器器后心电信号的时域波形和频谱图,如图5所示图5 经过布拉克曼窗低通滤波器器后心电信号的时域波形和频谱
13、图对比原始信号的时域波形图和频谱图可得通过低通滤波器后的心电信号波形图可以明显看出波形变得平滑，由工频干扰产生的毛刺被低通滤波器成功滤除.4、结果分析通过用matlab设计的IIR低通滤波器和FIR低通滤波器对带有噪声的心电信号进行处理和分析,消除了原心电信号中由工频干扰产生的毛刺现象,虽然仍存在一些误差，但基本上达到了预期效果。通过做课程设计熟悉并掌握了相关的matlab操作环境;进一步加深了对滤波器设计相关知识的理解，熟练了滤波器设计的方法和过程。5、参考文献：1 陈天华. 数字图像处理M。 北京. 清华大学出版社，20092 刘卫国. MALTAB程序设计与应用M. 北京。 高等教育出版
14、社，20083 程正兴。 小波分析算法与应用. 西安。 西安交通大学出版社,19984 程佩青. 数字信号处理. 北京。 清华大学出版社，20075 苏金明，王永利。 MALTAB应用指南M. 上册. 北京电子工业出版社，20046 夏良正. 数字图像处理(修订版）M. 南京. 东南大学出版社,19997 霍红涛，林小竹，何薇. 数字图像处理M. 北京。 北京理工大学出版社，20038 张开滋，郭继鸿，刘海洋。 临床心电信息学M. 长沙. 湖南科技出版社，20029 黄宝晨，朱怡然. 心电图基本知识J. 中国乡村医药杂志,2004(第7页）10 许原. 心电图解读心电图如何解读和诊断J. 中国
15、临床医生，2004(第5页) 附录：巴特沃斯数字低通滤波器：巴特沃斯滤波器是电子滤波器的一种。巴特沃斯滤波器的特点是通频带的频率响应曲线最平滑.巴特沃斯滤波器的特点是通频带内的频率响应曲线最大限度平坦，没有起伏,而在阻频带则逐渐下降为零。 在振幅的对数对角频率的波得图上,从某一边界角频率开始,振幅随着角频率的增加而逐步减少，趋向负无穷大. 巴特沃斯滤波器的振幅对角频率单调下降，并且也是唯一的无论阶数，振幅对角频率曲线都保持同样的形状的滤波器.只不过滤波器阶数越高,在阻频带振幅衰减速度越快。其他滤波器高阶的振幅对角频率图和低级数的振幅对角频率有不同的形状。窗函数法：窗函数法是设计FIR滤波器的最
16、主要方法之一，实际中遇到的离散时间信号总是有限长的，因此不可避免的要遇到数据截短的问题，在信号处理中，对离散序列的截短是通过序列与窗函数相乘来实现的.在信号处理中，窗函数是一种除在给定区间之外取值均为0的实函数。譬如:在给定区间内为常数而在区间外为0的窗函数被形象地称为矩形窗。任何函数与窗函数之积仍为窗函数，所以相乘的结果就像透过窗口“看”其他函数一样。窗函数在光谱分析、滤波器设计以及音频数据压缩等方面有广泛的应用。程序流程图：clear all；close all；%时域波形图以及频谱图b=load(C：UsersAdministratorDesktop课程设计101ecg。txt);t=b
17、(：，1）;时间c=b（：，2);%幅值figure（1);subplot（311）；plot(t，c）;title（原始心电信号的时域波形图）；xlabel（时间t/s）;ylabel(幅值/A）;n=3600;m=(fft(c，n）；fs=360;%采样频率f=fs/n*（0:n-1)；定位坐标subplot(312);plot（f,abs（m）；title（原始心电信号的频谱图）；xlabel（频率/HZ）;ylabel（幅值/db)；axis(0，360，0，150);subplot(313）;plot(f,angle（m)）；title（原始心电信号的相频图）;xlabel（频率/H
18、z）;ylabel(相角/rad)；axis（0，360,-5，5);% 低通IIR滤波器;wp=2*pi10/fs; 通带数字频率转换成模拟频率ws=2*pi15/fs; %通带截至数字频率转换成模拟频率rp=2； %通带最大衰减rs=20; 阻带最小衰减N,wc=buttord（wp,ws,rp，rs,s)； %确定最小阶数N和频率参数WcBz，Az=butter(N，wc）;得巴特沃斯归一化低通原型H，w=freqz（Bz，Az); 生成频率响应参数f1=w/pi*fs/2； %采样频率转换成模拟采样频率y1=filter（Bz，Az，c）；使用filter函数对信号进行滤波figure
19、（2）；subplot（211）; plot(f1，angle（H)）; xlabel(频率/Hz);ylabel（幅度）；title(低通滤波器相频特性）； subplot(212）;plot(f1,abs(H）；xlabel(频率/Hz）；ylabel(幅度);title（低通滤波器幅频特性)； 低通IIR滤波后图形figure(3）；subplot(311)；plot(t，y1)；title（滤波后时域波形)；xlabel(时间t/s)；ylabel（幅值/A）；subplot(312);plot(f,abs(fft(y1）；title（滤波后心电信号的频谱图)；xlabel(频率/Hz
20、);ylabel(幅值/db）;axis(0，360,0,150)；subplot(313）;plot（f，angle（fft(y1）)；axis（0，360，-5，5);title(滤波后心电信号的相频特性);xlabel（频率/Hz）;ylabel（相角/rad）；低通FIR滤波器N=50; 定义窗函数的长度wc=0.3；window=blackman(N)；根据N的值产生一个布拉克曼窗windowhn=fir1(N-1,wc，window）;%可以指定窗函数向量 window。如果缺省 window参数，则 fir1默认为 hamming窗。y2=filter（hn，1,b(：,2)）；
21、使用filter函数对信号进行滤波figure（4）;freqz(hn，1）； 低通滤波器滤波后图形figure(5)；subplot（311）；plot(t，y2)；title（滤波后时域波形)；xlabel(时间t/s);ylabel(幅值/A);subplot(312）;plot(f,abs（fft（y2);title(滤波后心电信号的频谱图)；xlabel(频率/Hz）;ylabel(幅值/db)；axis（0，360，0,150）；subplot（313）；plot（f，angle(fft（y2）);title(滤波后心电信号的相频特性）；xlabel（频率/Hz）;ylabel(相角/rad)；axis（0，360，-5,5）；
