基于某MATLAB的心电信号分析报告
计算机信息处理课程设计说明书题目：基于MATLAB的心电信号分析
学院（系）：
年级专业：
学号：
学生姓名：
指导教师：
燕山大学课程设计（论文）任务书
院（系）：基层教学单位：
说明：此表一式四份，学生、指导教师、基层教学单位、系部各一份。
2014年12月 01日
摘要
心电信号是人们认识最早、研究最早的人体生理电信号之一。目前心电检测已经成为重要的医疗检测手段，但是心电信号的相关试验及研究依然是医学工作者和生物医学工程人员的重要议题。
信号处理的基本概念和分析方法已应用于许多不同领域和学科中，尤其是数字计算机的出现和大规模集成技术的高度发展，有力地推动了数字信号处理技术的发展和应用。心电信号是人类最早研究并应用于医学临床的生物电信号之一，它比其他生物电信号便易于检测，并具有较直观的规律性，对某些疾病尤其是心血管疾病的诊断具有重要意义。它属于随机信号的一种，用数字信号处理的方法和Matlab软件对其进行分析后，可以得到许多有用的信息，对于诊断疾病有非常重要的参考价值。
本课题基于matlab对心电信号做了简单的初步分析。直接采用Matlab 语言编程的静态仿真方式、对输入的原始心电信号，进行线性插值处理，并通过matlab语言编程设计对其进行时域和频域的波形频谱分析，根据具体设计要求完成程序编写、调试及功能测试，得出一定的结论。
关键词： matlab 心电信号线性插值频谱分析
目录
一：课题的目的及意义 (1)
二：设计内容与步骤 (1)
1．心电信号的读取 (1)
2．对原始心电信号做线形插值 (3)
3．设计滤波器 (5)
4．对心电信号做频谱分析 (6)
三：总结 (7)
四：附录 (8)
五：参考文献 (12)
一：课题目的及意义
心电信号是人类最早研究并应用于医学临床的生物信号之一，它比其它生物电信号更易于检测，并且具有较直观的规律性，因而心电图分析技术促进了医学的发展。
然而，心电图自动诊断还未广泛应用于临床，从国内外的心电图机检测分析来看，自动分析精度还达不到可以替代医生的水平，仅可以为临床医生提供辅助信息。其主要原因是心电波形的识别不准，并且心电图诊断标准不统一。因此，探索新的方法以提高波形识别的准确率，寻找适合计算机实现又具诊断价值的诊断标准，是改进心电图自动诊断效果，扩大其应用范围的根本途径。如何把心电信号的特征更加精确的提取出来进行自动分析，判断出其异常的类型成了亟待解决的焦点问题。本课题通过matlab语言编程，对原始心电信号进行一定的分析处理。
二：设计内容与步骤
1.心电信号读取
美国麻省理工学院提供的MIT-BIH数据库是一个权威性的国际心电图检测标准库，近年来应用广泛，为我国的医学工程界所重视。MIT-BIH数据库共有48个病例，每个病例数据长30min，总计约有116000多个心拍，包含有正常心拍和各种异常心拍，内容丰富完整。
为了读取简单方便，采用其txt格式的数据文件作为我们的源心电信号数据。利用matlab提供的文件textread或textscan函数，读取txt数据文件中的信号，并且还原实际波形。
Txt格式的心电数据文件如下：
其中文件的第一列为采样时间，第二列是在以MLII这种导联方式所得到的采样数据，第三列式以V5这种导联方式所得到的采样数据，全文件记录了约为10s的心电数据，3600个采样数据，每一行数据之间用Tab符分隔。由于数据文件中后两列数据是对同一种心电信号进行不同的导联方式所得到的采样数据，所以可以只采用其中的一种采样数据，摒弃另外一种，即可完成对此心电信号的分析。全部的心电文件记录时间约为10s，共计12个左右周期的心电信号。
实际设计心电信号数据文件时应注意：
（1）数据文件的前两行为解释说明文字，不是真正的信号数据，读取信号程序要能够自动忽略前两行文字，只读取真正的数字信号数据（严禁自己手动删除源心电数据文件中的前两行数据，必须通过程序来实现忽略前两行文字的目的）。
（2）利用数组函数分别将文件的前两个列分别读入一个一维数组。
（3）最后利用已经转为数值的分别代表心电信号时间和幅值的两个一维
数组，图形化还原原始心电信号波形，在此推荐利用利用plot(x,y)函数对数据做图形化显示。
2.对原始心电信号做线性插值
由于原始心电信号数据并不是由等间隔采样而得到的，也就是说原始的
心电数据并不是均匀的，但是后面肯定要用到数字滤波器对心电信号做处理，而MATLAB 中用数字滤波器被处理的数据应该都是等间隔的，否则处理后会
出现偏差，因此首先对原始心电信号要做的就是做一个线性插值处理，使其
变为等间隔的数据。根据心电信号的特点, 把时间分隔成0.001s 。添加的幅值点采用一次线性插值。对二维数据进行插值,相连幅值间数据的插值根据时间
进行。
其中：
1--=?i i t t t 001.0/t N ?=1--=?i i A A A 001.01+=-j j t t N
A A A j j /1?+=-1,,,,1,3,2,1,,3,2,11111===-==----j i A A t t N j arraysize i i j i j
时数组 依次排列，即得到了插值后等间隔的新数据。
根据上文中提到的插值公式，以此为原理，设计matlab 程序，对心电信
号数据做线性插值处理。插值完以后的数据应该是时间均匀的、以0.001秒为
间隔的。此步骤的实现主要是基于matlab 中的数组操作函数来实现。 其中一种插值方法的思路是：将第一步中读取的心电信号数据的时间数据和幅值数据分别存放在一个一维数组中。然后利用for 循环结构把所有数据依次读取进来。判断时间数据数组中前后两个相邻的数据间隔是否为0.001s ，如果是则判断下一对相邻两个数据；如果间隔大于0.001s 则进行一维插值做处
理。
注意：对时间数据做插值的同时一定不要忘记对幅值数据同样做插值处
理，时间数据和幅值数据一定是相互对应的。
j j A t ,
3.设计滤波器
一般正常人的心电信号频率在0.7～100HZ范围内，幅度为10μv(胎儿)～5mv(成人)。人体心电信号微弱，信噪比小，因此，在采集心电信号时，易受到仪器、人体活动等因素的影响，而且所采集的心电信号常伴有干扰。采集心电数据时，由于人的说话呼吸，常常会混有约为0.1Hz到0.25Hz频段的干扰，对于这些低频干扰，可以让信号通过一个高频滤波器，低截止频率设置为0.25,来滤波低频信号，对于高频信号干扰，可以让信号再通过一个低频滤波器，其中截止频率设置为99Hz。
4.对心电信号做频谱分析
利用MATLAB软件对处理前后的心电信号编程显示其频谱，分析比较滤波前后的频谱，得出结论。
如果分析频谱，滤波效果不明显，则需变动滤波器参数指标，重新设计滤波器。通过频谱分析，多次试验确定最合适的滤波器。
三：总结
通过一周的计算机信号处理的课程设计，我先在图书馆里查找了相关的书籍，如MATLAB类的编程书籍，各类数据处理类的书籍，即丰富了自己的知识范围，又对与自己所学的知识有了更深的了解和认识，同时也对它的应用有了一个大体的认识。这样将会更加激励我好好学习相关的知识，不断的将所学的知识用于实践。于实践中牢牢的掌握它。
在设计的过程中，我也认识到了自己所学知识的不足。这也让我再次认识到知识是无尽的，只有不断的充实自己、完善自己的知识理论体系，才能够更好的胜任自己以后的工作。设计过程中知识的不足也让我更加坚定了终身学习的决心。
在设计的过程中，我也得到了很多同学和老师的帮助。这也加强了我与
其他同学合作的能力。查找资料的过程中我也增强自己学习的能力，这些都将让我在以后的学习、生活和工作中受益匪浅。
总之，对于这样的课程设计活动，我收获了很多东西，也将使我在以后的学习、工作中更加轻松和积极。这也正是参加这次活动的目的和意义。
四：附录（程序）
function [t,Xn]=duquexinhao1(w)
fid=fopen(w);
C=textscan(fid,'%8c %f %*f','headerlines',2);%去除前两行
fclose(fid);
a=C{2};
b=C{1};
k=length(b);
for i=1:k
c(i)=strread(b(i,:),'%*s %f','delimiter',':');
end
c=c';
d=[c,a];
t=d(:,1); %时间
Xn=d(:,2); %幅度
%线性插值
function [t3,Xn3]=xianxingchazhi(t,Xn)
m=max(t);
t3=0:0.001:m;
t3=t3';
Xn3=interp1(t,Xn,t3);
%保存插值前的信号
function baocun1(t,Xn)
fid = fopen('t.txt','wt');
fprintf(fid,'%g\n',t);
fclose(fid);
fid = fopen('Xn.txt','wt');
fprintf(fid,'%g\n',Xn);
fclose(fid);
%保存插值后的信号
function baocun2(t1,Xn1)
fid = fopen('t1.txt','wt');
fprintf(fid,'%g\n',t1);
fclose(fid);
fid = fopen('Xn1.txt','wt');
fprintf(fid,'%g\n',Xn1);
fclose(fid);
%画初始信号和即插值后信号频谱
function keshehuatu(t,Xn,t1,Xn1)
f=1000;
T=1/f;
m=1:length(Xn);
k1=length(Xn1);
m1=1:k1;
q=f*m/length(Xn);
q1=f*m1/k1;
subplot(2,2,1)
plot(t,Xn)
title('初始信号时域波形')
subplot(2,2,2)
Y=fft(Xn);
plot(q,abs(Y))
title('初始信号频谱')
subplot(2,2,3)
axis([0,1000,0,1000])
plot(t1,Xn1)
title('插值信号时域波形')
Y1=fft(Xn1);
subplot(2,2,4)
axis([0,1000,0,5000])
plot(q1,abs(Y1))
title('插值信号频谱')
%低通滤波器
function [H,f]=kesheditonglvboqi(wp,ws,Rp,As,Xn1) T=0.001;
f=1/T;
[N,Wc]=buttord(wp,ws,Rp,As,'s');
[b,a]=butter(N,Wc,'s');
f=(0:length(Xn1)-1)*f/length(Xn1);w=f*2*pi;
H=freqs(b,a,w);
%高通滤波器
function [H,f]=keshegaotonglvboqi(wp,ws,Rp,As,Xn1)
T=0.001;
fs=1/T;
[N,Wc]=buttord(wp,ws,Rp,As,'s');
[b,a]=butter(N,Wc,'high','s');
f=(0:length(Xn1)-1)*fs/length(Xn1);w=f*2*pi;
H=freqs(b,a,w);
%带阻滤波器
function [H,f]=keshedaizulvboqi(wp,ws,p,s,Xn1)
T=0.001;
f=1/T;
[N,Wc]=buttord(wp,ws,p,s,'s');
[b,a]=butter(N,Wc,'stop','s');
f=(0:length(Xn1)-1)*f/length(Xn1);w=f*2*pi;
H=freqs(b,a,w);
主函数如下
将信号通过低通、高通、带阻滤波器程序
[t,Xn]=duquexinhao1('117.txt');
baocun1(t,Xn) %保存读取信号
[t1,Xn1]=xianxingchazhi(t,Xn);
baocun2(t1,Xn1)%保存插值后信号
xy=[t1,Xn1]; %仿真输入二维数组
figure(1)
keshehuatu(t,Xn,t1,Xn1) %画原始信号和插值后信号波形和频谱wp=90*2*pi; %低通滤波器滤波
ws=99*2*pi;
p=1;
s=35;
[H1,f]=kesheditonglvboqi(wp,ws,p,s,Xn1);
wp=4*2*pi; %高通滤波器滤波
ws=0.25*2*pi;
p=1;
s=35;
[H2,f]=keshegaotonglvboqi(wp,ws,p,s,Xn1);
wp=[58,62]*2*pi; %带阻滤波器
ws=[59.9,60.1]*2*pi;
[H3,f]=keshedaizulvboqi(wp,ws,p,s,Xn1);
H=abs(H1).*abs(H2).*abs(H3); %低通和高通和带阻组合的滤波器Y=H'.*abs(fft(Xn1)); %经过滤波后心电信号频谱
y=ifft(Y); %滤波后心电信号时域波形
figure(2)
subplot(2,2,1)
plot(f,abs(H1))
axis([0,150,0,1.5])
title('低通滤波器')
subplot(2,2,2)
plot(f,abs(H2))
axis([0,50,0,1.5])
title('高通滤波器')
subplot(2,2,3)
plot(f,abs(H3))
axis([0,150,0,1.5])
title('带阻滤波器')
subplot(2,2,4)
plot(f,abs(H))
axis([0,100,0,1.5])
title('组合后滤波器')
figure(3)
plot(f,abs(Y))
axis([0,100,0,80])
title('滤波后心电信号频谱')
figure(4)
subplot(2,1,1)
plot(t1,Xn1)
title('滤波前信号')
subplot(2,1,2)
plot(t1,y)
title('滤波后信号')
五：参考文献
[1] 丁玉美.数字信号处理（第二版）.西安电子科技大学出版社，2001
[2] 吴大正. 信号与线性系统分析（第四版）. 高等教育出版社，2005,8
[3] 谢嘉奎. 电子线路--线性部分(第四版). 高等教育出版社，2003,2
[4] 陈后金. 信号分析与处理实验. 高等教育出版社，2006，8