FIR数字滤波器的(海明)窗函数法设计.doc

课程设计说明书 NO.15

FIR数字滤波器的（海明）窗函数法设计
1．课程设计目的
（1）熟悉并掌握MATLAB中有关声音（wave）录制、播放、存储和读取的函数。
（2）加深对FIR数字滤波器设计的理解，并用窗函数法进行FIR数字滤波器的设计。
（3）将设计出来的FIR数字滤波器利用MATLAB进行仿真。
（4）对一段音频文件进行加入噪声处理，对带有噪声的文件进行滤波处理。
2.设计方案论证
2.1 Matlab语言概述
MATLAB是一种以矩阵运算为基础的交互式程序语言，专门针对科学、工程计算及绘图的需求。随着版本的不断升级，内容不断扩充，功能更加强大，从而被广泛应用于仿真技术、自动控制和数字信号处理领域。
l 此高级语言可用于技术计算
l 此开发环境可对代码、文件和数据进行管理
l 交互式工具可以按迭代的方式探查、设计及求解问题
l 数学函数可用于线性代数、统计、傅立叶分析、筛选、优化以及数积分等
l 二维和三维图形函数可用于可视化数据
l 各种工具可用于构建自定义的图形用户界面
l 各种函数可将基于MATLAB的算法与外部应用程序和语言（如 C、C++、Fortran、Java、COM 以及 Microsoft Excel）集成
l 不支持大写输入，内核仅仅支持小写
2.2声音处理
语音是人类获取信息的重要来源和利用信息的重要手段。语音信号处理是一门发展十分迅速，应用非常广泛的前沿交叉学科，同时又是一门跨学科的综合性应用研究领域和新兴技术。声音是一种模拟信号，而计算机只能处理数字信息0和1。因此，首先要把模拟的声音信号变成计算机能够识别和处理的数字信号，这个过程称为数字化，也叫“模数转换”。在计算机对数字化后的声音信号处理完后，得到的依然是数字信号。必须把数字声音信号转变成模拟声音信号，然后再输出到扬声器，这个过程称为“数模转换”。
2.3数字滤波器的介绍
数字滤波器(digital filter)是由数字乘法器、加法器和延时单元组成的一种装置。其功能是对输入离散信号的数字代码进行运算处理，以达到改变信号频谱的目的。数字滤波器是一个离散时间系统（按预定的算法，将输入离散时间信号转换为所要求的输出离散时间信号的特定功能装置）。应用数字滤波器处理模拟信号时，首先须对输入模拟信号进行限带、抽样和模数转换。数字滤波器输入信号的抽样率应大于被处理信号带宽的两倍，其频率响应具有以抽样频率为间隔的周期重复特性，且以折叠频率即1／2抽样频率点呈镜像对称。为得到模拟信号，数字滤波器处理的输出数字信号须经数模转换、平滑。数字滤波器具有高精度、高可靠性、可程控改变特性或复用、便于集成等优点。
2.4 FIR滤波器基本结构
　　FIR滤波器的数学表达式为：

　　式中：N为FIR滤波器的抽头数；x(n)为第n时刻的输入样本；h(i)为FIR滤波器第i级抽头系数。
　　其相应的z变换为：

　　式中：z-i为N-1阶多项式。
　　在自适应处理、数据通信等领域中往往要求信号在传输过程中不能有明显的相位失真，FIR滤波器可以做到线性相位满足此要求。F1R滤波器实质上是一个分节的延迟线，把每一节的输出加权累加，得到滤波器的输出。对于FIR滤波器的单位脉冲响应h(i)只要满足以下2个条件之一，则为线性相位滤波器。

　　线性相位的FIR滤波器具有中心对称的特性，其对称中心在N／2处。

　　(2)由性能指标确定窗函数w(n)和窗口长度N
　　(3)求得实际滤波器的单位脉冲响应h(n)

　　(4)检验滤波器性能。
设计常用的窗函数有矩形窗、汉宁窗、海明窗、凯撒窗等。其中：海明窗的旁瓣峰值小于主瓣峰值的1％，99．963％的能量集中在主瓣内．且通过海明窗设计的FIR滤波器在较少的阶数下可以得到较小通带纹波，非常适合工程设计。
2.5利用Windows进行语音信号的采集
利用windows下的录音机，按照（开始—>程序—>附件—>娱乐—>录音机，文件—>属性—>立即转换—>8000KHz，8位，单声道）的顺序操作，录制一段自己的语音，录制时间为5秒，如图1，图2所示，将自己录好的语音文件保存为“zf.wav”。

图1 选择windows下的录音机

