基于MATLAB的心电信号分析心电信号分析(自己做的)带程序 带图片

1、基于MATLAB的心电信号分析心电信号分析(自己做的)带程序 带图片 河北工业大学信息工程学院 信号与线性系统课程设计 基于MATLAB的心电信号分析 摘要： 本课题设计了一个简单的心电信号分析系统。直接采用Matlab语言编程的静态仿真方式、采用Simulink进行动态建模和仿真的方式，对输入的原始心电信号，进行线性插值处理，并通过matlab语言编程设计对其进行时域和频域的波形频谱分析，根据具体设计要求完成系统的程序编写、调试及功能测试。得出一定的结论。 关键字：matlab、心电信号提取、线性插值、滤波、simulink仿真。 一、课题目的及意义 心电信号是人类最早研究并应用于医学临床的

2、生物信号之一，它比其它生物电信号更易于检测，并且具有较直观的规律性，因而心电图分析技术促进了医学的发展。 然而，心电图自动诊断还未广泛应用于临床，从国内外的心电图机检测分析来看，自动分析精度还达不到可以替代医生的水平，仅可以为临床医生提供辅助信息。其主要原因是心电波形的识别不准，并且心电图诊断标准不统一。因此，探索新的方法以提高波形识别的准确率，寻找适合计算机实现又具诊断价值的诊断标准，是改进心电图自动诊断效果，扩大其应用范围的根本途径。如何把心电信号的特征更加精确的提取出来进行自动分析，判断出其异常的类型成了亟待解决的焦点问题。本课题通过matlab语言编程，对原始心电信号进行一定的分析处理

3、。 二、课题任务及要求 1、必做部分 （1）利用Matlab对MIT-BIH数据库提供的数字心电信号进行读取，并还原（2）对原始心电信号做线性插值 （3）对处理前后的心电信号分别做频谱分析 利用Matlab软件对处理前后的心电信号编程显示其频谱，分析比较滤 （4） Simulink仿真 根据前面的设计，进行基于Simulink的动态仿真设计。实现心电信号 实际波形。 波前后的频谱，得出结论。 的分析和处理。 2、选作部分 1 河北工业大学信息工程学院 信号与线性系统课程设计 （1）只截取大约2.5s，三个周期左右，大约800个采样数据进行分析。 （2）60Hz工频陷波器设计 三、设计技术指标

4、四、设计方案论证 1、必做部分 2、选作部分 2 河北工业大学信息工程学院 信号与线性系统课程设计 五、设计内容及结果分析 1、基于matlab编写的程序如下： %读取心电信号并转化成数组形式 function t,Xn=duquexinhao1(w) fid=fopen(w); C=textscan(fid,%8c %f %*f,headerlines,2);%去除前两行 fclose(fid); a=C2; b=C1; k=length(b); for i=1:k c(i)=strread(b(i,:),%*s %f,delimiter,:); end c=c; d=c,a; t=d(:,

5、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,%gn,t); fclose(fid); fid = fopen(Xn.txt,wt); fprintf(fid,%gn,Xn); fclose(fid); %保存插值后的信号 function baocun2(t1,Xn1)

6、 fid = fopen(t1.txt,wt); fprintf(fid,%gn,t1); fclose(fid); fid = fopen(Xn1.txt,wt); fprintf(fid,%gn,Xn1); fclose(fid); %画初始信号和即插值后信号频谱 3 河北工业大学信息工程学院 信号与线性系统课程设计 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

7、(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=butto

8、rd(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); %带阻滤波器 functi

9、on H,f=keshedaizulvboqi(wp,ws,p,s,Xn1) T=0.001; 4 河北工业大学信息工程学院 信号与线性系统课程设计 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); 主函数如下 （1）、将信号通过低通、高通、带阻滤波器程序 t,Xn=duquexinhao1(117.txt); baocun1(t,Xn) %保存读取信号 t1,Xn1=xianxingchazhi(t,Xn)

10、; 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

11、*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,a

12、bs(H3) axis(0,150,0,1.5) title(带阻滤波器) 5 河北工业大学信息工程学院 信号与线性系统课程设计 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(滤波后信号) 所出图形如下 6 河北工业大学信息工程学院 信号与线

