吹管音乐滤波去噪——使用脉冲响应不变法设计的巴特沃斯滤波器.doc

吹管音乐滤波去噪-使用脉冲响应不变法设计的巴特沃斯滤波器第页共19页 吹管音乐滤波去噪使用脉冲响应不变法设计的巴特沃斯滤波器学生姓名： 指导老师： 摘 要 本课程设计主要内容是设计利用脉冲响应不变设计法选设计一个IIR巴特沃斯滤波器，对一段含噪吹管音乐信号进行滤波去噪处理并根据滤波前后的波形和频谱分析滤波性能。本课程设计仿真平台为MATLAB7.0，开发工具是M语言编程。首先在网上下载一段吹管音乐信号，并人为加入一单频噪声，然后对信号进行频谱分析以确定所加噪声频率为3100Hz，并设计滤波器进行滤波去噪处理，最后比较滤波前后的波形和频谱基本相同。由分析结果可知，滤波器后的吹管音乐信号与原始信号基本一致，即设计的IIR巴特沃斯滤波器能够去除信号中所加单频噪声，达到了设计目的。 关键词 课程设计；滤波去噪；IIR滤波器；脉冲响应不变法；巴特沃斯；MATLAB1 引 言本课程设计主要解决在含噪情况下对吹管音乐信号的滤波去噪处理，处理时采用的是利用脉冲响应不变设计法设计的IIR巴特沃斯滤波器1。1.1 课程设计目的课程设计是教学的最后一个步骤，课程设计有利于基础知识的理解，我们掌握了基础知识和基本技能，但是要接触才能真正理解课程的深入部分；还有利于逻辑思维的锻炼，在许多常规学科的日常教学中，我们不难发现这样一个现象，不少学生的思维常常处于混乱的状态，写起作文来前言不搭后语，解起数学题来步骤混乱，这些都是缺乏思维训练的结果，所以我们可以通过实践来分析问题、解决问题、预测目标等目的；同时也有利于与其他学科的整合，例如我们这次的课程设计就要运用MATLAB软件的帮助才能实现；最重要的有利于治学态度的培养，在课程设计中，我们可能经常犯很多小错误，可能要通过好几次的反复修改、调试才能成功，但这种现象会随着学校的深入而慢慢改观。这当中就有一个严谨治学、一丝不苟的科学精神的培养。数字信号处理课程设计是在学生完成数字信号处理和MATLAB的结合后的基本实验以后开设的。本课程设计的目的是为了让学生综合数字信号处理和MATLAB并实现一个较为完整的小型滤波系统。这一点与验证性的基本实验有本质性的区别。开设课程设计环节的主要目的是通过系统设计、软件仿真、程序安排与调试、写实习报告等步骤，使学生初步掌握工程设计的具体步骤和方法，提高分析问题和解决问题的能力，提高实际应用水平。1.2 课程设计的要求（1）滤波器指标必须符合工程设计。（2）设计完后应检查其频率响应曲线是否满足指标。（3）处理结果和分析结论应该一致，而且应符合理论。（4）独立完成课程设计并按要求编写课程设计报告。1.3 设计平台MATLAB是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平2。MATLAB作为一种科学计算的高级语言之所以受欢迎，就是因为它有丰富的函数资源和工具箱资源，编程人员可以根据自己的需要选择函数，而无需再去编写大量繁琐的程序代码，从而减轻了编程人员的工作负担，被称为第四代编程语言的MATLAB最大的特点就是简洁开放的程序代码和直观实用的开发环境。MATLAB和Mathematica、Maple并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域3。2 设计原理在网上采集一段吹管音乐信号，绘制波形并观察其频谱，给定相应技术指标，用脉冲响应不变法设计一个满足指标的巴特沃斯IIR滤波器，对该吹管音乐信号进行滤波去噪处理，比较滤波前后的波形和频谱并进行分析。2.1 IIR滤波器IIR滤波器具有无限长脉冲响应，因此能够与模拟滤波器相匹敌；一般来说，所有的模拟滤波器都有无限长脉冲响应。因此，IIR滤波器设计的基本方法是利用复值映射将大家熟知的模拟滤波器变换为数字滤波器。这一方法的优势在于各种模拟滤波器设计（AFD）表格和映射在文献中普遍都能获得。这个基本方法称为A/D（模拟-数字）滤波器变换。然而，AFD表格仅对低通滤波器适用，而同时也想要得到设计其他频率选择性滤波器（高通、带通、带阻等等）。为此，需要对低通滤波器实行频带变换这些变换也是复值映射，在各种文献中也能得到。这种IIR滤波器设计的基本方法存在两种途径：设计模拟低通滤波器途径1：实行滤波器变换SZ实行频带变换S S期望的IIR滤波器途径2：设计模拟低通滤波器实行频带变换ZZ实行滤波器变换SZ 期望的IIR滤波器在MATLAB中采用第1种途径设计IIR滤波器。这些MATLAB函数的直接使用并没有给出有关任何设计方法的细节。因此，本课程设计将研究第2种途经，因为它涉及数字域的频带变换。这样，在这种IIR滤波器设计方法将按下列步骤进行：（1）设计模拟低通滤波器。（2）研究并实行滤波器变换以得到数字低通滤波器。（3）研究并实行频带变换以便从数字低通滤波器得到其他数字滤波器。这些途经存在的主要问题是在IIR滤波器的相位特性上没有一点控制能力，因此IIR滤波器设计将仅作为幅度设计对待4。2.2 巴特沃斯滤波器这个滤波器是用这个性质表征的：它的幅度响应在通带和阻带都是平坦的。一个N阶低通滤波器的幅度平方响应给出为 (2.1)式中N是滤波器的阶，（rad/s）是截止频率。幅度平方响应的图如下图2-1所示。图2-1 巴特沃兹filter 振幅平方函数图中横坐标为，纵坐标为，交点处为。从这张图2-1可以看出下面几个性质：（1） =0，等于=1，对全部N。（2） =， =，对全部N，这意味着在有3dB衰减。（3） 是的单调下降的函数。（4） 随N向一个理想低通滤波器趋近。（5） 在=0是最大平坦，因为在这里所有阶的导数存在且等于零。为了确定系统函数,现将(2.1)式变换得到 (2.2)由（2.2）式分母多项式的根（或的极点）给出为 (2.3)式（2.3）的一种解释是：（1）总共有2N个极点，它们均匀分布在半径为的图上，相隔/N弧度。（2）对N为奇数，极点给出为 。（3）对N为偶数，极点给出为 。（4）极点对j轴是对称分布的。（5）极点永远不会落在虚轴上，且仅当N为奇数时才会落在实轴上5。通过选取在左半平面的极点就能给出一个稳定和因果的，并且能将写成 (2.4) 对设计滤波器这个过程的设计步骤是：（1） 选取T，这里设置T=1。并确定模拟频率如下： 和 (2.5) （2） 利用设参数，和设计一个模拟滤波器。这一步要用巴特沃斯滤波器为原型来完成。（3） 利用部分分式展开，将展开为： (2.6) （4） 现在将模拟极点变换为数字极点，得到数字滤波器如下式： (2.7) 并作化简得出作为z-1有理函数的。这样滤波器就可设计好了。2.3 脉冲响应不变法在这种设计方法中，我们想要这个数字滤波器的脉冲响应看起来与一个频率选择性模拟滤波器的冲激响应时“相似”的。为此，以某个采样间隔T对采样得到,即 (2.8)参数T要选成以使得的形状被它的样本“捕获”住。因为这是一种采样运算，所以模拟和数字频率由下式联系： 或 (2.9)由于是在单位圆上，是在虚轴上，所以有下面从s平面到z平面的变换 (2.10)系统函数和是经由频域混叠公式联系的：在（2.10）公式的映射关系下，复平面的变换如图2-2所示6。图2-2 脉冲响应不变法中的复平面映射3. 设计步骤3.1设计流程图开始吹管音乐信号滤波去噪使用脉冲不变响应法设计的巴特沃斯滤波器的设计流程如图3-1所示：吹管音乐信号的采集（wavread函数），加入单频噪声，画加噪前后时域图进行快速傅里叶变换，并且画加噪前后频谱图设定滤波器性能指标，单频噪声频率fn=3100，通带截止频率fp=fn-400=2700，阻带截止频率fc=fn-50=3050，通带波纹Rp=3，阻带波纹As=16.用脉冲响应不变法设计巴特沃斯滤波器验证滤波器是否符合指标并进行频谱分析设计好的滤波器进行滤波处理，回放滤波后的吹管音乐比较滤波前后吹管音乐信号的波形及频谱 结束图3-1 脉冲响应不变法巴特沃斯滤波器对吹管信号去噪流程图3.2吹管音乐信号的采集在网上下载一段吹管音乐信号，属性设置为8000Hz，8位单声道的音频格式。然后在Matlab软件平台下，利用函数wavread对语音信号进行采样，函数为x,fs,bits=wavread(e:dizi.wav)，记住采样频率和采样点数，如图3-2：图3-2 吹管音乐信号设置3.3吹管音乐信号的频谱分析画出加噪前后信号的时域波形和频域波形，程序截图3-3如下：图3-3程序截图图3-4 加噪前后时域和频域波形对比从图中可以看出时域波形加噪前后变化不是很明显，但变化是肯定的。频域波形加入单频噪声之后，在噪声频率为fn=3100Hz处很明显加入了噪声，频谱变化显而易见。3.4滤波器设计将数字滤波器的设计指标设为通带截止频率fn=3100Hz,fp=fn-400=2700Hz,阻带频率fc=fn-50=3050Hz,通带波纹Rp=3，阻带波纹As=16dB,要想要通过首先设计一个等效的模拟滤波器，然后再将它映射为所期望的数字滤波器来确定H(z)。根据前面原理和步骤，按要求用脉冲响应不变法设计好了巴特沃斯滤波器。程序截图3-5如下所示：图3-5 程序截图设计得到的滤波器其幅度响应，相位响应，脉冲响应如图3-6下所示： wp/pians = 0.6750 ws/pians = 0.7625图3-6 设计得到的巴特沃斯滤波器的幅度、相位和脉冲响应由图3-6还可以看出，最小阻带衰减大约是16dB，比16dB稍微大一点，设计的巴特沃斯滤波器达到指标As=16dB的要求。下图3-7为程序运行时的窗口显示，可以看出此巴特沃斯滤波器的阶数为16阶。图3-7 程序运行截图3.5 信号滤波处理前面已经用脉冲响应不变法设计好了我们要的巴特沃斯滤波器，接着就对吹管信号进行滤波处理，看自己设计的巴特沃斯滤波器有没有对吹管信号进行了滤波处理。所以就用filter函数进行滤波，即y_fil=filter(b,a,y)。我们将滤波前后的时域波形进行比较，并对其进行快速傅里叶变换，Y_fil=fft(y_fil);Y_fil=Y_fil(1:N/2)，对比前后的频谱，具体分析设计的滤波器是否达到设计要求。如图3-9为滤波前后的吹管音乐信号时域和频域波形对比，程序截图如图3-8所示：图3-8 程序截图图3-9 吹管音乐脉冲响应不变法巴特沃斯滤波器滤波前后时域和频域波形对比图3.6 结果分析我们先采集吹管音乐信号，再按照步骤用脉冲响应不变法设计巴特沃斯滤波器。由图3-7可知，设计的阻带截止频率wc，纵坐标差不多刚好在As处，所以设计的滤波器达到要求。我们观察到图3-9滤波前后吹管音乐信号的波形对比图，加噪之后的时域波形是有变化的，只是不明显。可能是因为我们采集的吹管信号噪声不是很大；但是可以看到在频域波形中，很明显地反应出设计的滤波器滤掉了我们采集的吹管音乐信号中的噪声。比较滤波前后的波形和频谱并进行分析，根据结果和学过的理论得出运用脉冲响应不变法设计的巴特沃斯滤波器达到了设计要求。与不同信源（如拉弦音乐信号），但相同滤波方法的同学比较可以得出虽然信源不同，拉弦音乐的幅度比较小，且变化幅度相对于吹管音乐比较小，但都不能完全滤除这两种音乐的噪声，滤波效果基本相同，如图3-10。与相同信源但不同滤波方法（如频率采样法FIR滤波器）的同学比较可以得出频率采样法FIR滤波器的滤波效果更好，如图3-11。图3-10拉弦音乐脉冲响应不变法巴特沃斯滤波器滤波前后时域和频域波形对比图图3-11吹管音乐频率采样法FIR滤波器滤波前后时域和频域波形对比图4.出现的问题及解决方法在这次的课程设计中我们设计的滤波器对吹管信号的滤噪处理或多或少会因为我们的疏忽或者实践能力不够还是会遇到问题的，当遇到问题时我们不但运用自己学到的知识还得到同学和老师的帮助将其解决，不断完善和修改。设计中出现的问题及解决方法如下：（1）在采集吹管音乐信号时，从网上下载下来的吹管音乐是立体音，需要自己重新设置一下，下载的吹管音乐音频格式应该设置为PCM 8000kHz，8位，单声道。如果不设置，频谱为蓝色和绿色混叠。下载音乐的时候有疑问，如果吹管音乐不是纯吹管的音乐，而是混有其他类型的音乐，会不会对比较的结果有影响。（2）由于刚开始设计出的巴特沃斯滤波器无法滤除加噪信号频谱中的噪声，滤除噪音后的吹管音乐感觉噪声更大了。向老师询问后才知道原来是因为fp=fn-400，fc=fn-50，定义通带和阻带截止频率的时候不合理，还有设置T=1，而不是T =1/fs，经更改后，就能滤波了，滤波后音乐听起来基本上与原来的吹管音乐无异。（3）还要在滤波器性能曲线的wp,ws处画一根竖线，这样更方便看出结果，其中wp,ws处线的确定还需计算出wp/pi,ws/pi的值。将wp/pi，ws/pi处以及As处的线画出来，更容易看出滤波器是否已经符合指标。（4）由滤波后的频谱与滤波前的频谱对比可知，单频噪声没有完全滤除，但是基本上已经滤除。更改巴特沃斯滤波器的参数和其他设置也不能完全滤除噪声，总体上已经满足要求。5.结束语这次DSP课程设计，不但巩固了以前学的知识，还学习了新内容，收获颇丰。原来滤去语音中的噪声我们运用自己学到的知识能做到，并更加巩固了MATLAB软件的使用。熟练掌握MATLAB软件的使用是非常重要的，整个课程设计过程中，我们都要用到它。在这个课程设计中我们须利用MATLAB设计滤波器。这次做的滤波器要滤去语音信号中的噪声，觉得很有意思，做完后很有成就感。其实用MATLAB软件做实验是要细心的，因为很多的语法和常量变量的定义我们都要仔细，一个不小心看错了或者输入不认真是容易出错误。在发现错误时，问了同学也问了老师。当发现错误原因时，很多次都是因为不细心造成。经过思考并仔细检查，最后得以编译成功。在做实验时一定要小心谨慎，也许一个小小的语法错误和常量变量的定义的错误就造成整个程序出现问题，得不到所需的波形和频谱图，导致实验结果不正确，或者得不出结果。经过这次的DSP课程设计，将自己学到的理论知识运用到实际中，提高了自己的动手能力和思维能力。在课程设计中发现自己的不足，所以在今后的学习中要更加注重实践，学习好专业知识并要能运用到实际。这次的设计中老师提出的几个问题，在老师的指导和同学的帮助下解决了，当然也少不了自己的努力和思考。将自己平时学到的东西能运用到实际中，让理论和实际得以结合还是很不错的。我在课程设计中找到了动手的乐趣和思考的快乐，有成就感。我从这次的课程设计中体会到，什么事都要我们真真正正用心去做了，才会使自己更加成长，没有学习就不可能有实践的机会，没有自己的实践就不会有所突破。最后，在此衷心地感谢老师和同学对我的帮助，也感谢学校给我们机会，让我们能够将自己学到的知识运用到实际中。参考文献1 张圣勤MATLAB7.0实用教程M北京：机械工程出版社，2006 2 张威MATLAB基础与编程入门M第2版西安：西安电子科技大学出版社，20083 (美)维纳K英格尔 ，（美）约翰G普罗克斯（著）刘树棠（译）数字信号处理（MATLAB版）M第2版西安：西安交通大学出版社，2008 4 程佩青数字信号处理教程M北京：清华大学出版社，20025 王大伦，王志新，王康数字信号处理M北京：清华大学出版社，20106 张立材，王民数字信号处理M北京：人民邮电出版社，2008附录1：吹管音乐信号滤波去噪设计源程序清单% 程序名称：dizi.m% 程序功能：采用脉冲响应不变法设计的巴特沃斯滤波器，对吹管音乐滤波去噪。% 程序作者：陈莹莹% 最后修改日期：2012-03-02%程序分为三段，以下为第一段。第一段程序主要是产生加噪的吹管音乐信号，并绘制加噪前后吹管音乐信号的时间波形图和频谱图(图3-4)。x,fs,bits=wavread(e:dizi.wav); % 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放N=length(x); % 计算信号x的长度t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率fn=3100; % 单频噪声频率x=x(:,1); y=x+0.1*sin(fn*2*pi*t); %形成加干扰信号sound(y,fs,bits); % 应该可以明显听出有尖锐的单频啸叫声X=abs(fft(x); % 对原始吹管信号进行fft变换,取幅度谱Y=abs(fft(y); % 对加噪吹管信号进行fft变换,取幅度谱subplot(2,2,1);plot(t,x); %画原始吹管信号图xlabel(时间（单位： s）);ylabel(幅度); title(原始吹管信号) ; %标明横纵坐标，标题subplot(2,2,3);plot(t,y); %画加噪的原始吹管信号图xlabel(时间（单位： s）);ylabel(幅度); title(加入单频干扰后的吹管信号) ; %标明横纵坐标，标题deltaf=fs/N; % 计算频谱的谱线间隔f=0:deltaf:fs/2-deltaf; % 计算频谱频率范围X=X(1:N/2); Y=Y(1:N/2); % 截取前半部分频谱subplot(2,2,2);plot(f,X); %画原始吹管信号幅度谱xlabel(频率（单位： Hz）);ylabel(幅度谱); title(原始吹管信号幅度谱图) ; %标明横纵坐标，标题subplot(2,2,4);plot(f,Y); %画加噪吹管信号幅度谱xlabel(频率（单位： Hz）);ylabel(幅度谱); title(加入单频干扰后的吹管信号幅度谱图) ; %标明横纵坐标，标题%以下为第二段。第二段程序主要是根据吹管音乐信号的频谱用脉冲响应不变法设计合格的巴特沃斯滤波器(图3-6)，并进行滤波处理。T=1; % Set T=1,定义采样间隔fn=3100; % 单频噪声频率fp=fn-400;fc=fn-50; %定义通带和阻带截止频率Rp=3;As=16; % 定义通带波纹和阻带衰减wp=fp/fs*2*pi;ws=fc/fs*2*pi; %计算对应的数字频率OmegaP=wp/T;OmegaS=ws/T; %截止频率线性变换Ripple=10(-Rp/20);Attn=10(-As/20); %把dB转换成绝对形式cs,ds=afd_butt(OmegaP,OmegaS,Rp,As); % 计算滤波器系统函数分子分母系数b,a=imp_invr(cs,ds,T); % 脉冲响应不变法得到数字滤波器系统函数分子分母系数db,mag,pha,grd,w=freqz_m(b,a); %验证滤波器是否达到指定性能delta=1,zeros(1,99);ha=filter(b,a,delta); % 计算脉冲响应Subplot(221);plot(w/pi,db) %画滤波器幅度特性（dB)xlabel(w/pi);ylabel(dB); title(滤波器幅度响应图) ;grid %标明横纵坐标，标题Subplot(222);plot(w/pi,mag) %画滤波器幅度特性（线性)xlabel(w/pi);ylabel(幅度mag); title(滤波器幅度响应图) ;grid %标明横纵坐标，标题Subplot(223);plot(w/pi,pha) %画滤波器相位特性xlabel(w/pi);ylabel(相位pha); title(滤波器相位响应图) ;grid %标明横纵坐标，标题Subplot(224);plot(ha) %画脉冲响应xlabel(n);ylabel(h(n); title(滤波器脉冲响应图) ;grid %标明横纵坐标，标题y_fil=filter(b,a,y); % IIR滤波器对信号signal进行滤波处理，其中b,a为上面设计好的滤波器参数。sound (y_fil,fs,bits); %滤波后吹管音乐回放%以下为第三段。第三段程序主要是比较滤波前后吹管音乐信号的波形及频谱(图3-9)。x,fs,bits=wavread(e:dizi.wav); % 输入参数为文件的全路径和文件名，N=length(x); % 计算信号x的长度t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率X=abs(fft(x);