语音信号滤波去噪——使用flattopwin窗设计FIR滤波器

1、 黄珺 吹管滤波去噪使用 FLATTOPWIN 窗设计的 FIR 滤波器 第0页，共22 页吹管音乐滤波去噪吹管音乐滤波去噪使用使用FLATTOPWIN 窗设计的窗设计的 FIR 滤波器滤波器学生姓名：黄珺学生姓名：黄珺 指导老师：高明指导老师：高明摘摘 要要 本课程设计主要使用 Flattopwin 窗设计的 FIR 滤波器对吹管音乐信号进行滤波去噪。课程设计的的平台为 MATLAB。采集一段吹管音乐，绘制波形并观察其频谱，给定相应技术指标，用窗函数法设计一个满足指标的 FIR 滤波器，对该吹管音乐信号进行滤波去噪处理。根据滤波后的时域图和原始音乐信号时域图的比较，以及滤波后信号的频谱图和原

2、始吹管音乐频谱图的比较，最后回放滤波后吹管音乐，滤波后的吹管音乐信号与原始吹管音乐信号一样清晰，成功地实现了滤波，达到了设计的要求。关键词关键词 课程设计；滤波去噪；FIR 滤波器；Flattopwin 窗；MATLAB 黄珺 吹管滤波去噪使用 FLATTOPWIN 窗设计的 FIR 滤波器 第1页，共22 页1 引引 言言本课程设计是采用 Flattopwin 窗设计的 FIR 滤波器对吹管音乐进行滤波去噪。通过课程设计了解 FIR 滤波器设计的原理和步骤，掌握用 Matlab 语言设计滤波器的方法，了解 DSP 对 FIR 滤波器的设计及编程方法。通过观察吹管音乐滤波前后的时域波形的比较，

3、加深对滤波器作用的理解。通过对比滤波前后波形图的比较和放滤波前后吹管音乐的对比 ，可以看出滤波器对有用信号无失真放大具有重大意义。1.1 课程设计目的课程设计目的熟悉 Matlab 语言环境，掌握 Matlab 语言的编程规则，利用 Flattopwin 窗函数设计法来设计符合要求的 FIR 滤波器来实现吹管音乐的滤波去噪。并绘制滤波前后的时域波形和频谱图。根据图形分析判断滤波器设计的正确性。通过本次课程设计熟悉利用 Flattopwin 窗函数法设计 FIR 滤波器的过程。增强自己独立解决问题的能力，提高自己的动手能力。加深对理论知识联系实际问题的理解。为以后的工作奠定坚实的基础。1.2 课

4、程设计要求课程设计要求录制一段吹管音乐，绘制观察波形及频谱图。根据 Flattopwin 的性能指标合理设计 FIR 滤波器，对吹管音乐加入干扰，再用滤波器对干扰语音信号进行滤波去噪，比较滤波前后的频谱图并进行分析。再回放吹管音乐对比原吹管音乐。看滤波器是否对吹管音乐进行了滤波去噪。1.3 课程设计平台课程设计平台MATLAB 名字由 MATrix 和 LABoratory 两词的前 3 个字母组合而成。20 世纪 70年代后期，时任美国新墨西哥大学计算机科学系主任的 Cleve Moler 教授出于减轻学生编程负担的动机，为学生设计了一组调用 LINPACK 和 EISPACK 库程序的“通

5、俗易用”的接口，此即用 FORTRAN 编写的萌芽状态的 MATLAB1。 MATLAB 软件包括五大通用功能：数值计算功能（Nemeric）；符号运算功能 黄珺 吹管滤波去噪使用 FLATTOPWIN 窗设计的 FIR 滤波器 第2页，共22 页（Symbolic）；数据可视化功能（Graphic）；数据图形文字统一处理功能（Notebook）和建模仿真可视化功能（Simulink）。其中，符号运算功能的实现是通过请求 MAPLE 内核计算并将结果返回到 MATLAB 命令窗口。该软件有三大特点：一是功能强大；二是界面友善、语言自然；三是开放性强。目前，Mathworks 公司已推出30 多

