数字滤波器设计及在心电信号滤波中的应用.docx
《数字滤波器设计及在心电信号滤波中的应用.docx》由会员分享，可在线阅读，更多相关《数字滤波器设计及在心电信号滤波中的应用.docx（22页珍藏版）》请在冰豆网上搜索。
数字滤波器设计及在心电信号滤波中的应用
课程设计报告
课程名称数字信号处理课程设计
课题名称数字滤波器设计及在心电信号滤波中的应用
专业通信工程
班级
学号
姓名
指导教师彭祯张鏖烽郭芳
2012年9月10日
湖南工程学院
课程设计任务书
课程名称数字信号处理课程设计
课题名称数字滤波器设计及在心电信号滤波中的应用
专业通信工程
班级
学号
姓名
指导教师彭祯张鏖烽郭芳
审批
任务书下达日期2012年9月1日
任务完成日期2012年9月10日
一、课程设计的性质与目的
《数字信号处理》课程是通信专业的一门重要专业基础课，是信息的数字化处理、存储和应用的基础。
通过该课程的课程设计实践，使学生对信号与信息的采集、处理、传输、显示、存储、分析和应用等有一个系统的掌握和理解；巩固和运用在《数字信号处理》课程中所学的理论知识和实验技能，掌握数字信号处理的基础理论和处理方法，提高分析和解决信号与信息处理相关问题的能力，为以后的工作和学习打下基础。
数字滤波器是一种用来过滤时间离散信号的数字系统，通过对抽样数据进行数学处理来达到频域滤波的目的。
根据其单位冲激响应函数的时域特性可分为两类：
无限冲激响应（IIR）滤波器和有限冲激响应（FIR）滤波器。
二、课程设计题目
数字滤波器设计及在心电信号滤波中的应用，其内容如下：
1、心电信号采集
心电信号作为心脏电活动在人体体表的表现，信号一般比较微弱，幅度在10μV~5mV，频率为0.05~100Hz。
在心电信号的采集、放大、检测及记录过程中，有来自外界的各种干扰。
记录一段时间内的人体心电信号波形，要求长度不小于10秒，并对记录的信号进行数字化，保存为数据文件；这里，请同学们使用美国的MIT/BIH心电原始数据，由实验老师给出一定长度的的心电原始数据，数据保存在文件“a01.txt～a10.txt”中，在MATLAB中通过如下语句读取：
%从当前路径下的a01.txt文件读取心电原始数据到变量a01中，a01为二维数据，第一列%为心电信号时间，第二列为心电信号幅度。
2、心电信号分析
使用MATLAB绘出数字化后的心电信号的时域波形和频谱图。
根据频谱图求出其带宽，并说明心电信号的基本特征。
3、含噪心电信号合成
在MATLAB软件平台下，给原始的心电信号叠加上噪声或干扰，干扰类型分为如下几种：
（1）白噪声；
（2）工频干扰（50Hz）；（3）谐波干扰（二次、三次谐波为主，分别为100Hz、150Hz）；（4）其它干扰，可设置为低频、高频、带限噪声，或冲激干扰。
绘出叠加噪声后的心电信号时域和频谱图，在视觉上与原始心电信号图形对比，绘出其时域波形差，分析频域基本特征变化。
4、数字滤波器设计及滤波
给定滤波器的规一化性能指标（参考指标,实际中依据每个同学所叠加噪声情况而定）例如：
通带截止频率wp=0.25*pi,阻通带截止频率ws=0.3*pi;通带最大衰减Rp=1dB;阻带最小衰减Rs=15dB，每个题目至少设计出5个用不同方法的不同类型滤波器。
采用双线性变换法与脉冲响应不变法，分别利用不同的原型低通滤波器（Butterworth型与切比雪夫I型）来设计各型IIR滤波器（低通、高通、带通、带阻中的至少3种类型），绘出滤波器的频域响应；并用这些数字滤波器对含噪心电信号分别进行滤波处理，比较不同方法下设计出来的数字滤波器的滤波效果，并从理论上进行分析（或解释）。
5、心电信号波形观察、频谱观察
对滤波后的心电信号观察其时域、频域特征变化。
绘出滤波后、滤波前、加噪后三个心电信号的差值波形，观察相互间的差异性；同时，分析频谱变化。
学生也可选用持续时间更长的心电原始数据，加上干扰后按上述要求设计滤波器。
三、课程设计要求
1、在一周内学生须上机16小时以上，程序调试完后，须由指导老师在机器上检查运行结果，经教师认可后的源程序可通过打印机输出，并请教师在程序清单上签字。
2、课程设计报告内容和格式：
设计题目，设计的详细步骤，设计过程中的结果、图形等，设计总结。
3、每组每人必须独立完成，成绩的考核按设计结果、答辩成绩及课程设计报告来综合评定。
成绩分为优、良、中、及格、不及格五级分评定。
4、指导教师：
彭祯，张鏖烽，郭芳。
四、设计进度安排
通信工程1001/1002:
1周周一上午，E-412,任务讲解与布置，学生分组选题，查找相关资料，准备课程设计，学生上机，按任务要求进行课程设计；分组选题；
1周周二上午，E-412,学生上机，按任务要求进行课程设计；分组选题；
1周周三上午，E-412,学生上机，按任务要求进行课程设计；分组选题；
1周周四上午，E-412,学生任务完成，答辩并提交课程设计报告。
附：
课程设计报告装订顺序：
封面、任务书、目录、正文、评分、附件（A4大小的图纸及程序清单）。
正文的格式:
一级标题用3号黑体,二级标题用四号宋体加粗,正文用小四号宋体;行距为22。
正文的内容:
一、课题的主要功能；二、课题的功能模块的划分；三、主要功能的实现；四、程序调试；五、总结；六、附件（所有程序的原代码，要求对程序写出必要的注释）；七、评分表。
目录
1课程设计的目的1
2课程设计的原理1
2.1Butterworth低通数字滤波器的设计1
2.2切比雪夫I型数字低通滤波器2
2.3IIR数字滤波器的性质..................................................................................................2
3课程设计设计步骤及结果分析4
3.1心电数据的导入4
3.2绘出心电信号的时域图和频谱图4
3.3加入噪声干扰5
3.3.1白噪声………………………………………………………………………………..5
3.3.2工频干扰（50Hz）………………………………………………………………….7
3.3.3带限chirp噪声……………………………………………………………………...8
3.4滤波器的设计9
3.4.1Butterwort型低通数字滤波器……………………………………………………...9
3.4.2切比雪夫I型数字低通滤波器………………………………………………….....12
3.4.3带阻滤波器的设计………………………………………………………………....15
4心得体会17
1课程设计的目的
数字滤波器是指输入,输出均为数字信号，通过数值运算处理改变输入信号所含频率成分的相对比例，或者滤除某些频率成分的数字器件或程序。
因此，数字滤波的概念和模拟滤波相同，只是信号的形式和实现滤波方法不同。
正因为数字滤波通过数值运算实现滤波，所以数字滤波器处理精度高，稳定，体积小，重量轻，灵活，不存在阻抗匹配问题，可以实现模拟滤波器无法实现的特殊滤波功能。
希望学生运用《数字信号处理》课程中所学的理论知识和实验技能，基本掌握数字信号处理的基础理论和处理方法，提高分析和解决信号与信息处理相关问题的能力，为以后的工作和学习打下基础。
2课程设计的原理
2.1Butterworth低通数字滤波器的设计
巴特沃斯低通滤波器的平方幅度响应为
其中，n为滤波器的阶数，
为低通滤波器的截止频率。
该滤波器具有一些特殊的性质：
①对所有的n,都有当
时,
；
②对所有的n,都有当
时,
；
③
是的单调递减函数,即不会出现幅度响应的起伏；
④当
时,巴特沃斯滤波器趋向于理想的低通滤波器；
⑤在
处平方幅度响应的各级导数均存在且等于0，因此
在该点上取得最大值，且具有最大平坦特性。
图1展示了2阶、4阶、8阶巴特沃斯低通滤波器的幅频特性。
可见阶数n越高，其幅频特性越好，低频检测信号保真度越高，过渡带变窄，即衰减加剧，但半功率点不变。
图1巴特沃斯低通滤波器的幅频特性
2.2切比雪夫I型数字低通滤波器
（1）确定数字低通滤波器的技术指标：
通带截止频率ωp、通带衰减p、阻带截止频率ωs、阻带衰减s
切比雪夫滤波器的振幅平方特性如图2所示：
图2切比雪夫滤波器的振幅平方特性
（2）将数字低通滤波器的技术指标转换成模拟低通滤波器的技术指标。
如果采用脉冲响不变法，边界频率的转换关系为:
如果采用双线性变换法，边界频率的转换关系为
（3）按照模拟低通滤波器的技术指标设计模拟低通滤波器。
（4）利用双线性变换法将模拟滤波器Ha（s），从s平面转换到z平面，得到数字低通滤波器系统函数H（z）。
（5）数字低通技术指标为：
ωp=0.4πrad,p=1dB;ωs=0.5πrad,s=40Db
（6）模拟低通的技术指标为：
归一化截止角频率wp=2pi*Fs/Ft;ws=2pi*Fs/Ft
（7）利用模拟切比雪夫滤波器设计数字滤波器。
通带截止频率为：
wp=0.4*pi;阻带截止频率为：
ws=0.5*pi;通带最大衰减为：
Rp=1;阻带最大衰减为：
As=15;设定周期为1s；模拟低通滤波器的生成：
[b,a]=cheby1（n,1,Wn,'low','s'）;满足设计指标的最小阶数和截止频率：
Wn[n,Wn]=cheb1ord（OmegaP,OmegaS,1,40,'s'）。
最后实现输入输出、幅频特性、
相频特性的图形。
2.3IIR数字滤波器的性质
无限长冲激响应（IIR）数字滤波器是数字滤器的一种，数字滤波器根据其冲激响应函数的时域特性,还包括有限长冲激响应（FIR）数字滤波器。
IIR数字滤波器的特征是：
①具有无限持续时间冲激响应；②需要用递归模型来实现，这可以从其差分方程
得出，也可以从其系数函数为：
得出。
数字巴特沃思滤波器属于IIR滤波器，该类滤波器具有特定的性质和设计方法。
目前比较成熟的IIR数字滤波器设计方法有两种：
1）直接法目前所用的方法主要是：
零极点累试法、频域幅度平方误差最小法和时域单位脉冲响应逼近法。
直接法的最大优点在于可以设计任意幅频特征的滤波器。
2）间接法，目前所用的方法主要是：
冲激响应不变法、阶跃响应不变法和双线性法。
它们都是借助于已经成熟的现有低通滤波器原型进行设计，即对数字低通数字滤波器，先将数字低通滤波器的技术指标按希望的设计方法转换为模拟低通滤波器的技术指标，再按指定的模拟低通滤波器的类型设计模拟滤波器H（s），然后，将模拟滤波器的系统函数H（s）从s平面转换到z平面，得到数字低通滤波器的系统函数H（z）；如所设计的数字滤波器为高通、带通或带阻滤波器，则可借助模拟滤波器的频带变换转换为低通模拟滤波器。
由于直接法设计巴特沃思滤波器相对复杂，在不需要任意幅频特征的情况下，一般采用问接法，同时由于冲激响应不变法和脉冲响应不变法，从s平面转换到z平面的映射为多值映射，容易造成频谱混叠，故而本文采用不会产生频谱混叠的双线性变换法。
3课程设计步骤及结果分析
3.1心电数据的导入
将老师给的心电信号原始数据存于桌面，然后我选用桌面上的第一组数据a=load（'C:
\DocumentsandSettings\hnie\桌面\心电信号数据\a01.txt'）;
3.2绘出心电信号的时域图和频谱图
将导入的两行数据分别用t,b来替换，然后通过调用plot函数来画出时域图，然后通过对1000个心电数据的幅值进行FFT运算，再次调用plot函数来绘出频域图，具体设计如下：
figure
（1）;
subplot（2,1,1）;
t=a（1:
1000,1）;b=a（1:
1000,2）;
plot（t,b）;
title（'原始波形图'）;xlabel（'时间（s）'）;ylabel（'幅值（A）'）;
y1=fft（a（:
2）,1000）;
f1=100*（0:
999）/1000;
subplot（2,1,2）;
plot（f1,abs（y1））;
title（'原始频谱图'）;xlabel（'频率（Hz）'）;ylabel（'幅度（dB）'）;图
图1.1原始时域和频谱图
通过导入的心电信号数据发现在其频谱图上的0~20Hz和80~100Hz之间的幅值比较大，而在30~70Hz之间的幅值相对较小。
3.3加入噪声干扰
这里我加入的噪声是：
白噪声，50Hz的工频噪声，带限chirp噪声。
3.3.1白噪声
通过用s来代表加入白噪声后的信号，并进行数字滤波器的频率响应，对s中的1000个频率点调用plot函数画出加入白噪声后的时域图，再对ws/pi,abshs调用plot函数画出加入白噪声后的频谱图，具体操作如下：
q=0.8*rand（1000,1）;
s=a（:
2）+q;
[hs,ws]=freqz（s,1,1024）;
abshs=abs（hs）;
figure
（2）;
subplot（2,1,1）;
plot（s（1:
1000））;
title（'加入白噪声后的时域图'）;xlabel（‘时间（s）'）;ylabel（'幅值（A）'）;
subplot（2,1,2）;
plot（ws/pi,abshs）;
title（'加入白噪声后的频谱图'）;xlabel（'幅度'）;ylabel（'Hz'）;
图1.2加入白噪声后的时域图和频谱图
通过观察加入白噪声后的时域图和频域图，将它与未加入白噪声进行比较，可以发现频谱图在0Hz时的幅度增加的很大，而且又在没有谱线的频率上竟然出现了频谱，这是由于白噪声在所有频率上都有频率造成的。
3.3.2工频干扰（50Hz）
用x表示加入工频干扰后的信号，再对工频信号的1000个频率点进行FFT运算，在进行相关的运算后，通过调用plot函数直接绘出加入工频干扰后的时域图和频谱图。
具体步骤如下：
x2=[sin（2*pi*50*t）];
t=0:
0.001:
0.001*（1000-1）;
x1=a（:
2）;
x=x1+x2;
y2=fft（x2,1000）;
f2=100*（0:
999）/1000;
figure（3）;subplot（2,1,1）;plot（t,x）;
title（'加入工频干扰后的时域图'）;xlabel（'时间（s）'）;ylabel（'幅值（A）'）;
subplot（2,1,2）;plot（f2,abs（y2））;
title（'加入工频干扰后的频谱图'）;xlabel（'幅度'）;ylabel（'Hz'）;
图1.3加入工频干扰后的时域图和频谱图
比较图1.1和图1.3，发现加噪声后的幅值有比较微小的变化，而在频谱图在38Hz和63Hz附近增值很大，这是由于50Hz的工频噪声造成的，而在其他频率范围内也有比较明显的变化。
3.3.3带限chirp噪声
同样用k表示加入带限chirp噪声后的信号，再用freqz进行滤波器的频率响应，调用plot对k画出加入带限噪声后的时域和频域图，具体描述如下：
p=0.5*chirp（a（:
1）,0,a（1000,1）,200）;
k=x1+p;
[hc,wc]=freqz（k,1,1024）;
abshc=abs（hc）;
figure（4）;
subplot（2,1,1）;plot（k（1:
1000））;
title（'加入带限噪声后的时域图'）;xlabel（'时间（s）'）;ylabel（'幅值（A）'）;
subplot（2,1,2）;
plot（wc/pi,abshc）;
title（'加入带限噪声后的频谱图'）;xlabel（'幅度'）;ylabel（'Hz'）;
图1.4加入带限噪声后的时域图和频谱图
加入带限chirp噪声干扰后，发现加噪之前时域图是均匀分布的，而加噪后，则变为前面稀疏后面密集的情况了，对于频谱图而言，加噪前是中间凹两边凸，加噪后是前凸后平，在个别点上幅值增加很多，这是由于加入chirp噪声的结果。
3.4滤波器的设计
3.4.1Butterworth型低通数字滤波器
用wp和ws表示分别将通带，阻带截止频率的角频率表示，在分别计算阶数n1和截止频率Wn,再设计低通Butterworth型模拟滤波器，然后采用双线性法将模拟滤波器系数变为数字滤波器系数，画出滤波器频谱图，调用filter实现对工频干扰的滤波，用plot函数画出滤除工频干扰后的时域图和滤除白噪声后的频谱图。
具体操作如下：
figure（5）;
fs=100;
f1=5;f2=10;
wp=（f1/fs）*2*pi;
ws=（f2/fs）*2*pi;
Omegap=2*fs*tan（wp/2）;
Omegas=2*fs*tan（ws/2）;
[n1,Wn]=buttord（Omegap,Omegas,1,50,'s'）;
[b,a]=butter（n1,Wn,'s'）;
[bz,az]=bilinear（b,a,fs）;
freqz（bz,az,512,fs）;
y=filter（bz,az,x2）;
figure（6）;
subplot（2,1,1）
plot（y）;
title（'滤除工频干扰后的时域图'）;
xlabel（'时间（s）'）;ylabel（'幅值（A）'）;
y3=fft（y,1000）;
f1=100*（0:
999）/1000;
subplot（2,1,2）;
plot（f1,abs（y3））;
title（'滤除白噪声后的频谱图'）;xlabel（'频率（Hz）'）;ylabel（'幅值（dB）'）;
图1.5Butterworth型低通数字滤波器频率响应图
图1.6滤除单频正弦波后的时域图和滤除噪声后的频谱图
由于信号处于频段的低频部分，而工频信号的频谱在整个上是呈对称分布的，与源信号的频段分布是相似的，采用低通滤波器将噪声信号的高频部分滤掉，由图1.3和1.6会发现加入工频干扰后的源信号，滤波后时域图的幅值比原来变小了，滤波后的频谱图频率在0~10Hz和90~100Hz之间的幅值变化较大，而在10~90Hz的区间，频谱图基本趋于直线。
3.4.2切比雪夫I型数字低通滤波器
用Wp1,Wp2,Ws1,Ws2表示分别将通带，阻带截止频率的角频率表示,算出频带宽带，计算阶数n1和截止频率WN，再设计切比雪夫I型模拟滤波器，采用双线性法将模拟滤波器系数变为数字滤波器系数，画出切比雪夫I型数字滤波器的频率响应，调用filter实现对白噪声的滤波，再最后调用plot函数画出滤除白噪声后的时域图和频域图。
具体过程如下：
figure（7）;
fs=100;
f11=10;f12=25;f21=5;f22=30;
Wp1=（f11/fs）*2*pi;
Ws1=（f21/fs）*2*pi;
Wp2=（f12/fs）*2*pi;
Ws2=（f22/fs）*2*pi;
Omegap1=2*fs*tan（Wp1/2）;Omegap2=2*fs*tan（Wp2/2）;
Omegas1=2*fs*tan（Ws1/2）;Omegas2=2*fs*tan（Wp2/2）;
BW=Omegap2-Omegap1;
W0=Omegap1*Omegap2;W00=sqrt（W0）;
WP=1;WS=WP*（W0^2-Ws1^2）/（Ws1*BW）;
[n1,WN]=buttord（WP,WS,1,50,'s'）;
[B,A]=cheby1（n1,1,WN,'s'）;
[BT,AT]=lp2bp（B,A,W00,BW）;
[num,den]=bilinear（BT,AT,0.5）;
freqz（num,den,64）;
y=filter（num,den,s）;
figure（8）;
subplot（2,1,1）;
plot（y）;
title（'滤除白噪声后的时域图'）;
xlabel（'时间（s）'）;ylabel（'幅值（A）'）;
s3=fft（y,1000）;
f1=100*（0:
999）/1000;
subplot（2,1,2）;
plot（f1,abs（s3））;
title（'滤除白噪声后的频谱图'）;xlabel（'频率（Hz）'）;ylabel（'幅值（dB）'）;
图1.7切比雪夫1型数字滤波器的频率响应图
图1.8滤除白噪声后的时域图和频域图
由于白噪声在整个频段上都存在，对与源信号共存的低频信号用选频滤波器是无法滤除的，故而采用低通滤波器将噪声的高频部分去掉。
比较图1.2和图1.7，它们在时域图上差不多是一致的，但在频谱图中的图1.2在0Hz上的幅值是很搞的，在滤波后的就明显变小了。
3.4.3带阻滤波器的设计
设计带阻滤波器，再画出滤波器的频率响应图，然后用filter实现对工频干扰信号的滤波，调用plot函数分别画出滤除工频干扰后的时域图和频谱图。
具体操作如下：
figure（9）;
d=fir1（1000,[0.1570.17],'stop'）;
freqz（d,512）;
y4=filter（d,1,x2）;
figure（10）;
subplot（2,1,1）;plot（y4）;
title（'滤除工频干扰后的时域图'）;
xlabel（'时间（s）'）;ylabel（'幅值（A）'）;
subplot（2,1,2）;
y5=fft（y4,1000）;
f1=100*（0:
999）/1000;
plot（f1,abs（y5））;
title（'滤除工频干扰后的频谱图'）;xlabel（'频率（Hz）'）;ylabel（'幅值（dB）'）;
图1.9带阻滤波器的频率响应图
图1.10滤除工频干扰后的时域图和频谱图
带阻滤波器滤除工频干扰在28Hz和74Hz的地方最为显著，带阻滤波器由于工频干扰信号是对称的，而且与源信号的高频部分在相同的频段上，因此加噪后高频部分的频段上的幅值会更高，故我们选择带阻滤波器来滤除噪声干扰的作用。
4心得体会
通过本次课程设计，使我基本懂得了《数字信号处理》的基础理论知识和设计方法，同时在课设时遇到了不少的问题，通过查阅相关书籍或请教老师同学获得了不少的帮助，在此期间也学到了不少的知识，提升了知识层面的容量，扩宽了设计的道路。
由于在课设中也遇到了关于MATLAB的相关问题，但也比较轻松的解决了，同时也学习到了在实践时的不同知识，真的让我获益良多。
课程设计中需要自我的耐心和解决问题的勇气，正好让我体验了一把，没有足够的耐心只会感觉到它的枯燥乏味，而没有体会到实践带来的动手能力的提高，这是在以后工作中不可缺少的能力之一，在此感谢辛勤的老师和帮助过我的同学，希望我们一起提高和进步。
评分表
课题名称：
项目
评价
设计方案的合理性与创造性
设计与调试结果
设计说明书的质量
答辩陈述与回答问题情况
课程设计周表现情况
综合成绩
教师签名：
日期：
升级会员 