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

弹拨音乐滤波去噪使用凯塞窗 设计的 FIR 滤波器 学生姓名：唐柯学生姓名：唐柯 指导老师：胡双红指导老师：胡双红 摘 要 本课程设计主要是通过使用凯塞窗设计一个 FIR 滤波器以对弹拨进行 滤波去噪处理。本设计首先从网下下载一段弹拨音乐，依据对该信号的频谱分 析，给定相关指标。以 MATLAB 软件为平台，采用凯塞窗设计满足指标的 FIR 滤 波器，以该弹拨音乐进行滤波去噪处理。通过对比滤波前后的波形图，深入了 解滤波器的基本方法。通过程序调试及完善，该设计基本满足设计要求。 关键词 滤波去噪；FIR 滤波器；凯塞窗函数；MATLAB 1 1 引引 言言 数字滤波器是一种用来过滤时间离散信号的数字系统,它是通过对抽样数据 进行数学处理来达到频域滤波的目的。随着现代通信的数字化,数字滤波器变得 更加重要。数字滤波器的种类很多，但总的来说可以分成两大类，一类是经典 滤波器，另一类可称为现代滤波器。从滤波特性方面考虑，数字滤波器可分成 数字高通、数字低通、数字带通和数字带阻等滤波器。从实现方法上考虑，将 滤波器分成两种，一种称为无限脉冲响应滤波器，简称 IIR（Infinite Impulse Response）滤波器，另一种称为 FIR（Finite Impulse Response）滤波器1。设计 FIR 数字滤波器的方法有窗函数法、频率采样法和等波纹最佳逼近法等。 1.11.1 课程设计目的课程设计目的 数字信号处理（Digital Signal Processing，简称 DSP）是一门涉及许多学科 而又广泛应用于许多领域的新兴学科。20 世纪 60 年代以来，随着计算机和信 息技术的飞速发展，数字信号处理技术应运而生并得到迅速的发展。在过去的 二十多年时间里，数字信号处理已经在通信等领域得到极为广泛的应用。数字 信号处理是利用计算机或专用处理设备，以数字形式对信号进行采集、变换、 滤波、估值、增强、压缩、识别等处理，以得到符合人们需要的信号形式。 在本次课程设计中，最主要的设计是设计 FIR 滤波器，FIR 滤波器的设计方法 主要分为两类：第一类是基于逼近理想滤波器器特性的方法包括窗函数法、频 率采样法、和等波纹最佳逼近法；第二类是最优设计法。 本次的课程设计主要采用的是第一类设计方法，是利用凯塞窗函数法设计 FIR 滤波器对一段弹拨音乐进行滤波去噪，通过这一过程，对滤波前后波形进 行对比分析得到结论。此课程设计比较简单，主要是将书本中的知识运用到现 实中，并且根据自己对设计题目的理解，运用软件编写出程序实现这一设计， 也是我们对数字信号处理的原理进行验证的一个过程。对此，也可以加深我们 对所学知识的理解，培养我们的动手能力。 1.21.2 课程设计的要求课程设计的要求 （1）通过利用各种不同的开发工具实现模拟信号数字化、信道编解码、基带数 字信号编解码、数字信号的调制解调和弹拨音乐的滤波去噪等课题，掌握数字 信号的分析方法和处理方法。 （2）按要求编写课程设计报告书，能正确阐述设计和实验结果等等。 （3）通过课程设计培养学生严谨的科学态度，认真的工作作风和团队协作精神。 （4）在老师的指导下，要求每个学生独立完成课程设计的全部内容。 1.31.3 设计平台设计平台 课程设计的主要设计平台式 MATLAB 7.0。如下图 1.1 所示：MATLAB 是 矩阵实验室（Matrix Laboratory）的简称，是美国 MathWorks 公司出品的商业 数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计 算语言和交互式环境，主要包括 MATLAB 和 Simulink 两大部分。MATLAB 和 Mathematica、Maple 并称为三大数学软件。它在数学类科技应用软件中在数 值计算方面首屈一指。MATLAB 可以进行矩阵运算、绘制函数和数据、实现算 法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制 设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。 图 1.1 MATLAB 7.0 的设计平台 MATLAB 的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式 十分相似，故用 MATLAB 来解算问题要比用 C，FORTRAN 等语言完成相同的 事情简捷得多，并且 MathWork 也吸收了像 Maple 等软件的优点,使 MATLAB 成为一个强大的数学软件。在新的版本中也加入了对 C，FORTRAN，C+ ，JAVA 的支持。可以直接调用,用户也可以将自己编写的实用程序导入到 MATLAB 函数库中方便自己以后调用，此外许多的 MATLAB 爱好者都编写了 一些经典的程序，用户可以直接进行下载就可以用2。 2 2 设计原理设计原理 2.12.1 FIRFIR 滤波器滤波器 FIR(Finite Impulse Response)滤波器：有限长单位冲激响应滤波器，是数字 信号处理系统中最基本的元件，它可以在保证任意幅频特性的同时具有严格的 线性相频特性，同时其单位抽样响应是有限长的，因而滤波器是稳定的系统。 因此，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.22.2 窗口设计法窗口设计法 窗函数设计法的基本思想是用 FIRDF 逼近洗完的滤波特性。设希望逼近的 滤波器的频率响应函数为，其单位脉冲响应为表示。为了设计简 单方便，通常选择为具有片段常数特性的理想滤波器。因此是无 限长非因果序列，不能直接作为 FIRDF 的单位脉冲响应。窗函数设计法就是截 取为有限长的一段因果序列，并用合适的窗函数进行加权作为 FIRDF 的 单位脉冲响应 h(n)。 用窗函数法设计 FIRDF 的具体设计步骤如下： （1）构造希望逼近的频率响应函数。 （2）求出 。 （3）加窗得到 FIRDF 的单位脉冲响应 h(n)。h(n)=w(n)。式中，w(n)称为 窗函数，其长度为。如果要求设计第一类线性相位 FIRDF，则要求 h(n)关于N (N-1)/2 点偶对称。而关于 n=t 点偶对称，所以要求 t=(N-1)/2。同时要求 w(n)关于(N-1)/2 点偶对称。 用窗函数法设计第一类线性相位 FIRDF 的步骤： （1）选择窗函数类型和长度，写出才窗函数的表达式。 根据阻带最小衰减选择窗函数的的类型，再根据过渡带宽度确定所选窗函 数的长度。用窗函数法设计的 FIRDF 通带波纹幅度近似等于阻带波纹幅度。N 一般阻带最小衰减达到 40dB 以上，则通带最大衰减就小于 0.1dB。所以用窗函 数法设计 FIRDF 时，通常只考虑阻带最小衰减就可以了。 （2）构造希望逼近的频率响应函数。 根据设计需要，一般选择线性相位理想滤波器（理想低通、理想高通、理想带 通、理想带阻） 。理想滤波器的截止频率近似为最终设计的 FIRDF 的过渡带 c 中心频率，幅度函数衰减一半（约-6dB） 。所以一般取，和 分别为通带边界频率和阻带边界频率。 （3）计算。 （4）加窗得到设计结果：h(n)=w(n)。 3 3 设计步骤设计步骤 3.13.1 设计流程设计流程 用麦克风采集一段音乐将格式改为.wav 加入单频噪声 对音乐信号进行频谱分析，画出时域和频域波形 图 用凯塞窗设计 FIR 滤波器 画出其频率响应 用 FIR 滤波器对音乐信号进行滤波 画出音乐信号滤波前后波形并且进行比较分析 开始 结束 图 3.1 设计流程 3.23.2 采集弹拨音乐采集弹拨音乐 在网络中找到自己想要的歌曲文件，并下载下来，即可对音乐文件进行文 件格式转换及声音属性的变化。在转换过程中，将音乐格式设置为 wav 格式。 然后再点击文件/另存为，即可将该音乐文件存在 MATLAB 文件下的 work 文件 夹里。然后在 Matlab 软件平台下，利用函数 wavread 对音乐信号进行采样，记 住采样频率和采样点数。通过 wavread 函数的使用，让我们很快理解了采样频 率、采样位数等概念。采集完成后在信号中加入一个单频噪声。 对应程序如下： x,fs,bits=wavread(Music.wav); sound(x,fs,bits); N=length(x); % 计算信号 x 的长度 fn=1900; % 单频噪声频率 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:length(X)/2); Y=Y(1:length(Y)/2); %截取前半部分 deltaf=fs/2/length(X); %计算频谱的谱线间隔 f=0:deltaf:fs/2-deltaf; %计算频谱频率范围 用 plot 函数画出原始语音信号及加噪语音信号的时域图和频谱图，如图 3.2 所 示： 图 3.2 频谱分析图 3.33.3 滤波器设计滤波器设计 在该滤波器的设计中，我们给出该滤波器的性能指标如下： fpd=1800; fsd=1850; fsu=1950; fpu=2000;Rp=1; As=100; 截止频率也可以任意自选，在单频噪声干扰附近即可。 我设计滤波器过程对应的程序如下： 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(As-7.95)/(2.285* dw)+1)+1; % 计算凯塞窗设计该滤波器时需要的 阶数 n=0:M-1; beta=0.1102*(As-8.7); % 定义时间范围 w_kaiser= kaiser(M,beta); % 产生 M 阶的凯塞窗 hd_bs=ideal_lp(wcd,M)+ideal_lp(pi,M)-ideal_lp(wcu,M); % 调用自编函数计算理想带阻滤波器的脉冲响应 h_bs=w_kaiser.*hd_bs; % 用窗口法计算实际滤波器脉冲响应 db,mag,pha,grd,w=freqz_m(h_bs,1); % 调用自编函数计算滤波器的频率特性 在将滤波器设计好，频率特性求出来之后，我编写如下程序，来画出该带阻滤 波器的特性图。通过该程序得到的图如下图 3.3.1 所示。 图 3.3 滤波器的特性图 3.43.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.53.5 结果分析结果分析 在将加噪信号滤波之后，我们将滤波前后弹拨音乐的波形及频谱图相互比 较。在同一张大图里分别绘制原始信号 x，加噪信号 y，滤波去噪信号 y_fil 的 时域波形和频谱，以便比较和分析。经过这段程序画出来的三个信号的时域波 形和频谱图如下图 3.4 所示。 图 3.4 三个信号的时域波形和频谱图比较 从图 3.5.1 中我们可以看出，原信号与滤波去噪信号的时域图基本相似，只 有边缘部分有点差异；原信号与滤波去噪信号的频谱图波形也大致相似。通过 观察可以看到，加噪信号的时域图中大部分都被加入的噪声给遮盖了，加噪信 号的频谱图中，我们可以很明显地看到与原信号频谱图相比，它在 1900Hz 左 右处有一个尖脉冲，而滤波去噪信号的频谱图中该尖脉冲已经消失，波形大致 与原图相似，可见滤波去噪效果基本不错。 在将三个信号的时域波形和频谱图比较之后，我们还要通过回放去滤波去 噪弹拨音乐，来跟原信号相比，以检验滤波器的效果。在 Matlab 中，函数 sound 可以对声音进行回放。其调用格式为：sound (x，fs，bits)。我用 sound(y_fil,fs,bits)语句回放该滤波去噪信号，便可以感觉到滤波后的弹拨音乐 与原信号差不多，但仍有一点点变化。 4 4 出现的问题和解决方法出现的问题和解决方法 在这前面一周的时间里，我们每天的上午都是待在学校的机房上机，同时， 老师跟我们进行讲解，然后，当我们有不懂的问题时可以随时向老师请教。在 老师没有进行讲解时，我总以为这次的课程设计会有不小的难度，可是当老师 跟我们讲解了之后，我却发现，我们的课程设计其实也不是很难，只要我们细 心一点，并有耐心，就能解决我们所遇到的问题。 首先，我按照课程设计的要求在网络上找到一段音乐，并对其进行格式转 换及属性变换，但是当我真正制作完成之后并存盘后，运行 MATLAB 后，才 感觉音乐文件过于长，造成出现的频谱图很模糊，根本看不清楚，后面看了一 下任务指导书及问了同学之后，截取其中的一小段，效果就好了许多。 在此次课程设计中，虽然没有预计的那么难，可还是遇到了一些小小的问 题。庆幸的是，能够在老师和同学的帮助下解决了。我觉得最重要的就是要有 耐心，要细致，要不懂就问，这样，只要自己坚持，一切问题都会解决的。 5 5 结束语结束语 在本次课程设计中，让我知道了以前所学信号与系统与数字信号处 理的结合，并在实际运用中设计滤波器的过程。课程设计不仅要求对滤波器 理论的研究，更重要的是实际设计中遇到的问题。 因为有了这次课程设计，我不得不对其设计原理进行更深一层次的理解， 对书中原来学到的理论，仅知道了其表面，而不知其原因。在设计中也使我对 一些概念有了更深刻的认识。例如：在指标方面，我混淆了模拟指标和数字指 标的概念，经过老师的点拨，自己更加明确，而且记忆深刻。还有在课程设计 中每一次的数据输入都有其重要意义，用 MATLAB 编译程序时，可以根据滤 波器指标的要求实时知道对滤波器的影响。例如，编好程序后，调试成功，任 意改变输入阻带或者通带衰减，可以看到输出波形的变化，改变截止频率 wc， 同样可以看到输出波形的变化。由此，对理论的理解就更加简单方便，而且记 忆力深刻。 除此之外，对程序的编译不是一蹴而就的，而是经过多次的编译与调试。 以前用 MATLAB 就是简单的输入，可以说都不是自己的劳动成果，但这次不 一样，课程设计没有别人给你编好，而是自己写，出错率就大大提高了，但这 是过程，学习就是在过程中进行的，经过自己几天的脑力劳动，再加上同学们 的帮助，不仅对读程序有了很大提高，而且自己的编译水平也上了一个新台阶， 更加熟系了 MATLAB 的应用，也对其中的函数有了大概的了解，对其中一些 函数也相当熟练，滤波器设计中用到了 ceil（）、freqz（）等很重要的函数， 虽然可以用 fir1（）、fir2（）可以直接调用，但那样就达不到对真正理论设计 过程的 理解和运用。 完成整个设计过程后，学到的东西已经不仅仅上面的那些东西，还有就是 同学之间的共同努力和探讨和设计过程中的每一个细节，也许每一个细节的错 误就可能导致结局的失败，所以我认为这次收获最大的莫过于静心，学习不能 急，一定要冷静，心无旁骛，不放过任何一个细节，就能带来凯旋的消息。 在此，感谢老师一学年以来对我们不仅学习上的关心，还有生活中的照顾，我 也不会辜负老师的期望，继续努力深造。 参考文献参考文献 1 程佩青数字信号处理教程M北京：清华大学出版社，2002 2 刘敏，魏玲Matlab 通信仿真与应用M北京：国防工业出版社， 2001 3 Sanjit K.Mitra 著，孙洪，余翔宇 译，数字信号处理实验指导书，电子 工业出版社，2005 年。 4 Edward W. Kamen, Bonnie S.Heck 编，信号与系统基础应用Web 和 MATLAB（第二版），科学出版社，2002 年。 5 张贤达 编, 现代信号处理（第二版），清华大学出版社，2002 年。 附录 1：弹拨音乐滤波去噪用凯塞窗设计的滤波器源程序清单 x,fs,bits=wavread(Music.wav); sound(x,fs,bits); N=length(x); fn=1900; 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); X=X(1:length(X)/2); Y=Y(1:length(Y)/2); deltaf=fs/N; f=0:deltaf:fs/2-deltaf; subplot(2,2,1); plot(t,x); grid on; xlabel(时间（单位：s）); ylabel(幅值); title( 原始音乐信号); axis(0,4,-1,1); subplot(2,2,2); plot(f,X); grid on; xlabel(频率（单位：Hz）); ylabel(幅度谱); title(音乐信号幅度谱图); axis(0,4000,0,5000); subplot(2,2,3); plot(t,y);grid on; xlabel(时间（单位：s）); ylabel(幅值); title( 加入单频干扰后音乐信号); axis(0,4,-1,1); subplot(2,2,4); plot(f,Y);grid on; xlabel(频率（单位：Hz）); ylabel(幅度谱); title(加入干扰后的音乐信号幅度谱图); axis(0,4000,0,5000); 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; ws