报告.docx
《报告.docx》由会员分享，可在线阅读，更多相关《报告.docx（59页珍藏版）》请在冰豆网上搜索。
报告
信号与线性系统课程设计报告
（课题二心电信号分析系统的设计与仿真）
成绩：
指导教师：
刘翠响
日期：
2016年1月7日
心电信号分析系统的设计与仿真
摘要：
心脏病发病率呈日渐上升趋势，是威胁全人类健康的首要疾病之一。
心电信号分析是判断心脏疾病最直接最有效的手段，研究心电信号的处理和分析系统是非常有意义的。
本文针对心电信号特征,采用MATLAB和LabVIEW两种软件实现心电信号分析系统的设计与仿真，主要实现对心电信号进行去噪和预处理。
MATLAB的设计包括对心电信号进行提取、线性插值、滤波和分析，实现simulink仿真，并且利用GUI做出系统的图形用户界面设计。
LabVIEW的设计包括在完成上述要求的基础上使用选项卡等控件做出简洁友好的人机交互界面。
关键词：
心电信号，预处理，MATLAB，LabVIEW，人机交互界面
1课程设计的目的及意义
1.1课程设计的目的
1．了解MATLAB软件的特点和使用方法，熟悉基于Simulink的动态建模和仿真的步骤和过程；
2．了解LabVIEW虚拟仪器软件的特点和使用方法，熟悉采用LabVIEW进行信号分析、系统设计及仿真的方法。
3．了解人体心电信号的时域特征和频谱特征；
4．通过设计具体的滤波器进一步加深对滤波器性能的理解；
5．掌握数字心电信号的分析方法，学会系统设计与软件仿真方法；
6．通过本课题的训练，培养学生运用所学知识分析和解决实际问题的能力。
1.2课程设计的意义
心电信号是组成人体心脏的肌细胞电活动在体表的综合表现,在一定程度上客观反映了心脏各部位的生理状况,因此,对心电信号的研究处理具有实用价值和重要意义。
但由于心电信号自身特征的复杂性,在对其研究处理中仍存在一些有待于进一步深入和解决的问题。
本课题通过使用Matlab和LabVIEW软件，实现简单的心电信号分析系统的设计及仿真，完成对心电信号的预处理及简单分析。
2课程设计任务及要求
2.1基于Matlab的简单心电信号分析系统设计
设计一个简单的心电信号分析系统。
其基本功能包括：
输入原始心电信号，对其做一定的数字信号处理，进行时域显示、分析及频谱分析。
采用Matlab软件设计相关程序。
对基于Matlab软件的程序设计，要求分别采用两种方式进行仿真，即直接采用Matlab语言编程的静态系统仿真方式、采用Simulink进行动态建模仿真的方式。
根据心电信号的具体特性参数设计系统各功能模块的源程序，进行调试。
具体要求如下：
1.对原始数字心电信号进行读取，由数字信号数据绘制出其时域波形并加以分析。
2.对数字信号数据做一次线性插值，使其成为均匀数字信号，以便后面的信号分析。
3.根据心电信号的频域特征（自己查阅相关资料），设计相应的滤波器去除噪声。
4.绘制进行信号处理前后的频谱，做频谱分析，得出相关结论。
5.使用GUI进行系统的图形用户界面设计，（包含以上功能）。
2.2基于LabVIEW虚拟仪器的简单心电信号分析系统设计
1.进行心电信号的频谱分析，根据心电信号的频域特征（自己查阅相关资料），设计相应的滤波器去除噪声。
要求给出系统的前面板和框图，并记录仿真结果。
2.根据心电信号的特征，针对系统进行功能拓展，记录仿真结果，并进行相应的分析。
3设计方案及论证
3.1基于Matlab的简单心电信号分析系统设计
3.1.1设计流程图
图1-1基于MATLAB的心电信号分析系统设计流程图
3.1.2设计原理
1．心电信号的读取：
文件的第一列为采样时间，后两列是不同的导联方式所得到的采样数据，所以只读取前两列数据，并且忽略前两行解释说明文字，分别存放在两个数组中。
2．插值：
把时间分隔成0.001s，添加的幅值点采用一次线性插值。
对二维数据进行插值,相连幅值间数据的插值根据时间进行。
对时间数据做插值的同时对幅值数据同样做插值处理，时间数据和幅值数相互对应。
3．根据心电信号的频域特征，设计相应滤波器：
一般正常人的心电信号频率在0.7～100HZ范围内，幅度为10（胎儿）～5mV（成人），采集心电数据时，由于人的说话呼吸，常常会混有约为0.1Hz到0.25Hz频段的干扰。
方案一：
根据指导书中的要求，设置一个巴特沃斯模拟低通滤波器滤除高频信号干扰，设置一个巴特沃斯模拟高通滤波器来滤除低频信号干扰。
低频滤波器的设计指标为：
通带截止频率为92Hz，阻带起始频率为99Hz，通带最大衰减为1dB，阻带最小衰减为30dB。
高频滤波器的设计指标为：
通带起始频率为0.7Hz，阻带起始频率为0.25Hz,通带最大衰减为1dB，阻带最小衰减为30dB。
方案二：
根据心电信号频谱范围，考虑到过渡带的衰减，我们设计了一个巴特沃斯数字带通滤波器并通过多次试验确定最合适的设计指标，滤波器的阻带范围为0~0.25Hz和170Hz及以上，通带范围为1.7Hz~90Hz。
将99Hz设在大概3dB截止频率处，这样能有效的滤除低频和高频信号干扰，又能避免有用的频带被滤除。
4．去除50Hz工频干扰：
由于电子设备采集到的信号经常会混有电源线干扰。
电源线干扰是以50Hz为中心的窄带噪声，带宽小于1Hz。
可以通过合适的带阻滤波器滤除电源线干扰。
方案一：
从零极点分布对系统频率特性的影响出发,设计一个数字陷波器。
系统函数完全由它的全部零点和极点来确定。
对于因果稳定系统，极点必须全部在单位圆内。
极点位置主要影响频响的峰值位置及尖锐程度，极点越靠近单位圆，峰值越高越尖锐；零点位置主要影响频响的谷点位置及形状，零点越靠近单位圆，谷值越接近零。
对进行陷波,则取零点。
但零点引起的凹谷会影响滤波器的通带范围，所以配置相应的极点来抵消这种影响，为了保证系统的因果稳定性，极点必须在单位圆内。
取极点,应小于1且接近于1。
于是,所设计的陷波器的系统函数为
对50Hz进行陷波，则），其中是采样频率，本实验中取，所以为，通过试验当取0.95时，滤波效果最好。
方案二：
根据工频干扰的频带特征，设计一个巴特沃斯数字带阻滤波器，技术指标为：
通带范围为0~45Hz和55Hz及以上，阻带范围为49.5Hz~50.5Hz，通带最大衰减2dB，阻带最小衰减30dB。
5．截取原始信号的约3个周期：
根据实验时选用的具体心电信号数据，截取心电信号中不正常的部分约3个周期。
截取相应的时间段，再根据截取部分在时间数组中的位置，截取相应的信号幅值。
6．二次采样：
根据采样频率及信号时域范围计算出采样间隔，根据采样间隔进行采样，就得到二次采样后的数据。
所用公式如下：
其中T是时域二次采样间隔，是二次采样频率，是时域范围，N是时域二次采样点数，是第一次采样点数，A是二次采样时数据数组元素间隔点数。
3.1.3仿真步骤、结果及分析
1．原始信号和插值后信号的时域频域分析
图1-2原始心电信号及插值后信号的波形和频谱
分析：
原始心电信号的频谱主要分布在0~100Hz,插值后相邻两点间的幅值之差减小了，相当于波形变平缓。
对应的频谱中可以看出，心电信号的低频分量幅度增加，高频分量幅度相对减少，能量重新分配。
2.相应的巴特沃斯低通、高通滤波器（其中以带通滤波器作为分析系统）
图1-3巴特沃斯模拟低通、高通滤波器的幅频特性
图1-4巴特沃斯带通滤波器的幅频特性
级联型系统的常数及分子、分母的系数：
b0=
3.6999e-004
B=
1.00002.00221.0022
1.00001.99780.9978
1.0000-1.99770.9977
1.0000-2.00231.0023
1.00002.00001.0000
1.0000-2.00001.0000
A=
1.0000-1.02800.2742
1.0000-1.13310.4123
1.0000-1.38550.7361
1.0000-1.97880.9790
1.0000-1.99910.9992
1.0000-1.98610.9861
分析：
采用级联型系统，因为级联型调整零极点方便，有利于控制频率响应，运算误差的累积相对直接型小。
该系统是一个6阶的因果稳定系统，可以实现带通滤波功能。
系统的信号流程图如下：
图1-5系统信号流图
图1-6系统的冲激响应和阶跃响应
分析：
在系统的时域特性中，冲激响应最终趋近于0，阶跃响应最终趋近于常数-0.38，说明系统是稳定的。
图1-7系统的零极点图
分析：
根据零极点图能定性分析系统函数的幅频响应。
在z=1处，部分零点与极点抵消，剩下的零点使系统函数的频响在ω=0处幅度为零，其它极点决定系统函数的频响的峰值位置，极点越靠近单位圆峰值越高越尖锐，在z=-1处，只有零点，使系统函数的频响在ω=π处幅度为零。
3．心电信号滤噪
图1-8通过组合滤波器和带通滤波器滤噪后的信号波形和频谱图
分析：
根据指导书要求设计的滤波器，最终有效的滤除了0.25Hz以下及99Hz以上的干扰，但靠近99Hz的部分有用信号也被滤除了。
自行设计的带通滤波器将99Hz大概设在3dB截止频率处，既滤除了高频干扰，又保留了有用的信号。
所以巴特沃斯数字带通滤波器的滤波效果较好。
4．滤除50Hz工频干扰
图1-9通过设置零极点设计的50Hz陷波器的幅频特性
图1-1050Hz陷波器的系统零极点图
图1-11巴特沃斯带阻滤波器的幅频特性
图1-12巴特沃斯带阻滤波器的零极点图
图1-13滤除50Hz工频干扰后的信号时域波形及频谱
分析：
采用设置零极点方法设计的陷波器只在60Hz处滤噪，阶数低，只有2阶，易实现，但是过渡带较宽，会影响其他频率分量的幅度。
巴特沃斯带阻滤波器设计为3阶，过渡带窄，避免了对其它频率分量的影响，但阻带有一定范围，会滤掉60Hz附近的部分频率分量。
滤噪后，心电信号的时域波形更平滑。
5．截取信号的3个周期
图1-14截取3个周期后的原始及插值之后的心电信号波形及频谱
图1-153个周期的心电信号滤噪后的波形及频谱
分析：
截取了心电信号中不正常的约3个周期的部分后，信号的时域范围缩小了，由可知，谱分辨率减小。
滤噪后，信号波形更平滑。
6．二次采样
图1-16二次采样后的信号波形及频谱
分析：
取第二次采样频率为200Hz，二次采样后谱分辨率不变，谱分析范围缩小。
3.1.4GUI设计及实现
图1-17GUI初始界面
点击“进入”，进入心电信号处理界面，如下：
图1-18心电信号分析系统的使用界面
导入数据后，点击“线性插值”、“截取三个周期”，选择巴特沃斯带通滤波器和50Hz陷波器分别进行滤波后，得到如下图形：
图1-19心电信号分析系统的GUI演示
点击“保存图片”，保存成功后，点击“退出使用”，就可以退出GUI界面。
保存的图片如下：
图1-20GUI界面中保存的心电信号图片
3.2基于simulink的简单心电信号分析系统设计
3.2.1仿真步骤及参数设置
1．数字带通滤波器和50Hz陷波器
图2-1建模方框图
图2-2数字带通滤波器参数设定及幅频特性
图2-350Hz陷波器参数设定及幅频特性
2．数字低通滤波器+数字高通滤波器+50Hz陷波器
图2-4数字低通+数字高通+50Hz陷波器建模方框图
图2-5数字低通滤波器参数设定及幅频特性
图2-6数字高通滤波器参数设定及幅频特性
3.Simulink线性插值模块
图2-7线性插值模块的建模方框图
3.2.2仿真结果及分析
图2-7例231.txt插值后的波形及频谱
1.数字带通滤波器和50Hz陷波器
图2-8滤波后的波形图频谱图
2．数字低通滤波器+数字高通滤波器+60Hz陷波器
图2-9滤波之后的波形图和频谱图
3．Simulink线性插值之后的波形图频谱图无明显变化。
分析：
经过多次修改参数，设计参数为：
ws1=0.25Hzwp1=1.7Hzwp2=90Hzws2=170HzRp=1Rs=30的数字带通滤波器滤波效果最为明显。
利用simulink进行动态仿真时我们设计并加上了50Hz工频陷波器，该陷波器是一个数字带阻滤波器设计参数为:
wp1=48Hzws1=49Hzws2=51Hzwp2=53HzRp=1Rs=30可以有效地滤除电源线干扰。
要求显示原始波形和频谱以及滤完波之后的波形和频谱，可以加上AveragingPowerSpectralDensity显示模块可同时显示波形和频谱。
3.3基于LabVIEW虚拟仪器的简单心电信号分析系统
3.3.1设计原理
1．心电信号的读取
本部分内容在以子vi的形式在主程序中调用。
图3-1读取心电信号.vi的程序流程图
图3-2读取心电信号.vi的程序框图
2．心电信号的线性插值处理
得到双精度型的时间和幅值序列之后，使用索引函数得到时间序列的第一个和最后一个数据，再使用斜坡信号函数得到间隔为0.001的时间序列，最后使用一维插值函数对幅度序列做线性插值。
图3-3心电信号的线性插值程序框图
3．设计相应的数字滤波器
利用枚举输入控件，while循环结构和条件结构实现对带通滤波器和带阻滤波器种类和参数的实时改变。
图3-4设计不同数字滤波器的程序框图
图3-5选择不同的数字滤波器的前面板
4.频谱分析
使用XY图显示控件对原始信号，插值信号显示。
使用波形图显示控件显示滤波前和滤波后的时域波形图和频谱图。
图3-6使用XY图显示控件显示原始信号和插值信号的程序框图
图3-7XY图显示原始信号和插值后信号前面板
图3-8使用波形图显示控件显示时域波形图和频谱图的前面板
5.使用选项卡控件设计简洁友好的人机交互界面
图3-9开始运行后的主界面
图3-10导入心电信号后的主界面
图3-11导入心电信号后的滤波前界面
图3-12导入心电信号的滤波后界面
3.3.2仿真调试步骤
1．导入心电信号文件
选择想要分析的心电信号文件，然后导入，并把插值后的幅值序列保存。
在导入文件后，程序运行，得到初始的时间幅值序列、插值后的时间幅值序列，插值前后的心电图以及滤波前后的波形图频谱图。
2．改变数字滤波器参数
对于带通滤波器和带阻滤波器，分别设计了Butterworth、Chebyshev、InverseChebyshev三种滤波器。
对于Butterworth滤波器可以更改滤波器阶数和3dB截止频率，对于Chebyshev和InverseChebyshev滤波器可以设置滤波器阶数、3dB截止频率和波纹衰减。
3.3.3分析
相同阶数下，Chebyshev和InverseChebyshev滤波器要比Butterworth滤波器的过渡带更窄，滤波效果更好；对于Chebyshev滤波器通带波纹衰减越小，滤波器性能越好；对于InverseChebyshev滤波器阻带波纹衰减越大，滤波器性能越好。
参考文献:
[1]王璨，章佳荣编著.LabVIEW2011程序设计与案例分析[M].北京：
北京航空航天大学出版社，2013
[2]罗华飞.MATLABGUI设计学习手记.北京航空航天大学出版社，2009.8
[3]丁玉美.数字信号处理（第二版）.西安电子科技大学出版社，2001
[4]吴大正.信号与线性系统分析（第四版）.高等教育出版社，2005.8
[5]丁亦农.Simulink与信号处理.北京航空航天大学出版社，2010.8
[6]王立会，潘东明.一种消除心电信号中工频干扰的陷波器设计.2007.7
附录：
1．心电信号的读取和插值程序（main.m文件）：
clear,clc
[t,xn]=duqu（'231.txt'）;%读取心电信号数据
subplot（2,2,1）,plot（t,xn）
xlabel（'t/s'）;ylabel（'xn'）;title（'原始心电信号波形'）
Fs=1000;
Xk=fft（xn,length（xn））;
n=0:
length（xn）-1;
f=Fs*n/length（xn）;
subplot（2,2,2）,plot（f,abs（Xk））
xlabel（'f/Hz'）;ylabel（'X（k）'）;title（'原始心电信号频谱'）
axis（[-10,150,0,1500]）
[t1,xn1]=chazhi（t,xn）;%对心电信号数据进行插值
baocun（t1,xn1）;
subplot（2,2,3）,plot（t1,xn1）
xlabel（'t/s'）;ylabel（'xn1'）;title（'插值后心电信号波形'）
Xk1=fft（xn1,length（xn1））;
n1=0:
length（xn1）-1;
f1=Fs*n1/length（xn1）;
subplot（2,2,4）,plot（f1,abs（Xk1））
xlabel（'f/Hz'）;ylabel（'X（k）1'）;title（'插值后心电信号频谱'）
axis（[-10,150,0,1500]）
%simulink数据输入端
simchazhi=[t1,xn1];
心电信号读取函数程序如下（duqu.m文件）：
function[t,Xn]=duqu（w）
fid=fopen（w）;%打开文档w
C=textscan（fid,'%8c%f%*f','headerlines',2）;%读出前两列数据且除去前两行
fclose（fid）;%关闭文件
a=C{2};
b=C{1};
k=length（b）;
fori=1:
k
c（i）=strread（b（i,:
）,'%*s%f','delimiter',':
'）;%以浮点数保留冒号之后的数据
end
c=c';
d=[c,a];
t=d（:
1）;
Xn=d（:
2）;
信号线性插值函数程序如下（chazhi.m文件）：
function[t2,Xn2]=chazhi1（t,Xn）
n=0;
y=0;
t=t.*1000;m=length（t）;
fori=1:
m
e（i）=round（t（i））;
end
fori=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;
forj=1:
N
z（y+j,1）=e（i）+（j-1）*1;
z（y+j,2）=Xn（i）+（j-1）*A
n=n+1;
j=j+1;
end
y=n;
end
i=i+1;
end
z（y+1,2）=Xn（i）;
z（y+1,1）=t（i）;
t2=z（:
1）;
t2=t2./1000;
Xn2=z（:
2）;
数据保存函数程序如下（baocun.m文件）：
functionbaocun（t,Xn）
fid=fopen（'t.txt','wt'）;
fprintf（fid,'%g\n',t）;%将t数组写入文档t.txt中
fclose（fid）;
fid=fopen（'Xn.txt','wt'）;
fprintf（fid,'%g\n',Xn）;%将t数组写入文档t.txt中
fclose（fid）;
2．滤波器的设计及使用
2.1巴特沃斯低通滤波器和高通滤波器（lowhigh.m文件）
%根据要求中给的设计指标，设计模拟巴特沃斯低通滤波器和高通滤波器
clear,clc
[t,xn]=duqu（'231.txt'）;%读取心电信号数据
[t1,xn1]=chazhi（t,xn）;%对心电信号数据进行插值
Fs=1000;
wp=92*2*pi;ws=99*2*pi;Rp=1;Rs=30;
[N,wc]=buttord（wp,ws,Rp,Rs,'s'）;
[b,a]=butter（N,wc,'s'）;
n=0:
length（xn1）-1;
f=n*Fs/length（xn1）;
w=f*2*pi;
h=freqs（b,a,w）;
subplot（2,1,1）,plot（f,abs（h））
xlabel（'f/Hz'）;ylabel（'H（z）'）;title（'巴特沃斯模拟低通滤波器'）
axis（[0,110,0,1.2]）
wph=0.7*2*pi;wsh=0.25*2*pi;Rp=1;Rs=30;
[Nh,wch]=buttord（wph,wsh,Rp,Rs,'s'）;
[bh,ah]=butter（Nh,wch,'high','s'）;
hh=freqs（bh,ah,w）;
subplot（2,1,2）,plot（f,abs（hh））
xlabel（'f/Hz'）;ylabel（'H（z）'）;title（'巴特沃斯模拟高通滤波器'）
axis（[0,10,0,1.2]）
2.2巴特沃斯带通滤波器（bandpass.m文件）
%巴特沃斯数字带通滤波器的分析
clear,clc
wp=[1.7*2*pi/1000,90*2*pi/1000];ws=[0.25*2*pi/1000,70*2*pi/1000];Rp=1;Rs=30;
[N,wc]=buttord（wp/pi,ws/pi,Rp,Rs）;
N
[b,a]=butter（N,wc）;
[h,w]=freqz（b,a）;
plot（w/max（w）,abs（h））;%幅频特性
xlabel（'w/pi'）;ylabel（'H（k）'）;title（'系统幅频特性'）;
axis（[-0.05,1,0,1.2]）;
%直接型转换成级联型
[b0,B,A]=dir2cas（b,a）
%分析系统的时域特性
y1=dimpulse（b,a,100）;%冲击响应
figure
subplot（2,1,1）,plot（y1）;
xlabel（'时间'）;ylabel（'幅值'）;title（'系统冲激响应'）;grid
axis（[0,100,-0.2,0.4]）;
y2=dstep（b,a,100）;%阶跃响应
subplot（2,1,2）,plot（y2）;
xlabel（'时间'）;ylabel（'幅值'）;title（'系统阶跃响应'）;grid
%分析系统零极点
figure
lingjdt（a,b）;
其中绘制零极点图函数程序如下（lingjdt.m文件）：
functionlingjdt（A,B）
p=roots（A）;%poles
q=roots（B）;%zeros
p=p';
q=q';
x=max（abs（[p,q]））;
x=x+0.1;
y=x;
clf
holdon
axis（[-xx-yy]）;
w=0:
pi/300:
2*pi;
t=exp（i*w）;
plot（t）;%画圆
axis（'square'）;
plot（[-xx],[00]）;
plot（[00],[-yy]）;
text（0.1,x,'jIM[z]'）;
text（y,1/10,'Re[z]'）;
plot（real（p）,imag（p）,'x'）;%poles
plot（real（q）,imag（q）,'o'）;%zeros
ti
升级会员 