课程设计--语音信号滤波去噪使用BOHMAN窗设计的频率采样型FIR滤波器.doc

使用BOHMAN窗设计的频率采样型FIR滤波器 第页 共23页语音信号滤波去噪使用BOHMAN窗设计的频率采样型FIR滤波器学生姓名： 指导老师：摘 要 本课程设计主要内容是，利用窗口设计法设计频率采样型FIR滤波器并对加噪语音信号进行滤波去噪处理。仿真平台为MATLAB7.0，开发工具是M语言编程。在课程设计中，首先利用Windows下的录音机工具录制一段格式为.wav的语音信号，然后在MATLAB中对语音信号进行加噪，并绘制原始语音信号和加噪语音信号的时域和频域波形，进行频谱分析以确定所加噪声频率，再利用Bohman窗设计FIR滤波器，并检测是否达到指标，最后使用滤波器对信号进行滤波去噪处理，并通过对比原始信号，加噪信号，滤波去噪信号的时域和频域波形，或回放语音信号，检测是否设计成功。通过程序调试及完善，本课程设计滤波后的语音信号与原始语音信号基本一致，即设计的滤波器能够从含噪信号中滤除单频噪声，还原原始信号，达到了设计目的。关键词 课程设计；滤波去噪；FIR滤波器；BOHMAN窗；MATLAB1 引 言 本课程设计主要解决在含噪情况下如何对语音信号进行滤波去噪处理。设计采用BOHMAN窗设计的频率采样型FIR滤波器对一段含噪语音信号进行滤波去噪处理。数字滤波器是一种用来过滤时间离散信号的数字系统，它是通过对抽样数据进行数字处理来达到频域滤波的目的。数字信号处理与模拟信号处理是信号处理的子集。数字信号处理的目的是对真实世界的连续模拟信号进行测量或滤波。随着现代通信的数字化，数字滤波器变得更加重要。数字滤波器的种类很多，但总的来说可以分成两大类，一类是经典滤波器，另一类可称为现代滤波器。从滤波特性方面考虑，数字滤波器可分为数字高通、数字低通、数字带通和数字帯阻等滤波器。从实现方法上考虑，将滤波器分成两种，一种称为无限脉冲响应滤波器，简称IIR滤波器，另一种为FIR滤波器。设计FIR数字滤波器的方法有窗口设计法、频率采样法和等波纹最佳逼近法等。11.1 课程设计目的 数字信号处理课程设计是数字信号处理课程的重要实践性环节，是通信工程一次较全面的工程师能力训练。(1)通过课程设计，加深对数字信号处理理论方面的认识；(2)熟悉MALTAB工具软件的运用，学会用计算机采集语音信号，以及对语音信号进行时域和频域的分析；(3)掌握数字信号处理中FIR滤波器的设计方法和步骤；(4)掌握使用MATLAB对信号分析和处理的方法；(5)提高实际运用、综合运用能力以及分析问题、解决问题能力；培养严肃认真、实事求是的科学态度，理论联系实际的工作作风和辩证思维能力。1.2 课程设计的要求(1)滤波器指标必须符合工程实际。(2)设计完后应检查其频率响应曲线是否满足指标。(3)处理结果和分析结论应该一致，而且应符合理论。(4)独立完成课程设计并按要求编写课程设计报告书。1.3 设计平台 本课程设计使用的设计平台为MATLAB7.0。MATLAB是美国MathWorks公司出品的商业数学软件，主要包括MATLAB和Simulink两大部分。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言的编辑模式，代表了当今国际科学计算软件的先进水平。MATLAB在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。MATLAB的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB来解算问题要比用其他语言完成相同的事情简捷得多，并且MATLAB也吸收了像Maple等软件的优点，使MATLAB成为一个强大的数学软件。在新的版本中也加入了对C，FORTRAN，C+，JAVA的支持。可以直接调用,用户也可以将自己编写的实用程序导入到MATLAB函数库中方便自己以后调用。22 设计原理2.1 FIR滤波器 FIR(Finite Impulse Response)滤波器：有限长单位冲激响应滤波器，是数字信号处理系统中最基本的元件，它可以在保证任意幅频特性的同时具有严格的线性相频特性，同时其单位抽样响应是有限长的，因而滤波器是稳定的系统。因此，FIR滤波器在通信、图像处理、模式识别等领域都有着广泛的应用。FIR滤波器相对于IIR滤波器具有易实现线性相位，系统总是稳定的，允许设计多通带（或多阻带）滤波器的优点，FIR滤波器主要有4种结构，直接型，级联型，线性相位型，频率采样型，与IIR结构相比，结构是相对简单的。FIR系统滤波器的缺点在于它的性能不如同样阶数的IIR滤波器，不过由于数字计算硬件的飞速发展，这一点已经不成为问题。再加上引入计算机辅助设计，FIR滤波器的设计也得到极大的简化。3 FIR滤波器的设计方法有多种，如窗口设计法，频率采样法，最优等波纹设计法等。随着MATLAB软件尤其是MATLAB的信号处理箱工具的不断完善，不仅数字滤波器的计算机辅助设计有了可能，而且还可以使设计达到最优化。 有限长单位冲激响应（FIR）滤波器有以下特点：(1) 系统的单位冲激响应h (n)在有限个n值处不为零(2) 系统函数H(z)在|z|0处收敛，极点全部在z = 0处（因果系统）(3) 结构上主要是非递归结构，没有输出到输入的反馈，但有些结构中（例如频率抽样结构）也包含有反馈的递归部分。设FIR滤波器的单位冲激响应h (n)为一个N点序列，0 n N 1，则滤波器的系统函数为H(z)=h(n)*z-n (2-1) 就是说，它有（N1）阶极点在z = 0处，有（N1）个零点位于有限z平面的任何位置。FIR滤波器的设计分三步完成：(1)技术要求：在设计滤波器之前，必须要有某些技术要求。这些技术要求是由用途决定的。(2)近似：一旦技术要求确定之后，就要用已学过的各种概念和数学提供一种滤波器的表述，它接近于所给出的一组技术要求。这一步是属于滤波器设计的范畴。(3)实现：上面一步的结果是一个滤波器的表述，它可能是一个差分方程的形式，或者是某一系统函数H(z)，或者是某一脉冲响应h(n)。依据这个表述要用硬件实现这个滤波器，或者在一台计算机上通过软件实现。42.2 窗口设计法窗口设计的基本思想是要选取某一种合适的理想频率选择性滤波器（这种滤波器总是有一个非因果，无限长的脉冲响应），然后将它的脉冲响应截断（或加窗）以得到一个线性相位和因果的FIR滤波器。因此，这种方法的重点在于选择某种恰当的窗函数和一种适合的理想滤波器。现在用Hd(ejw)代表一理想频率选择性滤波器，它在整个通带内有单位幅度增益和线性相位特性，而阻带内具有零响应。一理想带宽为wc的LPF由下式给出为 1e-jw , |w|wc Hd(ejw)= (2-2) 0 , wc x,fs,bits=wavread(F:yuyin.wav); % 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。 sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放 N=length(x); % 计算信号x的长度 fn=9000; % 单频噪声频率，此参数可改 t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率 x=x(:,1); % 将双声道转为单声道 y=x+0.1*sin(fn*2*pi*t); % 加噪声 sound(y,fs,bits); % 应该可以明显听出有尖锐的单频啸叫声 fs % 输出采样频率输出fs=22050Hz大于被采样信号频率的两倍，采样之后的数字信号完整地保留了原始信号中的信息。然后对语音号进行快速傅里叶变换，得到信号的频谱特性。对应程序如下：X=abs(fft(x); Y=abs(fft(y); % 对原始信号和加噪信号进行fft变换,取幅度谱X=X(1:N/2); Y=Y(1:N/2); % 截取前半部分deltaf=fs/N; % 计算频谱的谱线间隔f=0:deltaf:fs/2-deltaf; % 计算频谱频率范围 最后画出加噪前后语音信号的时域和频域波形，便于分析设计滤波器。如图3-3所示。通过图3-5，可以明显观察到，在时域波形图中，加噪后的语音信号明显被噪声遮盖了许多，在频域图中，在频率为9000Hz处，有一个尖脉冲，也就是噪声。图3-5 加噪前后语音信号的时域和频域波形3.3 滤波器设计 本课程设计利用BOHMAN窗实现FIR滤波器的设计。由于加入的单频噪声频率fn为9000Hz，根据加噪前后语言信号的幅度谱图分析可知，语言信号的有效带宽0Hz5000Hz左右，所以可设计一低通滤波器来滤除噪声。由于BOHMAN窗的过渡带宽定为，能提供的最小阻带衰减为51.5dB。所以阻带衰减As设置为一个适当小于51.5dB的参数。通带波纹Rp设置为1dB。设置的低通滤波器的性能指标如下：fp=4000Hz， fo=8800Hz， Rp=1dB， As=45dB 其中fp为通带截止频率，fo为阻带截止频率， Rp为通带波纹，As为阻带衰减。确定好滤波器的参数后，先计算出滤波器的阶数M，然后调用自编ideal_lp函数计算理想低通滤波器的脉冲响应hd，再调用bohmanwin(M)函数产生M阶的BOHMAN窗，并计算出实际滤波器脉冲响应，最后调用自编freqz_m函数计算出滤波器的频率特性。最后计算出真正的阻带衰减As，和真正的通带波纹Rp，检测设计的滤波器是否符合要求。 得出的M =24，Rp =0.0276，As = 49。所以滤波器的长度M=24，真正的通带波纹Rp是0.0276dB，小于设置的指标1dB，真正的阻带衰减As是49dB，大于设置的指标45dB,小于最小阻带衰减51.5dB。所以设计的低通滤波器满足要求。最后画出该低通滤波器的频谱特性图。如图3-6所示。图3-6滤波器的频谱图 滤波器设计过程对应的程序如下： Rp=1;As=45;fs=22050; fp=4000;fo=8800; % 低通滤波器设计指标 wp=fp/fs*2*pi;ws=fo/fs*2*pi;tr_width=ws-wp; % 将Hz为单位的模拟频率换算为rad为单位的数字频率 M=ceil(10*pi/tr_width)+1 % 计算BOHMAN窗设计该滤波器时需要的阶数 n=0:1:M-1; % 定义时间范围 wc=(ws+wp)/2; hd=ideal_lp(wc,M); % 调用自编函数计算理想低通滤波器的脉冲响应 w_bohmanwin=(bohmanwin(M); % 产生M阶的BOHMAN窗 h=hd.*w_bohmanwin; % 计算实际滤波器脉冲响应 db,mag,pha,grd,w=freqz_m(h,1); % 调用自编函数计算滤波器的频率特性 delta_w=2*pi/1000; Rp=-(min(db(1:1:wp/delta_w+1) % 真正的通带波纹 As=-round(max(db(ws/delta_w+1:1:501) % 真正的阻带衰减3.4 信号滤波处理将设计好的FIR低通滤波器滤波器对信号进行滤波去噪处理。对应程序如下： y_fil=filter(h,1,y);% 用设计好的滤波器对y进行滤波 Y_fil=fft(y_fil);Y_fil=Y_fil(1:N/2); % 计算频谱取前一半3.5 结果分析为了分析滤波器是否设计成功，首先用MATLAB绘制出原始信号，加噪信号，滤波去噪信号的时域波形和频谱，如图3-7所示。通过观察和分析图像，滤波后的语音信号时域波形与原始语音信号时域波形基本一致，位于幅度谱9000Hz处的单频噪声被滤除，与原始语音信号幅度谱基本一致。图3-7 原始信号，加噪信号，滤波去噪信号的时域波形和频谱然后通过回放原始语音信号，加噪语音信号，去噪语音信号，对比发现去噪语音信号与原始语音信号基本一致，即滤波器滤除了噪声的干扰。对应程序如下： sound(x,fs,bits); % 播放原始语音信号 sound(y,fs,bits); % 播放加噪语音信号 sound (y_fil,fs,bits); % 播放去噪语音信号通过对比分析可得出结论，本课程设计滤波后的语音信号与原始语音信号基本一致，即设计的滤波器能够从含噪信号中滤除单频噪声，还原了原始信号，达到了设计目的。3.6 绘制滤波器结构图本课程设计为频率采样型FIR滤波器，先求出滤波器的系数，对应语句如下： C,B,A=dir2fs(h) % 计算滤波器系数得出滤波器的系数为：C = 1.9987 1.9987 1.9965 1.9947 1.9577 1.6597 0.9734 0.3059 0.0344 0.0046 0.0021 0.99960.0000B = -0.9914 0.9914 0.9659 -0.9659 -0.9239 0.9239 0.8660 -0.8660 -0.7934 0.7934 0.7071 -0.7071 -0.6088 0.6088 0.5000 -0.5000 -0.3827 0.3827 0.2588 -0.2588 -0.1305 0.1305A = 1.0000 -1.9319 1.0000 1.0000 -1.7321 1.0000 1.0000 -1.4142 1.0000 1.0000 -1.0000 1.0000 1.0000 -0.5176 1.0000 1.0000 -0.0000 1.0000 1.0000 0.5176 1.0000 1.0000 1.0000 1.0000 1.0000 1.4142 1.0000 1.0000 1.7321 1.0000 1.0000 1.9319 1.0000 1.0000 -1.0000 01.0000 1.0000 0 绘制的频率采样型FIR滤波器如图3-8所示。 1.99871.9987 1.9965 1.9947 1.95771.6597 0.97340.3059 0.0344 0.0046 0.0021 0.9996 -1z-1z-1z-1z-1z-1z-1z-1z-1z-1z-1z-1z-1z-1z-1z-1z-1z-1z-1z-1z-1z-1z-1 z-1 1.9319 -11.7321-11.4142-11-10.5176-10-10.5176-1-1-1-1.4142-1-1.7321-1-1.9319 -1 1 -0.9914 0.9914 0.9659 - 0.9659 -0.9239 0.9239 0.8660 -0.8660 -0.7934 0.7934 0.7071 -0.7071 -0.6088 0.6088 0.5000 -0.5000 -0.3827 0.3827 0.2588 -0.2588 -0.1305 0.1305-z-24 y(n) x(n)图3-8 频率采样型FIR滤波器结构图4 出现的问题及解决方法最开始遇到的问题主要是对理论知识不够熟悉。由于数字信号处理是上个学期学习的课程，所以理论知识对我来说已经有点陌生了，最开始拿到课程设计题目时不知道从哪下手，在复习并掌握了FIR滤波器，窗口设计法相关的知识点后，我对本课程设计有了一定了理解。对语音信号加噪时遇到的问题：开始我将单频噪声fn设置为2000Hz，绘制出加噪前后语音信号的时域和频域波形后，发现语言信号的有效带宽从0Hz5000Hz左右噪声处于语音信号的有效频带范围内，要滤除噪声就要设置一个带阻滤波器，带阻滤波器相对于低通滤波器来说较为复杂，经过多次调试后，将fn设置为9000Hz，所以本课程设计中滤波器设置为低通滤波器。由于坐标设置不适当，最初显示的加噪前后语音信号的时域和频域波形对比不是很明显，对坐标进行调整后，通过对比加噪前后的频域图，可以明显观察到在频率为9000Hz处，有一个尖脉冲，也就是噪声。设计滤波器时遇到的问题：由于本课程设计采用的是BOHMAN窗，教材中并没有介绍，所以在设计时对我BOHMAN窗的性能并不了解，课程设计提供的资料中，BOHMAN窗的过渡带宽为，最小阻带衰减为51.5dB。但是运行程序后得出的真正的通带波纹Rp和真正的阻带衰减As是总是不达标，滤波器不符合设计要求，老师建议将过渡带宽增大，经过多次调试后，最终将过渡带宽定为，设计的低通滤波器也达到了指标。最初得到的BOHMAN窗的阶数M为456，由于要画出滤波器的结构图，所以阶数过高。老师建议增大滤波器的过渡带宽，或者对录制的语音信号的属性进行适当的调整。最后将语音信号的属性调整为22.050kHz。运行程序后最终BOHMAN窗的阶数M减小为24。为了尽量减少BOHMAM窗的阶数，减小M的值，低通滤波器的过渡带宽参数设置的过大，加噪语音信号通过低通滤波器后不能将噪声完全滤除，经过多次调试后将参数设置为了合适的值，能完全滤除噪声。绘制频率采样型FIR滤波器时遇到的问题：虽然M值调至为了较小的值，但是对于绘制滤波器的结构图来说还是较大，而且我对word绘图比较陌生，所以绘制滤波器时花费了较长的时间。程序中存在的问题：程序中有一些警告存在，多次修改调试程序还是没有解决，最后寻求了老师的帮助，老师解释到如果不影响整个程序的运行和滤波器的设计效果可以不多做考虑。编写报告时出现的问题：写报告时，主要是没有认真检查，存在较多的格式问题，报告内容也不够完善，经过老师的检查后，查出了较多的漏洞，对报告进行修正后，报告的内容变的更为充实和准确。5 结束语两周的课程设计说长不长说短不短，如果好好利用，我们不仅可以巩固学过的知识，也能拓宽自己的视野，提升自己的能力。从一开始的无从下手，到最后的自我提升，无可否认，这是对我们自身能力的重新认识和提升挑战。数字信号处理是通信、电子、信息领域中最重要的专业课之一，数字信号处理对于我们通信工程的学生来说是一门很重要的科目，不仅仅因为它涵盖了通信工程的许多理论知识，是考研的科目，而且其实践性很强。数字信号处理的目的是对真实世界的连续模拟信号进行测量或滤波。因此在进行数字信号处理之前需要将信号从模拟域转换到数字域，这通常通过模数转换器实现。而数字信号处理的输出经常也要变换到模拟域，这是通过数模转换器实现的。数字信号处理的算法需要利用计算机或专用处理设备如数字信号处理器（DSP）和专用集成电路（ASIC）等。数字信号处理技术及设备具有灵活、精确、抗干扰强、设备尺寸小、造价低、速度快等突出优点，这些都是模拟信号处理技术与设备所无法比拟的。数字信号处理的核心算法是离散傅里叶变换(DFT)，是DFT使信号在数字域和频域都实现了离散化，从而可以用通用计算机处理离散信号。而使数字信号处理从理论走向实用的是快速傅里叶变换(FFT)，FFT的出现大大减少了DFT的运算量，使实时的数字信号处理成为可能、极大促进了该学科的发展。数字信号处理课程设计也让我受益匪浅，通过此次课程设计，让我更深入的了解了语音信号滤波去噪，窗口设计法设置FIR滤波器的相关知识，注意到很多以前没有注意到的细节，学到了很多在书本上所没有学到过的知识。从一个小模块设计到具有一定功能的完整模块设计，从单个小问题的判断调试到整体化模块的分析思维，本次课程设计不但使我对数字信号处理了解得更透彻，也使我认识到理论当与实际相联系时才能使它的指导性充分地体现出来。当然，在设计过程中或多或少的遇到了一些问题，但在老师和同学的帮助下找出了原因所在。更培养了我利用自己的知识解决问题的能力。让我意识到将理论与实际结合起来的重要性，只有理论知识是远远不够的，只有把所学的理论知识与实践相结合起来，从理论中得出结论，才能将知识为我所用，从而运用到实际生活中，提高自己的实际动手能力和独立思考的能力。这次课程设计也我意识到遇到问题不要一味的依赖于老师，还可以通过上网，去图书馆查找资料解决问题，这样不仅提高了自己解决问题的能力，也更容易掌握知识。俗话说柴多火焰高，如果只凭自己的知识，我会浪费许多的时间和精力，可是滴水成海，在同学和老师的帮助下，我不仅把学到的知识运用的课程设计中，还体会到设计成功后带来的成就感。在本次课程设计中，特别感谢各位老师在本次课程设计中对我的细心指导，让我顺利完成了课程设计。参考文献1 高西全，丁玉美，阔永红数字信号处理原理、实现及应用北京：电子工业出版社，20092 张圣勤MATLAB7.0实用教材北京：机械工程出版社，20083 绿波杰能FIR滤波器百度文库：/view/2908166.htm，2013-2-274 张立材, 王民数字信号处理北京：人民邮电出版社，20085 维纳K英格尔，约翰G普罗克斯刘树棠数字信号处理(MATLAB版)西安：西安交通大学出版社，20086 赵红怡数字信号处理及其MATLAB实现北京：化学工业出版社，2002附录：语音信号滤波去噪设计源程序清单% 程序名称：voicenoisefilter.m% 程序功能：采用基于BOHMAN窗的窗口设计法，设计FIR滤波器对含噪语音进行滤波去噪处理。% 程序作者：王羚% 最后修改日期：2013-3-5%（1）对原始语音信号加噪x,fs,bits=wavread(F:yuyin.wav); % 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放N=length(x); % 计算信号x的长度fn=9000; % 单频噪声频率，此参数可改t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率x=x(:,1); % 将双声道转为单声道y=x+0.1*sin(fn*2*pi*t); % 加噪声sound(y,fs,bits); % 应该可以明显听出有尖锐的单频啸叫声fs % 输出采样频率%（2）画出加噪前后语音信号的时域和频域波形X=abs(fft(x); Y=abs(fft(y); % 对原始信号和加噪信号进行fft变换,取幅度谱X=X(1:N/2); Y=Y(1:N/2); % 截取前半部分deltaf=fs/N; % 计算频谱的谱线间隔f=0:deltaf:fs/2-deltaf; % 计算频谱频率范围magX=abs(X);magY=abs(Y);subplot(2,2,1);plot(t,x);gridxlabel(时间（单位： s）);ylabel(幅度);title(原始语音信号)subplot(2,2,2);plot(f,magX);gridxlabel(频率（单位： Hz）);ylabel(幅度谱);title(语音信号幅度谱图)subplot(2,2,3);plot(t,y);gridxlabel(时间（单位： s）);ylabel(幅度);title(加入单频干扰后的语音信号)subplot(2,2,4);plot(f,magY);gridxlabel(频率（单位： Hz）);ylabel(幅度谱);title(加入干扰后的语音信号幅度谱图)% （3）设置低通滤波器Rp=1;As=45;fs=22050; fp=4000;fo=8800; % 低通滤波器设计指标wp=fp/fs*2*pi;ws=fo/fs*2*pi;tr_width=ws-wp; % 将Hz为单位的模拟频率换算为rad为单位的数字频率M=ceil(10*pi/tr_width)+1 % 计算BOHMAN窗设计该滤波器时需要的阶数n=0:1:M-1; % 定义时间范围wc=(ws+wp)/2;hd=ideal_lp(wc,M); % 调用自编函数计算理想低通滤波器的脉冲响应w_bohmanwin=(bohmanwin(M); % 产生M阶的BOHMAN窗h=hd.*w_bohmanwin; % 计算实际滤波器脉冲响应db,mag,pha,grd,w=freqz_m(h,1); % 调用自编函数计算滤波器的频率特性delta_w=2*pi/1000;Rp=-(min(db(1:1:wp/delta_w+1) % 真正的通带波纹As=-round(max(db(ws/del