音乐信号滤波去噪用凯塞窗设计的FIR滤波器.docx

音乐信号滤波去噪用凯塞窗设计的FIR滤波器.docx

《音乐信号滤波去噪用凯塞窗设计的FIR滤波器.docx》由会员分享，可在线阅读，更多相关《音乐信号滤波去噪用凯塞窗设计的FIR滤波器.docx（18页珍藏版）》请在冰豆网上搜索。

音乐信号滤波去噪用凯塞窗设计的FIR滤波器

音乐信号滤波去噪—用凯塞窗设计的FIR滤波器

学生姓名：

指导老师：

摘要本课程设计主要是通过使用凯塞窗设计一个FIR滤波器以对音乐信号进行滤波去噪处理。

本设计首先通过麦克风采集一段音乐信号，依据对该信号的频谱分析，给定相关指标。

以MATLAB软件为平台，采用凯塞窗设计满足指标的FIR滤波器，以该音乐信号进行滤波去噪处理。

通过对比滤波前后的波形图，深入了解滤波器的基本方法。

通过程序调试及完善，该设计基本满足设计要求。

关键词滤波去噪；FIR滤波器；凯塞窗函数；MATLAB

1引言

数字滤波器是一种用来过滤时间离散信号的数字系统,它是通过对抽样数据进行数学处理来达到频域滤波的目的。

随着现代通信的数字化,数字滤波器变得更加重要。

数字滤波器的种类很多，但总的来说可以分成两大类，一类是经典滤波器，另一类可称为现代滤波器。

从滤波特性方面考虑，数字滤波器可分成数字高通、数字低通、数字带通和数字带阻等滤波器。

从实现方法上考虑，将滤波器分成两种，一种称为无限脉冲响应滤波器，简称IIR（InfiniteImpulseResponse）滤波器，另一种称为FIR（FiniteImpulseResponse）滤波器[1]。

设计FIR数字滤波器的方法有窗函数法、频率采样法和等波纹最佳逼近法等。

1.1课程设计目的

数字信号处理（DigitalSignalProcessing，简称DSP）是一门涉及许多学科而又广泛应用于许多领域的新兴学科。

20世纪60年代以来，随着计算机和信息技术的飞速发展，数字信号处理技术应运而生并得到迅速的发展。

在过去的二十多年时间里，数字信号处理已经在通信等领域得到极为广泛的应用。

数字信号处理是利用计算机或专用处理设备，以数字形式对信号进行采集、变换、滤波、估值、增强、压缩、识别等处理，以得到符合人们需要的信号形式。

在本次课程设计中，最主要的设计是设计FIR滤波器，FIR滤波器的设计方法主要分为两类：

第一类是基于逼近理想滤波器器特性的方法包括窗函数法、频率采样法、和等波纹最佳逼近法；第二类是最优设计法。

本次的课程设计主要采用的是第一类设计方法，是利用汉宁窗函数法设计FIR滤波器对一段语音进行滤波去噪，通过这一过程，对滤波前后波形进行对比分析得到结论。

此课程设计比较简单，主要是将书本中的知识运用到现实中，并且根据自己对设计题目的理解，运用软件编写出程序实现这一设计，也是我们对数字信号处理的原理进行验证的一个过程。

对此，也可以加深我们对所学知识的理解，培养我们的动手能力。

1.2课程设计的要求

（1）通过利用各种不同的开发工具实现模拟信号数字化、信道编解码、基带数字信号编解码、数字信号的调制解调和语音信号的滤波去噪等课题，掌握数字信号的分析方法和处理方法。

（2）按要求编写课程设计报告书，能正确阐述设计和实验结果等等。

（3）通过课程设计培养学生严谨的科学态度，认真的工作作风和团队协作精神。

（4）在老师的指导下，要求每个学生独立完成课程设计的全部内容。

1.3设计平台

课程设计的主要设计平台式MATLAB7.0。

如下图1.1所示：

MATLAB是矩阵实验室（MatrixLaboratory）的简称，是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。

MATLAB和Mathematica、Maple并称为三大数学软件。

它在数学类科技应用软件中在数值计算方面首屈一指。

MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。

图1.1MATLAB7.0的设计平台

MATLAB的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB来解算问题要比用C，FORTRAN等语言完成相同的事情简捷得多，并且MathWork也吸收了像Maple等软件的优点,使MATLAB成为一个强大的数学软件。

在新的版本中也加入了对C，FORTRAN，C++，JAVA的支持。

