音乐信号滤波去噪——使用CHEB窗设计的线性相位型的FIR滤波器.doc

音乐信号滤波去噪使用CHEB窗设计的线性相位型的FIR滤波器 学生姓名：胡国庆 指导老师：高明 摘 要 本课程设计主要使用CHEB窗设计的FIR滤波器对音乐信号进行滤波去噪。课程设计的的平台为MATLAB。采集一段音乐信号，绘制波形并观察其频谱，给定相应技术指标，用窗函数法设计一个满足指标的FIR滤波器，对该音乐信号进行滤波去噪处理。根据滤波后的时域图和原始音乐信号时域图的比较，以及滤波后信号的频谱图和原始音乐信号频谱图的比较，最后回放滤波后音乐信号，滤波后的音乐信号与原始音乐信号一样清晰，成功地实现了滤波，达到了设计的要求。 关键词 FIR滤波器；音乐信号；滤波去噪；CHEB窗设计；MATLAB 1 引 言本课程设计是采用CHEB窗设计的FIR滤波器对语音信号进行滤波去噪。通过课程设计了解FIR滤波器设计的原理和步骤，掌握用Matlab语言设计滤波器的方法，。通过观察音乐信号滤波前后的时域波形的比较，加深对滤波器作用的理解。通过对比滤波前后波形图的比较和放滤波前后音乐信号的对比 ，可以看出滤波器对有用信号无失真放大具有重大意义。1.1 课程设计目的熟悉Matlab语言环境，掌握Matlab语言的编程规则，利用CHEB窗函数设计法来设计符合要求的FIR滤波器来实现音乐信号的滤波去噪。并绘制滤波前后的时域波形和频谱图。根据图形分析判断滤波器设计的正确性。通过本次课程设计熟悉利用 CHEB窗函数法设计FIR滤波器的过程。增强自己独立解决问题的能力，提高自己的动手能力。加深对理论知识联系实际问题的理解。为以后的工作奠定坚实的基础。1.2 课程设计要求下载一段音乐信号，绘制观察波形及频谱图。根据CHEB的性能指标合理设计FIR滤波器，对音乐信号加入干扰，再用滤波器对干扰音乐信号进行滤波去噪，比较滤波前后的频谱图并进行分析。再回放音乐信号对比原音乐信号。看滤波器是否对音乐信号进行了滤波去噪。1.3课程设计平台MATLAB在信号与系统中的应用主要包括符号运算和数值计算仿真分析。由于信号与系统课程的许多内容都是基于公式演算，而MATLAB借助符号数学工具箱提供的符号运算功能能基本满足信号与系统课程的需求。例如，解微分方程、傅里叶正反变换、拉普拉斯正反变换、z正反变换等。MATLAB在信号与系统中的另一主要应用是数值计算与仿真分析，主要包括函数波形绘制、函数运算、冲激响应与阶跃响应仿真分析、信号的时域分析、信号的频谱分析、系统的S域分析、零极点图绘制等内容。数值计算仿真分析可以帮助学生更深入理解信号与系统的理论知识，并为将来使用MATLAB进行信号处理领域的各种分析和实际应用打下基础。2 设计原理2.1 FIR滤波器FIR滤波器：有限长单位冲激响应滤波器，是数字信号系统中最基本的元件，它可以在保证任意幅频特性的同时具有严格的线性相频特性，同时其单位抽样响应是有限长的，因而滤波器是稳定的系统。因此，FIR滤波器在通信、图像处理、模式识别等领域都有着广泛的应用。 FIR数字滤波器设计的基本步骤如下： （1）确定指标 在设计一个滤波器之前，必须首先根据工程实际的需要确定滤波器的技术指标。在很多实际应用中，数字滤波器常常被用来实现选频操作。因此，指标的形式一般在频域中给出幅度和相位响应。幅度指标主要以两种方式给出。第一种是绝对指标。它提供对幅度响应函数的要求，一般应用于FIR滤波器的设计。第二种指标是相对指标。它以分贝值的形式给出要求。（2）逼近确定了技术指标后，就可以建立一个目标的数字滤波器模型。通常用理想的数字滤波器模型。之后，利用数字滤波器的设计方法，设计出一个实际滤波器模型来逼近给定的目标。 （3）性能分析和计算机仿真上两步的结果是得到以差分或冲激响应描述的滤波器。根据这个描述就可以分析其频率特性和相位特性，以验证设计结果是否满足指标要求；或者利用计算机仿真实现设计的滤波器，再分析滤波结果来判断。 FIR滤波器的设计问题实际上是确定能满足所要求的转移序列或脉冲响应的常数的问题，设计方法主要有窗函数、频率响应法和等波纹最佳逼近法等。2.2 窗口设计法窗口设计法是一种通过截断和计权的方法使无限长非因果序列成为有限长脉冲响应序列的设计方法。通常在设计滤波器之前，应该先根据具体的工程应用确定滤波器的技术指标。在大多数实际应用中，数字滤波器常常被用来实现选频操作，所以指标的形式一般为在频域中以分贝值给出的相对幅度响应和相位响应。 窗口设计法步骤如下：（1）根据过渡带宽及阻带衰减要求，选择窗函数的类型并估计窗口长度N。窗函数的类型可根据最小阻带衰减AS独立选择。（2）根据待求滤波器的理想频率响应求出理想单位脉冲响应hd(n)。（3）由性能指标确定窗函数W(n)和长度N。（4）求得实际滤波器的单位脉冲响应h(n)， h(n)即为所设计FIR滤波器系数向量b(n)。 其性能如表2-1所示：表2-1 常见窗函数性能表2.3 切比雪夫窗切比雪夫又称Dolph-Chebyshev窗，它是由一个切比雪夫多项式在单位圆上作N点等间隔抽样，然后再作DFT反变化得到的。当给定窗的长度N，最大旁瓣衰减r（对于给定的旁瓣幅度.06hev低于主瓣的rdB分贝数。则其3dB带宽B的表达式为： （2.1）CHEB窗函数的时域幅度与频域幅度特性曲线的MATLAB实现的曲线图如图2-1所示。图2-1 CHEB窗函数的时域幅度与频域幅度特性曲线3 设计步骤3.1 设计流程图根据设计的要求，首先从网上下载一段音乐信号，对音乐信号加入噪声干扰，再利用CHEB窗设计合理的FIR滤波器。最后用滤波器对干扰后的音乐信号进行滤波去噪。具体设计流程图如下图3-1所示：图3-1 设计流程图3.2 录制音乐信号从网上下载一段音乐信号，并命名为“yinyue.wav”。其采样率fs为8000HZ,在Matlab平台下调用wavread函数可采集到下载的音乐信号。具体调用如下： x,fs,bits=wavread(e:yinyue.wav); % 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数 sound(x,fs,bits); %播放音乐信号 N=length(x);%计算音乐信号的长度 fn=3100; % 单频噪声频率 t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率 x=x; y=x+0. 1*sin(fn*2*pi*t); plot(t,x); xlabel(t ins);ylabel(x); % 画原始音乐信号的时域波形图 sound(y,fs,bits); % 应该可以明显听出有尖锐的单频啸叫声运行程序后得到的波形图如图3-2所示：图3-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; % 计算频谱频率范围 subplot(2,2,1);plot(t,x);xlabel(时间(t);ylabel(幅度); title(原始音乐信号); axis(0,15,-2,2)subplot(2,2,2);plot(f,X);xlabel(频率(f);ylabel(幅度谱);title(原始音乐信号幅度谱);axis(0,5000,0,1000); subplot(2,2,3);plot(t,y);xlabel(时间(t);ylabel(幅度);title(加干扰后的音乐信号); subplot(2,2,4);plot(f,Y);xlabel(频率(f);ylabel(幅度谱);title(加干扰后的音乐信号幅度谱);axis(0,5000,0,1000);运行程序后我们所得的波形图如图3-3所示：图3-3 音乐信号加入单频噪声前后的时域和频域波形图 wavwrite(y,fs,bits,e:zaoyin.wav); %运行写出噪音文件。由上图3-3我们可以看到语音信号加入单频噪声后的时域波形比未加之前在幅度范围内有了明显的增加，在频谱方面我们可以看到除了在加了噪声后的频谱图上的3100hz有个明显的冲激外，其余地方与未加时的频谱一模一样，这现象表现在语音播放时我们可以听到一声尖锐的噪声。3.3 滤波器设计滤波器的设计就是要找到一组能满足特定滤波要求的系数向量a和b，其中b为系统函数的分子系数，a为系统函数分母系数。在本次的课程设计中我所采用的就是利用CHED窗函数来设计FIR滤波器。设计的程序如下：fpd=2000;fsd=3000;fsu=3200;fpu=4200;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为单位的数字频率ws1=fsd/fs*2*pi;ws2=fsu/fs*2*pi;M=ceil(15.2*pi/dw)+1; % 计算CHEB窗设计该滤波器时需要的阶数n=0:M-1; % 定义时间范围w_cheb=(CHEBWIN(M); % 产生M阶的CHEB窗 hd_bs=IDEAL_LP(wcd,M)+ IDEAL_LP (pi,M)- IDEAL_LP (wcu,M); % 调用自编函数计算理想带阻滤波器的脉冲响应h_bs=w_cheb.*hd_bs; % 用窗口法计算实际滤波器脉冲响应db,mag,pha,grd,w=freqz_m(h_bs,1); % 调用自编函数计算滤波器的频率特性subplot(2,2,1);plot(w/pi,db);title(滤波器幅度响应图);xlabel(w/pi);ylabel(db);grid on;axis(0,1,-150,20);line(0,1.5,-Rp,-Rp,color,r,linestyle,-,LineWidth,2); line(0,1.5,-As,-As,color,r,linestyle,-,LineWidth,2); line(ws1/pi,ws1/pi,-140,20,color,r,linestyle,-,LineWidth,2); line(ws2/pi,ws2/pi,-140,20,color,r,linestyle,-,LineWidth,2);xlabel(w/pi);ylabel(db);title(以db为单位的滤波器幅度响应图);subplot(2,2,2);plot(w,mag);title(滤波器幅度响应图);xlabel(w/pi);ylabel(幅度mag); grid on;subplot(2,2,3);plot(w,pha);title(滤波器相位响应图);xlabel(w/pi);ylabel(相位pha);grid on; subplot(2,2,4);stem(n,h_bs);title(滤波器脉冲响图);xlabel(n);ylabel(h(n);axis(0,150,0,1);wp1=fpd/fs*2*pi;wp2=fpu/fs*2*pi;dww=2*pi/1000;RP=-min(db(1:1:ceil(wp1/dww)+1,ceil(wp2/dww)+1:1:end)AS=-max(db(ws1/dww+1:ws2/dww+1)图3-4 FIR滤波器的频率响应上图3-4为用CHEB窗函数法设计出的FIR滤波器图，我们可以看出，阻带最大衰减为-120dB,FIR滤波器的主瓣宽度很小,这样可以使过渡带很陡，旁瓣相对于主瓣也比较小。3.4信号滤波处理滤波器设计完成后，在MATLAB平台上用函数filter实现滤波，滤波程序如下：y_fil=filter(h_bs,1,y);% 用设计好的滤波器对y进行滤波Y_fil=fft(y_fil);Y_fil=Y_fil(1:N/2); % 计算频谱取前一半 subplot(3,2,1);plot(t,x);xlabel(时间(t);ylabel(幅度); title(原始音乐信号);axis(0,15,-2,2)subplot(3,2,2);plot(f,X);xlabel(频率(f);ylabel(幅度谱);title(原始音乐信号幅度谱);axis(0,5000,0,1000); subplot(3,2,3);plot(t,y);xlabel(时间(t);ylabel(幅度);title(加干扰后的音乐信号); subplot(3,2,4);plot(f,Y);xlabel(频率(f);ylabel(幅度谱);title(加干扰后的音乐信号幅度谱);axis(0,5000,0,1000); subplot(3,2,5);plot(t,y);xlabel(时间(t);ylabel(幅度); title(滤波后音乐信号); subplot(3,2,6);plot(f,Y_fil);xlabel(频率(f);ylabel(幅度谱);title(滤波后音乐信号幅度谱);axis(0,5000,0,1000);wavwrite(y_fil,fs,bits,e:quzao.wav);得到的滤波前后音乐信号的时域波形图和频谱图对比图如图3-5所示：图3-5 滤波前后的波形对比由上图3-5我们可以看出，加噪声的后的语音信号经过FIR滤波器的滤噪处理，时域和频域图都几乎完全一样，这说明噪声被完全滤掉，同时也说明FIR滤波器设计很理想，能满足所需要求。3.5结果分析在MATLAB中，经过sound(x,fs,bits)函数，对经过CHEB窗设计的FIR滤波器之后的音乐信号进行回放，可以听出滤波之后的信号跟原始信号一样清晰，完全滤除掉了噪声的干扰，通过下面语句来进行语音信号回放比较：sound(y_fil,fs,bits);%播放滤波后的音乐信号sound(x,fs,bits);%所得结果证明了用CHEB窗设计的FIR滤波器和音乐信号去噪设计是成功的。3.6滤波器结构设计本课程设计的滤波器类型为线性相位型FIR滤波器，在MATLAB中输入滤波器参数h_bs得到了相关系数以后，将滤波器结构用visio画好，填上系数。图3-6 相关系数图3-7 滤波器结构4出现的问题及解决方法在本次课程设计中我遇到的问题如下：1、不知道如何下载符合要求长度和格式的音乐信号。2、对利用CHEB窗函数设计FIR滤波器的设计步骤很生疏。3、在采用CHEB窗函数设计的FIR滤波器时得不到理想的滤波器。4、绘制滤波器结构时不知道As和Rp是否达标。5、绘制滤波器时因阶数较大导致工作量较大。针对以上问题，相应的解决方案如下：1、通过百度学习，然后借助酷狗音乐的铃声制作和格式，最后利用windows自带的录音软件转化其采样率fs为8000HZ的单频音乐信号。2、自己上网查阅资料，向同学请教，同时参考胡老师整理的模板，慢慢掌握利用CHEB窗函数设计FIR滤波器的方法和步骤。3、请教老师，老师帮我检查了之前的程序代码，没有发现错误后调整fpd、fpu等参数，通过不断设置参数的值，最终达到最理想的值。4、向同学寻求帮助，并参考胡老师整理的常见窗函数性能表，确定了As和Rp的理想值，从而顺利完成滤波器结构的绘制。5、通过调整fpd和fpu的值，减少了阶数，从而减少了工作量。5 结束语通过本次课程设计，总体上我对所学的数字信号处理有了更深的了解，对MATLAB的使用也更为熟练，同时提升我的自学能力。在课程设计之前，我们虽然学过数字信号处理，但我对大部分内容理解不够，对其运用更不了解，但通过本次课程设计，我对所学的课本知识有了更深的认识和理解，并且学会了运用。在做课设的过程中，胡老师总结的模板和程序对我而言有很大的帮助，减少了我的工作量，提升了效率，在此衷心感谢胡老师为我们认真总结的模板。在课设过程中我遇到的问题比较多，也比较杂，主要是通过百度和问同学的方式解决。最后，在团结合作能力方面，虽然每个人有每个人的题目，但是我们还是充分发挥了我们的团结合作能力，找到资料的同学帮没找到的同学找，有什么错误的大家一起讨论，一起找资料解决，特别是题目相似的同学，大家有资料都一起研究，而且还有我们的指导老师胡老师的精心指导，所以我们合作得很愉快。在这次设计过程中，既体现出了自己单独设计的能力以及综合运用知识的能力，又体会了学以致用、突出自己劳动成果的喜悦心情，从中也能发现自己平时学习的不足和薄弱环节，从而加以弥补。同时，也再次体会到了团结合作的快乐。在此我衷心感谢各位老师及同学给我的帮助，这是我能完成课程设计的主要支柱。通过此次课程设计我学到了更多的东西，学会了知识的实际运用，提升了自学的能力。参考文献1张圣勤.MATLAB7.0实用教程M.北京：机械工程出版社，2006 2 程佩青.数字信号处理.第四版.北京：清华大学出版社,2013年3 陈后金.数字信号处理.第三版.北京：高等教育业出版社，2004年4 刘敏,魏玲.MATLAB通信仿真与应用.第二版.北京：国防工业出版社，2001年5百科ROBOT,zivenwong.窗函数.百度百科,/view/2327611.htm?fr=ala0_1：2017-02-27附录1：用CHEB窗设计FIR滤波器程序清单%程序名称：FIR.M%程序功能：采用CHEB窗设计FIR滤波器的方法，给语音出去噪声。%程序作者：胡国庆%最后修改日期：2017-2-27x,fs,bits=wavread(e:yinyue.wav); % 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数sound(x,fs,bits); %播放音乐信号 N=length(x);%计算音乐信号的长度 fn=3100; % 单频噪声频率 t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率 x=x; y=x+0.1*sin(fn*2*pi*t); plot(t,x); xlabel(t ins);ylabel(x); % 画原始音乐信号的时域波形图 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; % 计算频谱频率范围 subplot(2,2,1);plot(t,x);xlabel(时间(t);ylabel(幅度); title(原始音乐信号); axis(0,15,-2,2)subplot(2,2,2);plot(f,X);xlabel(频率(f);ylabel(幅度谱);title(原始音乐信号幅度谱);axis(0,5000,0,1000); subplot(2,2,3);plot(t,y);xlabel(时间(t);ylabel(幅度);title(加干扰后的音乐信号); subplot(2,2,4);plot(f,Y);xlabel(频率(f);ylabel(幅度谱);title(加干扰后的音乐信号幅度谱);axis(0,5000,0,1000);fpd=2000;fsd=3000;fsu=3200;fpu=4200;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为单位的数字频率ws1=fsd/fs*2*pi;ws2=fsu/fs*2*pi;M=ceil(15.2*pi/dw)+1; % 计算CHEB窗设计该滤波器时需要的阶数n=0:M-1; % 定义时间范围w_cheb=(CHEBWIN(M); % 产生M阶的CHEB窗 hd_bs=IDEAL_LP(wcd,M)+ IDEAL_LP (pi,M)- IDEAL_LP (wcu,M); % 调用自编函数计算理想带阻滤波器的脉冲响应h_bs=w_cheb.*hd_bs; % 用窗口法计算实际滤波器脉冲响应db,mag,pha,grd,w=freqz_m(h_bs,1); % 调用自编函数计算滤波器的频率特性subplot(2,2,1);plot(w/pi,db);title(滤波器幅度响应图);xlabel(w/pi);ylabel(db);grid on;axis(0,1,-150,20);line(0,1.5,-Rp,-Rp,color,r,linestyle,-,LineWidth,2); line(0,1.5,-As,-As,color,r,linestyle,-,LineWidth,2); line(ws1/pi,ws1/pi,-140,20,color,r,linestyle,-,LineWidth,2); line(ws2/pi,ws2/pi,-140,20,color,r,linestyle,-,LineWidth,2);xlabel(w/pi);ylabel(db);ti