数字信号处理数字滤波器设计及在心电信号滤波中的应用的课程设计.docx

数字信号处理数字滤波器设计及在心电信号滤波中的应用的课程设计.docx

《数字信号处理数字滤波器设计及在心电信号滤波中的应用的课程设计.docx》由会员分享，可在线阅读，更多相关《数字信号处理数字滤波器设计及在心电信号滤波中的应用的课程设计.docx（28页珍藏版）》请在冰豆网上搜索。

数字信号处理数字滤波器设计及在心电信号滤波中的应用的课程设计

课程设计报告

课程名称数字信号处理课程设计

课题名称数字滤波器设计及在心电信号滤波中的应用

专业通信工程

班级1081

学号201013120103

姓名刘献文

指导教师彭祯谭小兰

2012年9月10日

湖南工程学院

课程设计任务书

课程名称数字信号处理课程设计

课题名称数字滤波器设计及在心电信号滤波中的应用

专业通信工程

班级1081

学号201013120103

姓名刘献文

指导教师彭祯谭小兰

审批

任务书下达日期2012年9月1日

任务完成日期2012年9月10日

《数字信号处理》课程设计指导

一、课程设计的性质与目的

《数字信号处理》课程是通信专业的一门重要专业基础课，是信息的数字化处理、存储和应用的基础。

通过该课程的课程设计实践，使学生对信号与信息的采集、处理、传输、显示、存储、分析和应用等有一个系统的掌握和理解；巩固和运用在《数字信号处理》课程中所学的理论知识和实验技能，掌握数字信号处理的基础理论和处理方法，提高分析和解决信号与信息处理相关问题的能力，为以后的工作和学习打下基础。

数字滤波器是一种用来过滤时间离散信号的数字系统，通过对抽样数据进行数学处理来达到频域滤波的目的。

根据其单位冲激响应函数的时域特性可分为两类：

无限冲激响应（IIR）滤波器和有限冲激响应（FIR）滤波器。

二、课程设计题目

课程设计题目分为2个方向，每班选做方向一的同学约占2/3,选做方向二的约占1/3。

方向一：

数字滤波器设计及在语音信号分析中的应用。

方向二：

数字滤波器设计及在心电信号滤波中的应用。

步骤：

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个用不同方法的不同类型滤波器。

题目

（1）：

采用窗函数法与等波纹法分别设计各型FIR滤波器（低通、高通、带通、带阻中的至少3种类型）来对叠加干扰前后的心电信号进行滤波处理，绘出滤波器的频域响应，绘出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；在相同的性能指标下比较各方法的滤波效果，并从理论上进行分析（或解释）。

题目

（2）：

采用双线性变换法与脉冲响应不变法，分别利用不同的原型低通滤波器（Butterworth型与切比雪夫I型）来设计各型IIR滤波器（低通、高通、带通、带阻中的至少3种类型），绘出滤波器的频域响应；并用这些数字滤波器对含噪心电信号分别进行滤波处理，比较不同方法下设计出来的数字滤波器的滤波效果，并从理论上进行分析（或解释）。

5、心电信号波形观察、频谱观察

对滤波后的心电信号观察其时域、频域特征变化。

绘出滤波后、滤波前、加噪后三个心电信号的差值波形，观察相互间的差异性；同时，分析频谱变化。

学生也可选用持续时间更长的心电原始数据，加上干扰后按上述要求设计滤波器。

三、课程设计要求

1、在一周内学生须上机16小时以上，程序调试完后，须由指导老师在机器上检查运行结果，经教师认可后的源程序可通过打印机输出，并请教师在程序清单上签字。

2、课程设计报告内容和格式：

设计题目，设计的详细步骤，设计过程中的结果、图形等，设计总结。

3、每组每人必须独立完成，成绩的考核按设计结果、答辩成绩及课程设计报告来综合评定。

成绩分为优、良、中、及格、不及格五级分评定。

4、指导教师：

彭祯，谭小兰。

四、设计进度安排

通信工程1081:

1周周一下午，E-412,任务讲解与布置，学生分组选题，查找相关资料，准备课程设计，学生上机，按任务要求进行课程设计；分组选题；

