数字信号处理课程设计.docx
《数字信号处理课程设计.docx》由会员分享，可在线阅读，更多相关《数字信号处理课程设计.docx（23页珍藏版）》请在冰豆网上搜索。
数字信号处理课程设计
课程设计报告
课程名称数字信号处理课程设计
课题名称数字滤波器设计及在语音信号分析中的应用
专业通信工程
班级
学号
姓名
指导教师张细政
2011年7月4日
湖南工程学院
课程设计任务书
课程名称数字信号处理
课题数字滤波器设计及语
音信号分析中的应用
专业班级通信工程0981
学生姓名
学号
指导老师张细政
审批
任务书下达日期2011年6月21日
任务完成日期2011年7月4日
目录
一、课程设计的性质与目的4
二、课程设计题目4
方向一：
数字滤波器设计及在语音信号分析中的应用4
方向二：
数字滤波器设计及在心电信号滤波中的应用5
三、课程设计要求6
四、设计进度安排7
五、课题的主要功能8
5.1问题的描述8
5.2功能要求8
六、设计原理8
6.1．用窗函数法设计FIR滤波器9
6.2.用双线性变换法设计IIR数字滤波器10
6.3.用脉冲响应不变法设计IIR数字滤波器12
七、课题设计13
7.1、信号的采集15
7.2、信号分析15
7.3、含噪信号合成15
7.4、数字滤波器15
八、设计体会22
九、参考文献23
课程设计评分表24
一、课程设计的性质与目的
《数字信号处理》课程是通信专业的一门重要专业基础课，是信息的数字化处理、存储和应用的基础。
通过该课程的课程设计实践，使学生对信号与信息的采集、处理、传输、显示、存储、分析和应用等有一个系统的掌握和理解；巩固和运用在《数字信号处理》课程中所学的理论知识和实验技能，掌握数字信号处理的基础理论和处理方法，提高分析和解决信号与信息处理相关问题的能力，为以后的工作和学习打下基础。
数字滤波器是一种用来过滤时间离散信号的数字系统，通过对抽样数据进行数学处理来达到频域滤波的目的。
根据其单位冲激响应函数的时域特性可分为两类：
无限冲激响应（IIR）滤波器和有限冲激响应（FIR）滤波器。
二、课程设计题目
课程设计题目分为2个方向，选做方向一的同学约占2/3,选做方向二的约占1/3。
方向一：
数字滤波器设计及在语音信号分析中的应用。
步骤：
1、语音信号采集
录制一段课程设计学生的语音信号并保存为文件，要求长度不小于10秒，并对录制的信号进行采样；录制时可以使用Windows自带的录音机，或者使用其它专业的录音软件，录制时需要配备录音硬件（如麦克风），为便于比较，需要在安静、干扰小的环境下录音。
2、语音信号分析
使用MATLAB绘出采样后的语音信号的时域波形和频谱图。
根据频谱图求出其带宽，并说明语音信号的采样频率不能低于多少赫兹。
3、含噪语音信号合成
在MATLAB软件平台下，给原始的语音信号叠加上噪声，噪声类型分为如下几种：
（1）白噪声；
（2）单频噪色（正弦干扰）；（3）多频噪声（多正弦干扰）；（4）其它干扰，可设置为低频、高频、带限噪声，或Chirp干扰、冲激干扰。
绘出叠加噪声后的语音信号时域和频谱图，在视觉上与原始语音信号图形对比，也可通过Windows播放软件从听觉上进行对比，分析并体会含噪语音信号频谱和时域波形的改变。
4、数字滤波器设计及滤波,完成以下题目中的一个
给定滤波器的规一化性能指标（参考指标,实际中依据每个同学所叠加噪声情况而定）例如：
通带截止频率wp=0.25*pi,阻通带截止频率ws=0.3*pi;通带最大衰减Rp=1dB;阻带最小衰减Rs=15dB，每个题目至少设计出5个用不同方法的不同类型滤波器。
题目
（1）：
采用窗函数法与等波纹法分别设计各型FIR滤波器（低通、高通、带通、带阻中的至少3种类型）来对叠加噪声前后的语音信号进行滤波处理，绘出滤波器的频域响应，绘出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；在相同的性能指标下比较各方法的滤波效果，并从理论上进行分析（或解释）。
题目
（2）：
采用双线性变换法与脉冲响应不变法，分别利用不同的原型低通滤波器（Butterworth型与切比雪夫I型）来设计各型IIR滤波器（低通、高通、带通、带阻中的至少3种类型），绘出滤波器的频域响应；并用这些数字滤波器对含噪语音信号分别进行滤波处理，比较不同方法下设计出来的数字滤波器的滤波效果，并从理论上进行分析（或解释）。
5、回放语音信号
对滤波后的语音信号进行回放，感觉滤波前后语音信号的变化。
6、设计一个语音信号分析与处理系统界面（选作）
利用MATLAB的界面设计功能，设计一个MATLAB环境下的人机交互界面，完成上述滤波器的设计及滤波过程，要求能够接收输入的设计参数，查看设计结果及滤波结果。
方向二：
数字滤波器设计及在心电信号滤波中的应用。
步骤：
1、心电信号采集
心电信号作为心脏电活动在人体体表的表现，信号一般比较微弱，幅度在10μV~5mV，频率为0.05~100Hz。
在心电信号的采集、放大、检测及记录过程中，有来自外界的各种干扰。
记录一段时间内的人体心电信号波形，要求长度不小于10秒，并对记录的信号进行数字化，保存为数据文件；这里，请同学们使用美国的MIT/BIH心电原始数据，由实验老师给出一定长度的的心电原始数据，数据保存在文件“a01.txt～a10.txt”中，在MATLAB中通过如下语句读取：
%从当前路径下的a01.txt文件读取心电原始数据到变量a01中，a01为二维数据，第一列%为心电信号时间，第二列为心电信号幅度。
load（‘a01.txt’）;
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、心电信号波形观察、频谱观察
对滤波后的心电信号观察其时域、频域特征变化。
绘出滤波后、滤波前、加噪后三个心电信号的差值波形，观察相互间的差异性；同时，分析频谱变化。
学生也可选用持续时间更长的心电原始数据，加上干扰后按上述要求设计滤波器。
3、课程设计要求
1、在两周内学生须上机30小时以上，程序调试完后，须由指导老师在机器上检查运行结果，经教师认可后的源程序可通过打印机输出，并请教师在程序清单上签字。
2、课程设计报告内容和格式：
设计题目，设计的详细步骤，设计过程中的结果、图形等，设计总结。
3、每组每人必须独立完成，成绩的考核按设计结果、答辩成绩及课程设计报告来综合评定。
成绩分为优、良、中、及格、不及格五级分评定。
4、指导教师：
张细政等。
四、设计进度安排
17周周二，E-512,学生上机，按任务要求进行课程设计；17周周三，E-512,学生上机，按任务要求进行课程设计；17周周四，E-512,学生上机，按任务要求进行课程设计；18周周四上午，E-512，完成课程设计任务，答辩并提交课程设计报告。
五、课题的主要功能
5.1问题的描述
自己录取一段语音信号，在MATLAB平台下，进行时域和频域的分析，然后给原始的语音信号叠加上噪声，噪声类型分为如下几种：
（1）白噪声；
（2）单频噪色（正弦干扰）；（3）多频噪声（多正弦干扰）；（4）其它干扰，可设置为低频、高频、带限噪声，或Chirp干扰、冲激干扰。
同样进行时域和频域的分析，在根据这些图谱，分析数据，然后设计不同滤波器，滤去噪音。
绘出叠加噪声后的语音信号时域和频谱图，在视觉上与原始语音信号图形对比，也可通过Windows播放软件从听觉上进行对比，分析并体会含噪语音信号频谱和时域波形的改变。
5.2功能要求
给定滤波器的规一化性能指标（参考指标,实际中依据每个同学所叠加噪声情况而定）例如：
通带截止频率wp=0.25*pi,阻通带截止频率ws=0.3*pi;通带最大衰减Rp=1dB;阻带最小衰减Rs=15dB，每个题目至少设计出5个用不同方法的不同类型滤波器。
题目
（1）：
采用窗函数法与等波纹法分别设计各型FIR滤波器（低通、高通、带通、带阻中的至少3种类型）来对叠加噪声前后的语音信号进行滤波处理，绘出滤波器的频域响应，绘出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；在相同的性能指标下比较各方法的滤波效果，并从理论上进行分析（或解释）。
题目
（2）：
采用双线性变换法与脉冲响应不变法，分别利用不同的原型低通滤波器（Butterworth型与切比雪夫I型）来设计各型IIR滤波器（低通、高通、带通、带阻中的至少3种类型），绘出滤波器的频域响应；并用这些数字滤波器对含噪语音信号分别进行滤波处理，比较不同方法下设计出来的数字滤波器的滤波效果，并从理论上进行分析（或解释）。
6、设计原理
6.1、滤波器的概念
滤波器是一种可通过或阻止某种信号频率的电路。
数字滤波器，事指输入、输出均为数字信号，通过一定运算关系改变输入信号所包含频率成分的相对比例或者滤除某种频率成分的器件。
根据滤波器的功能分类，滤波器可以分成四类：
低通、高通、带通和带阻滤波器。
数字和模拟滤波器一样。
低通滤波器,只允许0～f2的频率成分通过，而大于f2的频率成分衰减为零；高通滤波器与低通滤波器相反，它只允许f1～∞的频率成分通过，而f6.2．用窗函数法设计FIR滤波器
根据过渡带宽及阻带衰减要求，选择窗函数的类型并估计窗口长度N（或阶数M=N-1），窗函数类型可根据最小阻带衰减As独立选择，因为窗口长度N对最小阻带衰减As没有影响，在确定窗函数类型以后，可根据过渡带宽小于给定指标确定所拟用的窗函数的窗口长度N，设待求滤波器的过渡带宽为Δw，它与窗口长度N近似成反比，窗函数类型确定后，其计算公式也确定了，不过这些公式是近似的，得出的窗口长度还要在计算中逐步修正，原则是在保证阻带衰减满足要求的情况下，尽量选择较小的N，在N和窗函数类型确定后，即可调用MATLAB中的窗函数求出窗函数wd（n）。
根据待求滤波器的理想频率响应求出理想单位脉冲响应hd（n），如果给出待求滤波器频率应为Hd，则理想的单位脉冲响应可以用下面的傅里叶反变换式求出：
在一般情况下，hd（n）是不能用封闭公式表示的，需要采用数值方法表示；从w=0到w=2π采样N点，采用离散傅里叶反变换（IDFT）即可求出。
用窗函数wd（n）将hd（n）截断，并进行加权处理，得到
如果要求线性相位特性，则h（n）还必须满足：
根据上式中的正、负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。
要根据所设计的滤波特性正确选择其中一类。
例如，要设计线性相位低通特性可选择h（n）=h（N-1-n）一类，而不能选h（n）=-h（N-1-n）一类。
验算技术指标是否满足要求，为了计算数字滤波器在频域中的特性，可调用freqz子程序，如果不满足要求，可根据具体情况，调整窗函数类型或长度，直到满足要求为止。
6.3.用双线性变换法设计IIR数字滤波器
脉冲响应不变法的主要缺点是产生频率响应的混叠失真。
这是因为从S平面到Ｚ平面是多值的映射关系所造成的。
为了克服这一缺点，可以采用非线性频率压缩方法，将整个频率轴上的频率范围压缩到-π/T～π/T之间，再用z=esT转换到Z平面上。
也就是说，第一步先将整个S平面压缩映射到S1平面的-π/T～π/T一条横带里；第二步再通过标准变换关系z=es1T将此横带变换到整个Z平面上去。
这样就使S平面与Z平面建立了一一对应的单值关系，消除了多值变换性，也就消除了频谱混叠现象，映射关系如图1所示。
图1双线性变换的映射关系
为了将S平面的整个虚轴jΩ压缩到S1平面jΩ1轴上的-π/T到π/T段上，可以通过以下的正切变换实现
（1）
式中,T仍是采样间隔。
当Ω1由-π/T经过0变化到π/T时，Ω由-∞经过0变化到+∞，也即映射了整个jΩ轴。
将式
（1）写成
将此关系解析延拓到整个S平面和S1平面，令jΩ=s，jΩ1=s1，则得
再将S1平面通过以下标准变换关系映射到Z平面
z=es1T
从而得到S平面和Z平面的单值映射关系为：
（2）
（3）
式
（2）与式（3）是S平面与Z平面之间的单值映射关系，这种变换都是两个线性函数之比，因此称为双线性变换
式
（1）与式
（2）的双线性变换符合映射变换应满足的两点要求。
首先,把z=ejω，可得
（4）
即S平面的虚轴映射到Z平面的单位圆。
其次，将s=σ+jΩ代入式（4），得
因此
由此看出，当σ<0时，|z|<1；当σ>0时，|z|>1。
也就是说，S平面的左半平面映射到Z平面的单位圆内，S平面的右半平面映射到Z平面的单位圆外，S平面的虚轴映射到Z平面的单位圆上。
因此，稳定的模拟滤波器经双线性变换后所得的数字滤波器也一定是稳定的。
双线性变换法优缺点
双线性变换法与脉冲响应不变法相比，其主要的优点是避免了频率响应的混叠现象。
这是因为S平面与Z平面是单值的一一对应关系。
S平面整个jΩ轴单值地对应于Z平面单位圆一周，即频率轴是单值变换关系。
这个关系如式（4）所示，重写如下：
上式表明，S平面上Ω与Z平面的ω成非线性的正切关系，如图2所示。
由图2看出，在零频率附近，模拟角频率Ω与数字频率ω之间的变换关系接近于线性关系；但当Ω进一步增加时，ω增长得越来越慢，最后当Ω→∞时，ω终止在折叠频率ω=π处，因而双线性变换就不会出现由于高频部分超过折叠频率而混淆到低频部分去的现象，从而消除了频率混叠现象。
图2双线性变换法的频率变换关系
但是双线性变换的这个特点是靠频率的严重非线性关系而得到的，如式（4）及图2所示。
由于这种频率之间的非线性变换关系，就产生了新的问题。
首先，一个线性相位的模拟滤波器经双线性变换后得到非线性相位的数字滤波器，不再保持原有的线性相位了；其次，这种非线性关系要求模拟滤波器的幅频响应必须是分段常数型的，即某一频率段的幅频响应近似等于某一常数（这正是一般典型的低通、高通、带通、带阻型滤波器的响应特性），不然变换所产生的数字滤波器幅频响应相对于原模拟滤波器的幅频响应会有畸变，如图3所示。
图3双线性变换法幅度和相位特性的非线性映射
对于分段常数的滤波器，双线性变换后，仍得到幅频特性为分段常数的滤波器，但是各个分段边缘的临界频率点产生了畸变，这种频率的畸变，可以通过频率的预畸来加以校正。
也就是将临界模拟频率事先加以畸变，然后经变换后正好映射到所需要的数字频率上。
7、课程设计主要步骤
7.1、信号采集
录制一段课程设计学生的语音信号并保存为文件，要求长度不小于10秒，并对录制的信号进行采样；录制时可以使用Windows自带的录音机，或者使用其它专业的录音软件，录制时需要配备录音硬件（如：
麦克风），为便于比较，需要在安静、干扰小的环境下录音。
7.2、信号分析
利用函数wavread（‘文件名’）对录制的语音信号采样，使用MATLAB绘出采样后的语音信号的时域波形和频谱图。
根据频谱图求出其带宽，并说明语音信号的采样频率不能低于多少赫兹。
[x1,fs,bits]=wavread（'d:
/g'）;%利用wavread进行采样，存入x1，采样频率保存为Fs
sound（x1,fs,bits）;%播放录音
y1=fft（x1,32768）;%对录音取32768个点进行快速傅立叶变换
figure
（1）
subplot（2,1,1）%2表示将图分为两行
plot（x1）;
title（'原始语音信号时域波形'）;
subplot（2,1,2）
plot（abs（y1））;
title（'原始语音信号频谱'）
7.3、含噪信号合成
在MATLAB软件平台下，给原始的心电信号叠加上噪声或干扰，干扰类型分为如下几种：
（1）白噪声；
（2）工频干扰（50Hz）；（3）谐波干扰（二次、三次谐波为主，分别为100Hz、150Hz）；（4）其它干扰，可设置为低频、高频、带限噪声，或冲激干扰。
绘出叠加噪声后的心电信号时域和频谱图，在视觉上与原始心电信号图形对比，绘出其时域波形差，分析频域基本特征变化。
（1）白噪声
x2=awgn（x1,10）;%利用awgn（x,40）进行加噪，性噪比为40分贝，加噪后的信号用x1保存
sound（x2,fs,bits）;%播放录音
y2=fft（x2,32768）;%取32768个点进行快速傅立叶变换
figure
（2）;
subplot（2,1,1）
plot（x2）;
title（'加高斯白噪声语音信号的时域波形'）
subplot（2,1,2）
plot（abs（y2））;
title（'加高斯白噪声语音信号的频谱'）
（2）单频噪色（正弦干扰）；
N=length（x1）-1;%求数组长度，也就是信号长度
t=0:
1/fs:
N/fs;
d=[0.8*sin（2*pi*5000*t）]';%加载单频正弦信号
x3=x1+d;
sound（x3,fs,bits）;%播放语音
figure（3）
subplot（2,1,1）
plot（x3）
title（'加单频正弦信号的语音信号时域波形'）
y3=fft（x3,32768）;
subplot（2,1,2）
plot（abs（y3））;
title（'加单频正弦信号的语音信号频谱'）
（3）多频噪声（多正弦干扰）；
N=length（x1）-1;%求数组的长度，也就是信号长度
t=0:
1/fs:
N/fs;%N为点数
d1=[0.08*sin（2*pi*20*t）]';%加载多频正弦信号
x4=x1（:
1）+d1+d;
sound（x3,fs,bits）;%播放语音
figure（4）
subplot（2,1,1）
plot（x4）
title（'加多频正弦信号的语音信号时域波形'）
y4=fft（x4,32768）;
subplot（2,1,2）
plot（abs（y4））;
title（'加多频正弦信号的语音信号频谱'）
（4）其它干扰
t=0:
length（x1）-1;
Ft=22050;
Fn=150000;%干扰正弦信号的频率
omiga=2*pi*Fn/Ft;
noise1=0.1*sin（omiga*t）;%加入正弦干扰信号，频率1.5khz
x6=x1（:
1）+noise1';
subplot（311）;
plot（t,x6）;
grid;
xlabel（'时间（t）'）;
title（'加噪声的信号波形'）;
fz1=fftshift（fft（x6,length（x1）））;
subplot（312）;
plot（abs（fz1））;%画出加入高频干扰之后的语音信号幅频图
xlabel（'频率'）;
title（'加噪声的信号频谱'）;
grid;
subplot（313）;
plot（angle（fz1））;
xlabel（'频率'）;
title（'加噪声的信号相位'）;
gridon;
sound（x6,Ft）;%播放加入干扰信号后的声音信号
7.4、数字滤波器设计及滤波
给定滤波器的性能指标，完成各种滤波器的设计。
并对叠加噪色前后的语音信号进行滤波处理，绘出滤波器的频域响应，绘出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；在相同的性能指标下比较各方法的滤波效果，并从理论上进行分析（或解释）。
利用窗函数法设计低通滤波器
fp=3500,fc=5000,rs=40;
wp=2*pi*fp/fs;ws=2*pi*fc/fs;
Bt=ws-wp;
alph=0.5842*（rs-21）^0.4+0.07886*（rs-21）;
M=ceil（（rs-8）/2.285/Bt）;
wc=（wp+ws）/2/pi;
hn=fir1（M,wc,kaiser（M+1,alph））;
X=conv（hn,x3）;
sound（X,fs,bits）;
X1=fft（X,32768）
subplot（211）;
plot（X）;title（'处理后的信号波形'）;
subplot（212）;
plot（abs（X1））;
title（'处理后的信号频谱'）;
用双线性变换法设计切比雪夫II数字滤波器
wp=2*pi*3500;ws=2*pi*5000;Rp=10;As=50;
[N1,wp1]=cheb2ord（wp,ws,Rp,As,'s'）;
[B1,A1]=cheby2（N1,As,wp1,'s'）;
fk=0:
4500/32768:
4500;wk=2*pi*fk;
figure（5）
Hk=freqs（B1,A1,wk）;
plot（fk,20*log10（abs（Hk）））;gridon
xlabel（'频率（kHz）'）;ylabel（'幅度（dB）'）;
figure（6）
[Bz2,Az2]=bilinear（B1,A1,fs）;
xx=filter（Bz2,Az2,x3）;
sound（xx,fs,bits）;
subplot（2,1,1）
plot（xx）;
title（'滤波后的信号时域波形'）
yy1=fft（xx,327680）;
subplot（2,1,2）
plot（abs（yy1））;
title（'滤波后的信号信号频谱'）
用双线性变换法设计切比雪夫I数字滤波器
wp=2*pi*3500;ws=2*pi*5000;Rp=10;As=50;
[N1,wp1]=cheb1ord（wp,ws,Rp,As,'s'）;
[B1,A1]=cheby1（N1,As,wp1,'
升级会员 