《吹管音乐滤波去噪：基于汉宁创的FIR滤波器.doc》由会员分享，可在线阅读，更多相关《吹管音乐滤波去噪：基于汉宁创的FIR滤波器.doc（20页珍藏版）》请在三一办公上搜索。
1、吹管音乐滤波去噪基于汉宁窗的FIR滤波器学生姓名： 指导老师： 摘 要 从网站上下载一段吹管乐器演奏音乐，利用CE软件对音乐进行编辑。绘制波形并观察其频谱特点，加入一个带外单频噪声，用汉宁窗设计一个满足指标的FIR滤波器，对该含噪音乐信号进行滤波去噪处理，比较滤波前后波形和频谱并进行分析，根据结果和学过的理论得出合理结论。与不同信源相同滤波方法的同学比较各种信源的特点，与相同信源不同滤波方法的同学比较各种滤波方法性能优劣。关键词 滤波去噪；FIR滤波器；汉宁窗；MATLAB1 引 言本课程设计主要针对一段吹管音乐信号，在进行加噪后，利用窗函数设计法选择汉宁窗设计的FIR滤波器，对加噪后的吹管音
2、乐信号进行滤波去噪处理，并对前后时域波形和频域波形进行对比分析的程序设计。1.1 课程设计目的本次课设中的主要目的是让学生在熟悉Matlab语言环境，掌握其语言编程规则的前提下，利用汉宁窗设计一个符合要求的FIR滤波器来实现音乐信号的滤波去噪，并绘制滤波前后的时域波形和频谱图。根据图形分析判断滤波器设计的正确性。通过本次课设，我们能够学会如何综合运用课堂上学会的理论知识，增强自己的动能力与联系实际的能力，为以后的工作奠定基础。1.2 课程设计的要求（1）滤波器指标必须符合工程实际。（2）设计完后应检查其频率响应曲线是否满足指标。（3）处理结果和分析结论应该一致，而且应符合理论。（4）独立完成课
3、程设计并按要求编写课程设计报告书。1.3 课程设计平台课程设计的主要设计平台是MATLAB 7.0。MATLAB 的名称源自 Matrix Laboratory ，它是美国MathWorks公司生产的一个为科学和工程计算专门设计的交互式大型软件，是一个可以完成各种精确计算和数据处理的、可视化的、强大的计算工具。它集图示和精确计算于一身，在应用数学、物理、化工、机电工程、医药、金融和其他需要进行复杂数值计算的领域得到广泛应用。它不仅是一个在各类工程设计中便于使用的计算工具，而且也是一个在数学、数值分析和工程计算等课程教学中的优秀的教学工具，在世界各地的高等院校中十分流行，在各类工业应用中更有不俗
4、的表现。MATLAB可以在几乎所有的PC机和大型计算机上运行，适用于Windows、UNIX等各种系统平台。1MATLAB软件包括五大通用功能：数值计算功能（Nemeric）；符号运算功能（Symbolic）；数据可视化功能（Graphic）；数据图形文字统一处理功能（Notebook）和建模仿真可视化功能（Simulink）。其中，符号运算功能的实现是通过请求MAPLE内核计算并将结果返回到MATLAB命令窗口。该软件有三大特点：一是功能强大；二是界面友善、语言自然；三是开放性强。2MATLAB在信号与系统中的应用主要包括符号运算和数值计算仿真分析。由于信号与系统课程的许多内容都是基于公式演
5、算，而MATLAB借助符号数学工具箱提供的符号运算功能能基本满足信号与系统课程的需求。例如，解微分方程、傅里叶正反变换、拉普拉斯正反变换、z正反变换等。MATLAB在信号与系统中的另一主要应用是数值计算与仿真分析，主要包括函数波形绘制、函数运算、冲激响应与阶跃响应仿真分析、信号的时域分析、信号的频谱分析、系统的S域分析、零极点图绘制等内容。数值计算仿真分析可以帮助学生更深入理解信号与系统的理论知识，并为将来使用MATLAB进行信号处理领域的各种分析和实际应用打下基础3。2 设计原理2.1 数字信号处理数字信号处理（Digital Signal Processing，简称DSP）是一门涉及许多学
6、科而又广泛应用于许多领域的新兴学科4。20世纪60年代以来，随着计算机和信息技术的飞速发展，数字信号处理技术应运而生并得到迅速的发展。在过去的二十多年时间里，数字信号处理已经在通信等领域得到极为广泛的应用。数字信号处理是利用计算机或专用处理设备，以数字形式对信号进行采集、变换、滤波、估值、增强、压缩、识别等处理，以得到符合人们需要的信号形式。52.2 FIR滤波器有限长单位脉冲响应数字滤波器（Finite Impulse Response Digital Filter,缩写FIRDF）简称FIR滤波器，是数字信号处理系统中最基本的原件，其最大优点是可以实现线性相位滤波，可以在保证任意幅频特性的
7、同时具有严格的线性相频特性，满足了在数字通信和图像传输与处理等应用场合对线性相位的要求。FIR滤波器是全零点滤波器，硬件和软件实现结构简单，因而是十分稳定的系统。6FIR滤波器的设计方法主要分为两类：第一类是基于逼近理想滤波器器特性的方法包括窗函数法、频率采样法、和等波纹最佳逼近法；第二类是最优设计法。本次课设采用的是第一类设计法中的窗函数法。设FIR滤波器的单位脉冲响应的长度为，则其频率响应函数为 （2-1）一般将表示成如下形式： （2-2）式中，是的实函数（可以去负值）。与前面的表示形式，即相比， 与不同。与 不同。为了区别于幅频响应函数和相频响应函数，称为幅频特性函数，称为相频特性函数。
8、第一类线性相位FIR滤波器的相位特性函数是的严格线性函数： （2-3）2.3 窗口设计法窗口设计法是一种通过截断和计权的方法使无限长非因果序列成为有限长脉冲响应序列的设计方法。通常在设计滤波器之前，应该先根据具体的工程应用确定滤波器的技术指标。在大多数实际应用中，数字滤波器常常被用来实现选频操作，所以指标的形式一般为在频域中以分贝值给出的相对幅度响应和相位响应。6窗口设计法基本步骤如下：（1）根据过渡带宽及阻带衰减要求，选择窗函数的类型并估计窗口长度N。窗函数的类型可根据最小阻带衰减AS独立选择。（2）根据待求滤波器的理想频率响应求出理想单位脉冲响应hd(n)。（3）由性能指标确定窗函数W(n
9、)和长度N。（4）求得实际滤波器的单位脉冲响应h(n)，h(n)即为所设计FIR滤波器系数向量b(n)。 常见的窗函数性能表如下表2-1所示。图2.1 常见窗函数性能表名称滤波器过渡带宽最小阻带衰减名称滤波器过渡带宽最小阻带衰减矩形1.8/M21dBPARZENWIN6.6/M56dB巴特利特6.1/M25dBFLATTOPWIN19.6/M108dB汉宁6.2/M44dBGAUSSWIN5.8/M60dB汉明6.6/M51dBBARTHANNWIN3.6/M40dB布莱克曼11/M74dBBLACKMANHARRIS16.1/M109dBBOHMANWIN5.8/M51.5dBCHEBWIN
10、15.2/M113dBNUTTALLWIN15.4/M108dBTUKEYWIN2.4/M22dB2.4 汉宁窗（Hanning window）汉宁窗函数是余弦平方函数，又称之为升余弦函数，它的时域形式可以表为： （2-8）其中k=1,2,，k。它的频域幅度特性函数为： （2-9）其中为矩形窗函数的幅度频率特性函数。汉宁窗函数的最大旁瓣值比主瓣值低31dB，但是主瓣宽度比矩形窗函数的主瓣宽度增加了一倍，为。汉宁窗函数的时域幅度与频域幅度特性曲线的MATLAB实现的曲线图如图2-1所示。 图2-1 汉宁窗函数的时域幅度与频域幅度特性曲线3设计步骤3.1 设计流程图本课程设计主要是从网站上下载一段
11、吹管乐器演奏音乐，利用CE软件对音乐进行编辑。绘制波形并观察其频谱特点，加入一个带外单频噪声，用汉宁窗设计一个满足指标的FIR滤波器，对该含噪音乐信号进行滤波去噪处理，比较滤波前后波形和频谱并进行分析，根据结果和学过的理论得出合理结论。程序的设计流程图如下图3-1所示。 开始下载一段吹管音乐信号，用CE软件编辑格式为wav.加入单频干扰噪声对吹管音乐信号进行频谱分析，画出干扰前后的时域和频域波形图利用汉宁窗设计FIR滤波器对吹管音乐信号进行滤波比较滤波前后的时域波形和频谱图，并回放音乐信号，验证是否达到去噪效果是否达到去噪效果结束NOYES图3-1 程序设计流程图 3.2 编辑语音信号在网上下
12、载一段音乐，再利用CE软件将其转换成单声道的.格式文件，再将此.格式音乐控制在10秒内，以减少设计中的误差。然后在Matlab软件平台下，利用函数wavread对语音信号进行采样，记住采样频率和采样点数。CE软件操作界面如图3-2所示。 图3-2 CE软件操作界面3.3 语音加噪处理采集完成后在信号中加入一个单频噪声，绘制原音乐信号和加噪后的音乐信号的时域和频域的波形图。首先，输入原始音乐信号并播放一次。调用程序如下：x,fs,bits=wavread(h:2013DSPPurpleBambooTune.wav); % 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成
13、该波形文件时的采样率，bits是波形文件每样本的编码位数。sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放计算信号长度并加入噪声。调用程序如下：N=length(x); % 计算信号x的长度fn=2100; % 单频噪声频率t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率x=x(:,1);y=x+sin(fn*2*pi*t); % 加入一个单频噪声sound(y,fs,bits); % 可以明显听出有尖锐的单频啸叫声绘制原始音乐信号和加入噪声后的音乐信号的时域和频谱波形图。调用程序如下：X=abs(fft(x); Y=abs(fft(y);
14、% 对原始信号和加噪信号进行fft变换X=X(1:N/2); Y=Y(1:N/2); % 截取前半部分deltaf=fs/ N; % 计算频谱的谱线间隔f=0:deltaf:fs/2-deltaf; % 计算频谱频率范围figure(1);subplot(2,2,1);plot(t,x); % 布局为2*2的四个小图title(原始音乐信号);xlabel(时间(t);ylabel(幅度); %改变横纵坐标的范围axis(0,2,-1.5,1.5); %加上标题和横坐标名称grid on; % 加上网格subplot(2,2,2);plot(f,X);title(原始音乐信号频谱);xlabe
15、l(频率(f);ylabel(幅度谱);axis(0,3000,0,3000);grid on; subplot(2,2,3);plot(t,y);title(加入干扰后的音乐信号);xlabel(时间(t);ylabel(幅度);axis(0,2,-1.5,1.5);grid on;subplot(2,2,4);plot(f,Y);title(加入单频干扰后的音乐信号频谱);xlabel(频率(f);ylabel(幅度谱);axis(0,3000,0,3000);grid on;用绘图命令分别画出加噪前后的时域和频域波形，如下图3-3所示。图3-3 吹管音乐信号加入单频噪声前后的时域与频谱波
16、形图由上图可以看到，语音信号加入单频噪声后的时域波形比未加之前在幅度范围内有了明显的增加，在频谱方面可以看到除了在加了噪声之后的频谱图上的2100Hz出现一个明显的冲激信号外，其它地方均与未加时的原始吹管音乐信号频谱相同，这一现象表现在音乐播放时，可以听见一声尖锐的啸叫声。3.4 滤波器设计本次课程设计中主要应用汉宁窗设计出FIR滤波器。利用Matlab中的函数freqz画出各滤波器的频率响应，首先利用数字信号处理里面学过的知识，根据选定的参数，用汉宁窗函数法设计FIR数字滤波器，得到数字滤波器的参数b，a。其中b为系统函数的分子系数，a为系统函数分母系数。再调用freqz(b,a,512,f
17、s)即可得到该滤波器的频率响应。主程序如下：fpd=1450;fsd=1650;fsu=2250;fpu=2350; % FIR滤波器的上下截止频率Rp=1;As=37; % 带阻滤波器设计指标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.2*pi/dw)
18、+1; % 计算汉宁窗设计该滤波器时需要的阶数n=0:M-1; % 定义时间范围w_ham=hanning(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); % 调用自编函数计算滤波器的频率特性通过绘图工具可得出滤波器的波形图，如图3-4所示。图3-4 FIR滤波器的频率响应上图为用汉宁窗设计的FIR滤波器图，可以看出，阻带最
19、大衰减为-100dB，FIR滤波器的主瓣宽度很小，这样可以使过渡宽度很陡，旁瓣相对于主瓣也比较小。3.5 信号滤波处理用自己设计的各滤波器分别对采集的信号进行滤波，在Matlab中，FIR滤波器利用函数fftfilt对信号进行滤波，IIR滤波器利用函数filter对信号进行滤波，对语音信号进行滤波后，仔细对比滤波前和滤波后的语音信号图，得出结论。主程序如下：y_fil=filter(h_bs,1,y); % 用设计好的滤波器对y进行滤波Y_fil=fft(y_fil);Y_fil=Y_fil(1:N/2); % 计算频谱取前一半由绘图工具可以得出滤波前后的吹管音乐信号波形图、原始吹管音乐信号波
20、形图和加入噪声后的吹管音乐信号波形图，如图3-5所示。图3-5 滤波前后的波形图由上图可以看出，加噪后的吹管音乐信号经过FIR滤波器的滤噪处理，时域和频域图几乎相同，这说明噪声被完全滤掉，同时也说明FIR滤波器设计很理想，能满足课设要求。3.6 结果分析语音信号经过FIR滤波器的滤除噪声的处理，在Matlab中，函数sound可以对声音进行回放。其调用格式：sound (x，fs，bits)我们可以明显感觉滤波前后的声音有变化。声音中刺耳的声音没有了，几乎恢复成原始的声音，但较原始的声音更平滑一些。这说明用汉宁窗设计FIR滤波器滤掉了语音中的噪声同时，也把原始语音的很小的一部分也滤掉了，所以回
21、放语音的时候听起来比以前的更加平滑，说明这段程序设计是成功的。4出现的问题及解决方法 在设计课程设计时，出现了以下几个问题：1.、在编辑.wav音乐时，由于没有控制好所需音乐信号的时间，导致结果不理想。2、对各时域、频谱图的范围没有好的预计，导致出图时波形效果不理想。3、因为多次改动单频噪声频率值，所设计的滤波器的性能指标没有随之变化，出现了滤波上的错误。4、在确定单频噪声频率值后，由于不能很好的掌握其它参数的调试指标，导致多次调试都无法得到理想的滤波。5 结束语在此次课程设计中，我的任务是利用汉宁窗函数设计FIR滤波器，对吹管音乐信号进行滤波去噪处理。课设开始之前，我认真复习了有关窗函数特别
22、是汉宁窗的相关知识以及滤波器的设计方法，了解课设流程。在这两周里，我利用老师给出的模板，结合相关的专业知识，比较轻松的完成了课设的任务。不同于之前的理论课，虽然这次的课设内容并不是很难，但是仍然很考验我的动手实践能力。此间，我得到了以下收获：首先，在学习方面，虽然已经学习过了DSP课程，但是如何融洽的把理论和实际结合仍然是我需要面对的问题。其次，在matlab编程方面，由于老师给出了相关的程序资料，所以设计过程中并不算困难，可以说是顺利过关。最后，在团结合作方面，虽然每个人都单独分配了课题，但是总能找到与自己课题相关或者类似的同学，跟同学交流经验成为本次课设中的一个重要环节。在这次设计过程中，
23、既体现出了我自己单独设计的能力以及综合运用知识的能力，又让我体会到了学以致用、理论与实际贯通的喜悦，提高了我的团队协作能力。同时，我也从中也能发现自己平时学习的不足。在此，特别感谢指导老师高明，在您的指导下我成功完成了本次课设的任务,谢谢！参考文献1张威. MATLAB基础与编程入门. 西安：西安电子科技大学出版社，2008.12 张圣勤.MATLAB7.0实用教程.北京：机械工业出版社，20083 张志涌.精通MATLAB 6.5版M北京：北京航空航天大学出版社，20034 程佩青.数字信号处理教程.北京：清华大学出版社，20025维纳K英格尔，约翰G普罗克斯. 数字信号处理M. 西安：西安
24、交通大学出版社，2008.16张小虹.信号系统与数字信号处理M.第版.西安：西安电子科技出版社，2002附录1：语音信号滤波去噪设计源程序清单% 程序名称：DSPYFL.m% 程序功能：采用基于汉宁的窗口设计法，设计FIR滤波器对含噪语音进行滤波去噪处理。% 程序作者：余霏霖/% 最后修改日期：2013-3-8x,fs,bits=wavread(h:2013DSPPurpleBambooTune.wav); % 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。sound(x,fs,bits); % 按指定的采样
25、率和每样本编码位数回放N=length(x); % 计算信号x的长度fn=2100; % 单频噪声频率t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率x=x(:,1);y=x+sin(fn*2*pi*t); % 加入一个单频噪声sound(y,fs,bits); % 可以明显听出有尖锐的单频啸叫声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; % 计算频谱频率范围
26、figure(1)subplot(2,2,1);plot(t,x); % 布局为2*2的四个小图title(原始音乐信号);xlabel(时间(t);ylabel(幅度); %改变横纵坐标的范围axis(0,2,-1.5,1.5); %加上标题和横坐标名称grid on; % 加上网格subplot(2,2,2);plot(f,X);title(原始音乐信号频谱);xlabel(频率(f);ylabel(幅度谱);axis(0,3000,0,3000);grid on; subplot(2,2,3);plot(t,y);title(加入干扰后的音乐信号);xlabel(时间(t);ylabel
27、(幅度);axis(0,2,-1.5,1.5);grid on;subplot(2,2,4);plot(f,Y);title(加入单频干扰后的音乐信号频谱);xlabel(频率(f);ylabel(幅度谱);axis(0,3000,0,3000);grid on;fpd=1450;fsd=1650;fpu=2350;fsu=2250; % FIR滤波器的上下截止频率Rp=1;As=37; % 带阻滤波器设计指标fcd=(fpd+fsd)/2;fcu=(fpu+fsu)/2;df=min(fsd-fpd),(fpu-fsu); % 计算上下边带中心频率和频率间隔wcd=fcd/fs*2*pi;w
28、cu=fcu/fs*2*pi;dw=df/fs*2*pi; % 将Hz为单位的模拟频率换算为rad为单位的数字频率wsd=fsd/fs*2*pi;wsu=fsu/fs*2*pi;M=ceil(6.2*pi/dw)+1; % 计算汉宁窗设计该滤波器时需要的阶数n=0:M-1; % 定义时间范围w_ham=hanning(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,gr
29、d,w=freqz_m(h_bs,1); % 调用自编函数计算滤波器的频率特性figure(2)subplot(2,2,1);plot(w/pi,db);axis(0,0.4,-100,20);title(以db为单位的幅度特性);xlabel(w/pi);ylabel(db); grid on;subplot(2,2,2);plot(w/pi,mag);axis(0,0.4,-0.5,1.25);title(以线性为单位的幅度特性);xlabel(w/pi);ylabel(mag);grid on;subplot(2,2,3);plot(w,pha);title(滤波器相位响应图);xlab
30、el(w/pi);ylabel(相位（pha）);axis(0,3,-4,4);grid on;subplot(2,2,4);plot(h_bs);axis(0,800,-0.2,1);title(滤波器脉冲响应图);xlabel(n);ylabel(h(n);grid on; y_fil=filter(h_bs,1,y); % 用设计好的滤波器对y进行滤波Y_fil=fft(y_fil);Y_fil=Y_fil(1:N/2); % 计算频谱取前一半figure(3)subplot(3,2,1);plot(t,x); axis(0,2,-1.5,1.5);title(原始音乐信号时间x);xl
31、abel(时间（t）);ylabel(幅度);grid on;subplot(3,2,2);plot(f,X);axis(0,10000,0,1000);title(原始音乐信号幅度谱X);xlabel(频率（f）);ylabel(幅度);grid on;subplot(3,2,3);plot(t,y);axis(0,2,-1.5,1.5);title(加干扰音乐信号时间x1);xlabel(时间（t）);ylabel(幅度);grid on;subplot(3,2,4);plot(f,Y);axis(0,10000,0,1000);title(加干扰音乐信号幅度谱X1);xlabel(频率（
32、f）);ylabel(幅度);grid on;subplot(3,2,5);plot(t,y_fil);axis(0,2,-1.5,1.5);title(滤波后音乐信号时间y);xlabel(时间（t）);ylabel(幅度);grid on;subplot(3,2,6);plot(f,Y_fil);axis(0,10000,0,1000);title(滤波后音乐信号幅度谱Y);xlabel(频率（f）);ylabel(幅度);grid on;sound (y_fil,fs,bits);附录2：function hd = ideal_lp(wc,M);% 理想低通滤波器计算% -% hd =
33、ideal_lp(wc,M)% hd = 0 to M-1之间的理想脉冲响应% wc = 截止频率(弧度) % M = 理想滤波器的长度alpha = (M-1)/2;n = 0:1:(M-1);m = n - alpha + eps;hd = sin(wc*m) ./ (pi*m);附录3：function db,mag,pha,grd,w = freqz_m(b,a);% freqz 子程序的改进版本% -% db,mag,pha,grd,w = freqz_m(b,a);% db = 0 到pi弧度区间内的相对振幅(db)% mag = 0 到pi弧度区间内的绝对振幅% pha = 0
34、到pi弧度区间内的相位响应% grd = 0 到pi弧度区间内的群迟延% w = 0 到pi弧度区间内的501个频率样本向量% b = Ha(z)的分子多项式系数(对FIR b=h)% a = Ha(z)的分母多项式系数(对 FIR: a=1)% H,w = freqz(b,a,1000,whole);H = (H(1:1:501); w = (w(1:1:501);mag = abs(H);db = 20*log10(mag+eps)/max(mag);pha = angle(H);% pha = unwrap(angle(H);grd = grpdelay(b,a,w);% grd = diff(pha);% grd = grd(1) grd;% grd = 0 grd(1:1:500); grd; grd(2:1:501) 0;% grd = median(grd)*500/pi;