1周周二下午，E-412,学生上机，按任务要求进行课程设计；分组选题；

1周周三下午，E-412,学生上机，按任务要求进行课程设计；分组选题；

1周周四下午，E-412,学生任务完成，答辩并提交课程设计报告。

附：

课程设计报告装订顺序：

封面、任务书、目录、正文、评分、附件（A4大小的图纸及程序清单）。

正文的格式:

一级标题用3号黑体,二级标题用四号宋体加粗,正文用小四号宋体;行距为22。

正文的内容:

一、课题的主要功能；二、课题的功能模块的划分；三、主要功能的实现；四、程序调试；五、总结；六、附件（所有程序的原代码，要求对程序写出必要的注释）；七、评分表。

目录

1、课程设计的目的1

2、课程设计的原理1

2.1用窗函数法设计FIR滤波器1

2.2用巴特沃斯法设计IIR低通滤波器1

3、课程设计设计步骤及结果分析1

3.1心电信号采集1

3.2源数据的导入2

3.3绘出源心电信号的时域波形图和频谱图并进行分析2

3.4含噪心电信号合成3

3.5时域波形差5

3.6数字滤波器设计及滤波6

3.7心电信号时域波形及频谱的观察和分析12

4、心得体会18

5、附件：

（源代码）18

6、评分表22

1课程设计的目的

通过对常用数字滤波器的设计和实现，掌握数字信号处理的工作原理及设计方法；熟悉设计数字滤波器的方法，掌握利用数字滤波器对信号进行滤波的方法，掌握数字滤波器的计算机仿真，并能够对设计结果加以分析。

2课程设计的原理

2.1、用窗函数法设计FIR滤波器

根据过渡宽带及阻带衰减要求，选择窗函数的类型并估计窗口长度N，窗函数类型可根据最小阻带衰减As独立选择，因为窗口长度N对最小阻带衰减As没有影响，在确定窗函数类型以后，可根据过渡宽带小于给定指标确定所拟用的窗函数的窗口长度N，设待求滤波器的过渡宽带，它与窗口长度N近似成反比，窗函数类型确定后，其计算公式也确定了，不过这些公式是近似的，得出的窗口长度还要在计算中逐步修正，原则上在保证阻带衰减满足要求的情况下，尽量选择较小的N，在N和窗函数类型确定后，即可调用MATLAB中的窗函数求出窗函数。

2.2、用巴特沃斯法设计IIR低通滤波器

巴特沃斯低通滤波器的幅度平方函数中，N为滤波器的阶数。

在Ω=Ωc附近，随着Ω加大，幅度迅速下降。

幅度下降的速度与阶数N有关，N越大，通带愈平坦，过渡带愈窄，过渡带与阻带幅度下降的速度愈快，总的频响特性与理想低通滤波器的误差愈小。

3课程设计设计步骤及结果分析

3.1、心电信号采集

心电信号作为心脏电活动在人体体表的表现，信号一般比较微弱，幅度在10μV~5mV，频率为0.05~100Hz。

在心电信号的采集、放大、检测及记录过程中，有来自外界的各种干扰。

记录一段时间内的人体心电信号波形，要求长度不小于10秒，并对记录的信号进行数字化，保存为数据文件；这里，请同学们使用美国的MIT/BIH心电原始数据，由实验老师给出一定长度的的心电原始数据，数据保存在文件“a01.txt～a10.txt”中，在MATLAB中通过如下语句读取：

%从当前路径下的a01.txt文件读取心电原始数据到变量a01中，a01为二维数据，第一列%为心电信号时间，第二列为心电信号幅度。

3.2、源数据的导入

a=load（'F:

/心电信号数据/a16.txt'）;%选第十六组数据

3.3、绘出源心电信号的时域波形图和频谱图并进行分析

%时域波形图以及频谱图

a=load（'F:

/心电信号数据/a16.txt'）;

t=a（:

1）;b=a（:

2）;

figure

（1）;

subplot（211）;

plot（t,b）;

title（'心电信号的时域波形'）;

xlabel（'时间t/s'）;ylabel（'幅值/A'）;

n=1000;