图2 用8000Hz采样录音
2.6语音信号的分析
(1) 将上一步骤中保存下来的语音信号文件“zf.wav”复制到计算机装有Matlab软件的磁盘中相应Matlab目录中的“work”文件夹中:( C:\Program Files\MATLAB71\work)。
(2) 双击桌面上Matlab软件的快捷图标，打开Matlab软件。
(3) 在Matlab菜单栏中选择“FileànewàM-File”或是点击快捷按钮，打开m文件编辑器。
(4) 在m文件编辑器中输入相应的指令将自己的语音信号导入Matlab工作台。
(5) 画出原始语音信号s的波形，由于原始语音信号开始一段会是无用的语音信号，因此要截取掉，截取的一段语音信号为1至1+fs-1，即从1到7999，画出截取原始语音信号s1的波形，代码如下，波形如图3，所示：
close all
clear all
clc
[s, fs,bits]=wavread('C:\Program Files\MATLAB71\work\zf.wav');
s1=s(1:8000);
sound(s1,fs,bits);
figure(1);
subplot(211)
plot(s)
title('原始语音信号')
subplot(212)
plot(s1)
title('截短语音信号');
会得到下图：

图3原始语音信号和截短语音信号
(6) 对语音信号进行频谱分析，在Matlab中，利用函数FFT()对信号进行快速傅里叶变换，得到信号的频谱特性，如图4所示
wavwrite(s1,fs,'s1.wav'); %将被处理信号s1输出为语音文件“s1.wav”
S1=fft(s1);
figure(2)
subplot(311);
plot(s1);
title('截短预处理语音信号')
subplot(312)
plot(abs(S1))
title('预处理语音信号频谱');
subplot(313);
k=0:4000;
plot(k(1:4000)*1,abs(S1(1:4000)));
title('预处理语音信号单边带频谱')

图4 截短预处理语音信号
如图4所示，从右向左看，第一个较大的波峰所在的频率即为3db截止频率，第二个波峰所对应的频率为通带截止频率，在图中可以读出=610Hz；一般在3db截止频率右侧的波谷位置选择阻带截止频率=750Hz。通带截止频率、阻带截止频率数值的确定，就可以确定滤波器的基本指标。
图4里第二个图是信号的FFT结果，即是信号的实际频谱采样，本设计中信号的长度取的是L=8000点，所以，图中，的每两个相邻点之间的频率间隔大小，即频率分辨率： （Hz），所以，根据它的放大图，即图4中的第三个图的放大图，可以确定。
2.7滤波器的设计
2.7.1滤波器的参数设定
本次课设我设计的是一个线性FIR低通滤波器，利用的窗函数是hamming 窗，如上（图4）所述：
所以：
通带截止频率为 (rad/sample)
阻带截止频率为 (rad/sample)
(rad/sample)
由于海明窗过渡带满足：

求得滤波器阶数 189
=
N

94
2
1
=
-
=
N

（1） 给定所要求的频率响应函数

（2） 求单位采样响应

（3） 海明窗