可以直接调用,用户也可以将自己编写的实用程序导入到MATLAB函数库中方便自己以后调用，此外许多的MATLAB爱好者都编写了一些经典的程序，用户可以直接进行下载就可以用[2]。

2设计原理

2.1FIR滤波器

FIR（FiniteImpulseResponse）滤波器：

有限长单位冲激响应滤波器，是数字信号处理系统中最基本的元件，它可以在保证任意幅频特性的同时具有严格的线性相频特性，同时其单位抽样响应是有限长的，因而滤波器是稳定的系统。

因此，FIR滤波器在通信、图像处理、模式识别等领域都有着广泛的应用。

有限长单位冲激响应（FIR）滤波器有以下特点：

（1）系统的单位冲激响应h（n）在有限个n值处不为零；

（2）系统函数H（z）在|z|>0处收敛，极点全部在z=0处（因果系统）；

（3）结构上主要是非递归结构，没有输出到输入的反馈，但有些结构中（例如频率抽样结构）也包含有反馈的递归部分。

FIR滤波器的系统函数用下式表示：

。

H（n）就是FIR滤波器的单位脉冲响应。

FIR滤波器最重要的优点就是由于不存在系统极点，FIR滤波器是绝对稳定的系统。

相较于IIR滤波器,FIR滤波器有以下的优点：

（1）可以很容易地设计线性相位的滤波器。

线性相位滤波器延时输入信号，却并不扭曲其相位。

（2）实现简单。

在大多数DSP处理器，只需要对一个指令积习循环就可以完成FIR计算。

（3）适合于多采样率转换,它包括抽取（降低采样率），插值（增加采样率）操作。

无论是抽取或者插值，运用FIR滤波器可以省去一些计算，提高计算效率。

相反，如果使用IIR滤波器，每个输出都要逐一计算，不能省略,即使输出要丢弃。

（4）具有理想的数字特性。

在实际中，所有的DSP滤波器必须用有限精度（有限bit数目）实现，而在IIR滤波器中使用有限精度会产生很大的问题，由于采用的是反馈电路，因此IIR通常用非常少的bit实现，设计者就能解决更少的与非理想算术有关的问题。

（5）可以用小数实现。

不像IIR滤波器，FIR滤波器通常可能用小于1的系数来实现。

（如果需要，FIR滤波器的总的增益可以在输出调整）。

当使用定点DSP的时候，这也是一个考虑因素，它能使得实现更加地简单。

FIR滤波器的缺点在于它的性能不如同样阶数的IIR滤波器，不过由于数字计算硬件的飞速发展，这一点已经不成为问题。

再加上引入计算机辅助设计，FIR滤波器的设计也得到极大的简化。

基于上述原因，FIR滤波器比IIR滤波器的应用更广。

2.2窗口设计法

窗函数设计法的基本思想是用FIRDF逼近洗完的滤波特性。

设希望逼近的滤波器的频率响应函数为

，其单位脉冲响应为

表示。

为了设计简单方便，通常选择

为具有片段常数特性的理想滤波器。

因此

是无限长非因果序列，不能直接作为FIRDF的单位脉冲响应。

窗函数设计法就是截取

为有限长的一段因果序列，并用合适的窗函数进行加权作为FIRDF的单位脉冲响应h（n）。

用窗函数法设计FIRDF的具体设计步骤如下：

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

。

（2）求出

。

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

h（n）=

w（n）。

式中，w（n）称为窗函数，其长度为

。

如果要求设计第一类线性相位FIRDF，则要求h（n）关于（N-1）/2点偶对称。

而

关于n=t点偶对称，所以要求t=（N-1）/2。

同时要求w（n）关于（N-1）/2点偶对称。

用窗函数法设计第一类线性相位FIRDF的步骤：

（1）选择窗函数类型和长度，写出才窗函数

的表达式。

根据阻带最小衰减选择窗函数

的的类型，再根据过渡带宽度确定所选窗函数的长度

。

用窗函数法设计的FIRDF通带波纹幅度近似等于阻带波纹幅度。

一般阻带最小衰减达到40dB以上，则通带最大衰减就小于0.1dB。

所以用窗函数法设计FIRDF时，通常只考虑阻带最小衰减就可以了。

（2）构造希望逼近的频率响应函数

。

根据设计需要，一般选择线性相位理想滤波器（理想低通、理想高通、理想带通、理想带阻）。

理想滤波器的截止频率

近似为最终设计的FIRDF的过渡带中心频率，幅度函数衰减一半（约-6dB）。

所以一般取

，

和

分别为通带边界频率和阻带边界频率。

