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

河北工业大学信息工程学院 信号与线性系统课程设计 基于MATLAB的心电信号分析摘要： 本课题设计了一个简单的心电信号分析系统。直接采用Matlab语言编程的静态仿真方式、采用Simulink进行动态建模和仿真的方式，对输入的原始心电信号，进行线性插值处理，并通过matlab语言编程设计对其进行时域和频域的波形频谱分析，根据具体设计要求完成系统的程序编写、调试及功能测试。得出一定的结论。 关键字：matlab、心电信号提取、线性插值、滤波、simulink仿真。一、课题目的及意义心电信号是人类最早研究并应用于医学临床的生物信号之一，它比其它生物电信号更易于检测，并且具有较直观的规律性，因而心电图分析技术促进了医学的发展。然而，心电图自动诊断还未广泛应用于临床，从国内外的心电图机检测分析来看，自动分析精度还达不到可以替代医生的水平，仅可以为临床医生提供辅助信息。其主要原因是心电波形的识别不准，并且心电图诊断标准不统一。因此，探索新的方法以提高波形识别的准确率，寻找适合计算机实现又具诊断价值的诊断标准，是改进心电图自动诊断效果，扩大其应用范围的根本途径。如何把心电信号的特征更加精确的提取出来进行自动分析，判断出其异常的类型成了亟待解决的焦点问题。本课题通过matlab语言编程，对原始心电信号进行一定的分析处理。二、课题任务及要求1、必做部分（1）利用Matlab对MIT-BIH数据库提供的数字心电信号进行读取，并还原实际波形。（2）对原始心电信号做线性插值（3）对处理前后的心电信号分别做频谱分析 利用Matlab软件对处理前后的心电信号编程显示其频谱，分析比较滤 波前后的频谱，得出结论。（4） Simulink仿真 根据前面的设计，进行基于Simulink的动态仿真设计。实现心电信号 的分析和处理。2、 选作部分（1） 只截取大约2.5s，三个周期左右，大约800个采样数据进行分析。（2）60Hz工频陷波器设计三、设计技术指标四、设计方案论证1、必做部分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,:);endc=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,%gn,t); fclose(fid);fid = fopen(Xn.txt,wt);fprintf(fid,%gn,Xn); fclose(fid);%保存插值后的信号function baocun2(t1,Xn1)fid = fopen(t1.txt,wt);fprintf(fid,%gn,t1); fclose(fid);fid = fopen(Xn1.txt,wt);fprintf(fid,%gn,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);主函数如下（1） 、将信号通过低通、高通、带阻滤波器程序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(滤波后信号)所出图形如下结果分析：（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(H1) ; %带通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(带通滤波器)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)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(滤波后截取一部分信号) 所出图形如下结果分析：（3） 、将信号通过低通、高通组合成的带通滤波器程序t,Xn=duquexinhao1(117.txt);baocun1(t,Xn) %保存读取信号t1,Xn1=xianxingchazhi(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); %低通和高通组合的带通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,260)title(滤波后心电信号频谱)figure(3)subplot(2,1,1)plot(t1,Xn1)title(滤波前信号)subplot(2,1,2)plot(t1,y)title(滤波后信号) 所出图形如下结果分析：三种方案比较分析：（4）系统零极点分析（在此只以高通滤波器为例）%求高通滤波器的阶数及分子分母系数wp=0.1*pi;ws=0.25*pi;Rp=1;As=35;T=1;%数字指标OmegaP=(2/T)*tan(wp/2);%通带模拟频率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(digital filter Magnitude Response); axis(0,1,0,0.01)subplot(2,1,2);plot(w/pi,db);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）求上述高通滤波器的系统函数及其频谱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(magnitude);title(|H(jw)|);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 %-0.3674 %-0.2910 + 0.2156i %-0.2910 - 0.2156i图形如下系统函数及级联图：结果分析：（7）Simulink仿真：（在此只取第一种方案）图形如下：技术指标及结果分析：选作部分：1、 只截取大约2.5s，三个周期左右，大约800个采样数据进行分析程序如下：function t,Xn=duqu2(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,:);endc=c;d=c,a;%截取2.5s的心电信号for i=1:k if c(i)=2.5 e(i,:)=d(i,:); else break; endendt=e(:,1); %时间Xn=e(:,2); %幅度调用程序：图形如下：结果分析：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工频陷波器设计)图形如下：分析：课题总结附录一、参考文献1 北京迪阳正泰科技发展公司.综合通信实验系统信号与系统指导书（第二版）. 2006,62 丁玉美.数字信号处理（第二版）.西安电子科技大学出版社，20013 吴大正. 信号与线性系统分析（第四版）. 高等教育出版社，2005,84 谢嘉奎. 电子线路-线性部分(第四版). 高等教育出版社，2003,25 陈后金. 信号分析与处理实验. 高等教育出版社，2006，8二、 附录设计原理1.心电信号的读取txt格式的数据文件内容及格式如图1-1所示（以100.txt为例）。 图1-1 txt格式的心电数据文件其中文件的第一列为采样时间，第二列是在以MLII这种导联方式所得到的采样数据，第三列是以V5这种导联方式所得到的采样数据，全文件记录了约为10s的心电数据，3600个采样数据，每一行数据之间用Tab符分隔。由于数据文件中后两列数据是对同一种心电信号进行不同的导联方式所得到的采样数据，所以可以只采用其中的一种采样数据，摒弃另外一种，即可完成对此心电信号的分析。全部的心电文件记录时间约为10s，共计12个左右周期的心电信号。实际设计心电信号数据文件时应注意：2.心电信号的线性插值处理根据上文中提到的插值公式，以此为原理，设计Matlab程序，对心电信号数据做线性插值处理。插值完以后的数据应该是时间均匀的、以0.001秒为间隔的。此步骤的实现主要是基于Matlab中的数组操作函数来实现，建议一定首先熟悉并掌握Matlab中的所有数组操作函数的作用和操作方法。其中一种插值方法的思路是：将第一步中读取的心电信号数据的时间数据和幅值数据分别存放在一个一维数组中。然后利用for循环结构把所有数据依次读取进来。判断时间数据数组中前后两个相邻的数据间隔是否为0.001s，如果是则判断下一对相邻两个数据；如果间隔大于0.001s则进行一维插值处理。注意对时间数据做插值的同时一定不要忘记对幅值数据同样做插值处理，时间数据和幅值数据一定是相互对应的。3. 滤波器设计3.1 模拟滤波器设计原理（1）模拟巴特沃思滤波器原理 巴特沃斯滤波器具有单调下降的幅频特性：在小于截止频率的范围内，具有最平幅度的响应，而在后，幅频响应迅速下降。巴特沃思低通滤波器幅度平方函数为： (3-1)式中N为滤波器阶数，为3dB截止角频率。将幅度平方函数写成s的函数形式： (3-2) 该幅度平方函数有2N个等间隔分布在半径为的圆上的极点， 为了形成稳定的滤波器，取左半平面的N个极点构成，即： (3-3) 为使设计统一，将频率归一化，得到归一化极点，相应的归一化系统函数为： (3-4)多项式形式为： (3-5) （2）模拟切比雪夫滤波器原理 切比雪夫滤波器的幅频特性具有等波纹特性，有两种形式，在通带内等波纹、阻带单调的是型滤波器，在通带内单调、在阻带内等波