m=abs（fft（b,n））;

fs=100;

f=fs/n*（0:

n-1）;

subplot（212）;

plot（f,m）;

title（'心电信号的频谱图'）;

xlabel（'频率f/Hz'）;ylabel（'幅值/db'）;

3.4、含噪心电信号合成

1、加入单频正弦干扰

y=0.2*sin（2*40*pi*t）;

y1=y+b;

figure

（2）;

subplot（211）;

plot（t,y1）;

title（'加单频噪声时域波形'）;

xlabel（'时间t/s'）;ylabel（'幅值/A'）;

k=abs（fft（y1,n））;

subplot（212）;

plot（f,k）;

title（'加单频噪声频谱图'）;

xlabel（'频率/Hz'）;ylabel（'幅值'）

选择添加单频正弦信号。

通过转置后添加到原始信号b上，添加后的信号用y1表示。

对于变量、函数的操作是矩阵操作，两个信号相加必须长度一致。

2、继续加入白噪声

y2=awgn（y1,5）;

figure（3）;

subplot（211）;

plot（t,y2）;

title（'加单频,白噪声时域波形'）;

xlabel（'时间t/s'）;ylabel（'幅值/A'）;

h=abs（fft（y2,n））;

subplot（212）;

plot（f,h）;

title（'加单频,白噪声频谱图'）;

xlabel（'频率/Hz'）;ylabel（'幅值'）;

通过调用函数y2=awgn（y1,5）对已经加了单频噪声的信号进行加噪，加噪后的信号用y2保存。

对y2信号的分析，如图所示分别给出了它的时域波形图和频谱图，通过对比，时域波形变化不大，但频谱图有了大幅度的变化。

3.5、时域波形差

s=y2-b;

figure（4）;

subplot（111）;

plot（t,s）;

title（'时域波形差'）;

3.6、数字滤波器设计及滤波

1、设计低通FIR滤波器

%低通FIR滤波器

N=50;

wc=0.6;

window=blackman（N）;

hn=fir1（N-1,wc,window）;

lv=filter（hn,1,y2）;

freqz（hn,1）;

2、用低通FIR滤波器对单频噪声进行滤波

%低通滤波器滤波后图形

subplot（211）;

plot（t,lv）;

title（'滤波后时域波形'）;

xlabel（'时间t/s'）;ylabel（'幅值/A'）;

LV=abs（fft（lv））;

subplot（212）;

plot（f,LV）;

title（'滤波后频谱图'）;

xlabel（'频率/Hz'）;ylabel（'幅值'）

3、设计带阻FIR滤波器

%带阻FIR滤波器

wlp=0.2*pi;wls=0.4*pi;wus=1.2*pi;wup=0.6*pi;

B=wls-wlp;

M=ceil（12*pi/B）;

wp=[（wls+wlp）/2/pi,（wus+wup）/2/pi];

hn=fir1（M,wp,'stop',blackman（M+1））;

liu=filter（hn,1,y2）;

freqz（hn,1）;

4、用带阻FIR滤波器对单频噪声滤波后图形

%带阻FIR滤波器滤波后图形

subplot（211）;

plot（t,liu）;

title（'滤波后时域波形'）;

xlabel（'时间t/s'）;ylabel（'幅值/A'）;

LIU=abs（fft（liu））;

subplot（212）;

plot（f,LIU）;

title（'滤波后频谱图'）;

xlabel（'频率/Hz'）;ylabel（'幅值'）;

5、设计低通IIR滤波器

%低通IIR滤波器

wp=0.4*pi;ws=0.35*pi;rp=1;rs=10;

[N,wc]=buttord（wp,ws,rp,rs,'s'）;

[B,A]=butter（N,wc,'s'）;

[Bz,Az]=impinvar（B,A）;

wen=filter（Bz,Az,y2）;

freqz（Bz,Az）;

6、用低通IIR滤波器对单频噪声滤波后图形

%低通IIR滤波器滤波后图形

subplot（211）;

plot（t,wen）;

title（'滤波后时域波形'）;

xlabel（'时间t/s'）;ylabel（'幅值/A'）;

WEN=abs（fft（wen））;