（3）计算

。

（4）加窗得到设计结果：

h（n）=

w（n）。

3设计步骤

3.1设计流程

3.2采集音乐信号

在网络中找到自己想要的歌曲文件，并下载下来，即可对音乐文件进行文件格式转换及声音属性的变化。

在转换过程中，将音乐格式设置为wav格式，属性设置为4.410kHz，8位，单声道，128KB/秒，然后点击确定。

然后再点击文件/另存为，即可将该语音文件存在MATLAB文件下的work文件夹里。

然后在Matlab软件平台下，利用函数wavread对语音信号进行采样，记住采样频率和采样点数。

通过wavread函数的使用，让我们很快理解了采样频率、采样位数等概念。

采集完成后在信号中加入一个单频噪声。

对应程序如下：

[x,fs,bits]=wavread（‘yuntian.wav’）;%输入参数为文件的全路径和文件名（因为语音文件直接放在Matlab下的work文件夹里，所以直接使用文件名即可），输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。

sound（x,fs,bits）;%按指定的采样率和每样本编码位数回放所录语音文件

N=length（x）;%计算信号x的长度

fn=1900;%单频噪声频率，此参数可自己设计

t=0:

1/fs:

（N-1）/fs;%计算时间范围，样本数除以采样频率

x=x’;%将其转置

y=x+sin（fn*2*pi*t）;%在原信号上添加噪声

sound（y,fs,bits）;%应该可以明显听出有尖锐的单频啸叫声

3.3滤波器设计

在该滤波器的设计中，我们给出该滤波器的性能指标如下：

fpd=1800;fsd=1850;fsu=1950;fpu=2000;Rp=1;As=40;

截止频率也可以任意自选，在单频噪声干扰附近即可。

在这里，很重要的是通带截止db值的设置。

这个值一定要根据我们使用的设计滤波器的方法来设定。

因为我使用的是凯塞窗法，凯塞窗函数中，滤波器的过渡带宽为15.4

，最小阻带衰减为108db。

所以，一定要将通带截止db值设置的小于108，所以，我将其设置为100db。

在这里我是使用窗函数法设计上面要求的FIR滤波器。

在Matlab中，利用凯塞窗设计FIR滤波器，利用Matlab中的函数freqz画出该滤波器的频率响应。

首先，我们利用数字信号处理里面学过的知识，根据自己选定的参数，用指定的方法设计FIR滤波器，得到FIR滤波器的阶数M。

随后调用nuttallwin（M）函数产生M阶的凯塞窗。

然后，调用自编ideal_lp函数计算理想带阻滤波器的脉冲响应。

最后，再调用自编freqz_m（h_bs,1）函数即可计算得到该滤波器的频率特性。

我设计滤波器过程对应的程序如下：

fpd=1800;fsd=1850;fsu=1950;fpu=2000;Rp=1;As=40;%带阻滤波器设计指标

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;

M=ceil（6.6*pi/dw）+1;%计算凯塞窗设计该滤波器时需要的阶数

n=0:

M-1;%定义时间范围

w_ham=nuttallwin（M）;%产生M阶的凯塞窗

hd_bs=ideal_lp（wcd,M）+ideal_lp（pi,M）-ideal_lp（wcu,M）;

%调用自编函数计算理想带阻滤波器的脉冲响应

h_bs=w_ham'.*hd_bs;%用窗口法计算实际滤波器脉冲响应

[db,mag,pha,grd,w]=freqz_m（h_bs,1）;%调用自编函数计算滤波器的频率特性

在将滤波器设计好，频率特性求出来之后，我编写如下程序，来画出该带阻滤波器的特性图。

程序如下：

subplot（2,2,1）;plot（w,db）;gridon;

xlabel（'频率'）;ylabel（'db'）;

subplot（2,2,2）;plot（w,mag）;gridon;

xlabel（'频率'）;ylabel（'幅度'）;

subplot（2,2,3）;plot（w,pha）;gridon;

xlabel（'频率'）;ylabel（'相位'）;

subplot（2,2,4）;plot（h_bs）;gridon;

xlabel（'频率'）;ylabel（'滤波器脉冲响应'）;通过该程序得到的图如下图3.3.1所示。

图3.3.1滤波器的特性图

3.4信号滤波处理

在将滤波器设计好后，我们用自己设计的带阻滤波器对采集的语音信号进行滤波。

在Matlab中，FIR滤波器利用函数fftfilt对信号进行滤波。

我对信号进行滤波处理对应的程序如下：

