IIR滤波器的设计论文

IIR滤波器的设计

目 录

目 录 ............................... 错误！未定义书签。 一、简 介 ........................... 错误！未定义书签。

1.1 引 言 ......................................... 3 1.2 课程设计基本要求 .............................................................. 3 1.3 课程设计基本内容 .............................................................. 4 二、模拟滤波器的设计 .................................. 4

2.1 模拟滤波器的设计方法 .......................... 4 2.2 模拟原型滤波器级最小阶数的选择 ................ 6

2.2.1 巴特沃斯滤波器及最小阶数的选择 ........... 6 2.2.2低通原型滤波器的系统函数 ................ 11 2.2.3 椭圆滤波器及最小阶数的选择 .............. 12 2.2.4贝塞尔滤波器 .......................... 13

三、数字滤波器设计的总体方案 ......................... 13

3.1 数字流程图 ................................... 13 3.2 语音信号采集 ................................. 14 3.3 语音信号分析 .................................................................... 15 3.4 含噪语音信号合成 ............................................................ 15 3.5数字滤波器哦原理及步骤 ........................ 16

3.5.1相关原理 ............................... 16 3.5.2用窗函数数字滤波器设计的基本步骤 ........ 17

四、窗函数的设定 ..................................... 18

4.1窗函数的定义 .................................. 18 4.2 窗函数的分类 ................................. 19 4.3几种常用窗函数的性质和特点 .................... 19 五、利用窗函数设计给定参数的FIR滤波器 ............... 21

5.1指标的确定 .................................... 21 5.2凯塞窗函数的程序的设定 ................................................. 21 六、 滤波分析 ........................................ 22

6.1通过回放语音信号分析 .......................... 22 6.2通过信号滤波前后的波形图分析 .................. 22 七、总结 ............................................. 24

1

/ 26

一、简介

1.1 引言

数字信号处理的主要研究对象是数字信号，且是采用运算的方法达到处

理的目的的，因此，其实现方法，基本上分成两种实现方法，即软件和硬件实现方法。软件实现方法指的是按照原理和算法，自己编写程序或者采用现成的程序在通用计算机上实现，硬件实现指的是按照具体的要求和算法，设计硬件结构图，用乘法器加法器延时器、控制器、存储器以及输入输出接口部件实现的一种方法。显然前者灵活，只要改变程序中的有关参数，但是运算速度慢，一般达不到实时处理，因此，这种方法适合于科研和教学。后者运算速度快，可以达到实时处理要求，但是不灵活。目前DSP芯片已进入市场，且正在高速发展，速度高，体积小，性能优良，价格也在不断下降。可以说，用DSP芯片实现数字信号处理，正在变成工程技术领域的主要方法。用合适的DSP芯片，配有合适的芯片语言及任务要求的软件，来实现信号处理功能无疑是一种最佳的数字信号处理系统。本文仅使用软件实现。

MATLAB是由美国mathworks公司发布的主要面对科学计算、可视化

以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。

MATLAB 的应用范围非常广，包括信号和图像处理、通讯、控制系统

设计、测试和测量、财务建模和分析以及计算生物学等众多应用领域。附加的工具箱（单独提供的专用 MATLAB 函数集）扩展了 MATLAB 环境，以解决这些应用领域内特定类型的问题。 1.2课程设计基本要求

1. 熟悉离散信号和系统的时域特性。

2

/ 26

2. 掌握数字信号处理的基本概念，基本理论和基本方法。 3. 掌握序列快速傅里叶变换方法。

4. 学会MATLAB的使用，掌握MATLAB的程序设计方法。 5. 掌握利用MATLAB对语音信号进行频谱分析。

6. 掌握MATLAB设计FIR数字滤波器的方法和对信号进行滤波的方 1.3课程设计内容

1.利用Windows下的录音机或其他软件，录制一段自己的语音信号，时间控制在1s左右，并对录制的信号进行采样；

2.语音信号的频谱分析，画出采样后语音信号的时域波形和频谱图； 3.产生噪声信号并加到语音信号中，得到被污染的语音信号，并回放语音信号；

