窗函数设计滤波器.(DOC)

******************

实践教学

*******************

XXXX大学

计算机与通信学院

2014年春季学期

《信号处理》课程设计

题目：窗函数法设计FIR低通滤波器

专业班级：

姓名：

学号：

指导教师：

成绩：

摘要

在数字信号处理中．数字滤波器十分重要并已获得广泛应用，效字滤渡

器与模拟滤渡器比较，具有精度高、稳定、体积小、重量轻、灵活、不要求阻抗匹配以及实现模拟滤波器无法宴现的特殊滤波功能等优点．使数字滤波器成为信号处理领域的一十要研究课题。

本课程设计主要是对一段语音信号，加入噪声后，用汉宁窗设计出的FIR滤波器对加入噪声后的语音信号进行滤波去噪处理。在此次课程设计中，系统操作平台为Windows 7，程序设计的操作软件MATLAB 7.0。此课程设计首先是用音乐软件采集一段语音信号，加入噪声，然后采用汉宁窗函数法设计出FIR滤波器，再用设计出的滤波器对这段加噪后的语音信号进行滤波去噪，最后对前后时域和频域的波形图进行对比分析，从波形可以看出噪

声被完全滤除，达到了语音不失真的效果。

关键字滤波器 ;汉宁窗 ;去噪

III

目录

前言 (1)

第一章设计介绍 (2)

1.2 MATLAB在DSP中的运用 (2)

1.3 FIR数学滤波器介绍 (2)

1.4 窗函数法 (3)

1.5 汉宁窗 (5)

第二章基本原理 (5)

2.1 汉宁窗基本原理 (6)

2.2 窗函数法设计基本原理 (6)

第三章设计实现框图 (8)

第四章基于MATLAB的仿真程序及结果 (9)

4.1 读入语音数据并画出其时域与频域图 (9)

4.2 在MATLAB中对信号加入两个噪声 (9)

4.3 滤波器设计基本程序： (10)

参考文献： (15)

附录：程序源代码 (16)

致谢 (19)

前言

在当今数字信号处理中，滤波器十分重要并以获得广泛应用，它是去除信号中噪声的基本手段，因此，滤波运算时信号处理中的基本运算。滤波器设计问题也是数字信号处理中的基本问题。FIR 滤波器的设计方法主要是建立在对理想滤波器频率特性做某种近似的基础上的。这些近似方法有窗函数法、频率抽样法及最佳一致逼近法。

本设计以窗函数法设计FIR 数字滤波器，窗函数设计法的基本思想是用FIRDF 逼近希望的滤波特性。设希望逼近的滤波器的频率响应为()ω

j d e H ，其单

位脉冲响应用()n h d 表示。为了设计简单方便，通常选择

()ωj d e H 为具有片段常数特性的理想滤波器。因此()n h d 是无限长非因果序列，不能直接作为FIR 滤波器的单位脉冲响应。窗函数设计法就是截取

()n h d 为有限长的一段因果序

列，并用合适的窗口函数进行加权作为FIR 滤波器的单位脉冲响应()n h 。

第一章设计介绍

1.1 数字信号处理简介

数字信号处理（Digital Signal Processing，简称DSP）是一门涉及许多学科而又广泛应用于许多领域的新兴学科[2]。20世纪60年代以来，随着计算机和信息技术的飞速发展，数字信号处理技术应运而生并得到迅速的发展。在过去的二十多年时间里，数字信号处理已经在通信等领域得到极为广泛的应用。数字信号处理是利用计算机或专用处理设备，以数字形式对信号进行采集、变换、滤波、估值、增强、压缩、识别等处理，以得到符合人们需要的信号形式。

1.2 MATLAB在DSP中的运用