y_fil=fftfilt（h_bs,1,y）;%用设计好的滤波器对加噪信号y进行滤波

Y_fil=fft（y_fil）;Y_fil=Y_fil（1:

length（Y_fil）/2）;%计算频谱取前一半

3.5结果分析

在将加噪信号滤波之后，我们将滤波前后语音信号的波形及频谱图相互比较。

在同一张大图里分别绘制原始信号x，加噪信号y，滤波去噪信号y_fil的时域波形和频谱，以便比较和分析。

画这三个信号的比较图的对应的程序如下：

subplot（3,2,1）;plot（t,x）;gridon;xlabel（'时间'）;ylabel（'原始信号'）;

subplot（3,2,2）;plot（f,X）;gridon;xlabel（'时间'）;ylabel（'加入的噪声信号'）;

subplot（3,2,3）;plot（t,y）;gridon;xlabel（'频率'）;ylabel（'原始信号频谱'）;

subplot（3,2,4）;plot（f,Y）;gridon;xlabel（'频率'）;ylabel（'加入的噪声信号频谱'）;

subplot（3,2,5）;plot（t,y_fil）;gridon;xlabel（'时间'）;ylabel（'滤波之后的信号'）;

subplot（3,2,6）;plot（f,Y_fil）;gridon;xlabel（'频率'）;ylabel（'滤波之后的信号频谱'）;

经过这段程序画出来的三个信号的时域波形和频谱图如下图3.5.1所示。

图3.5.1三个信号的时域波形和频谱图比较

从图3.5.1中我们可以看出，原信号与滤波去噪信号的时域图基本相似，只有边缘部分有点差异；原信号与滤波去噪信号的频谱图波形也大致相似。

通过观察可以看到，加噪信号的时域图中大部分都被加入的噪声给遮盖了，加噪信号的频谱图中，我们可以很明显地看到与原信号频谱图相比，它在1900Hz左右处有一个尖脉冲，而滤波去噪信号的频谱图中该尖脉冲已经消失，波形大致与原图相似，可见滤波去噪效果基本不错。

在将三个信号的时域波形和频谱图比较之后，我们还要通过回放去滤波去噪语音信号，来跟原信号相比，以检验滤波器的效果。

在Matlab中，函数sound可以对声音进行回放。

其调用格式为：

sound（x，fs，bits）。

我用sound（y_fil,fs,bits）语句回放该滤波去噪信号，便可以感觉到滤波后的语音信号与原信号差不多，但仍有一点点变化。

4出现的问题及解决方法

原计划是准备二周的课程设计的，因为第二个礼拜的白天我们得去参加安博牛耳的一周免费实训，所以课程设计大部分内容都是在第一周完成的。

在这前面一周的时间里，我们每天的上午都是待在学校的机房上机，同时，老师跟我们进行讲解，然后，当我们有不懂的问题时可以随时向老师请教。

在老师没有进行讲解时，我总以为这次的课程设计会有不小的难度，可是当老师跟我们讲解了之后，我却发现，我们的课程设计其实也不是很难，只要我们细心一点，并有耐心，就能解决我们所遇到的问题。

首先，我按照课程设计的要求在网络上找到一段音乐，并对其进行格式转换及属性变换，但是当我真正制作完成之后并存盘后，运行MATLAB后，才感觉音乐文件过于长，造成出现的频谱图很模糊，根本看不清楚，后面看了一下任务指导书及问了同学之后，截取其中的一小段，效果就好了许多。

在此次课程设计中，虽然没有预计的那么难，可还是遇到了一些小小的问题。

庆幸的是，能够在老师和同学的帮助下解决了。

我觉得最重要的就是要有耐心，要细致，要不懂就问，这样，只要自己坚持，一切问题都会解决的。

5结束语

在本次课程设计中，让我知道了以前所学《信号与系统》与《数字信号处理》的结合，并在实际运用中设计滤波器的过程。

课程设计不仅要求对滤波器理论的研究，更重要的是实际设计中遇到的问题。

因为有了这次课程设计，我不得不对其设计原理进行更深一层次的理解，对书中原来学到的理论，仅知道了其表面，而不知其原因。

在设计中也使我对一些概念有了更深刻的认识。

例如：

在指标方面，我混淆了模拟指标和数字指标的概念，经过老师的点拨，自己更加明确，而且记忆深刻。

还有在课程设计中每一次的数据输入都有其重要意义，用MATLAB编译程序时，可以根据滤波器指标的要求实时知道对滤波器的影响。