subplot（212）;

plot（f,WEN）;

title（'滤波后频谱图'）;

xlabel（'频率/Hz'）;ylabel（'幅值'）;

3.7、心电信号时域波形及频谱的观察和分析

1、低通FIR滤波器滤波效果

subplot（331）;

plot（t,b）;

title（'心电信号的时域波形'）;

xlabel（'时间t/s'）;ylabel（'幅值/A'）;

subplot（332）;

plot（f,m）;

title（'心电信号的频谱图'）;

xlabel（'频率f/Hz'）;ylabel（'幅值/db'）;

subplot（333）;

plot（t,y1）;

title（'加单频噪声时域波形'）;

xlabel（'时间t/s'）;ylabel（'幅值/A'）

subplot（334）;

plot（f,k）;

title（'加单频噪声频谱图'）;

xlabel（'频率/Hz'）;ylabel（'幅值'）;

subplot（335）;

plot（t,y2）;

title（'加单频,白噪声时域波形'）;

xlabel（'时间t/s'）;ylabel（'幅值/A'）;

subplot（336）;

plot（f,h）;

title（'加单频,白噪声频谱图'）;

xlabel（'频率/Hz'）;ylabel（'幅值'）;

subplot（337）;

plot（t,lv）;

title（'滤单频波后时域波形'）;

xlabel（'时间t/s'）;ylabel（'幅值/A'）;

subplot（338）;

plot（f,LV）;

title（'滤单频波后频谱图'）;

xlabel（'频率/Hz'）;ylabel（'幅值'）;

2、带阻FIR滤波器滤波效果

subplot（331）;

plot（t,b）;

title（'心电信号的时域波形'）;

xlabel（'时间t/s'）;ylabel（'幅值/A'）;

subplot（332）;

plot（f,m）;

title（'心电信号的频谱图'）;

xlabel（'频率f/Hz'）;ylabel（'幅值/db'）;

subplot（333）;

plot（t,y1）;

title（'加单频噪声时域波形'）;

xlabel（'时间t/s'）;ylabel（'幅值/A'）

subplot（334）;

plot（f,k）;

title（'加单频噪声频谱图'）;

xlabel（'频率/Hz'）;ylabel（'幅值'）;

subplot（335）;

plot（t,y2）;

title（'加单频,白噪声时域波形'）;

xlabel（'时间t/s'）;ylabel（'幅值/A'）;

subplot（336）;

plot（f,h）;

title（'加单频,白噪声频谱图'）;

xlabel（'频率/Hz'）;ylabel（'幅值'）;

subplot（337）;

plot（t,liu）;

title（'滤单频波后时域波形'）;

xlabel（'时间t/s'）;ylabel（'幅值/A'）;

subplot（338）;

plot（f,LIU）;

title（'滤单频波后频谱图'）;

xlabel（'频率/Hz'）;ylabel（'幅值'）;

3、低通IIR滤波器滤波效果

subplot（331）;

plot（t,b）;

title（'心电信号的时域波形'）;

xlabel（'时间t/s'）;ylabel（'幅值/A'）;

subplot（332）;

plot（f,m）;

title（'心电信号的频谱图'）;

xlabel（'频率f/Hz'）;ylabel（'幅值/db'）;

subplot（333）;

plot（t,y1）;

title（'加单频噪声时域波形'）;

xlabel（'时间t/s'）;ylabel（'幅值/A'）

subplot（334）;

plot（f,k）;

title（'加单频噪声频谱图'）;

xlabel（'频率/Hz'）;ylabel（'幅值'）;

subplot（335）;

plot（t,y2）;

title（'加单频,白噪声时域波形'）;

xlabel（'时间t/s'）;ylabel（'幅值/A'）;

subplot（336）;

plot（f,h）;

title（'加单频,白噪声频谱图'）;

xlabel（'频率/Hz'）;ylabel（'幅值'）;

subplot（337）;

plot（t,wen）;

title（'滤单频波后时域波形'）;

xlabel（'时间t/s'）;ylabel（'幅值/A'）;

subplot（338）;

plot（f,WEN）;

title（'滤单频波后频谱图'）;

