基于MATLAB的心电信号分析系统的设计与仿真

课题二基于MATLAB的心电信号阐发系统的设计与仿真令狐采学摘要：本文是利用MATLAB软件对美国麻省理工学院提供的MITBIH 数据库的122号心电信号病例进行阐发，利用MATLAB软件及simulink平台对122号心电信号的病例进行读取、插值、高通滤波、低通滤波等的处理。

将心电信号中的高频和低频的杂波进行滤除后对插值前后滤波前后的时域波形及频谱进行阐发。

同时也将滤波器的系统函数进行读取，阐发，画出滤波的信号流程图，也阐发各个系统及级联后的系统的冲击响应、幅频响应、相位响应和零极点图来判断系统的稳定性，并用MATLAB软件将图形画出，以便于以后的对系统进行阐发。

关键词：MATLAB，simulink，心电信号，滤波器1.课程设计的目的、意义：本设计课题主要研究数字心电信号的初步阐发及滤波器的应用。

通过完本钱课题的设计，拟主要达到以下几个目的：(1)了解MATLAB软件的特点和使用办法，熟悉基于Simulink 的静态建模和仿真的步调和过程；(2)了解人体心电信号的时域特征和频谱特征；(3)进一步了解数字信号的阐发办法；(4)通过应用具体的滤波器进一步加深对滤波器理解；(5)通过本课题的设计，培养学生运用所学知识阐发和解决实际问题的能力。

2 设计任务及技术指标：设计一个简单的心电信号阐发系统。

对输入的原始心电信号，进行一定的数字信号处理，进行频谱阐发。

采取Matlab语言设计，要求辨别采取两种方法进行仿真，即直接采取Matlab语言编程的静态仿真方法、采取Simulink进行静态建模和仿真的方法。

根据具体设计要求完成系统的法度编写、调试及功能测试。

2.1必做部分：2.1.1读取原始心电信号美国麻省理工学院提供的MITBIH数据库是一个权威性的国际心电图检测标准库，近年来应用广泛，为我国的医学工程界所重视。

MITBIH数据库共有48个病例，每个病例数据长30min，总计约有116000多个心拍，包含有正常心拍和各种异常心拍，内容丰富完整。

为了读取简双便利，采取其txt格式的数据文件作为我们的原心电信号数据。

利用Matlab提供的文件textread或textscan函数，读取txt数据文件中的信号，并且还原实际波形。

2.1.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∆+=- 其中i t 是第i 个数据时间点，Ai 是与之对应的数据，N 是两数据之间需要的插值数，A ∆是需要插值的两点数据差，，，，，，，，，1321321-==N j arraysize i 1111----==i j i j A A t t ，， 1=j i ，时数组j j A t ,依次排列,即获得了插值后等间隔的新数据。

2.1.3根据心电信号的频域特征，设计相应的低通和带通滤波器一般正常人的心电信号频率在0.7～100HZ 规模内，幅度为V 10μ(胎儿)～5mV(成人)。

人体心电信号微弱，信噪比小，因此，在收集心电信号时，易受到仪器、人体活动等因素的影响，并且所收集的心电信号常伴随干扰。

收集心电数据时，由于人的说话呼吸，经常会混有约为0.1Hz 到0.25Hz 频段的干扰，对这些低频干扰，可以让信号通过一个高频滤波器，低截止频率设置为0.25,来滤除低频信号，对高频信号干扰，可以让信号再通过一个低频滤波器，其中截止频率设置为99Hz 。

也可以直接应用带通滤波器设计。

（1）根据以上指标，设计模拟巴特沃斯（切比雪夫）低通、高通或带通滤波器，画出幅频特性（模拟滤波器幅频特性freqs ）。

（2）根据心电信号频谱规模设计一个3阶以上模拟滤波器对心电信号进行预滤波；（3）采取直接、级联或并联方法，实现该系统，并画出系统的信号流图；（4）阐发系统的时域特性（阶跃响应、冲击响应等），并用Matlab 绘出相关波形；（5）用Matlab阐发幅频特性，并绘出相关波形；（6）阐发系统函数零极点与幅频特性的关系。

2.1.4对处理前后的心电信号辨别做频谱阐发利用Matlab软件对处理前后的心电信号编程显示其频谱，阐发比较滤波前后的频谱，得出结论。

如果阐发频谱，滤波效果不明显，则需变动滤波器参数指标，重新设计滤波器。

通过频谱阐发，屡次试验确定最合适的滤波器。

2.1.5 Simulink仿真根据前面的设计，进行基于Simulink的静态仿真设计。

实现心电信号的阐发和处理。