4.污染信号的频谱分析，画出被污染的语音信号时域波形和频谱； 5.根据有关的频谱特性，设计FIR数字滤波器； 6.用自己设计的滤波器对被污染的语音信号进行滤波；

7.分析得到信号的频谱，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化； 8.回放语音信号。

二、数字滤波器的设计的总体方案

2.1模拟滤波器的设计

模拟滤波器的理论和设计方法已经发展的相当成熟，且有若干典型的模拟滤波器供我们选择，如巴特沃斯（Butterworth滤波器.切比雪夫（Chebyshev）滤波器等。这些工作的理论分析和设计方法在20世纪30年代就完成，然而烦琐.冗长的数字计算使它难以付诸实用。直到50年代，由于计算机技术的逐步成熟，求出大量设计参数和图表，这种方法才得到广泛应用。这些典型的滤波器各有特点：巴特沃斯滤波器具有单调下降的幅频特性；切比雪夫滤波器的幅频特性在通带或者阻带有波动发，可以提高选择性。这样根据具体要求可以选择不同类型的滤波器。

模拟滤波器按幅度特征可以分成低通、高通、带通和带阻滤波器。它们的理想幅度特性如图2.1所示，但我们设计滤波器时，总是先设计低通滤波器，再通过频率变换将低通滤波器转换成希望类型的滤波器

3

/ 26

低通 高通 带通

图2.1 模拟滤波器理想幅度特性

带阻

利用频率变换设计模拟滤波器的步骤为：

（1）给定模拟滤波器的性能指标，如截止频率0或上、下边界频率1,2等。 （2）确定滤波器阶数

（3）设计模拟低通原型滤波器。

（4）按频率变换设计模拟滤波器（低通、高通、带通、带阻）。

模拟低通滤波器的设计指标有p，p和s，其中p和s分别称为通带截止频率和阻带截止频率。p 是通带Ω（=0—p）中的最大衰减系数，s是阻带Ω≥s的最小衰减系数，p和s一般用dB表示。对于单调下降的幅度特性，可表示成：

p10lgHa(j0)22 （2.1）

Ha(jp))Ha(j0)2

p10lgHa(js))2 （2.2）

如果Ω=0处幅度已归一化为一，即Haj1，p和s表示为

p10lgHa(jp) （2.3） s10lgHa(js) （2.4）

以上技术指标用图2.2表示，图中c称为3dB 截止频率，因

4

/ 26

22

Hajc12

,-20Hajc3dB

1 0.707 

 jΩ

图2.2 低通滤波器的幅度特性

滤波器的技术指标给定以后，需要设计一个传输函数Has，希望其幅度平方函数

Haj满足给定的指标p和s，一般滤波器的单位冲激响应为实数，因此

Ha(j)Ha(a)Ha(s)|sj

 =Ha(j)Ha(j) (2.5)

2如果能由p，p，s，s求出Haj，那么就可以求出所需的Has，对于上面

2介绍的典型滤波器，其幅度平方函数有自己的表达式，可以直接引用。这里要说明的是HasHaj必须是稳定的。因此极点必须落在s平面的左半平面，相应的Has的极点落在右半平面。

2.2 模拟原型滤波器及最小阶数的选择

2.2.1 巴特沃斯滤波器及最小阶数的选择

巴特沃斯滤波器是最基本的逼近方法形式之一。它的幅频特性模平方为

2Ha(j)(1(12N)c)2 (2.6)

式中N是滤波器的阶数。当Ω=0时，Haj1；当Ω=c时，Haj5

12，

/ 26

c是3dB截止频率。

不同阶数N的巴特沃斯滤波器特性如图2.3所示，这一幅频特性具有下列特点： （1）最大平坦性：可以证明：在Ω=0点，它的前（2N-1）阶导数都等于0，这表明巴特沃斯滤波器在Ω=0附近一段范围内是非常平直的，它以原点的最大平坦性来逼近理想低通滤波器。“最平响应”即由此而来。

（2）通带，阻带下降的单调性。这种滤波器具有良好的相频特性。