xlabel（'频率/Hz'）;ylabel（'幅值'）;

4心得体会

通过这次的课程设计，我从中得到了许多经验和MATLAB软件设计的一些新思路；在学习MATLAB的过程中，我明白了“实践出真知”这句话的真谛，在书上看到的内容，如果不懂的话，我们可以把程序输入代码窗口，运行之后，看得出的结果，然后再加上书本上的讲解，就很好理解了，相反，如果我们只是想着这个程序或者这个函数有什么功能，肯定想很长时间都不明白。

通过这次课程设计，MATLAB的使用不仅仅只是看书就能够学会的，要自己在计算机上面动手操作，才能熟练的使用MATLAB软件，我还发现了MATLAB一个强大的功能是它为我们提供一套功能强大的绘图命令，这些命令可以根据输入的数据自动完成图形的绘制，可以完成对图形的加标号，加标题等操作。

总结一下这次课程设计，发现自己虽然在不仅在理论上没有掌握牢固，并且在实践的时候也遇到了问题，所以自己还是远远的不足，不管是在MATLAB的设计上，还是其他专业课上，在以后的一段学习时间里必须坚持自己思考，自己多动脑，多动手，这样才能脱离理论，让自己的学习更上一层楼。

5附件：

（源代码）

a=load（'F:

/心电信号数据/a16.txt'）;

t=a（:

1）;b=a（:

2）;

figure

（1）;

subplot（211）;

plot（t,b）;

title（'心电信号的时域波形'）;

xlabel（'时间t/s'）;ylabel（'幅值/A'）;

n=1000;

m=abs（fft（b,n））;

fs=100;

f=fs/n*（0:

n-1）;

subplot（212）;

plot（f,m）;

title（'心电信号的频谱图'）;

xlabel（'频率f/Hz'）;ylabel（'幅值/db'）;

y=0.2*sin（2*40*pi*t）;

y1=y+b;

figure

（2）;

subplot（211）;

plot（t,y1）;

title（'加单频噪声时域波形'）;

xlabel（'时间t/s'）;ylabel（'幅值/A'）;

k=abs（fft（y1,n））;

subplot（212）;

plot（f,k）;

title（'加单频噪声频谱图'）;

xlabel（'频率/Hz'）;ylabel（'幅值'）;

y2=awgn（b,5）;

figure（3）;

subplot（211）;

plot（t,y2）;

title（'加单频,白噪声时域波形'）;

xlabel（'时间t/s'）;ylabel（'幅值/A'）;

h=abs（fft（y2,n））;

subplot（212）;

plot（f,h）;

title（'加单频,白噪声频谱图'）;

xlabel（'频率/Hz'）;ylabel（'幅值'）;

s=y2-b;

figure（4）;

subplot（111）;

plot（t,s）;

title（'时域波形差'）;

%低通FIR滤波器

N=50;

wc=0.6;

window=blackman（N）;

hn=fir1（N-1,wc,window）;

lv=filter（hn,1,y2）;

freqz（hn,1）;

%低通滤波器滤波后图形

subplot（211）;

plot（t,lv）;

title（'滤波后时域波形'）;

xlabel（'时间t/s'）;ylabel（'幅值/A'）;

LV=abs（fft（lv））;

subplot（212）;

plot（f,LV）;

title（'滤波后频谱图'）;

xlabel（'频率/Hz'）;ylabel（'幅值'）;

%带阻FIR滤波器

wlp=0.2*pi;wls=0.4*pi;wus=1.2*pi;wup=0.6*pi;

B=wls-wlp;

M=ceil（12*pi/B）;

wp=[（wls+wlp）/2/pi,（wus+wup）/2/pi];

hn=fir1（M,wp,'stop',blackman（M+1））;

liu=filter（hn,1,y2）;

freqz（hn,1）;

%带阻FIR滤波器滤波后图形

subplot（211）;

plot（t,liu）;

title（'滤波后时域波形'）;

xlabel（'时间t/s'）;ylabel（'幅值/A'）;

LIU=abs（fft（liu））;

subplot（212）;

plot（f,LIU）;

title（'滤波后频谱图'）;

xlabel（'频