6、个应用工具箱。MATLAB 在线性代数、矩阵分析、数值及优化、数理统计和随机信号分析、电路与系统、系统动力学、信号和图像处理、控制理论分析和系统设计、过程控制、建模和仿真、通信系统、以及财政金融等众多领域的理论研究和工程设计中得到了广泛应用。 MATLAB 在信号与系统中的应用主要包括符号运算和数值计算仿真分析。由于信号与系统课程的许多内容都是基于公式演算，而 MATLAB 借助符号数学工具箱提供的符号运算功能能基本满足信号与系统课程的需求。例如，解微分方程、傅里叶正反变换、拉普拉斯正反变换、z 正反变换等。MATLAB 在信号与系统中的另一主要应用是数值计算与仿真分析，主要包括函数波形绘制、

7、函数运算、冲激响应与阶跃响应仿真分析、信号的时域分析、信号的频谱分析、系统的 S 域分析、零极点图绘制等内容。数值计算仿真分析可以帮助学生更深入理解信号与系统的理论知识，并为将来使用MATLAB 进行信号处理领域的各种分析和实际应用打下基础2。2 基本理论基本理论2.1 FIR 滤波器滤波器FIR 滤波器：有限长单位冲激响应滤波器，是数字信号系统中最基本的元件，它可以在保证任意幅频特性的同时具有严格的线性相频特性，同时其单位抽样响应是有限长的，因而滤波器是稳定的系统。因此，FIR 滤波器在通信、图像处理、模式识别等领域都有着广泛的应用。FIR 数字滤波器设计的基本步骤如下：（1）确定指标在设计

8、一个滤波器之前，必须首先根据工程实际的需要确定滤波器的技术指标。在很多实际应用中，数字滤波器常常被用来实现选频操作。因此，指标的形式一般在频域中给出幅度和相位响应。幅度指标主要以两种方式给出。第一种是绝对指标。它提供对幅度响应函数的要求，一般应用于 FIR 滤波器的设计。第二种指标是相对指标。 黄珺 吹管滤波去噪使用 FLATTOPWIN 窗设计的 FIR 滤波器 第3页，共22 页它以分贝值的形式给出要求。 （2）逼近 确定了技术指标后，就可以建立一个目标的数字滤波器模型。通常用理想的数字滤波器模型。之后，利用数字滤波器的设计方法，设计出一个实际滤波器模型来逼近给定的目标。 （3）性能分析和

9、计算机仿真上两步的结果是得到以差分或冲激响应描述的滤波器。根据这个描述就可以分析其频率特性和相位特性，以验证设计结果是否满足指标要求；或者利用计算机仿真实现设计的滤波器，再分析滤波结果来判断。FIR 滤波器的设计问题实际上是确定能满足所要求的转移序列或脉冲响应的常数的问题，设计方法主要有窗函数、频率响应法和等波纹最佳逼近法等。2.2 窗口设计法窗口设计法窗口设计法是一种通过截断和计权的方法使无限长非因果序列成为有限长脉冲响应序列的设计方法。通常在设计滤波器之前，应该先根据具体的工程应用确定滤波器的技术指标。在大多数实际应用中，数字滤波器常常被用来实现选频操作，所以指标的形式一般为在频域中以分贝

10、值给出的相对幅度响应和相位响应。窗口设计法步骤如下：（1）根据过渡带宽及阻带衰减要求，选择窗函数的类型并估计窗口长度 N。窗函数的类型可根据最小阻带衰减 AS 独立选择。（2）根据待求滤波器的理想频率响应求出理想单位脉冲响应 hd(n)。（3）由性能指标确定窗函数 W(n)和长度 N。（4）求得实际滤波器的单位脉冲响应 h(n)， h(n)即为所设计 FIR 滤波器系数向量b(n)。 ( )( ).( )dh nh n W n (2.1) 黄珺 吹管滤波去噪使用 FLATTOPWIN 窗设计的 FIR 滤波器 第4页，共22 页常见的窗函数性能表如下图2.1所示：名称滤波器过渡带宽最小阻带衰减

11、名称滤波器过渡带宽最小阻带衰减矩形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.5dBCHEBWIN15.2/M113dBNUTTALLWIN15.4/M108dBTUKEYWIN2.4/M22dB图2.1 常见窗函数性能表2.3 FLATTOPWIN 窗窗w=Flattopwin (L) 返回