（4） 滤波器的单位采样响应：

)
(
)]
94
cos(
46
.
0
54
.
0
)][
94
(
[
17
.
0
189
n
R
n
n
w
Sa
c
p
-
-
=

2.7.2滤波器的MATLAB仿真
在M文件中继续编写代码，把计算出来的参数带入代码中。代码如下：
%%%%%%%%%加噪声完成信号截取
[s, fs,bits]=wavread('C:\Program Files\MATLAB71\work\zf.wav');
s1=s(1:8000);
sound(s1,fs,bits);
figure(1);
subplot(211)
plot(s)
title('原始语音信号')
subplot(212)
plot(s1)
title('截短语音信号');
wavwrite(s1,fs,'s1.wav');
S1=fft(s1);
figure(2)
subplot(311);
plot(s1);
title('截短预处理语音信号')
subplot(312)
plot(abs(S1))
title('预处理语音信号频谱');
subplot(313);
k=0:4000;
plot(k(1:4000)*1,abs(S1(1:4000)));
title('预处理语音信号单边带频谱')
s2=awgn(s1,15); %%%%%%%%%%%%%%%%%%%%%%%%%%完成加噪！15db
wavwrite(s2,fs,'s2.wav');
figure(3);
subplot(211);
plot(s2);
title('加噪后语音信号');
subplot(212);
S2=fft(s2);
plot(abs(S1));
title('加噪后信号频谱');
figure(4)
subplot(211);
plot(s1);
title('语音信号');
subplot(212);
plot(s2);
title('加噪后语音信号');

%%%%%%%%%滤波器完成相关参数配置
wp=610*2*pi/8000;
wst=750*2*pi/8000;
wc=(wp+wst)/2;
N=ceil(3.3*2*pi/(wst-wp))+1;
r=(N-1)/2;
hn1=fir1(N-1,wc/pi,'low',hamming(N));%
s3=conv(s2,hn1);
wavwrite(s3,fs,'s3.wav');
S3=fft(s3);
figure(5)
freqz(hn1);
title('滤波器幅频特性与相频特性')
figure(6)
subplot('111')
plot(hn1);
title('滤波器系统函数');
figure(7)
subplot(211)
plot(s3)
title('滤波器处理之后信号图')
subplot(212);
plot(abs(S3));
title('滤波器处理之后频谱');
figure(8)
subplot(211)
plot(s2);
title('加噪后语音信号');
subplot(212);
plot(s3);
title('滤波器处理之后信号图');
%%%%%%%%%%%%%%求信噪比snr在Workspace中体现！
s4=conv(s1,hn1);
%p1=sum(s1.^2);
%p2=sum(s2.^2)-sum(s1.^2);
%SNR1=10*log10(p1/p2);
p1=sum(s4.^2)/8000;
p2=sum(s3.^2)/8000-sum(s4.^2)/8000;
SNR2=10*log10(p1/p2);
运行代码效果如图5所示。的噪声生成叠加指令为：awgn，所加的噪声为15 。

图5 加噪后语音信号和频谱

图6 语音信号和加噪后语音信号
图6为原信号和加噪后信号的波形对比图，所加噪声为计算机随机生成的高斯白噪声。

图7 滤波器幅频特性与相频特性
设计的滤波器是用单位采样响应h(n)表示的，可以利用带噪声语音与h(n)做时域卷积，即： 。在Matlab中，卷积运算可以用函数“conv( )”实现。

图8滤波器系统函数

图9 加噪后语音信号与滤波器处理之后信号图
3.设计结果与分析
3.1设计结果
3.1.1比较滤波前（含噪声信号的文件）和滤波后的语音信号效果。
图10利用hamming窗FIR低通滤波器滤波效果图
3.1.2滤波前（含噪声）的信号和滤波后信号的信噪比
利用“信噪比代码段”做出来的效果图，图中名称snr为最终算出来的信噪比，由于加入的噪声为15 ，而用MATLAB代码算出来的信噪比最终达到了效果。

图11 相关参数的截图，其中snr为信噪比
3.1.3信号采样的频率分辨率
频率分辨率计算公式，其中,N=8000
△f≈1 Hz
3.2结果分析
在运行过程中，每次程序运行后得到信噪比SNR2的数值都是不同的，其原因是在先前用awgn（）函数所加的高斯白噪声是系统随机产生的。从结果来看，可以看到滤波后的信噪比SNR2=15.838大于先前加入高斯白噪声的信号的信噪比SNR1=15db，证明滤波器各项参数设置正确，分别试听截取后语音信号s1、s2、s3，发现滤波后的信号s3相比加噪后s2清晰度有明显改善，原来的噪声基本被滤除，证明滤波器能够达到要求。但是在实际应用中并不会出现随机噪声的情况，信噪比会表现的比较稳定。相比于模拟信号，数字信号对于信噪比的要求不是十分严格。
4．设计体会
通过这次课程设计使我懂得了理论与实际相结合是很重要的，只有理论知识是远远不够的，只有把所学的理论知识与实践相结合起来，从理论中得出结论，才能真正为社会服务，从而提高自己的实际动手能力和独立思考的能力。在设计的过程中遇到问题，可以说得上是困难重重，然而经过同学们不断的查找资料、分析处理，使得我加强了对数字信号处理过程的认识，更加强了对程序的理解，同时在做课设的时候也分析了信号的流程，切熟知了原理图。同时在设计的过程中发现了自己的不足之处，对以前所学过的知识理解得不够深刻，掌握得不够牢固，比如说双线性过程不是很了解……通过这次课程设计之后，一定把以前所学过的知识重新温故。
这次课程设计终于顺利完成了，在设计中遇到了很多地方不是很理解，最后在王老师的辛勤指导下，终于游逆而解。同时，在王老师的身上我学得到很多实用的知识，在次我表示感谢！同时，对给过我帮助的所有同学和各位指导老师再次表示忠心的感谢！
5．参考文献
1] 黄建明，隋燕. 基于数字滤波器设计的讨论. 南开大学学报(自然科学版),2003,26(3):15-18
[2] 丁玉美，高西全.数字信号处理（第二版）. 西安:西安电子科技大学出版社,2000：151-317
[3] 陈怀琛，吴大正，高西全.MATLAB及在电子信息课程中的应用（第二版）.北京:电子工业出版社,2003：174-233
[4] 程佩青，数字信号处理教程.北京:清华大学出版社,2001:186-275