（3）3dB的不变性：随着N的增加，频带边缘下降越陡峭，越接近理想特性，但不管N是多少，幅频特性都通过-3dB点。当Ω≥c时，特性以20NdB/dec速度下降。

图2.3 不同阶数N的巴特沃斯滤波器特性

现根据式（2.6）求巴特沃斯滤波器的系统函数Ha（s）。令Ω=s/j，带入式（2.6）

Ha(j)2sj(jc)2N1 Ha(s)Ha(s)2N2Ns2Ns(jc)1()jc2N对应的极点：s2Njc0

22N

skjc1ce1j2k122N (2.7)

sk即为HasHas的极点，此极点分布有下列特点：

(1)HasHas的2N个极点以π/N为间隔均匀分布在半径为c的圆周上，这个圆称为巴特沃斯圆。

(2)所有极点以jΩ轴为对称轴成对称分布，jΩ轴上没有极点。

(3)当N为奇数时，有两个极点分布在sc的实轴上；N为偶函数时，实轴上没有

6

/ 26

极点。所有复数极点两两呈共轭对称分布。图2.4画出了N=3时的HasHas极点分布。全部零点位于s=∞处。

jΩ  c cδ

图2.4 N=3时Ha（s）Ha（-s）极点分布

为得到稳定的Has，取全部左半平面的极点。

NcHassskk1N （2.8）

当N为偶数时

HasNcN2NcN2 （2.9）

22k12sssss2cosskkcc22Nk1k1当N为奇数时

HasNc22k12s2cossccsc2N2k1为使用方便把式（2.9）和式（2.10）对c进行归一化处理，为此，分子分母各除以c，并令sNN12 （2.10）

s，s称为归一化复频率： cHas1（N为偶数） （2.11）

N22k1s2coss122Nk17

/ 26

Has1 （N为奇数）（2.12）

sk1N1222k12coss1s122N用归一化频率/c表示的频率特性称为原型滤波特性（Ωَ即归一化复频率sَ 的虚部）。对式（2.6）所示的低通巴特沃斯特性用Ωَ表示得到：

Haj212N1 （2.13）

称Haj为巴特沃斯低通原型滤波器幅频特性。在低通原型滤波频率特性上，截止频率

c=1。

若给出模拟低通滤波器的设计性能指标要求：通带边界频率p，阻带边界频率s，通带波纹Rp(dB)，阻带衰减RS(dB)，要确定butterworth ,，低通滤波器最小阶数N及截止频率C(3dB)。p，S，RS，RP的意义如图所示。 当

=P时，H(j)10RP20 即RP10lg H(j)，以截至频率c（幅值下降3dB）

2为1，化为相对为相对的相对频率由上式可写为1c10lg10RPcp2N1()c。 同理，当=c时，H(j)10RS201，

10lg10RPS)2N1(C 。 (10由此可见NcS2NRP/101)(10RP/101)p2log)10(s N

应向上取整cp2N(10RP/101)，

(10RS/101)再用MATLAB 编程计算滤波器最小阶数N和截止频率c。

sk就是切比雪夫滤波器HasHas的极点，给定N，c，ε即可求的2N个极点分

8

/ 26

布。由式（2.22）实部与虚部的正弦和余弦函数平方约束关系可以看出，此极点分布满足椭圆方程，其短轴和长轴分别为

11asinharcsinhc （2.23） N11bcosharcsinhcN图2.7画出了N=3时切比雪夫滤波器的极点分布。

jΩ a b ζ

图2.7

极点所在的椭圆可以和半径为a的圆和半径为b的圆联系起来，这两个圆分别称为巴特沃斯小圆和巴特沃斯大圆。N阶切比雪夫滤波器极点的纵坐标，而横坐标等于N阶巴特沃斯小圆极点的横坐标取左半平面的极点：

2k1asink2N k=1,2„，N （2.24） 2k1bcosk2n则切比雪夫滤波器的系统函数：

HasA(ssk1N （2.25）

k)cNjc其中skkjk，常数A=。因而切比雪夫滤波器的系统函数表示为： N12HasNc/2N1sskk1N （2.27）

9

/ 26

切比雪夫滤波器的截止角频率c不是像巴特沃斯滤波器中所规定的（-3dB）处角频率，而是通带边缘的频率。若波纹参数满足

10.5，可以求的-3dB处的角频率为 21113dBccosharcosh （2.28）

N将式（2.27）表示的Has对c归一化，得到切比雪夫I型 2.2.2低通原型滤波器的系统函数

HassN12N1 （2.29） N1aN1sa1sa0对不同的N，式（2.29）的分母多项式已制成表格，供设计参考。

和butterworth低通模拟滤波器设计一样，若给定性能指标要求：c，s，RP，Rs确定Chebyshev低通模拟滤波器最小阶数N和截止频率c（-3dB频率）。

2.2.2.1 Chbbyshev I型 由式HajA()221RS/20RP/10A10可得 故101221CN(/C) 式中，g阶数N可由下式求得Nlog10(gg21)log10(s(s)21))pp(A21)/2，截至