MATLAB在信号与系统中的应用主要包括符号运算和数值计算仿真分析。由于信号与系统课程的许多内容都是基于公式演算，而MATLAB借助符号数学工具箱提供的符号运算功能，能基本满足信号与系统课程的需求。例如解微分方程、傅里叶正反变换、拉普拉斯正反变换和z正反变换等。MATLAB在信号与系统中的另一主要应用是数值计算与仿真分析，主要包括函数波形绘制、函数运算、冲击响应与阶跃响应仿真分析、信号的时域分析、信号的频谱分析、系统的S域分析和零极点图绘制等内容。数值计算仿真分析可以帮助学生更深入地理解理论知识，并为将来使用MATLAB进行信号处理领域的各种分析和实际应用打下基础。

1.3 FIR数学滤波器介绍

数字滤波器对信号滤波的方法是：用数字计算机对数字信号进行处理，处理就是按照预先编制的程序进行计算。数字滤波器的原理如图所示，它的核心是数字信号处理器。在本设计中，要求学生掌握利用窗函数的方法设计一个FIR低通数字滤波器。数字滤波器是一个能够完成特定任务的离散时间系统．它可以利用有限精度算法来实现。当采用一个因果稳定的离散线性移不变系统的系鲈函数去逼近滤渡器所要求的性能指标时．由于系统函数有无

限长单位冲激响应ⅡI R1系统函数和有限长单位冲激响应系统函数两种，相应地数字滤渡器也就有无限长单位冲激响应滤渡器和有限长单位冲激响应滤渡器两种剐。FIR 滤渡器的设计问题实质上是确定能满足所要求的转移序列或脉冲响应的常数的问题．设计方法主要有窗函数法、频率采样法和等波纹最佳逼近法等。FIR 的突出优点是：系统总是稳定的、易于实现线性相位、允许设计多通带（或多阻带）滤波器，但与IIR 相比，在满足同样阻带衰减的情况下需要的阶数较高，滤波器的阶数越高，占用的运算时间越多，因此在满足指标要求的情况下应尽量减少滤波器的阶数。

1.4 窗函数法

FIR 数字滤波器设计最简单的方法是窗函数法，通常也称为傅立叶级数法。它是在时域进行的，因而必须由理想滤波器的频率响应)(jw d e H 推导出其单位冲激响应)(n h d ，在设计一个FIR 数字滤波器的单位冲激响应)(n h 去逼近)(n h d 。窗函数是一种用一定宽度窗函数截取无限长脉冲响应序列获取有限长脉冲响应序列的设计方法。丽其设计Fm 滤披器的基本思想：根据培定的滤渡器技术指标选取滤波器长度N 和窗函数“靠J ，使其具有最窄宽度的主瓣和最小的旁瓣。其核心是从给定的频率特性．通过加窗确定有限长单位脉冲响应序列h(n)即实际滤波器的系散向量．其是由理想滤渡器脉冲响应hd “n 与窗函敷函数相加相乘得到。

使用窗函数法设计时要满足以下两个条件： ○

1窗谱主瓣尽可能地窄，以获得较陡的过渡带； ○

2尽量减少窗谱的最大旁瓣的相对幅度，也就是使能量尽量集中于主瓣，减小峰肩和纹波，进行增加阻带的衰减。 窗函数的选择原则是：

[1]具有较低的旁瓣幅度，尤其是第一旁瓣的幅度； [2]旁瓣的幅度下降的速率要快，以利于增加阻带的衰减； [3]主瓣的宽度要窄，这样可以得到比较窄的过渡带。

通常上述的几点难以同时满足。实际中设计FIR 数字滤波器往往要求是线性相位的，因此要求)(n w 满足线性相位的条件，即要求w(n)满足：

)1()(n N w n w --= （1-1）

所以，窗函数不仅有截短的作用，而且能够起到平滑的作用，在很多领域得到了应用。

窗口设计法基本步骤如下：

（1）构造希望逼近的频率响应函数()

ωj d e H 。以低通线性相位FIRDF 设计为例，一般选择()

ωj d e H 为线性理想低通滤波器，即

()

???=-，

，0ωτω

j j d e e H πωωωω≤<≤c c （1-2）

（2）求出()n h d 。对()

ωj d e H 进行IFT 得到

()()

()[]()

τπτωωπππωω--==?-n n d e e H n h c n

j j d d sin 21 （1-3）

（3）加窗得到FIRDF 的单位脉冲响应()n h ，

()()()n w n h n h d = （1-4）

式中，()n w 称为窗口函数，其长度为N 。如果要求第一类线性相位FIRDF ，则要求()n h 关于()2/1-N 点偶对称。而()n h d 关于τ=n 点偶对称，所

()2/1-=N τ，同时要求()n w 关于()2/1-N 点偶对称。

常见的窗函数，可以分为以下主要类型：

（1）幂窗--采用时间变量某种幂次的函数，如矩形、三角形、梯形或其

它时间（t ）的高次幂；

（2）三角函数窗--应用三角函数，即正弦或余弦函数等组合成复合函

数，例如汉宁窗、海明窗等；

（3）指数窗--采用指数时间函数，例如高斯窗等。

各种类型的窗函数可达到的阻带最小衰减和过渡带宽度见下表中所示；

表2-1常见窗函数性能表

名称滤波器

过渡带

宽

最小阻

带衰减

名称

滤波器

过渡带

宽

最小阻带

衰减

矩形 1.8π/M21dB PARZENWIN 6.6π/M 56dB

巴特利特 6.1π/M25dB FLATTOPWIN 19.6π/

M

108dB

汉宁 6.2π/M44dB GAUSSWIN 5.8π/M 60dB 汉明 6.6π/M51dB BARTHANNWIN 3.6π/M 40dB

布莱克曼11π/M74dB 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

1.5 汉宁窗

汉宁窗（Hanning Window）又称升余弦窗，汉宁窗可以看作是3个矩形时间窗的频谱之和，或者说是3个()t

sin型函数之和，而括号中的两项相对于第一个谱窗向左、右各移动了π/T，从而使旁瓣互相抵消，消去高频干扰和漏能。可以看出，汉宁窗主瓣加宽并降低，旁瓣则显著减小，从减小泄漏观点出发，汉宁窗优于矩形窗．但汉宁窗主瓣加宽，相当于分析带宽加宽，频率分辨力下降。

第二章 基本原理

2.1 汉宁窗基本原理

汉宁窗函数是余弦平方函数，又称之为升余弦函数，它的时域形式可以表为：

))1