13、性系统课程设计 7 河北工业大学信息工程学院 信号与线性系统课程设计 结果分析： （2）、直接通过带通滤波器程序 t,Xn=duquexinhao1(117.txt); baocun1(t,Xn) %保存读取信号 t1,Xn1=xianxingchazhi(t,Xn); baocun2(t1,Xn1)%保存插值后信号 figure(1) keshehuatu(t,Xn,t1,Xn1) %画原始信号和插值后信号波形和频谱 wp=2,80*2*pi; ws=0.25,99*2*pi; p=1; s=35; H1,f=kesheditonglvboqi(wp,ws,p,s,Xn1); H=abs(H

14、1) ; %带通 Y=H.*abs(fft(Xn1);%经过滤波后心电信号频谱 y=ifft(Y); %滤波后心电信号时域波形 figure(2) subplot(1,2,1) plot(f,abs(H1) axis(0,200,0,1.5) title(带通滤波器) 8 河北工业大学信息工程学院 信号与线性系统课程设计 subplot(1,2,2) plot(f,abs(Y) axis(0,100,0,80) title(滤波后心电信号频谱) figure(3) subplot(2,2,1) plot(t1,Xn1) title(滤波前信号) subplot(2,2,2) plot(t1,y

15、) title(滤波后信号) subplot(2,2,3) plot(t1,Xn1) axis(0,1.5,-1.5,1.5) title(滤波前截取一部分信号) subplot(2,2,4) plot(t1,y) axis(0,1.5,-1.5,1.5) title(滤波后截取一部分信号) 所出图形如下 9 河北工业大学信息工程学院 信号与线性系统课程设计 结果分析： 10 河北工业大学信息工程学院 信号与线性系统课程设计 （3）、将信号通过低通、高通组合成的带通滤波器程序 t,Xn=duquexinhao1(117.txt); baocun1(t,Xn) %保存读取信号 t1,Xn1=xi

16、anxingchazhi(t,Xn); baocun2(t1,Xn1)%保存插值后信号 figure(1) keshehuatu(t,Xn,t1,Xn1) %画原始信号和插值后信号波形和频谱 xy=t1,Xn1; wp=0.52*2*pi; %低通滤波器滤波 ws=0.62*2*pi; p=1; s=35; H1,f=kesheditonglvboqi(wp,ws,p,s,Xn1); wp=0.10*2*pi; %高通滤波器滤波 ws=0.25*2*pi; p=1; s=35; H2,f=keshegaotonglvboqi(wp,ws,p,s,Xn1); H=abs(H1).*abs(H2)

17、; %低通和高通组合的带通 Y=H.*abs(fft(Xn1); %经过滤波后心电信号频谱 y=ifft(Y); %滤波后心电信号时域波形 figure(2) subplot(2,2,1) plot(f,abs(H1) axis(0,1,0,1.5) title(低通滤波器) subplot(2,2,2) plot(f,abs(H2) axis(0,1,0,1.5) title(高通滤波器) subplot(2,2,3) plot(f,abs(H) axis(0,1,0,1.5) title(组合 带通滤波器) subplot(2,2,4) plot(f,abs(Y) axis(0,1,0,2

18、60) title(滤波后心电信号频谱) figure(3) subplot(2,1,1) plot(t1,Xn1) title(滤波前信号) subplot(2,1,2) 11 河北工业大学信息工程学院 信号与线性系统课程设计 plot(t1,y) title(滤波后信号) 所出图形如下 12 河北工业大学信息工程学院 信号与线性系统课程设计 结果分析： 三种方案比较分析： （4）系统零极点分析（在此只以高通滤波器为例） %求高通滤波器的阶数及分子分母系数 wp=0.1*pi;ws=0.25*pi;Rp=1;As=35;T=1;%数字指标 OmegaP=(2/T)*tan(wp/2);%通带

19、模拟频率 OmegaS=(2/T)*tan(ws/2);%阻带模拟频率 cs,ds=afd_butt(OmegaP,OmegaS,Rp,As);%归一化巴特沃斯滤波器原型系统函数 N=ceil(log10(10(Rp/10)-1)/(10(As/10)-1)/(2*log10(wp/ws) OmegaC=wp/(10(Rp/10)-1)(1/(2*N); %求对应于N的3db截止频率； b,a=u_buttap(N,OmegaC);%去归一化巴特沃斯滤波器原型系统函数 db,mag,pha,w=freqz_m(b,a); subplot(2,1,1);plot(w/pi,mag); title

20、(digital filter Magnitude Response); axis(0,1,0,0.01) subplot(2,1,2);plot(w/pi,db); 13 河北工业大学信息工程学院 信号与线性系统课程设计 title(digital filter Magnitude in DB); axis(0,1,-30,5); %结果：N = 5 % Butterworth Filter Order= 5 %OmegaC = 0.3626 %b = 0.0063 %a = 1.0000 1.1734 0.6884 0.2496 0.0559 0.0063 %N=6 （5）求上述高通滤波器

21、的系统函数及其频谱 b=0.0063; a= 1.0000 1.1734 0.6884 0.2496 0.0559 0.0063; h=impz(b,a); %系统的单位取样响应 figure(1);plot(h) %画出单位取样响应 title(h(n) figure(2) fs=1000; H,f=freqz(b,a,256,fs); %求出系统的频率响应 mag=abs(H); %幅度响应 ph=angle(H); %相位响应 ph=ph*180/pi; subplot(2,1,1),plot(f,mag);grid %画出幅度响应 xlabel(frequency(Hz);ylabel

22、(magnitude); title(|H(jw)|); 14 河北工业大学信息工程学院 信号与线性系统课程设计 subplot(2,1,2);plot(f,ph);grid %画出相位响应 xlabel(frequency(Hz); ylabel(phase); title(相位); figure(3) zr=roots(b) %求出系统的零点 pk=roots(a) %求出系统的极点 zplane(b,a); %zplane函数画出零极点图 %结果：zr = Empty matrix: 0-by-1 %pk = -0.1120 + 0.3438i % -0.1120 - 0.3438i %

23、-0.3674 %-0.2910 + 0.2156i %-0.2910 - 0.2156i 图形如下 系统函数及级联图： 15 河北工业大学信息工程学院 信号与线性系统课程设计 结果分析： 16 河北工业大学信息工程学院 信号与线性系统课程设计 （7）Simulink仿真：（在此只取第一种方案） 图形如下： 17 河北工业大学信息工程学院 信号与线性系统课程设计 技术指标及结果分析： 选作部分： 1、 只截取大约2.5s，三个周期左右，大约800个采样数据进行分析 程序如下： function t,Xn=duqu2(w) fid=fopen(w); C=textscan(fid,%8c %f

24、%*f,headerlines,2); fclose(fid); a=C2; b=C1; k=length(b); for i=1:k c(i)=strread(b(i,:),%*s %f,delimiter,:); end c=c; d=c,a; %截取2.5s的心电信号 for i=1:k 18 河北工业大学信息工程学院 信号与线性系统课程设计 if c(i)<=2.5 e(i,:)=d(i,:); else break; end end t=e(:,1); %时间 Xn=e(:,2); %幅度 调用程序： 图形如下： 19 河北工业大学信息工程学院 信号与线性系统课程设计 20 河

25、北工业大学信息工程学院 信号与线性系统课程设计 结果分析： 2、60Hz工频陷波器设计： 程序如下： %60Hz工频陷波器设计 wp=58,62*2*pi; ws=59.9,60.1*2*pi; H3,f=keshedaizulvboqi(wp,ws,p,s,Xn1); plot(f,abs(H3) axis(0,150,0,1.5) title(60Hz工频陷波器设计) 图形如下： 21 河北工业大学信息工程学院 信号与线性系统课程设计 分析： 课题总结 22 河北工业大学信息工程学院 信号与线性系统课程设计 附录 一、参考文献 1 北京迪阳正泰科技发展公司.综合通信实验系统信号与系统指导书

26、 （第二版）. 2006,6 2 丁玉美.数字信号处理（第二版）.西安电子科技大学出版社，2001 3 吴大正. 信号与线性系统分析（第四版）. 高等教育出版社，2005,8 4 谢嘉奎. 电子线路-线性部分(第四版). 高等教育出版社，2003,2 5 陈后金. 信号分析与处理实验. 高等教育出版社，2006，8 二、 附录设计原理 1.心电信号的读取 txt格式的数据文件内容及格式如图1-1所示（以100.txt为例）。 图1-1 txt格式的心电数据文件 其中文件的第一列为采样时间，第二列是在以MLII这种导联方式所得到的采样数据，第三列是以V5这种导联方式所得到的采样数据，全文件记录了

27、约为10s的心电数据，3600个采样数据，每一行数据之间用Tab符分隔。 由于数据文件中后两列数据是对同一种心电信号进行不同的导联方式所得到的采样数据，所以可以只采用其中的一种采样数据，摒弃另外一种，即可完成对此心电信号的分析。全部的心电文件记录时间约为10s，共计12个左右周期的心电信号。 实际设计心电信号数据文件时应注意： 23 河北工业大学信息工程学院 信号与线性系统课程设计 2.心电信号的线性插值处理 根据上文中提到的插值公式，以此为原理，设计Matlab程序，对心电信号数据做线性插值处理。插值完以后的数据应该是时间均匀的、以0.001秒为间隔的。 此步骤的实现主要是基于Matlab中

28、的数组操作函数来实现，建议一定首先熟悉并掌握Matlab中的所有数组操作函数的作用和操作方法。 其中一种插值方法的思路是：将第一步中读取的心电信号数据的时间数据和幅值数据分别存放在一个一维数组中。然后利用for循环结构把所有数据依次读取进来。判断时间数据数组中前后两个相邻的数据间隔是否为0.001s，如果是则判断下一对相邻两个数据；如果间隔大于0.001s则进行一维插值处理。 注意对时间数据做插值的同时一定不要忘记对幅值数据同样做插值处理，时间数据和幅值数据一定是相互对应的。 3. 滤波器设计 3.1 模拟滤波器设计原理 （1）模拟巴特沃思滤波器原理 巴特沃斯滤波器具有单调下降的幅频特性：在小

29、于截止频率?c的范围内， 具有最平幅度的响应，而在?c后，幅频响应迅速下降。 巴特沃思低通滤波器幅度平方函数为： Ha(j?)2?1 (3-1) ?2N1?()?c 式中N为滤波器阶数，?c为3dB截止角频率。将幅度平方函数写成s的函数形式： (s?) Ha(s)Ha?1 (3-2) 2N1?)j?c 该幅度平方函数有2N个等间隔分布在半径为?c的圆上的极点 sk?ce12k?1j?(?)22N ，k?0,1,.2N?1 为了形成稳定的滤波器，取左半平面的N个极点构成Ha(s)，即： 24 河北工业大学信息工程学院 信号与线性系统课程设计 Ha(s)? Nc ?(s?s) (3-3) k k?

30、0 N?1 为使设计统一，将频率归一化，得到归一化极点pk?e化系统函数为： Ha(p)?12k?1j?(?)22N ，相应的归一 ?(p?p) (3-4) k k?0 N?1 多项式形式为： Ha(p)?(b0?b1p?.?pN) (3-5) （2）模拟切比雪夫滤波器原理 切比雪夫滤波器的幅频特性具有等波纹特性，有两种形式，在通带内等波纹、阻带单调的是?型滤波器，在通带内单调、在阻带内等波纹的是滤波器。以?型滤波器为例。 切比雪夫滤波器的幅度平方函数为： A2(?)?H(j?)2? a 1 2 1?2CN( )?p (3-6) 为小于1的正数，表示通带内幅度波动的程度。p称为通带截止频率。令=/p，称为对p的归一化频率。CN(x)为N阶切比雪夫多项式。幅度平方函数的极点是分布在bp为长半轴，ap为短半轴的椭圆上的点。同样取s平面左半平面的极点构成Ha(s)： Ha(s)? Np ?2 N?1 N?1 ?(s?si) (3