1、信号与线性系统课程设计报告课题二 心电信号分析系统的设计与仿真班级:通信 C114姓名：xxx学号：xxx成绩：指导教师：日期：2013.1.5摘要：本课题设计了一个简单的心电信号分析系统。直接采用 Matlab 语言编程的静态仿真方式、采用 Simulink 进行动态建模和仿真的方式，对输入的原始心电信号，进行线性插值处理，并通过 matlab 语言编程设计对其进行时域和频域的波形频谱分析，根据具体设计要求完成系统的程序编写、调试及功能测试。得出一定的结论。 关键字：matlab、心电信号提取、线性插值、滤波、simulink 仿真。本课题的目的本设计课题主要研究数字心电信号的初步分析方法及
2、滤波器的设计与应用。通过完成本课题的任务，拟主要达到以下几个目的：1了解 MATLAB 软件的特点和使用方法，熟悉基于 Simulink 的动态建模和仿真的步骤和过程；2. 了解 LabVIEW 虚拟仪器软件的特点和使用方法，熟悉采用 LabVIEW 进行信号分析、系统设计及仿真的方法。3了解人体心电信号的时域特征和频谱特征；4通过设计具体的滤波器进一步加深对滤波器性能的理解；5掌握数字心电信号的分析方法，学会系统设计与软件仿真方法；6通过本课题的训练，培养学生运用所学知识分析和解决实际问题的能力。2 设计任务及技术指标设计一个简单的心电信号分析系统。其基本功能包括：输入原始心电信号，对其做一
3、定的数字信号处理，进行时域显示、分析及频谱分析。采用 Matlab 软件（或 LabVIEW 虚拟仪器软件）设计相关程序。对基于 Matlab 软件的程序设计，要求分别采用两种方式进行仿真，即直接采用 Matlab 语言编程的静态系统仿真方式、采用 Simulink 进行动态建模仿真的方式。根据心电信号的具体特性参数设计系统各功能模块的源程序，进行调试。1.对原始数字心电信号进行读取，由数字信号数据绘制出其时域波形并加以分析。2.对数字信号数据做一次线性插值，使其成为均匀数字信号，以便后面的信号分析。3.根据心电信号的频域特征（自己查阅相关资料） ，设计相应的滤波器去除噪声。4.绘制进行信号处
4、理前后的频谱，做频谱分析，得出相关结论。5.对系统功能进行综合测试，整理数据，撰写设计报告。3 主要设备和软件1PC 机一台。2. MATLAB6.5 以上版本软件，一套。4 设计内容以及实验结果与分析读取文件的 function 函数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:kc(i)=strread(b(i,:),%*s %f,delimiter,:);endc=c;d=c,a;t=d
5、(:,1); Xn=d(:,2);保存文件 1function 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);保存文件 2function 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);插
6、值文件function t2,Xn2=chazhi1(t,Xn)n=0;y=0;t=t.*1000; m=length(t);for i=1:me(i)=round (t(i);endfor i=1:(length(t)-1)if(e(i+1)-e(i)=1)N=(e(i+1)-e(i)/1;A=(Xn(i+1)-Xn(i)/N;for j=1:Nz(y+j,1)=e(i)+(j-1)*1;z(y+j,2)=Xn(i)+(j-1)*An=n+1;j=j+1;endy=n;endi=i+1;endz(y+1,2)=Xn(i);z(y+1,1)=t(i);t2=z(:,1);t2=t2./1000
7、;Xn2=z(:,2);afd_butt 文件function b,a = afd_butt(Wp,Ws,Rp,As);if Wp = 0error(Passband edge must be larger than 0)endif Ws = Wperror(Stopband edge must be larger than Passband edge)endif (Rp = 0) | (As 0)error(PB ripple and/or SB attenuation ust be larger than 0)endN = ceil(log10(10(Rp/10)-1)/(10(As/10
8、)-1)/(2*log10(Wp/Ws);fprintf(n* Butterworth Filter Order = %2.0f n,N)OmegaC = Wp/(10(Rp/10)-1)(1/(2*N);b,a=u_buttap(N,OmegaC);freqz_m文件functiondb,mag,pha,grd,w=freqz_m(b,a)H,w=freqz(b,a,1000,whole); H=(H(1:1:501);w=(w(1:501); mag=abs(H); db=20*log10(mag+eps)/max(mag); pha=angle(H); grd=grpdelay(b,a,
9、w); u_buttap文件function b,a = u_buttap(N,Omegac);% Unnormalized Butterworth Analog Lowpass Filter Prototype% -% b,a = u_buttap(N,Omegac);% b = numerator polynomial coefficients of Ha(s)% a = denominator polynomial coefficients of Ha(s)% N = Order of the Butterworth Filter% Omegac = Cutoff frequency i
10、n radians/sec%z,p,k = buttap(N);p = p*Omegac;k = k*OmegacN;B = real(poly(z);b0 = k;b = k*B;a = real(poly(p);一.提取 txt 格 式 心 电 信 号 ：fid=fopen(1.txt);C=textscan(fid,%8c %f %*f,headerlines,2);fclose(fid);a=C2;b=C1;k=length(b);for i=1:kc(i)=strread(b(i,:),%*s %f,delimiter,:);endc=c;plot(c,a)二 .对 原 数 据 进
11、行 插 值 ：tict,Xn=duquexinhao1(1.txt)baocun1(t,Xn);figure(1)plot(t,Xn)t2,Xn2=chazhi1(t,Xn);baocun2(t2,Xn2);figure(2)plot(t2,Xn2)toc0 1 2 3 4 5 6 7 8 9 10-1-0.500.511.520 1 2 3 4 5 6 7 8 9 10-1-0.500.511.52三.插值前后波形对比：subplot(2,2,1) plot(t,Xn) title(初始信号时域波形 ) axis(0 2.5 -1 2) subplot(2,2,2) fs=1000; N=l
12、ength(Xn) n=1:N; f1=n*fs/N; Y1=fft(Xn); plot(f1,abs(Y1) title(初始信号频谱 ) axis(0 1000 0 200) subplot(2,2,3) plot(t2,Xn2) title(插值后信号时域波形 ) axis(0 2.5 -1 2) M=length(Xn2); m=1:M; f2=m*fs/M; Y2=fft(Xn2); subplot(2,2,4) plot(f2,abs(Y2) title(插值后信号频谱 ) axis(0 1000 0 200) tt=t2,Xn20 0.5 1 1.5 2 2.5-1012 位位位
13、位位位位位0 500 1000050100150200 位位位位位位0 0.5 1 1.5 2 2.5-1012 位位位位位位位位位0 500 1000050100150200 位位位位位位位四 .把 数 据 读 到 txt 中 ：fid = fopen(1.txt,wt); fprintf(fid,%gn,t); fclose(fid); fid = fopen(F.txt,wt); fprintf(fid,%gn,Xn2); fclose(fid); 五 .模 拟 低 通 滤 波 器 ：wp=50*2*pi;ws=99*2*pi;Rp=1;As=40;%设置滤波器参数N,wc=buttor
14、d(wp,ws,Rp,As,s) %计算滤波器阶数N 和3dB 截止频率wc 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); axis(0,200,-40,5);结 果N =8wc =349.7984B =1.0e+020 *Columns 1 through 6 0 0 0 0 0 0Columns 7 through 9 0 0 2.2415A =1.0e+020 *Columns 1 through
15、6 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000Columns 7 through 9 0.0002 0.0328 2.24150 20 40 60 80 100 120 140 160 180 200-40-35-30-25-20-15-10-505六 .模 拟 高 通 滤 波 器 ：wp=0.6*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;BH,AH=lp2hp(B0,A0,wph); h,w=freqs(
16、BH,AH); plot(w,20*log10(abs(h); axis(0,1,-80,5);结果N = 8wc =2.79330 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-80-70-60-50-40-30-20-100七 .滤 波 前 后 图 形 对 比 ：fid=fopen(1.txt);C=textscan(fid,%8c %f %*f,headerlines,2);fclose(fid);a=C2;b=C1;k=length(b);for i=1:kc(i)=strread(b(i,:),%*s %f,delimiter,:);endc=c;fig
17、ure(1)subplot(3,1,1); plot(c,a) wp=0.6*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=impinvar(B,A,T) db,mag,pha,w=freqz_m(b,a); y1=filter(b,a,Xn2); subplot(3,1,2);plot(y1); title(高通滤波后 ) wp1=2*pi*50;ws1=2*pi*99;Rp1=0.1;As1=40;T1=1000; OmegaP1=wp1/T1;OmegaS1=
18、ws1/T1; cs1,ds1=afd_butt(OmegaP1,OmegaS1,Rp1,As1); b1,a1=impinvar(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(Xn2); m=1:M; fs=1000; f2=m*fs/M; F1=fft(Xn2); Y1=fft(y1); Y2=fft(y2) figure(2) subplot(3,1,1) plot(f2,abs(F1) axis(0,1000
19、,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)0 1 2 3 4 5 6 7 8 9 10-10120 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000-202 位位位位位0 1000 2000 3000 4000 5000 6000 7000
20、8000 9000 10000-202 位位位位位0 100 200 300 400 500 600 700 800 900 10000100200位位位位位位9.9970 100 200 300 400 500 600 700 800 900 10000100200位位位位位位位位位9.9970 100 200 300 400 500 600 700 800 900 10000100200位位位位位位位位位9.997八 .系 统 函 数 的 级 联 、 冲 击 响 应 、 幅 度 响 应 、 相 位 响 应 ：figure(3) H1=impz(b,a); H2=impz(b1,a1); H
21、n=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(B,A,256,fs); %系统一的幅频特性曲线mag=abs(H); %幅度响应ph=angle(H); %相位响应ph=ph*180/pi; subplot(2
22、,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);grid title(h2幅度响应 ) subplot(2,2,4);plot(f,ph);grid title(h2相位
23、响应 ) figure(5) zr=roots(B) %系统一的零点图pk=roots(A) %系统一的极点图zplane(B,A); %zplane函数画出系统一的零极点图figure(6) zr1=roots(cs1) %系统二的零点图pk1=roots(ds1) %系统二的极点图zplane(cs1,ds1); %zplane函数画出系统一的零极点图0 2 4 6 8 10 12 14 16 18-0.500.51 位位位位位位位位位位位位0 20 40 60 80 100 120 140 160 180-0.200.2 位位位位位位位位位位位位0 20 40 60 80 100 120
24、 140 160 180 200-0.200.2 位位位位位位位位位位位位位0 200 400 60002468 h1位位位位0 200 400 600-150-100-500 h1位位位位 0 200 400 60005101520 h2位位位位0 200 400 600050100150200 h2位位位位-4 -3 -2 -1 0 1 2-3-2-10123Real PartImaginaryPart-1 -0.5 0 0.5 1-1-0.8-0.6-0.4-0.200.20.40.60.8110Real PartImaginaryPart9.Simulink 仿真以及仿真结果图高通滤波
25、器：低通滤波器：带阻滤波器：原信号波形：滤波后信号波形：滤波后滤波前简要分析：5 设计方案论证1.心 电 信 号 读 取 设 计 方 案 论 证 与 选 择美国麻省理工学院提供的 MIT-BIH 数据库是一个权威性的国际心电图检测标准库，近年来应用广泛，为我国的医学工程界所重视。MIT-BIH 数据库共有48 个病例，每个病例数据长 30min，总计约有 116000 多个心拍，包含有正常心拍和各种异常心拍，内容丰富完整。为了读取简单方便，采用其 txt 格式的数据文件作为我们的原心电信号数据。利用 Matlab 提供的文件 textread 或 textscan 函数，读取 txt 数据文件
26、中的信号，并且还原实际波形。 （或者利用 labVIEW 提供的文件 I/O 函数，读取 txt 数据文件中的信号，并且还原实际波形。 ）由于不熟悉 Labview 提供的文件 I/O 函数，所以采用老师提供的 Matlab 提供的文件 textscan。2.线性插值的设计方案与选择由于原始心电信号数据不是通过等间隔采样得到的，也就是说原始的心电数据并不是均匀的，而用 Matlab 中提供的数字滤波器处理数据时，要求数据是等间隔的。因此设计的系统首先应对原始心电信号做线性插值处理，使其变为等间隔的数字信号，否则直接处理后会出现偏差，根据心电信号的特点, 把时间分隔成 0.001s。添加的幅值点
27、采用一次线性插值。对二维数据进行插值 ,相连幅值间数据的插值根据时间进行,运算公式如下: ， ， ， ，1it01./tN1iA01.jt NAjj/1其中 是第 i 个数据时间点，A i 是与之对应的数据， N 是两数据之间需要的插值it数， 是需要插值的两点数据差， ， 3232jarysize， 时数组 依次排列,即得到了插值后等间隔的新11ijijt， j， jt数据。虽然 Labview 提供的插值法有很多，但是由于自身对 Labview 不熟悉，所以采用教师提供的线性插值法（即上述所列公式） 。3.滤波器的设计方案论证与选择一般正常人的心电信号频率在 0.7100HZ 范围内，幅度
28、为 (胎儿)V105mV(成人)。人体心电信号微弱，信噪比小，因此，在采集心电信号时，易受到仪器、人体活动等因素的影响，而且所采集的心电信号常伴有干扰。采集心电数据时，由于人的说话呼吸，常常会混有约为 0.1Hz 到 0.25Hz 频段的干扰，对于这些低频干扰，可以让信号通过一个高频滤波器，低截止频率设置为0.25，来滤除低频信号，对于高频信号干扰，可以让信号再通过一个低频滤波器，其中截止频率设置为 99Hz。所以我选择 butterworth 滤波器。6 总结7 参考文献1 北京迪阳正泰科技发展公司.综合通信实验系统信号与系统指导书（第二版）. 2006,62 丁玉美.数字信号处理（第二版）
29、.西安电子科技大学出版社，20013 吴大正. 信号与线性系统分析（第四版）. 高等教育出版社，2005,84 谢嘉奎. 电子线路-线性部分(第四版). 高等教育出版社，2003,25 陈后金. 信号分析与处理实验. 高等教育出版社，2006，86陈锡辉，张银鸿编著.LabVIEW 8.20 程序设计从入门到精通M.北京：清华大学出版社，2007.8 附录1.心电信号的读取txt 格式的数据文件内容及格式如图 1-1 所示（以 100.txt 为例） 。图 1-1 txt 格式的心电数据文件其中文件的第一列为采样时间，第二列是在以 MLII 这种导联方式所得到的采样数据，第三列是以 V5 这种
30、导联方式所得到的采样数据，全文件记录了约为 10s 的心电数据，3600 个采样数据，每一行数据之间用 Tab 符分隔。由于数据文件中后两列数据是对同一种心电信号进行不同的导联方式所得到的采样数据，所以可以只采用其中的一种采样数据，摒弃另外一种，即可完成对此心电信号的分析。全部的心电文件记录时间约为 10s，共计 12 个左右周期的心电信号。实际设计心电信号数据文件时应注意：2.心电信号的线性插值处理根据上文中提到的插值公式，以此为原理，设计 Matlab 程序，对心电信号数据做线性插值处理。插值完以后的数据应该是时间均匀的、以 0.001 秒为间隔 的。此步骤的实现主要是基于 Matlab
31、中的数组操作函数来实现，建议一定首先熟悉并掌握 Matlab 中的所有数组操作函数的作用和操作方法。其中一种插值方法的思路是：将第一步中读取的心电信号数据的时间数据和幅值数据分别存放在一个一维数组中。然后利用 for 循环结构把所有数据依次读取进来。判断时间数据数组中前后两个相邻的数据间隔是否为 0.001s，如果是则判断下一对相邻两个数据；如果间隔大于 0.001s 则进行一维插值处理。注意对时间数据做插值的同时一定不要忘记对幅值数据同样做插值处理，时间数据和幅值数据一定是相互对应的。3. 滤波器设计3.1 模拟滤波器设计原理（1）模拟巴特沃思滤波器原理巴特沃斯滤波器具有单调下降的幅频特性：
32、在小于截止频率 的范围内，c具有最平幅度的响应，而在 后，幅频响应迅速下降。c巴特沃思低通滤波器幅度平方函数为： (3-1)221()()aNcHj式中 N 为滤波器阶数， 为 3dB 截止角频率。将幅度平方函数写成 s 的函c数形式：(3-2)21()()a Ncssj该幅度平方函数有 2N 个等间隔分布在半径为 的圆上的极点， 为了形成稳定的滤波器，取左半平面的 N12() kjNkcse0,1.个极点构成 ，即：aHs(3-3)10()()NackkHss为使设计统一，将频率归一化，得到归一化极点 ，相应的归12()kjNkpe一化系统函数为： (3-4)10()()Nakkp多项式形式
33、为： 01().)NaHpbp(3-5)（2）模拟切比雪夫滤波器原理切比雪夫滤波器的幅频特性具有等波纹特性，有两种形式，在通带内等波纹、阻带单调的是 型滤波器，在通带内单调、在阻带内等波纹的是滤波器。以 型滤波器为例。切比雪夫滤波器的幅度平方函数为： (3-6) 2221()()aNpAHjC 为小于 1 的正数，表示通带内幅度波动的程度。p 称为通带截止频率。令 =/p，称为对 p 的归一化频率。 CN(x)为 N 阶切比雪夫多项式。幅度平方函数的极点是分布在 bp 为长半轴，ap 为短半轴的椭圆上的点。同样取 s平面左半平面的极点构成 ：()Hs(3-7)12()NNapiis进行归一化，
34、得到： 1()iip其中 ， 21(2)sincoskkpchjhiNN 1()Arsh3.2 模拟滤波器数字化原理将模拟滤波器转化为数字滤波器在工程上常用的有脉冲响应不变法和双线性变换法。脉冲响应不变法是一种时域上的转换方法，它是数字滤波器的单位取样响应在采样点上等于模拟滤波器的单位冲激响应，即：(3-8)()hnTa设模拟滤波器只有单阶极点，其系统函数为：(3-9)1()NiaiAHss对 进行拉氏反变换得到 ，对 进行等间隔采样，得到()aHsht()at，对 进行 Z 变换，得到数字滤波器系统函数：()hnT()hn(3-10)11()iNsTiAzez这种方法 s 和 z 的关系是：
35、 。该方法的优点是频率坐标变换时线性的切数字滤波器的单位脉冲响应完全模仿模拟滤波器的单位冲激响应，时域特性逼近好；缺点是会产生频谱混叠现象，适合低通、带通滤波器的设计，不适合高通、带阻滤波器的设计。 双线性变换法为了克服频谱混叠现象，采用非线性频率压缩方法，将整个频率轴上的频率范围压缩到 之间，再用 转换到 Z 平面上。/TsTze这种方法 s 和 z 的关系是： 。该方法克服了频谱混1(2)/)s叠现象，但带来了频率坐标变换的非线性： ，由模拟滤波器(2tan(/系统函数转换为数字滤波器系统函数公式为：(3-11)12()|azTHzs3.3 数字高通、带通、带阻滤波器的设计上述滤波器可以借助于模拟滤波器的频率变换设计一个所需类型的模拟滤波器， 再通过双线性变换法将其转换成所需类型的数字滤波器。首先确定所需类型数字滤波器的技术指标；然后将数字滤波器技术指标按照公式 转换成所需类型滤波器的模拟域技术指标；将所需类(2/)tan(/)T型滤波器的模拟域技术指标转换成低通滤波器技术指标；设计归一化模拟低通滤波器；去归一化得到模拟低通滤波器的系统函数；将模拟低通滤波器转换为所需类型的模拟滤波器；最后通过双线性变换法转换成所需类型的数字滤波器。9 课程设计心得体会