2cos(1(5.0)(+-=n k

k w π

（2-1） 其中1,2,,k n =…。它的频域幅度特性函数为：

1

()222()0.5()0.25[()()]11N j R R R W W W W e

N N ω

ππωωωω--??=+-++??--??

（2-2）

其中()R W ω为矩形窗函数的幅度频率特性函数。汉宁窗函数的最大旁瓣值比主瓣值低31dB ，但是主瓣宽度比矩形窗函数的主瓣宽度增加了一倍，为

8N

π

。汉宁窗函数的时域幅度与频域幅度特性曲线的MATLAB 实现的曲线图如图2-1所示。

图2.1 汉宁窗函数的时域幅度与频域幅度特性曲线

2.2 窗函数法设计基本原理

用窗函数设FIR 滤波器的基本思路：从时域出发设计 h(n)逼近理想hd(n)。设理想滤波器的单位响应在时域表达为hd(n),则Hd(n) ,一般是无限长的，且是非因果的，不能直接作为FIR 滤波器的单位脉冲响应。要想得到一个因果的有限长的滤波器单位抽样响应h(n)最直接的方法是先将hd(n)往右平移，再进行截断，即截取为有限长因果序列：h(n)=hd(n)w(n)，并用合适的窗函数进行加权作为FIR 滤波器的单位脉冲响应。按照线性相位滤波器的要求，线性相位FIR数字低通滤波器的单位抽样响应h(n)必须是偶对称的。对称中心必须等于滤波器的延时常数，即用矩形窗设计的FIR 低通滤波器，所设计滤波器的幅度函数在通带和阻带都呈现出振荡现象，最大波纹大约为幅度的9%，这个现象称为吉布斯（Gibbs）效应。为了消除吉布斯效应，一般采用其他类型的窗函数。MATLAB 设计 FIR 滤波器有多种方法和对应的函数。窗函数设计法不仅在数字滤波器的设计中占有重要的地位，同时可以用于功率谱的估计，从根本上讲，使用窗函数的目的就是消除由无限序列的截短而引起的Gibbs现象所带来的影响。

第三章设计实现框图

开始

剪切语音将语音格式转换为wav

加入噪声

对语音信号进行频谱分析

用汉宁窗设计FIR

频率响应

用FIR滤波器对信号滤波

画出音乐信号滤波前后波形并比较

结束

图3.1 程序设计流程图

第四章基于MATLAB的仿真程序及结果

4.1 读入语音数据并画出其时域与频域图

用音乐剪切工具剪切一段音乐信号，并以”.wav”格式保存，本设计剪切音乐片段为“1.wav”,然后将该片段读入MATLAB，画出该信号的时域和频域的波形图。再对该片段进行抽样，得到抽样后的预处理信号的时域和频域波形图。

读入音乐信号：

[s, fs,bits]=wavread('C:\Users\Administrator\Desktop\1.wav');

ss=s(:,1);

%sound(s,fs,bits)

作出该原信号的时域波形:

s1=s(17000:end);

figure(1);

subplot(321)

plot(s)

title('原始语音信号');

作出该原信号的频域波形:

S=fft(s);

subplot(322);

plot(abs(S));

title('原始语音信号频谱');

subplot(323)

plot(s1)

对信号进行抽样：

title('预处理截短语音信号');

wavwrite(s1,fs,'s1.wav');

S1=fft(s1);

subplot(324)

plot(abs(S1))

title('预处理语音信号频谱');

N=length(s);

fn=10^5;

t=0:1/fs:(N-1)/fs;

4.2 在MATLAB中对信号加入两个噪声

s2=ss'+0.22*sin(fn*3*pi*t)+0.22*sin(fn*1.5*pi*t)+0.22*sin(fn*2*pi* t);

ss2=s2+awgn(s2,15);

读入加噪后的音乐信号并作出加噪后的时域和频域波形图:

wavwrite(ss2,fs,'s2.wav');

subplot(325);

plot(ss2);

title('加噪后语音信号');

subplot(326);

S2=fft(ss2);

plot(abs(S2));

title('加噪后信号频谱');

sound(ss2,fs,bits);

图4.1 原始、预处理及加噪的信号的时域频域波形图截图中给出了原始语音信号的时域与频域的波形图，及第一和第二张图。中间两幅图给出了预处理的时域和频域的波形图。最后两幅图是加噪后的时域和频域的波形图。

4.3 滤波器设计基本程序：

wp=.2*pi;

ws=.3*pi;

wc=(wp+ws)/2;

N=ceil(3.1*2*pi/(ws-wp));

n=0:1:N-1;

r=(N-1)/2;

hn1=fir1(N-1,wc/pi,'low',hanning(N));

hd=ideal_hp1(wc,N);

w_han=(hanning(N))';

h=hd.*w_han;

figure(2)

subplot(211)

stem(n,hd)

title('理想单位脉冲响应hd(n)');

subplot(212)

stem(n,h);

title('实际单位脉冲响应hd(n)');

figure(3)

freqz(hn1);

title('滤波器幅频特性与相频特性');

图4.2理想单位脉冲与实际单位脉冲

该截图给出了理想单位脉冲与实际单位脉冲图。s3=conv(s2,hn1); wavwrite(s3,fs,'s3.wav');

S3=fft(s3);

sound(s3,fs,bit

图4.3 滤波器特性图

该图给出了滤波器参数的相关波形图。给出了滤波器的幅频和相频特性波形图。

滤波处理后的信号频谱图：

figure(4)

subplot(211);

plot(s3);

title('滤波器处理之后信号图');

subplot(212);

plot(abs(S3));

title('滤波器处理之后频谱')

图4.4 滤波后的时域与频域波形图

该图是滤波后的时域频域的波形图像，从图像中不难看出经过FIR滤波器处理后的图像在频谱上相对于加噪的波形图像要规则的多，而且从声音上分辨，基本去除了信号中加给的两种噪音信号，达到了滤波的效果。

第五章实验结论及总结

通过本次课程设计，我巩固了课堂中学习的理论知识，并能够用所学习的理论知识正确分析数字信号处理的基本问题，和解释数字信号处理的基本现象，掌握了用窗函数法设计FIR数字低通滤波器的正确方法和步骤，以及用MATLAB软件编写程序实现该滤波器的仿真，还学到了一些用所学知识解决实际问题的技巧。还体会到了与其他高级语言的程序设计相比，MATLAB环境下可以更方便、快捷地设计具有严格线性相位的FIR滤波器，节省大量的编程时间，提高编程效率，且参数的修改也十分方便，还可以进一步进行优化设计。另外，由于课本知识掌握的不牢固和MATLAB软件编程环境的不熟悉，在课程设计的过程中也遇到了一些难题，例如不知道用窗函数法设计FIR数字低通滤波器的详细方法和步骤，题目所给的参数不知道如何转换成低通滤波器的技术指标，还有此程序指令的MATLAB实现也让我们很头疼。最后通过认真看了几遍课本和其他有关书籍，网上查阅相关资料，才使得这些问题得以解决。

我们也都深深意识到了团队合作的重要性，任务的完成需要团队每个人的密切配合，合理的分工与合作，共同面对和解决设计中遇到的各种问题，才能值得此次课程设计逐渐趋于完善。

参考文献：

[1] 丁玉美. 数字信号处理[M]. 西安电子科技大学出版社，2003，3.

[2] 朱冰莲. 数字信号处理[M]. 电子工业出版社，2003，7.

[3] 程佩青.数字信号处理教程（第二版)[M].北京：清华大学出版社，2001.

[4] 韩纪庆，张磊，郑铁然.语音信号处理[M].北京：清华大学出版社，2004.

［5］ Sanjit K.Mitra 著，孙洪，余翔宇译，《数字信号处理实验指导书》，电子工业出版社，2005 年。

[6] 张葛祥，李娜. MATLAB仿真技术与应用.北京：清华大学出版社，2002年。

[7] 董长虹等. MATLAB信号处理与应用.北京：国防工业出版社，2005。

附录：程序源代码

[s, fs,bits]=wavread('C:\Users\Administrator\Desktop\1.wav');

ss=s(:,1);

%sound(s,fs,bits)

s1=s(17000:end);

figure(1);

subplot(321)

plot(s)

title('原始语音信号');

S=fft(s);

subplot(322);

plot(abs(S));

title('原始语音信号频谱');

subplot(323)

plot(s1)

title('预处理截短语音信号');

wavwrite(s1,fs,'s1.wav');

S1=fft(s1);

subplot(324)

plot(abs(S1))

title('预处理语音信号频谱');

N=length(s);

fn=10^5;

t=0:1/fs:(N-1)/fs;

s2=ss'+

0.22*sin(fn*3*pi*t)+0.22*sin(fn*1.5*pi*t)+0.22*sin(fn*2*pi*t);

ss2=s2+awgn(s2,15);

wavwrite(ss2,fs,'s2.wav');

subplot(325);

plot(ss2);

title('加噪后语音信号');

S2=fft(ss2);

plot(abs(S2));

title('加噪后信号频谱');

sound(ss2,fs,bits);

wp=.2*pi;

ws=.3*pi;

wc=(wp+ws)/2;

N=ceil(3.1*2*pi/(ws-wp));

n=0:1:N-1;

r=(N-1)/2;

hn1=fir1(N-1,wc/pi,'low',hanning(N)); hd=ideal_hp1(wc,N);

w_han=(hanning(N))';

h=hd.*w_han;

figure(2)

subplot(211)

stem(n,hd)

title('理想单位脉冲响应hd(n)'); subplot(212)

stem(n,h);

title('实际单位脉冲响应hd(n)'); figure(3)

freqz(hn1);

title('滤波器幅频特性与相频特性');

s3=conv(s2,hn1);

wavwrite(s3,fs,'s3.wav');

S3=fft(s3);

sound(s3,fs,bits);

figure(4)