频率cp由上面两式用Matlab 编程计算滤波器最小阶数N和截止频率

2.2.2.2 Chbbyshev II型

c

Chbbyshev II型通带内是平滑的，而阻带具有等波纹起伏特性。因此，在阶数N的计算公式上是相同的，而-3dB截止频率c则不同。

2.2.3 椭圆滤波器及最小阶数的选择

椭圆的模拟低通滤波器圆形的平方幅值响应函数为

10

/ 26

H(j)2A(2)112E2N(/C)

式中，为小于1的正书，表示波纹情况；c为截止频率；En(/c）为椭圆函数，定义为

222K当N为偶数（N=2m）时,EN()1k1mmK22

当N为奇数（N=2m+1）时,EN()1k1K222K 其中/c

椭圆模拟滤波器特点是：在通带和阻带内均具有等波纹起伏特性。何以上滤波器相比，相同的性能指标所需要的阶数最小。但频率响应应具有明显的非线性。由式

HajA(2)21 221EN(/C)RP/10滤波器的阶数可由下式确定101，A10RS/20 NK(k)K(1K12K(k1)K(1K)2，

/2pd k1 K(x)由上式计算滤波器的np 式中k222sA11xsin0最小阶数N和截止频率C。

2.2.4贝塞尔滤波器

贝塞尔模拟低通滤波器原型的特点是在零频时具有最平坦的群延迟，并在整个通带内延

(2N)!迟几乎不变。在零频时的群延迟为NN!21/N。由于这一特点，贝塞尔模拟滤波器通带内保

持信号形状不变。滤波器传递汉书具有下面形式

]H(S)k

(sp(1))(sp(2))(sp(n))11

/ 26

三、数字滤波器设计的方案 3.1流程框图

方案一：

在一个相对较安静的环境下，录下1s左右的wav声音信号，然后对声音进行采样，画出其时域波形和频谱图，利用程序编一个噪声信号加载在原声音信号里面，将这个被污染的语音信号通过baterworth带通滤波器，将滤波后的信号进行抽样再和原始信号进行比较，其流程图如下所示：

开 始 噪 声 录 音 频 域 波 形 噪声+录音 时 域 波 形

滤 波 （baterworth） over 方案二：

12

/ 26

我们知道在我们的现实生活中噪声信号是无处不在的，所以，声音信号中本来就存在着噪声，我们只需设计一个合适的滤波器将原始声音中的杂音滤除就可以了，其流程图如下所示：

添加噪声 录制语音信号并读取该音频信号 语音信号采样 对语音信号进行频谱分析

设计滤波器

对语音信号进行滤波

3.2、语音信号采集 输出各部分频谱图 录制一段课程设计学生的语音信号并保存为文件，要求长度控制在1秒，并对录制的信号进行采样；录制时可以使用Windows自带的录音机，或者使用其它专业的录音软件，录制时需要配备录音硬件（如麦克风），为便于比较，需要在安静、干扰小的环境下录音。

3.3、语音信号分析

通过用windows录音之后，将录音的文件导入到MATLAB，并使用MATLAB绘出采样后的语音信号的时域波形和频谱图。根据频谱图求出其带宽，并说明语音信号的采样频率不能低于多少赫兹。

语音信号采集的程序： fs=8000,bits=8,T=1,Ts=1/fs; N=T/Ts;

[x,fs,bits]=wavread('E:\\MATLAB7\\work\\yuyin.wav',[1024 6023]); x=x(:,1);

subplot(321);plot(x);

13

/ 26

sound(x,fs,bits);

title('语音信号时域波形图') y=fft(x,1024);

n=(fs/1024)*[1:1024]; subplot(322);

plot(n(1:512),abs(y(1:512))); title('语音信号频谱图');

3.4、含噪语音信号合成

在MATLAB软件平台下，给原始的语音信号叠加上噪声，噪声类型分为如下几种：（1）白噪声；（2）单频噪色（正弦干扰）；（3）多频噪声（多正弦干扰）；（4）其它干扰，可设置为低频、高频、带限噪声，或Chirp干扰、冲激干扰。绘出叠加噪声后的语音信号时域和频谱图，在视觉上与原始语音信号图形对比，也可通过Windows播放软件从听觉上进行对比，分析并体会含噪语音信号频谱和时域波形的改变。这里，我们加的是单频的正弦干扰。

噪声信号与语音信号合成的程序： dt=0.1*sin(2*pi*5000*(1:size(x))/fs); % dt1=fft(dt,1024);

% plot(n(1:512),abs(dt1(1:512))); x1=x+dt'; sound(x1,fs,bits); y1=fft(x1,1024); subplot(323);plot(x1); title('含噪声信号波形'); subplot(324);

plot(n(1:512),abs(y1(1:512))); title('含噪声信号频谱');

3.5、数字滤波器原理及步骤

给定滤波器的规一化性能指标（参考指标,实际中依据每个同学所叠加噪声情况而定）例如：通带截止频率wp=0.25*pi, 阻通带截止频率ws=0.3*pi; 通带最

14

/ 26

大衰减Rp=1 dB; 阻带最小衰减Rs=15 dB，可以采用窗函数法与等波纹法分别设计各型FIR滤波器（低通、高通、带通、带阻中的至少3种类型）来对叠加噪声前后的语音信号进行滤波处理，绘出滤波器的频域响应，绘出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；在相同的性能指标下比较各方法的滤波效果，并从理论上进行分析（或解释）。

3.5.1相关原理：

设计数字滤波器的任务就是寻求一个因果稳定的线性时不变系统，并使系统函数H(z)具有指定的频率特性。 数字滤波器从实现的网络结构或者从单位冲激响应分类，可以分成无限长单位冲激响应（IIR）数字滤波器和有限长单位冲激响应（FIR）数字滤波器。 数字滤波器频率响应的三个参数： （1） 幅度平方响应： （2） 相位响应 （3） 群时延响应

FIR数字滤波器：

FIR滤波器通常采用窗函数方法来设计。窗设计的基本思想是，首先选择一个适当的理想选频滤波器（它总是具有一个非因果，无限持续时间脉冲响应），然后截取（加窗）它的脉冲响应得到线性相位和因果FIR滤波器。因此这种方法

jH(e)表示理想的选的重点是选择一个合适的窗函数和理想滤波器。我我们用d频滤波器，它在通带上具有单位增益和线性相位，在阻带上具有零响应。一个带宽c的低通滤波器由下式给定：

1ej,||cHd(e)0,c||

j15

/ 26

为了从hd(n)得到一个FIR滤波器，必须同时在两边截取hd(n)。而要得到一个因果的线性相位滤波器，它的h(n)长度为N，必须有：

h(n),0nN1h(n)d0,其他

N12

这种操作叫做加窗，h(n)可以看作是hd(n)与窗函数w(n)的乘积：

h(n)hd(n)w(n)

其中

0nN1关于对称，w(n)0,其他

根据w(n)的不同定义，可以得到不同的窗结构。

jjjH(e)H(e)W(e)的周期卷d在频域中，因果FIR滤波器响应由和窗响应

积得到，即：

1H(e)Hd(e)W(e)2jjjjj()W(e)H(e)dd

常用的窗函数有矩形窗、巴特利特（BARTLETT）窗、汉宁（HANNING）窗、海明（HAMMING）窗、布莱克曼（BLACKMAN）窗、凯泽（KAISER）窗等。

3.5.2 用窗函数数字滤波器设计的基本步骤

（1）确定技术指标

在设计一个滤波器之前，必须首先根据工程实际的需要确定滤波器的技术指标。在很多实际应用中，数字滤波器常常被用来实现选频操作。因此，指标的形

16

/ 26

式一般在频域中给出幅度和相位响应。幅度指标主要以两种方式给出。第一种是绝对指标。它提供对幅度响应函数的要求，一般应用于FIR滤波器的设计。第二种指标是相对指标。它以分贝值的形式给出要求。在工程实际中，这种指标最受欢迎。对于相位响应指标形式，通常希望系统在通频带中人有线性相位。运用线性相位响应指标进行滤波器设计具有如下优点：①只包含实数算法，不涉及复数运算；②不存在延迟失真，只有固定数量的延迟；③长度为N的滤波器（阶数为N-1），计算量为N/2数量级。因此，本文中滤波器的设计就以线性相位FIR滤波器的设计为例。 （2）逼近

确定了技术指标后，就可以建立一个目标的数字滤波器模型。通常采用理想的数字滤波器模型。之后，利用数字滤波器的设计方法，设计出一个实际滤波器模型来逼近给定的目标。 （3）性能分析和计算机仿真

上两步的结果是得到以差分或系统函数或冲激响应描述的滤波器。根据这个描述就可以分析其频率特性和相位特性，以验证设计结果是否满足指标要求；或者利用计算机仿真实现设计的滤波器，再分析滤波结果来判断。 数字滤波 器技术指 标 指标参数变 换 模拟滤波器技术指 标 相应的模拟滤波器设计 模拟滤波器离散 化 数字滤波器

图1 数字滤波器设计步骤

四、 窗函数的设定

4.1 窗函数的定义

17

/ 26

数字信号处理的主要数学工具是博里叶变换．而傅里叶变换是研究整个时间域和频率域的关系。不过，当运用计算机实现工程测试信号处理时，不可能对无限长的信号进行测量和运算，而是取其有限的时间片段进行分析。做法是从信号中截取一个时间片段，然后用观察的信号时间片段进行周期延拓处理，得到虚拟的无限长的信号，然后就可以对信号进行傅里叶变换、相关分析等数学处理。 无线长的信号被截断以后，其频谱发生了畸变，原来集中在f0处的能量被分散到两个较宽的频带中去了，这种现象称之为频谱能量泄漏。信号截断以后产生的能量泄漏现象是必然的，因为窗函数w(t)是一个频带无限的函数，所以即使原信号x(t)是限带宽信号，而在截断以后也必然成为无限带宽的函数，即信号在频域的能量与分布被扩展了。又从采样定理可知，无论采样频率多高，只要信号一经截断，就不可避免地引起混叠，因此信号截断必然导致一些误差。 为了减少频谱能量泄漏，可采用不同的截取函数对信号进行截断，截断函数称为窗函数，简称为窗。泄漏与窗函数频谱的两侧旁瓣有关，如果两侧瓣的高度趋于零，而使能量相对集中在主瓣，就可以较为接近于真实的频谱，为此，在时间域中可采用不同的窗函数来截断信号。

4.2实际应用的窗函数，可分为以下主要类型：

4.2.1幂窗

采用时间变量某种幂次的函数，如矩形、三角形、梯形或其它时间（t）的高次幂；

4.2.2 三角函数窗

应用三角函数，即正弦或余弦函数等组合成复合函数，例如汉宁窗、海明窗等；

4.2.3 指数窗

采用指数时间函数，如 形式，例如高斯窗等。

4.3几种常用窗函数的性质和特点。

4.3.1 矩形窗

矩形窗属于时间变量的零次幂窗。矩形窗使用最多，习惯上不加窗就是使信号通过了矩形窗。这种窗的优点是主瓣比较集中，缺点是旁瓣较高，并有负旁瓣，导致变换中带进了高频干扰和泄漏，甚至出现负谱现象

18

/ 26

4.3.2 三角窗

三角窗亦称费杰（Fejer）窗，是幂窗的一次方形式。与矩形窗比较，主瓣，宽约等于矩形窗的两倍，但旁瓣小，而且无负旁瓣。 4.3.3汉宁（Hanning）窗

汉宁窗又称升余弦窗，汉宁窗可以看作是3个矩形时间窗的频谱之和，或者说是 3个 sine（t）型函数之和，而括号中的两项相对于第一个谱窗向左、右各移动了 π/T，从而使旁瓣互相抵消，消去高频干扰和漏能。可以看出，汉宁窗主瓣加宽并降低，旁瓣则显著减小，从减小泄漏观点出发，汉宁窗优于矩形窗．但汉宁窗主瓣加宽，相当于分析带宽加宽，频率分辨力下降。 4.3.4 海明（Hamming）窗

海明窗也是余弦窗的一种，又称改进的升余弦窗。海明窗与汉宁窗都是余弦窗，只是加权系数不同。海明窗加权的系数能使旁瓣达到更小。分析表明，海明窗的第一旁瓣衰减为一42dB．海明窗的频谱也是由3个矩形时窗的频谱合成，但其旁瓣衰减速度为20dB／（10oct），这比汉宁窗衰减速度慢。海明窗与汉宁窗都是很有用的窗函数。 4.3.5高斯窗

高斯窗是一种指数窗。高斯窗谱无负的旁瓣，第一旁瓣衰减达一55dB。高斯富谱的主瓣较宽，故而频率分辨力低．高斯窗函数常被用来截断一些非周期信号，如指数衰减信号等。

不同的窗函数对信号频谱的影响是不一样的，这主要是因为不同的窗函数，产生泄漏的大小不一样，频率分辨能力也不一样。信号的截断产生了能量泄漏，而用FFT算法计算频谱又产生了栅栏效应，从原理上讲这两种误差都是不能消除的，但是我们可以通过选择不同的窗函数对它们的影响进行抑制。图6.5是几种常用的窗函数的时域和频域波形，其中矩形窗主瓣窄，旁瓣大，频率识别精度最高，幅值识别精度最低；布莱克曼窗主瓣宽，旁瓣小，频率识别精度最低，但幅值识别精度最高。

对于窗函数的选择，应考虑被分析信号的性质与处理要求。如果仅要求精确读出主瓣频率，而不考虑幅值精度，则可选用主瓣宽度比较窄而便于分辨的矩形窗，

19

/ 26

例如测量物体的自振频率等；如果分析窄带信号，且有较强的干扰噪声，则应选用旁瓣幅度小的窗函数，如汉宁窗、三角窗等；对于随时间按指数衰减的函数，可采用指数窗来提高信噪比。 五、 利用窗函数设计给定参数的

FIR滤波器

5.1指标的确定

对于不同类型的滤波器，参数wp和ws有一些限制：对于低通滤波器，wpws； 首先，根据所录得音跟所加正弦噪声合成后的频谱观察，确定wp和ws的大致范围，设定好了wp和ws之后，变可以确定过度带的大小了，根据过度带的及阻带衰减的指标要求，来选择窗函数的类型。根据前面介绍的录音噪声信号跟源信号合成后的频谱图可以看出噪声频率大概在3KHZ左右，所以就可以设定fp,fs从而确定wp,ws。接着，按照阻带衰减选择窗函数的类型。

5.2 凯塞窗函数的程序的设定

滤波器的程序：

fp=2800,fc=3600,rs=90; %设置通带频率和阻带频率，单位是Hz 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; %理想滤波器通带截止频率

20

/ 26

hn=fir1(M,wc,kaiser(M+1,alph)); %理想低通滤波器脉冲响应并形成了凯塞函数

X=conv(hn,x); %过滤后的信号 sound(X,fs,bits); X1=fft(X,1024)

subplot(325);plot(X);title('处理后的信号波形');

subplot(326);plot(n(1:512),abs(X1(1:512)));; title('处理后的信号频谱')

六、滤波分析

6.1 通过回放语音信号分析

从以上种音乐滤波去噪设计，我们可以从信号滤波前后的波形图以及频谱图上看出变化。当然，也可以用sound()函数来播放滤波后的音乐，从听觉上直接感受音乐信号的变化，本次设计比较明显，滤波后的音乐比滤波前要清晰明亮很多，当然由于人耳听力的限制，有些情况下我们是很难听出异同的。

sound(x1,fs,bits); %回放含噪音乐 sound(X,fs,bits); %回放去噪声后音乐

6.2 通过信号滤波前后的波形图分析

滤波前后信号在时间域和频率域上的对照图如下所示： 6.2.1图5-1到5-3原信号以及滤波前后信号在时间域的对照图 信号时域图形如下所示：

21

/ 26

IIR数字带通滤波器幅频—相频特性

滤波前后信号波形对比如下：

滤波前信号波形

滤波前加噪声信号波形

22

/ 26

滤波后信号波形

滤波前后信号频谱对比如下：

滤波前信号频谱

滤波前加噪声信号频谱

滤波后信号频谱

七、总结

在本次课程设计的过程中，主要遇到的问题是原始音乐信号的选取。有些音乐信号在MATLAB中运行时出错，开始以为自己的操作问题，在跟同学商讨多次更换音乐信号后得到了解决。其次就是在滤波器参数的选取上，开始时不清楚怎么选取，又在课本中找了好久才找到。

学习的过程是相互讨论共同进步的，多多讨论课题中遇到的问题，可以巩固我们的知识掌握能力，增加熟练运用度。

从本次课程设计的中心来看，课题是希望将数字信号处理技术应用于某一实际领域，这里就是指对音乐的处理。作为存储于计算机中的音乐信号，其本身就是离散化了的向量，我们只需将这些离散的量提取出来，就可以对其进行处理了。

23

/ 26

在这里，用到了处理数字信号的强有力工具MATLAB，通过MATLAB里几个命令函数的调用，很轻易的在实际化音乐与数字信号的理论之间搭了一座桥。课题的特色在于它将音乐看作了一个向量，于是音乐数字化了，则可以完全利用数字信号处理的知识来解决。我们可以像给一般信号做频谱分析一样，来给音乐信号做频谱分析，也可以较容易的用数字滤波器来对音乐进行滤波处理。改变参数，理论结合实际，分析各参数对图形的影响，从而加深对各个参数的理解。在完成这次课程设计过程中学到了许多东西，进一步理解了滤波器设计方法和各参数意义，通过分析信号时域和频域的关系等，加深了对滤波性能的理解，而且学会了使用Matlab一些基本函数，增加了进一步学习Matlab软件的兴趣。同时，通过本次课程设计，锻炼了我的动手能力，和提高了我分析问题，解决问题的能力。

附录

源程序：

fs=8000,bits=8,T=1,Ts=1/fs; N=T/Ts;

[x,fs,bits]=wavread('E:\\MATLAB7\\work\\shengyin.wav'); x=x(:,1);

subplot(321);plot(x); sound(x,fs,bits); title('时域波形图') y=fft(x,1024);

n=(fs/1024)*[1:1024]; subplot(322);

plot(n(1:512),abs(y(1:512))); title('频谱图');

dt=0.01*sin(2*pi*7000*(1:size(x))/fs); dt1=fft(dt,1024);

% plot(n(1:512),abs(dt1(1:512)));

24

/ 26

x1=x+dt'; sound(x1,fs,bits); y1=fft(x1,1024); subplot(323); plot(x1);

title('污染信号波形'); subplot(324);

plot(n(1:1024),abs(y1(1:1024))); axis([0,12000,0,0.4]); title('污染信号频谱'); fp=5000,fc=6600,rs=100; 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,x); sound(X,fs,bits); X1=fft(X,1024) subplot(325);

plot(X);title('处理后的信号波形'); subplot(326);

plot(n(1:512),abs(X1(1:512)));; title('处理后的信号频谱');

运行结果：

25

/ 26

时域波形图 频谱图

污染信号波形 污染信号频谱

处理后的信号波形 处理后的信号频谱

26

/ 26