给出系统的基于Simulink的静态建模和仿真的系统方框图，同时记录系统的各个输出点的波形和频谱图。

2.2选作部分：2.2.1减少阐发数据的工作量试验只截取年夜约2.5s，三个周期左右，年夜约800个采样数据进行阐发；2.2.2 50Hz工频陷波器设计由于电子设备收集到的信号经常会混有电源线干扰。

电源线干扰是以50 Hz为中心的窄带噪声，带宽小于1Hz。

设计相应的带阻滤波器滤除电源线干扰，并对处理后的信号做频谱阐发。

3 设计计划论证：4 设计内容（法度清单附带图片）：4.1MATLAB法度清单及图片:4.1.1 提取txt格式心电信号：fid=fopen('122.txt');C=textscan(fid,'%8c %f %*f','headerlines',2); fclose(fid);a=C{1};y=C{2};k=length(a)for i=1:kc(i)=strread(a(i,:),'%*s %f','delimiter',':'); endx=c';plot(x,y)4.1.2 对原始心电信号进行线性插值t=0:0.001:2.5;F=interp1(x,y,t);F=F';t=t';plot(t,F)4.1.3 把数据读到txt中fid = fopen('t.txt','wt');fprintf(fid,'%g\n',t);fclose(fid);fid = fopen('F.txt','wt');fprintf(fid,'%g\n',F);fclose(fid);4.1.4插值前后波形比较subplot(2,2,1)plot(x,y)title('初始信号时域波形') axis([0 2.5 2 1])subplot(2,2,2)fs=1000;N=length(y)n=1:N;f1=n*fs/N;Y1=fft(y);plot(f1,abs(Y1))title('初始信号频谱')axis([0 1000 0 200]) subplot(2,2,3)plot(t,F)title('差值后信号时域波形') axis([0 2.5 2 1])M=length(F);m=1:M;f2=m*fs/M;Y2=fft(F);subplot(2,2,4)plot(f2,abs(Y2))title('插值后信号频谱')axis([0 1000 0 200])4.1.5模拟低通滤波器：wp=60*2*pi;ws=99*2*pi;Rp=1;As=40; [N,wc]=buttord(wp,ws,Rp,As,'s')[B,A]=butter(N,wc,'s')k=0:511;fk=0:1000/512:1000;wk=2*pi*fk;Hk=freqs(B,A,wk);plot(fk,20*log10(abs(Hk)));grid onN =11wc =409.2596B =1.0e+028 *Columns 1 through 100 0 0 0 0 0 0 0 0 0 0 5.3949Columns 11 through 12A =1.0e+028 *Columns 1 through 100.0000 0.0000 0.0000 0.0000 0.00000.0000 0.0000 0.0000 0.0000 0.0008Columns 11 through 120.0926 5.39494.1.6模拟高通滤波器：wp=0.7*2*pi;ws=0.25*2*pi;Rp=0.1;As=40;[N,wc]=buttord(wp,ws,Rp,As,'s')[B0,A0]=butter(N,wc,'s');wph=2*pi*0.25;hk=freqs(B0,A0,wph);[BH,AH]=lp2hp(B,A,wph);[h,w]=freqs(BH,AH);plot(w,20*log10(abs(h)));axis([0,1,80,5]);N =7wc =3.0327B0 =1.0e+003 *0 0 0 0 0 0 0 2.3595A0 =1.0e+003 *0.0010 0.0136 0.0929 0.4070 1.23432.59053.4964 2.35954.1.7滤波前后图形比较：t=0:0.001:9.997;F=interp1(x,y,t);F=F';t=t';figure(1)subplot(3,1,1);plot(1000*t,F);wp=0.7*2*pi;ws=0.25*2*pi;Rp=0.1;As=40;T=1; [N,wc]=buttord(wp,ws,Rp,As,'s')[B,A]=butter(N,wc,'s');[b,a]=imp_invr(B,A,T)[db,mag,pha,w]=freqz_m(b,a);y1=filter(b,a,F);subplot(3,1,2);plot(y1);title('高通滤波后')wp1=2*pi*60;ws1=2*pi*99;Rp1=0.1;As1=40;T1=1000; OmegaP1=wp1/T1;OmegaS1=ws1/T1;[cs1,ds1]=afd_butt(OmegaP1,OmegaS1,Rp1,As1); [b1,a1]=imp_invr(cs1,ds1,T)[db1,mag1,pha1,w1]=freqz_m(b1,a1);y2=filter(b1,a1,y1);subplot(3,1,3);plot(y2);title('低通滤波后')M=length(F);m=1:M;fs=1000;f2=m*fs/M;F1=fft(F);Y1=fft(y1);Y2=fft(y2)figure(2)subplot(3,1,1)plot(f2,abs(F1))axis([0,1000,0,200])title('原始信号频谱_{9.997}')subplot(3,1,2)plot(f2,abs(Y1))axis([0,1000,0,200])title('高通滤波后信号频谱_{9.997}')subplot(3,1,3)plot(f2,abs(Y2))axis([0,1000,0,200])title('低通滤波后信号频谱_{9.997}')4.1.8系统函数的级联、冲击响应、幅度响应、相位响应figure(3)H1=impz(b,a);H2=impz(b1,a1);Hn=conv(H1,H2);%将系统一的系统函数与系统二的函数相级联subplot(3,1,1);plot(H1);title('高通系统函数冲击响应曲线'); %系统函数的冲击响应曲线subplot(3,1,2);plot(H2);title('低通系统函数冲击响应曲线');subplot(3,1,3);plot(Hn);title('级联后系统函数冲击响应曲线');figure(4)fs=1000;[H,f]=freqz(cs,ds,256,fs); %系统一的幅频特性曲线mag=abs(H); %幅度响应ph=angle(H); %相位响应ph=ph*180/pi;subplot(2,2,1),plot(f,mag);grid %h1的幅度响应title('h1幅度响应')%xlabel('frequency(Hz)');ylabel('magnitude');subplot(2,2,2);plot(f,ph);grid %h1的相位响应title('h1相位响应 ')fs=1000;[H,f]=freqz(BH,AH,256,fs);mag=abs(H);ph=angle(H);ph=ph*180/pi;subplot(2,2,3),plot(f,mag);gridtitle('h2幅度响应')subplot(2,2,4);plot(f,ph);gridtitle('h2相位响应')figure(6)zr=roots(B) %系统一的零点图pk=roots(A) %系统一的极点图zplane(B,A); %zplane函数画出系统一的零极点图figure(7)zr1=roots(cs1) %系统二的零点图pk1=roots(ds1) %系统二的极点图zplane(cs1,ds1); %zplane函数画出系统一的零极点图4.1.9截取到2.5s对截取的部分进行滤涉及频谱阐发t1=0:0.001:2.5;F0=interp1(x,y,t1);F0=F0';t1=t1';figure(8)subplot(3,1,1);plot(1000*t1,F0);wp=0.7*2*pi;ws=0.25*2*pi;Rp=0.1;As=40;T=1;[N,wc]=buttord(wp,ws,Rp,As,'s')[B,A]=butter(N,wc,'s');[b,a]=imp_invr(B,A,T)[db,mag,pha,w]=freqz_m(b,a);y11=filter(b,a,F0);subplot(3,1,2);plot(y11);title('高通滤波后_{2.5}')wp1=2*pi*60;ws1=2*pi*99;Rp1=0.1;As1=40;T1=1000; OmegaP1=wp1/T1;OmegaS1=ws1/T1;[cs1,ds1]=afd_butt(OmegaP1,OmegaS1,Rp1,As1); [b1,a1]=imp_invr(cs1,ds1,T)[db1,mag1,pha1,w1]=freqz_m(b1,a1);y21=filter(b1,a1,y11);subplot(3,1,3);plot(y21);title('低通滤波后_{2.5}')M=length(F0);m=1:M;fs=1000;f2=m*fs/M;F01=fft(F0);Y11=fft(y11);Y21=fft(y21)figure(9)subplot(3,1,1)plot(f2,abs(F01))axis([0,1000,0,200])title('原始信号频谱_{2.5}')subplot(3,1,2)plot(f2,abs(Y11))axis([0,1000,0,200])title('高通滤波后信号频谱_{2.5}')subplot(3,1,3)plot(f2,abs(Y21))axis([0,1000,0,200])title('低通滤波后信号频谱_{2.5}')4.2．Simulink仿真参数设置及波形比较：4.2.1滤波器参数设置4.2.2滤波前后的图形进行比较原始波形低通滤波后高通滤波后5 实验结果与阐发：6 总结7.实验收获及心得体会参考文献：[1] 北京迪阳正泰科技成长公司.综合通信实验系统——信号与系统指导书（第二版）. ,6[2] 丁玉美.数字信号处理（第三版）.西安电子科技年夜学出版社，[3] 吴年夜正. 信号与线性系统阐发（第四版）. 高等教育出版社，,8[4] 谢嘉奎. 电子线路—非线性部分(第五版). 高等教育出版社，,2[5] 陈后金. 信号阐发与处理实验. 高等教育出版社，，8。