12、L-点 Flattopwin 窗口中列向量。Flattopwin 窗的滤波器的过渡带宽为 19.6/M，最小阻带衰减 108db。定义式：定义式：)8cos()6cos()4cos()2cos()(43210NnaNnaNnaNnaan时间波形和幅度谱：时间波形和幅度谱： 黄珺 吹管滤波去噪使用 FLATTOPWIN 窗设计的 FIR 滤波器 第5页，共22 页图 2.2 时间波形图 2.3 幅度谱 黄珺 吹管滤波去噪使用 FLATTOPWIN 窗设计的 FIR 滤波器 第6页，共22 页 3 3 设计步骤设计步骤3.1 设计流程图设计流程图根据设计的要求，首先从网上下载一段音乐信号，对音乐信

13、号加入噪声干扰，再利用 Flattopwin 窗设计合理的 FIR 滤波器。最后用滤波器对干扰后的语音信号进行滤波去噪。具体设计流程图如下图 3.1 所示： NO YES 图 3.1 设计流程图 3.2 录制音乐信号录制音乐信号在网上下载一段音乐，再利用视频转换器将其转换成*.格式，且为单声3mpwav开始下载一段吹管音乐信号，命名为“yinyue.wav” ，绘制出其时域波形和频谱图。加入噪声干扰，绘制干扰后的时域波形和频谱图。利用 Flattopwin 窗设计合理的 FIR 滤波器。将干扰后的吹管音乐信号通过自己设计的滤波器，进行滤波去噪。比较滤波前后的时域波形和频谱图，并回放音乐信号，验

14、证是否达到去噪效果。结束是否达到去噪 效果？ 黄珺 吹管滤波去噪使用 FLATTOPWIN 窗设计的 FIR 滤波器 第7页，共22 页道，再将此.格式音乐控制在 5 秒内，以减少设计中的误差。然后在 Matlab 软件平wav台下，利用函数 wavread 对语音信号进行采样，记住采样频率和采样点数，采集完成后在信号中加入一个单频噪声，绘制原音乐信号和加噪后的音乐信号的时域和频域的波形图。具体调用如下：x,fs,bits=wavread(g:yinyue.wav); % 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs 是生成该波形文件时的采样率，bits 是波形文件每样

15、本的编码位数。 sound(x,fs,bits);%播放音乐信号 N=length(x);%计算音乐信号的长度 fn=2200; % 单频噪声频率 t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率 x=x(:,1); y=x+sin(fn*2*pi*t); plot(t,x); xlabel(t ins);ylabel(x);%画原始音乐信号的时域波形图sound(y,fs,bits);%播放干扰后的音乐信号，应该可以明显听出有尖锐的单频啸叫声运行程序后得到的波形图如图 3.2 所示： 黄珺 吹管滤波去噪使用 FLATTOPWIN 窗设计的 FIR 滤波器 第8页，共

16、22 页图 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; % 计算频谱频率范围figure(2) subplot(2,2,1);plot(t,x);xlabel(时间(t);ylabel(幅度); title(原始音

17、乐信号); subplot(2,2,2);plot(f,X);xlabel(频率(f);ylabel(幅度谱);title(原始音乐信号幅度谱);axis(0,8000,0,600); subplot(2,2,3);plot(t,y);xlabel(时间(t);ylabel(幅度);title(加干扰后的音乐信号); subplot(2,2,4);plot(f,Y);xlabel(频率(f);ylabel(幅度谱);title(加干扰后的音乐信号幅度谱);axis(0,8000,0,600);运行程序后我们所得的波形图如图 3.3 所示： 黄珺 吹管滤波去噪使用 FLATTOPWIN 窗设计的

18、 FIR 滤波器 第9页，共22 页图 3.3 干扰前后信号时域图与频谱图的比较3.3 滤波器设计滤波器设计滤波器的设计就是要找到一组能满足特定滤波要求的系数向量 a 和 b，其中 b 为系统函数的分子系数，a 为系统函数分母系数。在本次的课程设计中我所采用的就是利用Flattopwin 窗函数来设计 FIR 滤波器。 具体调用如下： fpd=2100;fsd=2150;fsu=2250;fpu=2300;Rp=1;As=20;% 带阻滤波器设计指标 fcd=(fpd+fsd)/2;fcu=(fpu+fsu)/2;df=min(fsd-fpd),(fpu-fsu);% 计算上下边带中心频率，和

19、频率间隔 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.1*pi/dw)+1; % 计算窗设计该滤波器时需要的阶数 n=0:M-1; % 定义时间范围 w_par=(Flattopwin (M); % 产生 M 阶的 Flattopwin 窗 hd_bs=ideal_lp(wcd,M)+ideal_lp(pi,M)-ideal_lp(wcu,M);% 调用自编函数计算理想带阻滤波器的脉冲响应 h_bs

20、=w_par.*hd_bs; % 用窗口法计算实际滤波器脉冲响应 db,mag,pha,grd,w=freqz_m(h_bs,1);% 调用自编函数计算滤波器的频率特性figure(3)subplot(2,2,1);plot(w,db);title(滤波器幅度响应图);xlabel(w/pi);ylabel(db);axis(0,0.5,-30,10); line(0,0.5,-As,-As,color,r,linestyle,-,LineWidth,2); line(0,0.5,-Rp,-Rp,color,r,linestyle,-,LineWidth,2); line(wsd,wsd,-3

21、0,10,color,r,linestyle,-,LineWidth,2); line(wsu,wsu,-30,10,color,r,linestyle,-,LineWidth,2); subplot(2,2,2);plot(w,mag);title(滤波器幅度响应图);xlabel(w/pi);ylabel(幅度 mag);axis(0,1,-0.5,1.5); subplot(2,2,3);plot(w,pha);title(滤波器相位响应图);xlabel(w/pi);ylabel(相位 pha); 黄珺 吹管滤波去噪使用 FLATTOPWIN 窗设计的 FIR 滤波器 第10页，共22

22、 页axis(0,3,-4,4); subplot(2,2,4);stem(n,h_bs);title(滤波器脉冲响应图);xlabel(n);ylabel(h(n);axis(0,3000,0,1);程序中用到的二个函数定义函数 IDEAL_LP.M 定义：function hd = ideal_lp(wc,M);% 理想低通滤波器计算% -% hd = ideal_lp(wc,M)% hd = 0 to M-1 之间的理想脉冲响应% wc = 截止频率(弧度) % M = 理想滤波器的长度%alpha = (M-1)/2;n = 0:1:(M-1);m = n - alpha + eps;

23、hd = sin(wc*m) ./ (pi*m);函数 FREQZ_M.M 定义：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 到 pi 弧度区间内的相位响应% grd = 0 到 pi 弧度区间内的群迟延% w = 0 到 pi 弧度区间内的 501 个频率样本向量% b = Ha(z)的分子多项式系数(对 FIR b=h)%

24、a = Ha(z)的分母多项式系数(对 FIR: a=1) 黄珺 吹管滤波去噪使用 FLATTOPWIN 窗设计的 FIR 滤波器 第11页，共22 页%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(

25、1:1:500); grd; grd(2:1:501) 0;% grd = median(grd)*500/pi;得到的波形图如图 3.4 所示： 图 3.4 设计的滤波器图形 黄珺 吹管滤波去噪使用 FLATTOPWIN 窗设计的 FIR 滤波器 第12页，共22 页3.4 信号滤波处理信号滤波处理 滤波器设计完成后，在 MATLAB 平台上用函数 filter 实现滤波。 滤波程序如下： y_fil=filter(h_bs,1,y);% 用设计好的滤波器对 y 进行滤波 Y_fil=fft(y_fil);Y_fil=Y_fil(1:N/2); % 计算频谱取前一半figure(4) sub

26、plot(3,2,1);plot(t,x);xlabel(时间(t);ylabel(幅度); title(原始音乐信号); subplot(3,2,2);plot(f,X);xlabel(频率(f);ylabel(幅度谱);title(原始音乐信号幅度谱);axis(0,8000,0,600); subplot(3,2,3);plot(t,y);xlabel(时间(t);ylabel(幅度);title(加干扰后的音乐信号); subplot(3,2,4);plot(f,Y);xlabel(频率(f);ylabel(幅度谱);title(加干扰后的音乐信号幅度谱);axis(0,8000,0,

27、600); subplot(3,2,5);plot(t, y_fil);xlabel(时间(t);ylabel(幅度); title(滤波后音乐信号); subplot(3,2,6);plot(f,Y_fil);xlabel(频率(f);ylabel(幅度谱);title(滤波后音乐信号幅度谱);axis(0,8000,0,600);得到的滤波前后音乐信号的时域波形图和频谱图对比图如图 3.5 所示： 图 3.5 滤波前后音乐信号的时域波形图和频谱图对比图 黄珺 吹管滤波去噪使用 FLATTOPWIN 窗设计的 FIR 滤波器 第13页，共22 页3.5 结果分析结果分析在 MATLAB 中，

28、经过 sound(x,fs,bits)函数，对经过 Flattopwin 窗设计的 FIR 滤波器之后的音乐信号进行回放，可以听出滤波之后的音乐信号跟原始音乐信号一样清晰，完全滤除掉了噪声的干扰，通过下面语句来进行音乐信号回放比较： sound (y_fil,fs,bits); 播放滤波后的音乐信号 所得结果证明了用 Flattopwin 窗设计的 FIR 滤波器和音乐信号去噪设计是成功的。 4 出现的问题及解决方法在本次课程设计中我遇到的问题如下：1、 在录音时，没有将录音的属性改为 wav 的格式，当在 Matlab 软件平台上调用时，出现无法调用的提示，不知道如何将音乐信号截断成自己所需

29、的长度。2、对在开始录制音乐信号并将其导入 MATLAB 中的时候出现过错误原因是我所下载的音乐信号是双音频信号，利用 Flattopwin 窗函数设计 FIR 滤波器的设计步骤很生疏。 3、在在最后对加噪音乐信号进行滤波并绘图的时也出现了错误，绘制出最后的滤波信号的频谱图时纵坐标出现了负值，但是原信号的频谱图中却没有，采用 Flattopwin 窗函数设计的 FIR 滤波器时得不到理想的滤波器，因而信号的恢复不是特别理想。 针对以上问题，相应的解决方案如下： 1、通过与同学讨论发现，不能直接更改文件后缀，应用转换器将其转换为 wav 格式， 向同学学习如何截取音乐。 2、在老师的指导下我将音

30、乐信号变成了单声道信号，自己上网查阅资料，或向图书馆借阅资料，掌握利用 Flattopwin 窗函数设计 FIR 滤波器的方法和步骤。 3、对滤波的信号进行 fft 变换之后取绝对值，通过过不断设置参数的值，最终达到最理想的值，设计出理想的滤波器，使信号得到理想恢复。 黄珺 吹管滤波去噪使用 FLATTOPWIN 窗设计的 FIR 滤波器 第14页，共22 页 5 结束语 本次的课程设计，我的任务是利用 Flattopwin 窗函数设计 FIR 滤波器对吹管音乐信号滤波去噪。在本次课程设计之前，我对 Flattopwin 窗函数完全没有了解，因此在看到这个题目时，我是一头雾水。但是通过自己翻阅

31、资料和询问同学，我掌握了用Flattopwin 窗函数设计 FIR 滤波器的方法步骤，了解了窗函数的基本设计流程。经过这两周的课程设计，我学会了很多东西。 在此我要感谢学校为我们提供这次课程设计的机会，感谢老师对我的悉心指导，也感谢同学对我的帮助。这次的课程设计让我理论联系实际，不仅巩固了我们的理论知识，还提高了我的动手能力，经过这两周的课程设计经过为期两周的课程设计，这次课程设计比较的简单，根据老师给出的模板，再结合以前学过的知识，并在同学和老师的帮助下我顺利的完成了任务。不同于在教室里上的理论，这次的课程设计需要将我们平时所学习的知识运用到实践之中，将知识学以致用。因为是以所学理论为基础，

32、所以在课程设计的过程中，我又重温滤波器、窗函数等知识，更加熟悉了MATLAB 的操作。课程设计是我们运用所学知识，动手实践的一个很好的机会。它既可以帮助我们加深对所学知识的理解，又能提高我们运用知识，联系实际，动手实践的能力。而且在设计过程中可能用到我们没学过的知识，需要我们去查阅资料获取相关信息，这又提高了我们查找信息和学习新知识的能力。在实物的调试与检测过程中，又会遇到许多意想不到的问题，需要我们去分析原因和解决问题。也体会到真正的去独立地完成一件事情是很困难的，同学以及老师的帮助和提醒是必须的。通过这次课程设计，我拓宽了知识面，锻炼了实际操作能力，综合素质也得到了提高，进一步加深了了我们

33、对专业的认识和激发了我们对专业的兴趣。虽然课程设计结束了，但是我们的学习还没结束，对知识的进一步学习还需要继续，很开心成功地完成了这次设计。 黄珺 吹管滤波去噪使用 FLATTOPWIN 窗设计的 FIR 滤波器 第15页，共22 页 参考文献参考文献1张志涌精通 MATLAB 6.5 版M北京：北京航空航天大学出版社，20032约翰普罗克斯数字信号处理M西安：西安交通大学出版社，20093张小虹信号系统与数字信号处理M第 1 版西安：西安电子科技出版社，20024谢德芳数字信号处理M北京：科学出版社，20055郝建军数字通信M第 2 版北京：北京邮电大学出版社，20106张威MATLAB 基

34、础与编程入门M西安：西安电子科技大学出版社，2010 黄珺 吹管滤波去噪使用 FLATTOPWIN 窗设计的 FIR 滤波器 第16页，共22 页附录一： 音乐信号滤波去噪使用 FLARTOPWIN 滤波器%程序功能：在 Matlab 中，用窗口设计法设计 FIR 滤波器%程序作者：赵伯政%最后修改日期：2012-3-4%程序一： x,fs,bits=wavread(yinyue.wav);%读取音乐信号 sound(x,fs,bits);%播放音乐信号 N=length(x);%计算音乐信号的长度 fn=2200; % 单频噪声频率 t=0:1/fs:(N-1)/fs; % 计算时间范围，样

35、本数除以采样频率 x=x(:,1); y=x+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; % 计算频谱频率范围figure(2) subpl

36、ot(2,2,1);plot(t,x);xlabel(时间(t);ylabel(幅度); title(原始音乐信号); subplot(2,2,2);plot(f,X);xlabel(频率(f);ylabel(幅度谱);title(原始音乐信号幅度谱);axis(0,8000,0,600); subplot(2,2,3);plot(t,y);xlabel(时间(t);ylabel(幅度);title(加干扰后的音乐信号); 黄珺 吹管滤波去噪使用 FLATTOPWIN 窗设计的 FIR 滤波器 第17页，共22 页 subplot(2,2,4);plot(f,Y);xlabel(频率(f);y

37、label(幅度谱);title(加干扰后的音乐信号幅度谱);axis(0,8000,0,600);%程序三： fpd=2100;fsd=2150;fsu=2250;fpu=2300;Rp=1;As=20;% 带阻滤波器设计指标 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*

38、pi; M=ceil(6.1*pi/dw)+1; % 计算窗设计该滤波器时需要的阶数 n=0:M-1; % 定义时间范围 w_par=(Flattopwin (M); % 产生 M 阶的 Flattopwin 窗 hd_bs=ideal_lp(wcd,M)+ideal_lp(pi,M)-ideal_lp(wcu,M);% 调用自编函数计算理想带阻滤波器的脉冲响应 h_bs=w_par.*hd_bs; % 用窗口法计算实际滤波器脉冲响应 db,mag,pha,grd,w=freqz_m(h_bs,1);% 调用自编函数计算滤波器的频率特性figure(3)subplot(2,2,1);plot(

39、w,db);title(滤波器幅度响应图);xlabel(w/pi);ylabel(db);axis(0,0.5,-30,10); line(0,0.5,-As,-As,color,r,linestyle,-,LineWidth,2); line(0,0.5,-Rp,-Rp,color,r,linestyle,-,LineWidth,2); line(wsd,wsd,-30,10,color,r,linestyle,-,LineWidth,2); line(wsu,wsu,-30,10,color,r,linestyle,-,LineWidth,2); subplot(2,2,2);plot(

40、w,mag);title(滤波器幅度响应图);xlabel(w/pi);ylabel(幅度 mag);axis(0,1,-0.5,1.5); subplot(2,2,3);plot(w,pha);title(滤波器相位响应图);xlabel(w/pi);ylabel(相位 pha);axis(0,3,-4,4); subplot(2,2,4);stem(n,h_bs);title(滤波器脉冲响应图);xlabel(n);ylabel(h(n); 黄珺 吹管滤波去噪使用 FLATTOPWIN 窗设计的 FIR 滤波器 第18页，共22 页axis(0,3000,0,1);%程序四： y_fil=filter(h_bs,1,y);% 用设计好的滤波器对 y 进行滤波 Y_fil=fft(y_fil);Y_fil=Y_fil(1:N/2); % 计算频谱取前一半figure(4) subplot(3,2,1);plot(t,x);xlabel(时间(t);ylabel(幅度); title(原始音乐信号); s