例如，编好程序后，调试成功，任意改变输入阻带或者通带衰减，可以看到输出波形的变化，改变截止频率wc，同样可以看到输出波形的变化。

由此，对理论的理解就更加简单方便，而且记忆力深刻。

除此之外，对程序的编译不是一蹴而就的，而是经过多次的编译与调试。

以前用MATLAB就是简单的输入，可以说都不是自己的劳动成果，但这次不一样，课程设计没有别人给你编好，而是自己写，出错率就大大提高了，但这是过程，学习就是在过程中进行的，经过自己几天的脑力劳动，再加上同学们的帮助，不仅对读程序有了很大提高，而且自己的编译水平也上了一个新台阶，更加熟系了MATLAB的应用，也对其中的函数有了大概的了解，对其中一些函数也相当熟练，滤波器设计中用到了ceil（）、freqz（）等很重要的函数，虽然可以用fir1（）、fir2（）可以直接调用，但那样就达不到对真正理论设计过程的理解和运用。

完成整个设计过程后，学到的东西已经不仅仅上面的那些东西，还有就是同学之间的共同努力和探讨和设计过程中的每一个细节，也许每一个细节的错误就可能导致结局的失败，所以我认为这次收获最大的莫过于静心，学习不能急，一定要冷静，心无旁骛，不放过任何一个细节，就能带来凯旋的消息。

在此，感谢老师一学年以来对我们不仅学习上的关心，还有生活中的照顾，我也不会辜负老师的期望，继续努力深造。

参考文献

[1]程佩青．数字信号处理教程〔M〕．北京：

清华大学出版社，2002．

[2]刘敏，魏玲．Matlab通信仿真与应用〔M〕．北京：

国防工业出版社，2001．

［3］SanjitK.Mitra著，孙洪，余翔宇译，《数字信号处理实验指导书》，电子工业

出版社，2005年。

［4］EdwardW.Kamen,BonnieS.Heck编，《信号与系统基础——应用Web和MATLAB

（第二版）》，科学出版社，2002年。

［5］张贤达编,《现代信号处理（第二版）》，清华大学出版社，2002年。

附录1：

语音信号滤波去噪——用凯塞窗设计的滤波器源程序清单

//程序名称：

1.m

//程序功能：

在Matlab中，用凯塞窗法设计滤波器

//程序作者：

傅滔

//最后修改日期：

2010-3-7

[x,fs,bits]=wavread（'贝多芬的悲伤.wav'）;

sound（x,fs,bits）;%按指定的采样率和每样本编码位数回放

N=length（x）;%计算信号x的长度

fn=1900;

t=0:

1/fs:

（N-1）/fs;%计算时间范围，样本数除以采样频率

x=x（:

1）;

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:

length（X）/2）;Y=Y（1:

length（Y）/2）;%截取前半部分

deltaf=fs/2/length（X）;%计算频谱的谱线间隔

f=0:

deltaf:

fs/2-deltaf;%计算频谱频率范围

figure

（1）;

subplot（2,2,1）;plot（t,x）;gridon;

xlabel（'时间'）;ylabel（'原始信号'）;

subplot（2,2,2）;plot（f,X）;gridon;

xlabel（'时间'）;ylabel（'加入的噪声信号'）;

subplot（2,2,3）;plot（t,y）;gridon;

xlabel（'频率'）;ylabel（'原始信号频谱'）;

subplot（2,2,4）;plot（f,Y）;gridon;

xlabel（'频率'）;ylabel（'加入的噪声信号频谱'）;

fpd=1800;fsd=1850;fsu=1950;fpu=2000;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;

M=ceil（15.4*pi/dw）+1;%计算凯塞窗设计该滤波器时需要的阶数

n=0:

M-1;%定义时间范围

w_ham=nuttallwin（M）;%产生M阶的凯塞窗

hd_bs=ideal_lp（wcd,M）+ideal_lp（pi,M）-ideal_lp（wcu,M）;

%调用自编函数计算理想带阻滤波器的脉冲响应

h_bs=w_ham'.*hd_bs;%用窗口法计算实际滤波器脉冲响应

[db,mag,pha,grd,w]=freqz_m（h_bs,1）;%调用自编函数计算滤波器的频率特性

figure

（2）;

subplot（2,2,1）;plot（w,db）;gridon;

xlabel（'频率'）;ylabel（'db'）;

subplot（2,2,2）;plot（w,mag）;gridon;

xlabel（'频率'）;ylabel（'幅度'）;

subplot（2,2,3）;plot（w,pha）;gridon;

xlabel（'频率'）;ylabel（'相位