吹管音乐滤波去噪—使用BARTHANNWIN窗设计的FIR滤波器.doc

大学数字信号处理课程设计报告王二 校徽学 院 XXX 专 业 通信工程 班 级 XXX 学 号 XX 学生姓名 王二 指导教师 XXX 课程成绩 完成日期 2012年3月2日 课程设计成绩评定学 院 XXX 专 业 通信工程 班 级 XXX 学 号 XX 学生姓名 王二 指导教师 XXX 完成日期 XXXX 年 X月 X 日 指导教师对学生在课程设计中的评价评分项目优良中及格不及格课程设计中的创造性成果学生掌握课程内容的程度课程设计完成情况课程设计动手能力文字表达学习态度规范要求课程设计论文的质量指导教师对课程设计的评定意见综合成绩 指导教师签字 XXXX年X月 X日课程设计任务书 XXX 系 XXX 专业 课程名称数字信号处理课程设计时间XXXXXX学年第X学期XX周学生姓名王二指导老师XXX题 目吹管音乐滤波去噪使用BARTHANNWIN设计的FIR滤波器主要内容： 从网站上下载一段吹管乐器演奏音乐，绘制波形并观察其频谱特点，加入一个带外单频噪声，用BARTHANNWIN窗设计一个满足指标的FIR滤波器，对该含噪音乐信号进行滤波去噪处理，比较滤波前后的波形和频谱并进行分析，根据结果和学过的理论得出合理的结论。与不同信源相同滤波方法的同学比较各种信源的特点，与相同信源不同滤波方法的同学比较各种滤波方法性能的优劣。要求：（1）滤波器指标必须符合工程实际。（2）设计完后应检查其频率响应曲线是否满足指标。（3）处理结果和分析结论应该一致，而且应符合理论。（4）独立完成课程设计并按要求编写课程设计报告书。应当提交的文件：（1）课程设计学年论文。（2）课程设计附件（主要是源程序）。校徽 王二 吹管音乐滤波去噪使用BARTHANNWIN窗设计的FIR滤波器26共吹管音乐滤波去噪使用BARTHANNWIN窗设计的FIR滤波器学生姓名：王二 指导老师：XXX摘 要 本课程设计主要是利用Barthannwin窗设计满足指标的FIR滤波器对一段含噪吹管乐器音频进行滤波去噪处理。以Matlab为设计平台，利用wavread函数对音频信号进行采样，获取参数，对音频信号进行读取和加噪，根据相应技术指标，设计Barthannwin FIR滤波器，对该音频信号进行滤波去噪处理，还原原始信号。通过比较滤波前后的波形和频谱图并进行分析，回放滤波后音乐信号，滤波后音乐信号与原始音乐信号一样清晰。本课程设计成功地对吹管音乐音频信号进行了滤波去噪，初步完成了设计指标。关键词 课程设计；滤波去噪；FIR滤波器；Barthannwin窗；MATLAB目录1 引 言31.1课程设计目的31.2 课程设计的要求41.3 设计平台42 设计原理52.1 FIR滤波器52.2窗口设计法62.3 BARTHANNWIN窗83 设计步骤93.1设计流程图93.2录制语音信号103.3 滤波器设计113.4 信号滤波处理133.5 结果分析144 出现的问题及解决方法165 结束语17参考文献19附录1：吹管音乐音频信号提取和加噪源程序清单20附录2：使用BARTHANNWIN窗设计的FIR滤波器源程序清单22附录3：信号滤波处理源程序清单24附录4：freqz_m的定义源程序清单25附录5：ideal_lp的定义源程序清单261 引 言本课程设计是采用Barthannwin窗设计的FIR滤波器对含噪吹管音乐信号进行滤波去噪。通过课程设计了解FIR滤波器设计的原理和步骤，掌握用MATLAB语言设计滤波器的方法，了解DSP对FIR滤波器的设计及编程方法。通过观察音乐信号滤波前后的时域波形的比较，加深对滤波器作用的理解。通过对比滤波前后波形图的比较和放滤波前后音乐信号的对比，可以看出滤波器对有用信号无失真放大具有重大意义。1.1课程设计目的在本次课程设计中，最主要的设计是设计FIR滤波器，FIR滤波器的设计方法主要分为两类：第一类是基于逼近理想滤波器器特性的方法包括窗函数法、频率采样法、和等波纹最佳逼近法；第二类是最优设计法。本次的课程设计主要采用的是第一类设计方法，是利用Barthannwin窗函数法设计FIR滤波器对一段语音进行滤波去噪，并绘制滤波前后的时域波形形和频谱图，对滤波前后波形进行对比分析得到结论。通过这一过程能够熟悉MATLAB语言环境，掌握MATLAB语言的编程规则。通过本次课程设计熟悉利用Barthannwin窗函数法设计FIR滤波器的过程增强自己独立解决问题的能力，提高自己独立解决问题的动手能力同时加深对理论知识联系实际的理解，为以后的工作奠定坚实的基础。1.2 课程设计的要求下载一段吹管乐器演奏音乐，绘制观察波形及其频谱图特点，加入一个带外单频噪声，用Barthannwin窗设计一个满足指标的FIR滤波器，对该含噪音乐信号进行滤波去噪处理，比较滤波前后的波形和频谱图并进行分析。再回放音乐信号对比原音乐信号，看滤波器是否对音乐信号进行了滤波去噪。（1）滤波器指标必须符合工程实际。（2）设计完后应检查其频率响应曲线是否满足指标。（3）处理结果和分析结论应该一致，而且应符合理论。（4）独立完成课程设计并按要求编写课程设计报告书。1.3 设计平台本课程设计的主要设计平台式MATLAB7.0，MATLAB是矩阵实验室（Matrix Laboratory）的简称。MATLAB是由美国MathWorks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。MATLAB和Mathematica、Maple并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域1。MATLAB在信号与系统中的应用主要包括符号运算和数值计算仿真分析。由于信号与系统课程的许多内容都是基于公式演算，而MATLAB借助符号数学工具箱提供的符号运算功能能基本满足信号与系统课程的需求。例如，解微分方程、傅里叶正反变换、拉普拉斯正反变换、z正反变换等。MATLAB在信号与系统中的另一主要应用是数值计算与仿真分析，主要包括函数波形绘制、函数运算、冲激响应与阶跃响应仿真分析、信号的时域分析、信号的频谱分析、系统的S域分析、零极点图绘制等内容。数值计算仿真分析可以帮助学生更深入理解信号与系统的理论知识，并为将来使用MATLAB进行信号处理领域的各种分析和实际应用打下基础2。2 设计原理2.1 FIR滤波器数字滤波器根据其冲激响应函数的时域特性，可分为2种，即无限长冲激响应（IIR）滤波器和有限长冲激响应（FIR）滤波器。FIR(Finite Impulse Response)滤波器：有限长单位冲激响应滤波器，是数字信号系统中最基本的元件，它可以在保证任意幅频特性的同时具有严格的线性相频特性，同时其单位抽样响应是有限长的，因而滤波器是稳定的系统。因此，FIR滤波器在通信、图像处理、模式识别等领域都有着广泛的应用。FIR系统有自己突出的优点：（1）可以很容易地设计线性相位的滤波器。线性相位滤波器延时输入信号，却并不扭曲其相位。（2）实现简单。在大多数DSP处理器，只需要对一个指令积习循环就可以完成FIR计算。（3）适合于多采样率转换,它包括抽取(降低采样率)，插值(增加采样率)操作。无论是抽取或者插值，运用FIR滤波器可以省去一些计算，提高计算效率。相反，如果使用IIR滤波器，每个输出都要逐一计算，不能省略,即使输出要丢弃。（4）具有理想的数字特性。在实际中，所有的DSP滤波器必须用有限精度（有限bit数目）实现，而在IIR滤波器中使用有限精度会产生很大的问题，由于采用的是反馈电路，因此IIR通常用非常少的bit实现，设计者就能解决更少的与非理想算术有关的问题。（5）可以用小数实现。不像IIR滤波器，FIR滤波器通常可能用小于1的系数来实现。（如果需要，FIR滤波器的总的增益可以在输出调整）。当使用定点DSP的时候，这也是一个考虑因素，它能使得实现更加地简单3。FIR数字滤波器设计的基本步骤如下：（1）确定技术指标 在设计一个滤波器之前，必须首先根据工程实际的需要确定滤波器的技术指标。在很多实际应用中，数字滤波器常被用来实现选频操作。因此，指标的形式一般在频域中给出幅度和相位响应。幅度指标主要以2种方式给出。第一种是绝对指标。他提供对幅度响应函数的要求，一般应用于FIR滤波器的设计。第二种指标是相对指标。他以分贝值的形式给出要求。本文中滤波器的设计就以线性相位FIR滤波器的设计为例。（2）逼近 确定了技术指标后，就可以建立一个目标的数字滤波器模型（通常采用理想的数字滤波器模型）。之后，利用数字滤波器的设计方法（窗函数法、频率采样法等），设计出一个实际滤波器模型来逼近给定的目标。 （3）性能分析和计算机仿真 上两步的结果是得到以差分或系统函数或冲激响应描述的滤波器。根据这个描述就可以分析其频率特性和相位特性，以验证设计结果是否满足指标要求；或者利用计算机仿真实现设计的滤波器，再分析滤波结果来判断。FIR滤波器的设计问题实际上是确定能满足所要求的转移序列或脉冲响应的常数的问题，设计方法主要有窗函数、频率响应法和等波纹最佳逼近法等。FIR滤波器的缺点在于它的性能不如同样阶数的IIR滤波器，不过由于数字计算硬件的飞速发展，这一点已经不成为问题。再加上引入计算机辅助设计，FIR滤波器的设计也得到极大的简化。基于上述原因，FIR滤波器比IIR滤波器的应用更广4。2.2窗口设计法FIR滤波器的设计方法有许多种，如窗函数设计法、频率采样设计法和最优化设计法等。窗函数法是设计FIR数字滤波器的最简单的方法。它在设计FIR数字滤波器中有很重要的作用，正确地选择窗函数可以提高设计数字滤波器的性能，或者在满足设计要求的情况下，减小FIR数字滤波器的阶次。常用的窗函数有以下几种：矩形窗(Rectangular window)、三角窗(Triangular window)、汉宁窗(Hanning window)、汉明窗(Hamming window)、布拉克曼窗(Blackman window)等等。窗函数设计法的基本原理是用一定宽度窗函数截取无限脉冲响应序列获得有限长的脉冲响应序列，主要设计步骤为：（1）通过傅里叶逆变换获得理想滤波器的单位脉冲响应hd(n)。 (2-1)（2）由性能指标确定窗函数W(n)和窗口长度N。（3）求得实际滤波器的单位脉冲响应h(n)， h(n)即为所设计FIR滤波器系数向量b(n)。 (2-2)（4）检验滤波器性能常见的窗函数性能如表2-1所示：表2-1 常见窗函数性能表名称滤波器过渡带宽最小阻带衰减名称滤波器过渡带宽最小阻带衰减矩形1.8/M21dBPARZENWIN6.6/M56dB巴特利特6.1/M25dBFLATTOPWIN19.6/M108dB汉宁6.2/M44dBGAUSSWIN5.8/M60dB汉明6.6/M51dBBARTHANNWIN5.6/M40dB布莱克曼11/M74dBBLACKMANHARRIS16.1/M109dBBOHMANWIN5.8/M51.5dBCHEBWIN15.2/M113dBNUTTALLWIN15.4/M108dBTUKEYWIN2.4/M22dB2.3 BARTHANNWIN窗BARTHANNWIN窗是修正的BARTLETT-HANN窗，它的最小阻带衰减为40dB，过度带宽为5.6/M（M为滤波器阶数），公式给出如下5：wn=0.62-0.48+0.38cos(2，0nN 0 , 其余n(2-3)在MATLAB中的调用格式为：w=barthannwin(M)；相应的程序段：M=64;n=0:M-1;w_barth=barthannwin(M);hd_bs=ideal_lp(0.45*pi,M);h_bs=hd_bs.*w_barth;db,mag,pha,grd,w=freqz_m(h_bs,1); 用绘图命令画出的窗函数时域频域图如图2-1所示：图2-1 Barthannwin窗时域频域响应图3 设计步骤3.1设计流程图吹管音乐音频信号滤波去噪使用Bathannwin窗设计的FIR滤波器的设计流程如图3-1所示：开始吹管音乐音频信号的采集在MATLAB集成环境下用wavread函数将语音信号录入，并对语音信号加入噪声，画出语音信号加噪声前后的时域和频域波形图并进行比较。设定带阻滤波器的性能指标，下通带fpd=3000;下阻带fsd=3400;上阻带fsu=3600;上通带fpu=4000;通带波纹Rp=1;阻带波纹As=30 利用Barthannwin窗设计合理的FIR滤波器。将干扰后的音乐信号通过自己设计的滤波器，进行滤波去噪。比较滤波前后的时域波形和频谱图，并回放音乐信号，验证是否达到去噪效果。 是否达到去噪效果？ N Y 结束 图3-1 使用Bathannwin窗设计的FIR滤波器对吹管音乐音频信号去噪流程图3.2录制语音信号在网上下载一段吹管乐器演奏音乐，时间在2s左右，在转换过程中，将音乐格式设置为wav格式，属性设置为16000Hz，8位，单声道，128KB/秒，然后保存在目录e:下，命名为Taihuajiao.wav，最后在Matlab软件平台下调用函数wavread对音乐信号进行采样，记录采样频率和采样点数。具体调用如下： x,fs,bits=wavread(e:wang.wav); %读取音乐信号运行后得出fs=16000，bits=8。输出的第一个参数x是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放N=length(x); % 计算信号x的长度 fn=3500; % 单频噪声频率t=0:1/fs:(N-1)/fs; %计算时间范围，样本数除以采样频率 x=x(:,1);y=x+0.5*sin(fn*2*pi*t); %加入噪声函数sound(y,fs,bits); % 应该可以明显听出有尖锐的单频啸叫声X=fft(x); Y=fft(y); % 对原始信号和加噪信号进行fft变换,取幅度谱X=X(1:N/2); Y=Y(1:N/2); % 截取X的前半部分， 截取Y的前半部分deltaf=fs/N; % 计算频谱的谱线间隔f=0:deltaf:fs/2-deltaf; % 计算频谱频率范围然后用绘图命令分别画出加噪前后信号的时域和频域波形，运行程序后我们所得的波形图如图3-2所示：图3-2 加噪前后信号的时域图和频域比较图3.3 滤波器设计滤波器的设计就是要找到一组能满足特定滤波要求的系数向量a和b，其中b为系统函数的分子系数，a为系统函数分母系数。在本次的课程设计中我所采用的就是利用Barthannwin窗函数来设计FIR滤波器。在该滤波器的设计中，我们给出该滤波器的性能指标如下：fpd=3000; fsd=3400; fsu=3600; fpu=4000;Rp=1; As=30;截止频率也可以任意自选，在单频噪声干扰附近即可。在这里，很重要的是通带截止db值的设置。这个值根据我们使用的设计滤波器的方法来设定。因为使用的是Barthannwin窗法，Barthannwin窗函数中，滤波器的过渡带宽为5.6，最小阻带衰减为40db。所以，一定要将通带截止db值设置的小于40，所以，我将其设置为30db。用窗函数法，在MATLAB中利用Barthannwin窗设计FIR滤波器利用Matlab中的函数freqz画出滤波器的频率响应。利用数字信号处理里面学过的知识，根据自己选定的参数，用指定的方法设计FIR滤波器，得到FIR滤波器的阶数M。随后调用BARTHANNWIN(M)函数产生M阶的Barthannwin窗。然后，调用自编ideal_lp函数计算理想带阻滤波器的脉冲响应。最后，再调用自编freqz_m(h_bs,1)函数即可计算得到该滤波器的频率特性6。 fpd=3000;fsd=3400;fsu=3600;fpu=4000;Rp=1;As=30; % 带阻滤波器设计指标 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; wsd=fsd/fs*2*pi;wsu=fsu/fs*2*pi; % 将Hz为单位的模拟频率换算为rad为单位的数字频率 M=ceil(5.6*pi/dw)+1; % 计算BARTHANNWIN窗设计该滤波器时需要的阶数 n=0:M-1; % 定义时间范围 w_barth=BARTHANNWIN(M); % 产生M阶的BARTHANNWIN窗 hd_bs=ideal_lp(wcd,M)+ideal_lp(pi,M)-ideal_lp(wcu,M);% 调用自编函数计算理想带阻滤波器的脉冲响应 h_bs=w_barth.*hd_bs; % 用窗口法计算实际滤波器脉冲响应 db,mag,pha,grd,w=freqz_m(h_bs,1); % 调用自编函数计算滤波器的频率特性然后用绘图命令分别画出窗口法设计出来的滤波器的幅度特性和相位响应图，运行程序后所得到设计的滤波器的图形如图3-3所示：图3-3 滤波器的幅度特性以及相位响应图3.4 信号滤波处理滤波器设计完成后，在MATLAB中，FIR滤波器利用函数fftfilt对信号进行滤波。程序如下： y_fil=filter(h_bs,1,y); % 用设计好的滤波器对y进行滤波 Y_fil= abs(fft(y_fil);Y_fil=Y_fil(1:N/2)；% 计算频谱取前一半用绘图命令分别画出滤波前后音乐信号的时域波形图和频谱图对比图，如图3-4所示：图3-4 滤波前和滤波后的时域和频域波形对比图3.5 结果分析从图3-4中可以观察到，加噪信号的时域图中大部分都被加入的噪声给遮盖了，滤波后的双音频信号发生了衰减，滤波去噪信号与原信号的时域图基本相似，说明滤波器起到了滤波作用。加噪信号的频谱图中，我们可以很明显地看到与原信号频谱图相比，它在3500Hz左右处有一个尖脉冲，而滤波去噪信号的频谱图中该尖脉冲已经消失，波形大致与原图相似，可见滤波去噪效果基本不错。在将三个信号的时域波形和频谱图比较之后，我们还要通过回放去滤波去噪语音信号，来跟原信号相比，以检验滤波器的效果。在MATLAB中，经过sound(x,fs,bits)函数，对经过Barthannwin窗设计的FIR滤波器滤波去噪之后的音乐信号进行回放，可以听出滤波之后的信号跟原始信号一样清晰，完全滤除掉了噪声的干扰，通过下面语句来进行语音信号回放比较：sound(x,fs,bits); %显示原始双音频信号sound(y,fs,bits); %显示滤波后的双音频信号分别听原始语音和滤波后的语音信号，发现滤波后的语音信号噪声基本消失，同时原始信号强度稍有变化，基本达到了滤波的效果，所得结果证明了用Barthannwin窗设计的FIR滤波器和音乐信号去噪设计是成功的。4 出现的问题及解决方法本次课程设计中当遇到问题时我在运用自己所学的知识思考的同时还通过上网搜索资料，参考相关书籍以及询问同学和老师，集思广益，使问题尽量得到完善的解决，在整个设计中出现的主要问题及解决方法如下：（1）、起初对利用Barthannwin窗函数设计FIR滤波器的设计步骤很生疏，在采用Barthannwin窗函数设计的FIR滤波器时得不到理想的滤波器，因而信号的恢复不是特别理想，通过不断设置参数的值，最终设计出理想的滤波器，使信号得到较为理想的恢复；（2）、加入单频噪声时，单频噪声为余弦信号，其系数为1，该余弦信号是n行1列的矩阵，而我所录制的语音是1行n列的矩阵，所以需要将我所录制的语音信号进行转置，并且将余弦信号的系数改为0.5或者更小，否则程序会报错或者所画出来的图形将都是一条直线；（3）、在画滤波器性能波形的时候，要横坐标要计算w/pi的值，而不是直接取w的值；（4）、在对加入噪声后的语音用FIR滤波器进行语音去噪处理时，所编写的函数出现错误，绘制不出图形，分析错误主要是调用FIR滤波器时所用的fftfilt函数改为fitler函数。filter函数主要运用于IIR滤波器，所调用的格式为filter(b,a,x)，而fftfilt函数的调用为fftfilt(b,x)，所以在调用时应将a删去；（5）、绘制出滤波后的波形，发现FIR滤波器没有滤掉单频噪声。通过自己的仔细检查，是单频噪声的频率改动后，FIR滤波器的频率没有改动。所以单频噪声的频率也应该自己先定义，FIR滤波器的截止频率应该以单频噪声的频率为中心，这样重新运行后，结果正确；（6）、在做课程设计的文档时，发现流程图不会绘制。然后自己查了一些相关书籍得知，绘制流程图时一般用椭圆表示“开始”与“结束”，行动方案普通工作环节用矩形表示，问题判断或判定（审核/审批/评审）环节用菱形表示 ，箭头代表工作流方向；（7）、在写报告的时候，没有仔细注意报告版式、格式的要求，导致要修改的地方较多，后经仔细检查及修改最终达到报告的要求。5 结束语通过这次课程设计，我对其设计原理进行了更深一层次的理解，对书中原来学到的理论，有了用实践检验的机会，也有了更加确切的掌握，在设计中也使我对一些概念有了更深刻的理解。例如：在指标方面，我混淆了模拟指标和数字指标的概念，经过老师的点拨，自己更加明确，而且记忆深刻。还有在课程设计中每一次的数据输入都有其重要意义，用MATLAB编译程序时，可以根据滤波器指标的要求实时知道对滤波器的影响。例如，编好程序后，调试成功，任意改变输入阻带或者通带衰减，可以看到输出波形的变化，改变截止频率wc，同样可以看到输出波形的变化。由此，对理论的理解就更加简单方便，而且记忆力深刻。除此之外，由于老师只给出了部分的程序段模版，再加上设计滤波方法的不同以及信源的差异，对程序的编译不是一蹴而就的，而要经过多次的编译与调试，所以很多地方是自己，出错率也大大提高了，其中的错误原因很多次都是因为自己的不细心造成。但经过思考并仔细检查，最后得以编译成功。做设计时一定要小心谨慎，也许一个小小的语法错误和常量变量的定义的错误就造成整个程序出现问题，得不到所需的波形，导致实验结果不正确。但这就是过程，学习就是在过程中进行的，经过自己的不懈努力再加上同学们的帮助，不仅程序识读能力提高了很多，自己的编译水平也上了一个新台阶，更加熟系了MATLAB的应用，也对其中的函数如freqz、ideal函数等有了细致的了解。完成整个设计过程后，学到的东西已经不仅仅上面的那些东西，能将自己平时学到的东西能运用到实际中，能让理论和实际得以结合还是很不错的，整个课程设计的过程中也让我找到了动手的乐趣和思考的快乐，很有成就感。我从这次的课程设计中体会到，什么事都要我们真真正正用心去做了，才会使自己更加成长，没有学习就不可能有实践的机会，没有自己的实践就不会有所突破，希望这次的经历能让我们在以后的学习生活中不断成长。在这次的课程设计中，犯了些错误，也长了教训，发现自己还不够认真不够细心以前的知识学习的不是很牢固，但是这次课程设计，不但巩固了以前学的知识，还学习了新内容，收获颇丰。我认为这次收获最大的莫过于静心，学习不能急，一定要冷静，心无旁骛，不放过任何一个细节，就能取得十足的进步。本次课程设计在老师的悉心指导和严格要求下业已完成，从课题选择、方案论证到具体设计和调试，无不凝聚着老师的心血和汗水，在课程设计期间，也始终感受着导师的精心指导和无私的关怀，我受益匪浅。在此向老师表示深深的感谢和崇高的敬意。正是有了他们的悉心帮助和支持，才使我的课程设计顺利完成。参考文献1 张圣勤MATLAB7.0实用教程M 北京：机械工程出版社，20063 2 维纳K英格尔，约翰G普罗克斯（著）刘树棠（译）数字信号处理（MATLAB版）M，西安：交通大学出版社，200813 罗军辉，MALAB7.0在数字信号处理中的应用，北京：机械工业出版社，200554 张小虹信号系统与数字信号处理M第版西安：西安电子科技出版社，20025 谢德芳数字信号处理M北京：科学出版社，20056 陈怀琛，数字信号处理教程：MATLAB释义与实现，北京：电子工业出版社，200412附录1：吹管音乐音频信号提取和加噪源程序清单% 程序名称：吹管音乐音频信号提取和加噪% 程序功能：用MATLAB语言读出自己所录得双音频信号并对双音频信号加噪% 程序作者：王二% 最后修改日期：- x,fs,bits=wavread(e:Taihuajiao.wav);% 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。 sound(x,fs,bits);% 按指定的采样率和每样本编码位数回放 N=length(x); % 计算信号x的长度 fn=3500;% 单频噪声频率 t=0:1/fs:(N-1)/fs;% 计算时间范围，样本数除以采样频率 x=x(:,1);y=x+0.5*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);% 截取X的前半部分， 截取Y的前半部分 deltaf=fs/N;% 计算频谱的谱线间隔 f=0:deltaf:fs/2-deltaf;% 计算频谱频率范围 figure(1) subplot(2,2,1);plot(t,x);grid %画出加噪前后信号的时域和频域波形，布局为2*2子图，加上横纵坐标，网格和标题。 axis(0,2,-1.1,1.1) xlabel(t);ylabel(x); title(加噪声前的时域波形图); subplot(2,2,2);plot(f,X);grid axis(0,8000,0,4000) xlabel(f);ylabel(X);title(加噪声前的频域波形图); subplot(2,2,3);plot(t,y);grid axis(0,2,-2.1,2.1) xlabel(t);ylabel(y); title(加噪声后的时域波形图); subplot(2,2,4);plot(f,Y);grid axis(0,8000,0,6000) xlabel(f);ylabel(Y); title(加噪声后的频域波形图);附录2：使用BARTHANNWIN窗设计的FIR滤波器源程序清单% 程序名称：使用BARTHANNWIN窗设计的FIR滤波器% 程序功能：使用BARTHANNWIN窗设计FIR滤波器，观察滤波器的相位和幅度特性图形% 程序作者：王二% 最后修改日期：- fpd=3000;fsd=3400;fsu=3600;fpu=4000;Rp=1;As=30; % 带阻滤波器设计指标 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(5.6*pi/dw); % 计算BARTHANNWIN窗设计该滤波器时需要的阶数 n=0:M-1; % 定义时间范围 w_barth=BARTHANNWIN(M); % 产生M阶的BARTHANNWIN窗 hd_bs=ideal_lp(wcd,M)+ideal_lp(pi,M)-ideal_lp(wcu,M); % 调用自编函数计算理想带阻滤波器的脉冲响应 h_bs=w_barth.*hd_bs; % 用窗口法计算实际滤波器脉冲响应 db,mag,pha,grd,w=freqz_m(h_bs,1); % 调用自编函数计算滤波器的频率特性 figure(2) subplot(2,2,1);plot(w/pi,db);grid axis(0,0.6,-40,10) xlabel(w/pi);ylabel(db); title(以db为单位的幅度特性);a=0 3;b=-30 -30;line(a, b,Color,r,LineStyle,-);a=0 3;b=-1 -1;line(a, b,Color,r,LineStyle,-);a=wsd/pi wsd/pi;b=-40 10;line(a, b,Color,r,LineStyle,-);a=wsu/pi wsu/pi;b=-40 10;line(a, b,Color,r,LineStyle,-); subplot(2,2,2);plot(w/pi,mag);grid axis(0,0.6,0,1.2) xlabel(w/pi);ylabel(幅度mag); title(以线性为单位的幅度特性); subplot(2,2,3);plot(w/pi,pha);grid axis(0,1.2,-4,4) xlabel(w/pi);ylabel(相位phg); title(滤波器相位响应图);axis tight; subplot(2,2,4);stem(n,h_bs);grid axis(0,150,-0.3,0.7) xlabel(n);ylabel(h(n); title(滤波器响应图);axis tight;附录3：信号滤波处理源程序清单% 程序名称：信号滤波处理% 程序功能：对信号进行滤波处理得到滤波后的时域和频率波形图 % 程序作者：王二% 最后修改日期：-figure(3) y_fil=fftfilt(h_bs,y); % 用设计好的滤波器对y进行滤波 Y_fil= abs(fft(y_fil);Y_fil=Y_fil(1:N/2); % 计算频谱取前一半 subplot(3,2,1);plot(t,x);grid axis