课程设计报告正文311.docx
《课程设计报告正文311.docx》由会员分享，可在线阅读，更多相关《课程设计报告正文311.docx（23页珍藏版）》请在冰豆网上搜索。
课程设计报告正文311
语音信号滤波去噪
——基于FLATTOP窗的线性FIR滤波器
学生姓名：
徐康指导老师：
曹敦
摘要本课程设计主要内容是设计利用窗口设计法选择FLATTOP设计一个线性相位FIR滤波器[1]，对一段含噪语音信号进行滤波去噪处理并根据滤波前后的波形和频谱分析滤波性能。
本课程设计仿真平台为MATLAB7.0，开发工具是M语言编程。
首先利用WindowsXP系统自带的录音软件录制一段语音信号，并人为加入一单频噪声，然后对信号进行频谱分析以确定所加噪声频率，并设计滤波器进行滤波去噪处理，最后比较滤波前后的波形和频谱并进行分析。
由分析结果可知，滤波器后的语音信号与原始信号基本一致，即设计的FIR滤波器能够去除信号中所加单频噪声，达到了设计目的。
关键词滤波去噪；FIR滤波器；FLATTOP窗；MATLAB
1引言
本课程设计主要解决在含噪情况下对语音信号的滤波去噪处理，处理时采用的是利用窗口设计法选择FLATTOP设计的线性相位型FIR滤波器。
数字滤波器是对数字信号实现滤波的线性时不变系统。
数字滤波实质上是一种运算过程，实现对信号的运算处理。
输入数字信号（数字序列）通过特定的运算转变为输出的数字序列，因此，数字滤波器本质上是一个完成特定运算的数字计算过程，也可以理解为是一台计算机。
描述离散系统输出与输入关系的卷积和差分方程只是给数字信号滤波器提供运算规则，使其按照这个规则完成对输入数据的处理。
1.1课程设计目的
通过课程设计掌握熟悉Matlab编程环境，掌握带通数字滤波器的设计方法并能在实际中学会选择和应用；用FLATTOP窗设计满足指标的有限长脉冲响应数字滤波器（FIR）并且以图文的形式展示出来，比较器滤波前后的波形和频谱并进行分析，根据结果和学过的理论得出合理的结论。
对信号传输过程中的去燥过程有更为深入明确的理解。
1.2课程设计的要求
要求：
（1）滤波器指标必须符合工程实际。
（2）设计完后应检查其频率响应曲线是否满足指标。
（3）处理结果和分析结论应该一致，而且应符合理论。
（4）独立完成课程设计并按要求编写课程设计报告书。
2基本理论
2.1FIR滤波器的特点
FIR数字滤波器的设计方便而且简单，特别是采用了双线性变换法设计的数字滤波器没有频谱混叠的问题，效果较好。
在有些实际的应用场合，例如数据传输和图象处理等对滤波器的线性相位特性要求较为严格，所以在这种场合中，一般选用FIR数字滤波器。
FIR数字滤波器可设计得具有严格的线性相位，而且其幅度特性可以随意设计。
FIR滤波器的单位冲激响应h（n）是有限长序列，其系统函数的极点位于Z平面原点，因此FIR滤波器不存在稳定问题[2]。
另外，FIR滤波器还可以采用FFT方法实现其功能，从而大大提高了其效率。
因此，FIR数字滤波器日渐引起了人们的注意。
对于线性相位滤波器而言，通过滤波器的信号不存在相位失真，只存在固定单位的延迟。
而这一点在数字通信、图象处理和语音合成等许多应用中都极为的重要。
FIR数字滤波器具有线性相位特征的条件是它的单位冲激响应h（n）满足偶对称或者奇对称。
若h（n）满足偶对称，则其相位特性
为：
（1-1）
若h（n）满足奇对称，则其相位特性
为：
（1-2）
从频域上看h（n）是长度为N的实序列，所以h（n）的N点DFT满足共轭对称性，即：
H（k）=H*（N-k）而幅度函数则一定呈偶对称特性，即：
（1-3）
当h（n）为偶对称，信号通过该滤波器时，其通带就产生（N-1）/2点群延迟，滤波器的相频特性是相位过原点的一条直线，其斜率为[-（N-1）/2]，其幅频特性相当于对信号低通滤波；当h（n）为奇对称，信号通过该滤波器时，其通带不但产生（N-1）/2点群延迟，还对所有通带内各种频率的信号都有一个90°移相，这相当于信号先经过一个90°移相器，然后再作带通滤波。
FIR数字滤波器设计的基本步骤如下：
（1）确定指标
在设计一个滤波器之前，必须首先根据工程实际的需要确定滤波器的技术指标。
在很多实际应用中，数字滤波器常常被用来实现选频操作。
因此，指标的形式一般在频域中给出幅度和相位响应。
幅度指标主要以两种方式给出。
第一种是绝对指标。
它提供对幅度响应函数的要求，一般应用于FIR滤波器的设计。
第二种指标是相对指标。
它以分贝值的形式给出要求。
（2）逼近
确定了技术指标后，就可以建立一个目标的数字滤波器模型。
通常用理想的数字滤波器模型。
之后，利用数字滤波器的设计方法，设计出一个实际滤波器模型来逼近给定的目标。
（3）性能分析和计算机仿真
上两步的结果是得到以差分或冲激响应描述的滤波器。
根据这个描述就可以分析其频率特性和相位特性，以验证设计结果是否满足指标要求；或者利用计算机仿真实现设计的滤波器，再分析滤波结果来判断。
FIR滤波器的设计问题实际上是确定能满足所要求的转移序列或脉冲响应的常数的问题，设计方法主要有窗函数、频率响应法和等波纹最佳逼近法等。
2.2窗口设计法
窗口设计法是一种通过截断和计权的方法使无限长非因果序列成为有限长脉冲响应序列的设计方法。
通常在设计滤波器之前，应该先根据具体的工程应用确定滤波器的技术指标。
在大多数实际应用中，数字滤波器常常被用来实现选频操作，所以指标的形式一般为在频域中以分贝值给出的相对幅度响应和相位响应。
窗口设计法步骤如下：
（1）根据过渡带宽及阻带衰减要求，选择窗函数的类型并估计窗口长度N。
窗函数的类型可根据最小阻带衰减AS独立选择。
（2）根据待求滤波器的理想频率响应求出理想单位脉冲响应hd（n）。
（3）由性能指标确定窗函数W（n）和长度N。
（4）求得实际滤波器的单位脉冲响应h（n），h（n）即为所设计FIR滤波器系数向量b（n）。
（2.1）
常见的窗函数性能表如下2.1图所示：
名称
滤波器
过渡带宽
最小阻带衰减
名称
滤波器
过渡带宽
最小阻带衰减
矩形
1.8π/M
21dB
PARZENWIN
6.6π/M
56dB
巴特利特
6.1π/M
25dB
FLATTOPWIN
19.6π/M
108dB
汉宁
6.2π/M
44dB
GAUSSWIN
5.8π/M
60dB
汉明
6.6π/M
51dB
BARTHANNWIN
3.6π/M
40dB
布莱克曼
11π/M
74dB
BLACKMANHARRIS
16.1π/M
109dB
BOHMANWIN
5.8π/M
51.5dB
CHEBWIN
15.2π/M
113dB
NUTTALLWIN
15.4π/M
108dB
TUKEYWIN
2.4π/M
22dB
图2.1常见窗函数性能表
2.3FALTTOP窗
w=Flattopwin（L）返回L-点Flattopwin窗口中列向量。
Flattopwin窗的滤波器的过渡带宽为19.6π/M，最小阻带衰减108db。
定义为：
（2.2）
时间波形和幅度谱：
（2.3）
图2.2FLATTOP窗的时间波形和幅度谱
3设计步骤
3.1设计流程图
根据设计的要求，首先录制一段长约1～2s的自我介绍语音信号，对录制的语音信号加入噪声干扰，再利用Flattopwin窗设计合理的FIR滤波器。
最后用滤波器对干扰后的语音信号进行滤波去噪。
具体设计流程图如下图3.1所示：
无效
有效
图3.1设计流程图
3.2录制语音信号
在WindowsXP环境下，打开“开始——程序——附件——录音机”，在录音机程序的“文件——属性——立即转换”中将声音选定属性选择为“8.000kHz，8位，单声道，7kB/s”,
图3.2XP系统录音机设置
录制语音信号“大家好，我是徐康”，时间在1～2s左右，然后保存为“.wav”格式的音频文件
3.3语音信号的频谱分析
然后在Matlab软件平台下，利用函数wavread对语音信号进行采样，记住采样频率和采样点数。
通过wavread函数的使用，理解采样频率、采样位数等概念。
采集完成后在信号中加入一个单频噪声，设计的任务即为从含噪信号中滤除单频噪声，还原原始信号。
具体调用为新建“voice_noise.m”文件，代码如下:
[x,fs,bits]=wavread（'d:
\astonxk.wav'）;%%输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。
sound（x,fs,bits）;%按指定的采样率和每样本编码位数回放
N=length（x）;%计算信号x的长度
fn=2000;%单频噪声频率，2000Hz
t=0:
1/fs:
（N-1）/fs;%计算时间范围，样本数除以采样频率
x=x';y=x+0.1*sin（fn*2*pi*t）;%将原始语音信号与单频噪声混叠
sound（y,fs,bits）;%应该可以明显听出有尖锐的单频啸叫声
X=abs（fft（x））;Y=abs（fft（y））;%对原始信号和加噪信号进行fft变换,取幅度谱
X=X（1:
N/2）;Y=Y（1:
N/2）;%截取前半部分
deltaf=fs/N;%计算频谱的谱线间隔
f=0:
deltaf:
fs/2-deltaf;
原始语音信号经过加噪处理后，可以通过Matlab软件的plot函数分别对原始语音信号和加噪之后的混叠信号进行时域和频域的图形对比,获得最为直观的差异，具体程序如下：
subplot（2,2,1）;plot（t,x）
title（'原始语音信号'）;xlabel（'时间（单位：
s）'）;ylabel（'幅度'）;
axis（[0,2,-1,1]）;gridon
%title（'原始语音信号的时域图'）;xlabel（'t'）;ylabel（'x'）;
subplot（2,2,2）;plot（f,X）;axis（[0,3000,0,300]）
title（'语音信号幅度谱图'）;xlabel（'频率（单位：
Hz）'）;ylabel（'幅度谱'）;gridon
%title（'原始语音信号的频域图'）;xlabel（'f'）;ylabel（'X'）;
subplot（2,2,3）;plot（t,y）;
title（'加入单频干扰后的语音号'）;xlabel（'时间（单位：
s）'）;ylabel（'幅度'）;axis（[0,2,-1,1]）;
gridon
%title（'加噪干扰后的时域图'）;xlabel（'t'）;ylabel（'y'）;
subplot（2,2,4）;plot（f,Y）
title（'加入干扰后的语音信号幅度谱图'）;
xlabel（'频率（单位Hz）'）;ylabel（'幅度谱'）;axis（[0,3000,0,300]）;gridon
%title（'加噪干扰的频域图'）;xlabel（'f'）;ylabel（'Y'）;
运行程序后得到的波形图如图所示：
图3.3加噪前后信号时域图与频谱图的比较
3.4滤波器的设计
滤波器的设计就是要找到一组能满足特定滤波要求的系数向量a和b，其中b为系统函数的分子系数，a为系统函数分母系数。
在本次的课程设计中我所采用的就是利用Flattopwin窗函数来设计FIR滤波器。
因为加入的单频噪声的频率为2000Hz，HzH所以滤波器的性能指标定为：
fpd=200;fsd=1950;fsu=2050;fpu=3800;Rp=1;As=100;%带阻滤波器设计指标
fcd=（fpd+fsd）/2;fcu=（fpu+fsu）/2;df=min（（fsd-fpd）,（fpu-fsu））;
%计算上下边带中心频率，和频率间隔
wcd=fcd/fs*2*pi;wcu=fcu/fs*2*pi;dw=df/fs*2*pi;
%将Hz为单位的模拟频率换算为rad为单位的数字频率
wsd=fsd/fs*2*pi;wsu=fsu/fs*2*pi;
由于是基于FLATTOP窗设计的滤波器，参考常见窗函数性能表，可知，FLATTOP窗的过渡带宽为
，最小阻带衰减为108dB；所以：
M=ceil（19.6*pi/dw）+1;%计算窗设计该滤波器时需要的阶数
n=0:
M-1;%定义时间范围
w_par=（Flattopwin（M））;%产生M阶的Flattopwin窗
hd_bs=ideal_lp（wcd,M）+ideal_lp（pi,M）-ideal_lp（wcu,M）;%调用自编函数计算理想带阻滤波器的脉冲响应
h_bs=w_par'.*hd_bs;%用窗口法计算实际滤波器脉冲响应
[db,mag,pha,grd,w]=freqz_m（h_bs,1）;%调用自编函数计算滤波器的频率特性
figure（3）
subplot（2,2,1）;plot（w,db）;title（'滤波器幅度响应图'）;xlabel（'w/pi'）;ylabel（'相对振幅db'）;
line（[0,3.5],[-As,-As],'color','r','linestyle','--','LineWidth',1）;
line（[0,3.5],[-Rp,-Rp],'color','r','linestyle','--','LineWidth',1）;
line（[wsd,wsd],[-350,50],'color','r','linestyle','--','LineWidth',1）;
line（[wsu,wsu],[-350,50],'color','r','linestyle','--','LineWidth',1）;
subplot（2,2,2）;plot（w,mag）;title（'滤波器幅度响应图'）;xlabel（'w/pi'）;ylabel（'幅度mag'）;
subplot（2,2,3）;plot（w,pha）;title（'滤波器相位响应图'）;xlabel（'w/pi'）;ylabel（'相位pha'）;
axis（[0,3.5,-4,4]）;
subplot（2,2,4）;stem（n,h_bs）;title（'滤波器脉冲响应图'）;xlabel（'n'）;ylabel（'脉冲响应h（n）'）;
axis（[0,60,0,1]）;
图3.4设计的滤波器图形
3.5信号滤波处理
滤波器设计完成后，在MATLAB平台上用函数filter实现滤波。
新建“filter.m”文件，滤波程序如下：
y_fil=filter（h_bs,1,y）;%用设计好的滤波器对y进行滤波
Y_fil=fft（y_fil）;Y_fil=Y_fil（1:
N/2）;%计算频谱取前一半
figure（4）
suptitle（'语音信号滤波去噪--单频噪声频率为2k,'）;
subplot（3,2,1）;plot（t,x）;xlabel（'时间（t）'）;ylabel（'幅度'）;title（'原始语音信号'）;
subplot（3,2,2）;plot（f,X）;xlabel（'频率（f）'）;ylabel（'幅度谱'）;
title（'原始语音信号幅度谱'）;axis（[0,5000,0,300]）;
subplot（3,2,3）;plot（t,y）;xlabel（'时间（t）'）;ylabel（'幅度'）;title（'加干扰后的语音信号'）;
subplot（3,2,4）;plot（f,Y）;xlabel（'频率（f）'）;ylabel（'幅度谱'）;
title（'加干扰后的语音信号幅度谱'）;axis（[0,5000,0,300]）;
subplot（3,2,5）;plot（t,y_fil）;xlabel（'时间（t）'）;ylabel（'幅度'）;title（'滤波后语音信号'）;
subplot（3,2,6）;plot（f,Y_fil）;xlabel（'频率（f）'）;ylabel（'幅度谱'）;
title（'滤波后语音信号幅度谱'）;axis（[0,5000,0,300]）;
sound（x,fs,bits）;%播放原始语音信号
sound（y,fs,bits）;%播放原始语音与单频噪声混叠之后的信号
sound（y_fil,fs,bits）;%播放经过滤波器滤波之后的语音信号
得到的滤波前后音乐信号的时域波形图和频谱图对比图如图所示：
图3.5滤波前后语音信号的时域波形图和频谱图对比图
调用wavwrite函数导出原始语音与噪声的混叠音频文件和滤波之后的音频文件
wavwrite（y,fs,bits,'Befor_filter.wav'）;
wavwrite（y,fs,bits,'After_filter.wav'）;
3.6结果分析
由图3.3所示，一方面此滤波器的阻带截止频率为0.6左右，而此处第一幅小图里面的As对应的纵坐标同时也在16dB左右；另一方面在图3.4中示，滤波前的噪声频率大约在3000Hz左右，滤波后的噪声频率大约为350Hz，大致符合滤波前后公式db（As） =20*log10（（mag+eps）/max（mag））（eps，mag分别为滤波前后的噪声频率）中对滤波前后噪声频率的要求。
综合上述两个方面的分析，可以确定的是此滤波器符合设计要求。
在MATLAB中，经过sound（x,fs,bits）函数，对经过Flattopwin窗设计的FIR滤波器之后的音乐信号进行回放，可以听出滤波之后的音乐信号跟原始音乐信号一样清晰，完全滤除掉了噪声的干扰，通过下面语句来进行音乐信号回放比较：
>>sound（y_fil,fs,bits）;播放滤波后的音乐信号
所得结果证明了用Flattopwin窗设计的FIR滤波器和音乐信号去噪设计是成功的。
3.7绘制滤波器结构图
由前面程序运行的结果，直接可以在Matlab的“Workspace”里面查找到设计的FIR滤波器的脉冲相应序列h_bs，其值为：
2.2289e-006-5.8891e-0066.1146e-0070.00019324-0.000386840.00028061-0.001059-0.000332590.000338350.0024567-0.000945010.0037688-0.002341-0.000112260.00411090.0004012-0.0065299-0.0552170.050638-0.12170.260480.0510830.314950.314950.0510830.26048-0.12170.050638-0.055217-0.00652990.00040120.0041109-0.00011226-0.0023410.0037688-0.000945010.00245670.00033835-0.00033259-0.0010590.00028061-0.000386840.000193246.1146e-007-5.8891e-0062.2289e-006
系数计算好以后，将滤波器结构用visio或其他绘图软件画好，填好系数即可。
图3.6设计的线性相位FIR滤波器结构图
4异常处理
1．在录音时，原是在自己笔记本电脑的Windows7系统下使用系统自带的录音软件；录音的，结果录制的音频文件不是所要求的pcm格式单声道8k信号。
2．绘制原始语音信号与加噪之后的混叠信号的频谱图时，由于横纵坐标设置不当导致图像信息不清，难以识别。
3．在确定带通滤波器的上下边带各频率时，由于取值不当导致滤波器的阶高达几百甚至一千多，直接影响后续的Visio软件绘图。
针对以上问题，相应的解决方案如下：
1．经网上查询以后，了解到Win7系统自带的录音软件不具备设置音频文件属性，后在学校机房电脑的XP系统环境下录制了符合要求的音频文件。
2．增大过渡带宽以后将滤波器阶数降到了五十阶以下。
3．使用axis（）函数调整坐标后最终达到了理想的状态。
5结束语
本次的课程设计，我的任务是利用Flattopwin窗函数设计线性相位FIR滤波器对语音信号滤波去噪。
在本次课程设计之前，我对Flattopwin窗函数完全没有了解，因此在看到这个题目时，我是一头雾水。
但是通过自己翻阅资料和询问同学，我掌握了用Flattopwin窗函数设计FIR滤波器的方法步骤，了解了窗函数的基本设计流程。
经过这两周的课程设计，我学会了很多东西。
在此我要感谢学校为我们提供这次课程设计的机会，感谢老师对我的悉心指导，也感谢同学对我的帮助。
这次的课程设计让我理论联系实际，不仅巩固了我们的理论知识，还提高了我的动手能力，经过这两周的课程设计经过为期两周的课程设计，这次课程设计比较的简单，根据老师给出的模板，再结合以前学过的知识，并在同学和老师的帮助下我顺利的完成了任务。
不同于在教室里上的理论，这次的课程设计需要将我们平时所学习的知识运用到实践之中，将知识学以致用。
因为是以所学理论为基础，所以在课程设计的过程中，我又重温滤波器、窗函数等知识，更加熟悉了MATLAB的操作。
课程设计是我们运用所学知识，动手实践的一个很好的机会。
它既可以帮助我们加深对所学知识的理解，又能提高我们运用知识，联系实际，动手实践的能力。
而且在设计过程中可能用到我们没学过的知识，需要我们去查阅资料获取相关信息，这又提高了我们查找信息和学习新知识的能力。
在实物的调试与检测过程中，又会遇到许多意想不到的问题，需要我们去分析原因和解决问题。
也体会到真正的去独立地完成一件事情是很困难的，同学以及老师的帮助和提醒是必须的。
通过这次课程设计，我拓宽了知识面，锻炼了实际操作能力，综合素质也得到了提高，进一步加深了了我们对专业的认识和激发了我们对专业的兴趣。
虽然课程设计结束了，但是我们的学习还没结束，对知识的进一步学习还需要继续，很开心成功地完成了这次设计。
参考文献
[1]程佩青.数字信号处理教程（第四版）.清华大学出版社，2013-2-1
[2]林志源.FIR和IIR数字滤波器的设计及稳定性研究林志源.《广西师范大学》2009年硕士论文
附录1：
语音信号滤波去噪设计源程序清单
%程序名称：
voice_noise.m
%程序功能：
原始语音信号的频谱分析。
%程序作者：
徐康
/%最后修改日期：
2015-3-10
[x,fs,bits]=wavread（'d:
\astonxk.wav'）;
sound（x,fs,bits）;%按指定的采样率和每样本编码位数回放
N=length（x）;%计算信号x的长度
fn=2000;%单频噪声频率，此参数可改
t=0:
1/fs:
（N-1）/fs;%计算时间范围，样本数除以采样频率
x=x';y=x+0.1*sin（fn*2*pi*t）;
sound（y,fs,bits）;
X=abs（fft（x））;Y=abs（fft（y））;
升级会员 