弹拨音乐滤波去噪——使用GAUSSWIN设计的FIR滤波器

弹拨音乐滤波去噪

——使用GAUSSWIN设计的FIR滤波器学生姓名：季政府指导老师：胡双红

摘要本课程设计主要内容是设计利用窗口设计法选择GAUSSWIN设计一个FIR滤波器，对一段含噪弹拨音乐信号进行滤波去噪处理并根据滤波前后的波形和频谱分析滤波性能。本课程设计仿真平台为MATLAB7.0，开发工具是M语言编程。首先从网站上下载一段弹拨乐器演奏音乐，绘制波形并观察其频谱特点，然后加入一个带外单频噪声，用高斯窗设计一个满足指标的FIR滤波器，最后对该音乐信号进行滤波去噪处理，比较滤波前后波形和频谱并进行分析。由分析结果可知，滤波器后的弹拨音乐信号与原始信号基本一致，即设计的FIR滤波器能够去除信号中所加单频噪声，达到了设计目的。

关键词课程设计；滤波去噪；FIR滤波器；GAUSSWIN；MATLAB7.0

1 引言

本课程设计主要是将一段弹拨音乐信号加入噪声，然后用某种函数法设计出的FIR 滤波器对加入噪声后的弹拨音乐信号进行滤波去噪处理，处理时采用的是利用窗口设计法选择GAUSSWIN设计的FIR滤波器，通过课程设计了解FIR 滤波器设计的原理和步骤，掌握用Matlab 语言设计滤波器的方法，观察音乐信号滤波前后的时域波形的比较，加深对滤波器作用的理解[1]。

1.1课程设计的目的

数字信号处理（Digital Signal Proccessing，简称DSP）是一门涉及许多学科而广泛应用于许多领域的新型学科。20世纪60年代后，随着计算机和信息技术的飞速发展，数字信号处理技术应运而生并并得到迅速的发展。在过去的二十多年时间里，数字信号处理已经在通信领域得到极为广泛的应用。数字信号处理是利用计算机或专用处理设备，以数字形式对信号进行采集、变换、滤波、估值、增强、压缩、识别等处理，以得

到符合人们的需要的信号形式。

在本次课程设计中，最主要的设计是设计FIR滤波器，FIR滤波器的设计方法主要分为两类，第一类是基于逼近理想滤波器特性的方法，包括窗函数法、频域采样法和等波纹最佳逼近法；第二类是最优设计法。

本次的课程设计主要采用的是第一类设计方法，是利用GAUSSWIN法设计FIR滤波器对一段弹拨音乐信号进行滤波去噪，通过这一过程，对滤波前后波形进行对比分析得到结论。此课程设计比较简单，主要是将书本中的知识应用到现实中，并且根据自己对设计题目的理解，运用软件编写出程序实现这一设计，也是我们对数字信号处理的原理进行验证的一个过程。对此，也可以加深我们对所学知识的理解，培养我们的动手能力。

1.2课程设计的要求

（1）滤波器指标必须符合工程实际。

（2）设计完后应检查其频率响应曲线是否满足指标。

（3）处理结果和分析结论应该一致，而且应符合理论。

（4）独立完成课程设计并按要求编写课程设计报告书。

1.3设计平台MATLAB

MATLAB是由美国Math Works公司20世纪80年代中期推出的数学软件。MATLAB 是“Matric Laboratory”的缩写，意及“矩阵实验室”，优秀的数值计算能力和卓越的数据可视化能力使其很快在数学软件中脱颖而出。Matlab已经发展成为多学科、多种工作平台的功能强大的大型软件。在欧美的高校和研究机构中，MATLAB是一种非常流行的计算机语言，许多重要的学术刊物上发表的论文均是用MATLAB来分析计算以及绘制出各种图形。

MATLAB是一完整的并可扩展的计算机环境，是一种进行科学和工程计算的交互式程序语言。它的基本数据单元是不需要指定维数的矩阵，它可直接用于表达数学的算式和技术概念，而普通的高级语言只能对一个个具体的数据单元进行操作。因此，解决同样的数值计算问题，使用MATLAB要比使用Basic、Fortran和C语言等提高效率许多倍。许多人赞誉它为万能的数学“演算纸”。MATLAB采用开放式的环境，你可以读到它的算法，并能改变当前的函数或增添你自己编写的函数[2]。

MATLAB包含的内容非常丰富，功能强大，可以概括为以下几个方面：

(1)可以在多种操作系统下运行，如DOS、Windows 95/98/2000/2000/NT、Compaq Alpha、LinuxSun Solaris等。

(2)有超过500种的数学、统计、科学及工程方面的函数，使用简单快捷，并且有很强的用户自定义函数的能力。

(3)有强大的图形绘制和可视化功能，可以进行视觉数据处理和分析，进行图形、图像的显示及编辑，能够绘制二维、三维图形，使用户可以制作高质量的图形，从而写出图文并茂的文章。

(4)有从外部文件及外部硬件设备读入数据的能力。

(5)有丰富的工具箱〔toolbox〕。各个领域的专家学者将众多学科领域中常用的算法编写为一个个子程序，即m文件，这些m文件包含在一个个工具箱中。其工具箱可以分为两大类，即功能性工具箱和科学性工具箱。功能性工具箱主要用来扩充MATLAB 的符号计算、图形可视化、建模仿真、文字处理等功能以及与硬件实时交互的功能。学科性工具箱是按学科领域来分类的，如信号处理、控制、通信、神经网络图像处理、系统辨识、鲁棒控制、模糊逻辑、小波等工具箱。

MATLAB中的信号处理工具箱内容丰富，使用简便。在数字信号处理中常用的算法，如FFT，卷积，相关，滤波器设计，参数模型等，几乎都只用一条语句即可以调用。数字信号处理所常用的函数有波形的产生、滤波器的分析和设计、傅里叶变换、Z变换等[3]。

2 设计原理

2.1 FIR滤波器

滤波器根据其冲激响应函数的时域特性，可分为2种，即无限长冲激响应（IIR）滤波器和有限长冲激响应（FIR）滤波器。FIR和IIR的滤波原理都是进行卷积，就是对数入信号进行某种计算。FIR用处就在于对数字信号进行必要的处理，得到所需的输出信号。

FIR系统有自己突出的优点：系统总是稳定的；易实现线性相位；允许设计多通带（或多阻带）滤波器，后两项都是IIR系统不易实现的。FIR数字滤波器的设计方法有多种，如窗函数设计法、频率采样法和Chebyshev逼近法等。随着Matlab软件尤其是Matlab的信号处理工作箱的不断完善，不仅数字滤波器的计算机辅助设计有了可能，而

且还可以使设计达到最优化。

FIR 数字滤波器设计的基本步骤如下：

（1）确定技术指标

在设计一个滤波器之前，必须首先根据工程实际的需要确定滤波器的技术指标。在很多实际应用中，数字滤波器常被用来实现选频操作。因此，指标的形式一般在频域中给出幅度和相位响应。幅度指标主要以2种方式给出。第一种是绝对指标。他提供对幅度响应函数的要求，一般应用于FIR 滤波器的设计。第二种指标是相对指标。他以分贝值的形式给出要求。本文中滤波器的设计就以线性相位FIR 滤波器的设计为例。

（2）逼近

确定了技术指标后，就可以建立一个目标的数字滤波器模型（通常采用理想的数字滤波器模型）。之后，利用数字滤波器的设计方法（窗函数法、频率采样法等），设计出一个实际滤波器模型来逼近给定的目标。

（3）性能分析和计算机仿真

上两步的结果是得到以差分或系统函数或冲激响应描述的滤波器。根据这个描述就可以分析其频率特性和相位特性，以验证设计结果是否满足指标要求；或者利用计算机仿真实现设计的滤波器，再分析滤波结果来判断[4]。

2.2窗口设计法

窗口法设计的基本想法是要选取某一种合适的理想频率选择性滤波器（这种滤波器总是有一个非因果，无限长的脉冲响应），然后将它的脉冲响应截断（或加窗）以得到一个线性相位和因果的FIR 滤波器。因此，这种方法的重点在于选择某种恰当的窗函数和一个合适的理想滤波器。现用)(jw d e H 代表一理想频率选择性滤波器，它在整个通带内有单位幅度增益和线性相位特性，而阻带内有零响应[5]。

用窗口设计法基本步骤如下：

（1）构造希望逼近的频率响应函数()ωj d e H 。以低通线性相位FIRDF 设计为例，一般选择()

ωj d e H 为线性理想低通滤波器，即 ()???=-，

，0ωτωj j d e e H πωωωω≤<≤c c （2-1） 其中c w 也称为截止频率，α称为样本延迟。

（2）这个滤波器的脉冲响应应该具有无限长，求出()n h d 。对()ωj d e H 进行IFT 得到 )

()](sin[)(21

)(τπτπππ--==?-n n w dw e e H n h c jwn jw d d （2-2） 注意，)(n H d 是关于α对称的，这一点对于线性相位FIR 滤波器来说是有用的。

（3）为了从)(n H d 得到一个FIRl 滤波器必须在)(n H d 两边将它截断。为了得到一个长度为M 的因果且线性相位的FIR 滤波器)(n h ，就必须有

()?

??-≤≤=n 010)(，其他，M n n h n h d 和21-M =? （2-3） 这种运算叫做“加窗”，一般来说，()n h 可以当做是()n h d 和某一窗函数)(n w 相乘而得到的即

()()n w n h n h d =)( （2-4）

式中，()n w 称为窗口函数，其长度为N 。

窗口法设计的基本思想：对于给定的滤波器技术要求，选择滤波器长度M 和具有最窄主瓣宽度和尽可能小的旁瓣衰减的某个窗函数)(n w [6]。

常见的窗函数性能表如下表2-1所示。

表2-1 常见的窗函数性能表

2.3 GAUSSIAN 窗

高斯窗是一种指数窗。高斯窗谱无负的旁瓣，第一旁瓣衰减达-55dB 。高斯幅度谱的主瓣较宽，故而频率分辨力低。高斯窗函数常被用来截短一些非周期信号，如指数衰减信号等。

GAUSSWIN 的时域形式可以表示为：

??2)2/2/(211N N k e k w --=+α (2-5)

其中k 与α的取值范围必须满足0≤≤k N,≥α 2

以下是画高斯图和幅度响应图的主程序：

M=100;

n=0:M-1;

gauss=gausswin(M);

hd_bs=ideal_lp(0.5*pi,M);

h_bs=hd_bs.*gauss';

[db,mag,pha,grd,w]=freqz_m(h_bs,1);

高斯窗及其幅度响应图如图2-1所示。

00.10.20.30.40.5

0.60.70.80.91-100-50

w/pi d B 幅度响应0102030405060

708090100

n

w (n )高斯窗

图2-1高斯窗及其幅度响应图

3设计步骤

3.1设计流程图

滤波器设计流程图如图3-1所示。

图3-1 滤波器设计流程图

3.2 下载并截取弹拨音乐信号

从网上下载一段弹拨音乐，从中截取一段格式为.wav的语音信号，时间为4—5S，并加入噪声。在MALAB平台上，观察原始语音信号与加入噪声后的时域和频谱图。

原程序如下：

>>[x,fs,bits]=wavread('高山流水2.wav'); %读取音乐信号的数据

sound(x,fs,bits);

N=length(x); % 计算信号x的长度

fn=5000; % 单频噪声频率

t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率x=x(:,1);x=x';

y=x+ 0.1*sin(fn*2*pi*t); %加单频噪声

sound(y,fs,bits); % 应该可以明显听出有尖锐的单频啸叫声X=abs(fft(x));

Y=abs(fft(y)); % 对原始信号和加噪信号进行fft变换

X=X(1:length(X)/2);

Y=Y(1:length(Y)/2); %截取前半部分

deltaf=fs/2/length(X); %计算频谱的谱线间隔

f=0:deltaf:fs/2-deltaf; %计算频谱频率范围

figure (1);

原始信号与含噪信号的时域图和频域图如图3-2所示。

图3-2 原始信号和含噪的信号时域图和频谱图

3.3 滤波器设计

截取好原始信号，接下来的工作是设计一个GAUSSIAN 窗滤波器，

利用公式：

Rp=-(min(db(1:1:wp/delta_w+1))) （3-1）

As=-round(max(db(ws/delta_w+1:1:501))) （3-2）分别得到通带波纹及最小阻带衰减，并输入相应滤波器（高斯窗）的频率参数，接着用freqz求得滤波器的相应频率特性。

源程序如下：

fpd=4700;fsd=4850;fsu=5150;fpu=5300;Rp=1;As=60 %带阻滤波器设计指标

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(10.8*pi/dw)+1; %计算GAUSSWIN窗设计该滤波器时需要的阶数

n=0:M-1; %定义时间范围

gauss=gausswin(M);

hd_bs=ideal_lp(wcd,M)+ideal_lp(pi,M)-ideal_lp(wcu,M);% 调用自编函数计算理想带阻滤波器的脉冲响应

h_bs=gauss'.*hd_bs; % 用窗口法计算实际滤波器脉冲响应

[db,mag,pha,grd,w]=freqz_m(h_bs,1); %调用自编函数计算滤波器的频率特性figure (2);

滤波器幅度响应图、相位响应图、脉冲响应图如图3-3所示。

图3-3 滤波器幅度和相位响应图

3.4信号滤波处理

设计好滤波器后，我们要对语音信号进行滤波，对比滤波前和滤波后的语音信号。

原程序如下：

y_fil=fftfilt(h_bs,y); %用设计好的滤波器对y进行滤波

Y_fil=abs(fft(y_fil));

Y_fil=Y_fil(1:length(Y_fil)/2); %计算频谱取前一半

y_fil=fftfilt(h_bs,y); %用设计好的滤波器对y进行滤波

Y_fil=abs(fft(y_fil));Y_fil=Y_fil(1:length(Y_fil)/2); %计算频谱取前一半

figure (3);

弹拨音乐原始信号、加噪后信号与滤波后信号的时域、频域图如图3-4所示。

图3-4滤波前和滤波后的对比

3.5结果分析

通过观察上图，滤波后的语音信号发生了衰减，说明滤波器起到了滤波作用，同时通过频谱对比，可以看出滤波器滤掉了一部分频率范围内的信号。分别听原始语音和滤波后的语音信号，发现滤波后的语音信号噪声减小了，同时原始信号强度稍有减弱，基本达到了滤波的效果。

4出现的问题及解决方法

在本次课程设计中我们设计的滤波器对音乐信号的滤噪处理中，我遇到了许多细节疏忽和一些比较难处理的问题，当遇到问题时，我们不但运用自己学到的知识将其解决，还得到了同学老师的帮助，将其不断完善和修改。设计中出现的问题及解决方法如下：

1、在调用音乐时，用了语句[x,fs,bits]=wavread('高山流水2.wav');出现错误，后来得知我下载的高山流水高山流水音乐是MP3格式，不是wav格式，后来用千千静听软件将其改为wav格式之后就正确了。

2、在运用Matlab设计滤波器时，当编辑完前面两条程序时无法发出声音，后来发现我们应当把音乐信号wav文件放到Matlab的work文件夹中。

3、在对加入噪声后的音乐信号用FIR滤波器进行去噪处理时，所编写的函数出现错误，绘制不出图形，分析错误主要是因为调用FIR滤波器是所用的函数改为fftfilt函数。

4、在进行音乐信号的提取时没设采样频率，结果频率没达到要求，通过向老师询问，在windows自带录音机软件中改变音乐信号频率的属性，使时域图达到了要求。

5、绘制出滤波后的的波形，发现FIR滤波器没有滤掉单频噪声。通过自己仔细的检查，是单频噪声的频率改动后，FIR滤波器没有改动。所以单频噪声的频率也应该先定义，FIR滤波器的截止频率应该以单频信号的频率为中心，这样重新运行后，结果正确。

6、在采用gausswin 窗设计的FIR 滤波器时得不到理想的滤波器，把通带截止频率与阻带起始频率之间的差值设置的太小或者太大，后来通过适当的选择参数fp 和fs，绘制出来的图形效果比较明显，基本符合设计指标。这样通过MATLAB 运算出来的滤波器的阻带波纹达到了滤波器设计的要求，得到了比较理想的滤波器。

7、还要在滤波器性能曲线的wc处画一根竖线，这样更方便看出结果，其中wc处线的确定还需计算出wb/pi的值。

8、所有的时间波形横坐标都要化为时间，滤波前后频谱的横坐标应是频率，这样在观察同代截止频率和阻带截止频率时更加精确，误差较小。对利用gausswin函数设计FIR 滤波器的设计步骤很生疏,不知如何计算M 的值，自己上网查阅资料，或向图书馆借阅资料，掌握利用gausswin 窗函数设计FIR 滤波器的方法和步骤。通过同学的帮助，利用matlab 编写关于gausswin 的M 文件得到text5 函数，再通过 4 个不同参数代入text5 函数，计算出过渡带宽。

9、在对滤波器进行设计的过程中，由于程序较多，总会或多或少出现些错误，但是通过请教同学，在程序包中找到了freqz_m.m和ideal_lp.m这两个程序包，从而解决了对滤波器的设计问题。

10、最后当观察滤波前和滤波后的时域和频域波形对比图时，发现滤波后的信号频

谱仍有一小段噪声没有滤除，但是与未滤波前的噪声相比，这一小段噪声可以忽略，通过听取前后两段语音信号也没有出现任何问题。

5结束语

至此，本次课程设计就告一个段落了。这次课程设计中，是对上学期的专业课数字信号处理和MATLAB的一次复习，和通过实践加深了对所学知识的理解，所以这次课程设计对每位同学来说都是受益颇多的。

学习的过程是相互讨论共同进步的，多多讨论课题中遇到的问题，可以巩固我们的知识掌握能力，增加熟练运用度。

从本次课程设计的中心来看，课题是希望将数字信号处理技术应用于某一实际领域，这里就是指对音乐的处理。作为存储于计算机中的音乐信号，其本身就是离散化了的向量，我们只需将这些离散的量提取出来，就可以对其进行处理了。

在这里，用到了处理数字信号的强有力工具MATLAB，通过MATLAB里几个命令函数的调用，很轻易的在实际化音乐与数字信号的理论之间搭了一座桥。课题的特色在于它将音乐看作了一个向量，于是音乐数字化了，则可以完全利用数字信号处理的知识来解决。我们可以像给一般信号做频谱分析一样，来给音乐信号做频谱分析，也可以较容易的用数字滤波器来对音乐进行滤波处理。改变参数，理论结合实际，分析各参数对图形的影响，从而加深对各个参数的理解。在完成这次课程设计过程中学到了许多东西，进一步理解了滤波器设计方法和各参数意义，通过分析信号时域和频域的关系等，加深了对滤波性能的理解，而且学会了使用Matlab一些基本函数，增加了进一步学习Matlab软件的兴趣。同时，通过本次课程设计，锻炼了我的动手能力，和提高了我分析问题，解决问题的能力。

同学们在一起思考问题，通过亲自动手做一个FIR滤波器，发现了许多潜在的问题，而这些问题是在平常被我们所忽略的，甚至认为是不成问题的。让我再次感受到我们应当把所学知识和实践相结合，才能够超越水平，取得更大的发展。在这次课程设计中，胡双红老师给了我们很大的帮助和启发，随时随地为我们答疑，引导我们思考，而不是一味的接受。在这里向所有所有老师和同学给予的热心支持和帮助表示感谢，谢谢你们。

参考文献

[1] 程佩青．数字信号处理教程［M］．北京：清华大学出版社，2002

[2] 薛年喜主编．MATLAB 在数字信号处理中的应用．北京：清华大学出版社

[3] 维纳·K·恩格尔，约翰·G·普罗克斯.《数字信号处理》.西安交通大学出版社，2002

[4] 陈其宗.数字信号处理频谱计算与滤波器设计[M].北京：电子工业出版社，2002

[5] 刘敏，魏玲．Matlab通信仿真与应用[M]．北京：国防工业出版社，2001

[6] 百科ROBOT,zivenwong.窗函数.百度百科：

附录1：弹拨音乐信号滤波去噪设计源程序清单

% 程序名称：voicenoisefilter.m

% 程序功能：采用基于gausswin的窗口设计法，设计FIR滤波器对含噪语音进行滤波去噪处理。

% 程序作者：季政府

/% 最后修改日期：2012-3-2

[x,fs,bits]=wavread('高山流水2.wav');

sound(x,fs,bits);

N=length(x); % 计算信号x的长度

fn=5000; % 单频噪声频率

t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率

x=x(:,1);

x=x';

y=x+ 0.1*sin(fn*2*pi*t); %加单频噪声

sound(y,fs,bits); % 应该可以明显听出有尖锐的单频啸叫声

X=abs(fft(x)); Y=abs(fft(y)); % 对原始信号和加噪信号进行fft变换

X=X(1:length(X)/2); Y=Y(1:length(Y)/2); %截取前半部分

deltaf=fs/2/length(X); %计算频谱的谱线间隔

f=0:deltaf:fs/2-deltaf; %计算频谱频率范围

figure (1);

subplot(2,2,1);plot(t,x);grid on;

xlabel('时间（单位：s）');ylabel('幅值');title(' 原始语音信号');axis tight;

subplot(2,2,2);plot(f,X);grid on;

xlabel('频率（单位：Hz）');ylabel('幅度谱');title('语音信号幅度谱图');axis tight;

subplot(2,2,3);plot(t,y);grid on;

xlabel('时间（单位：s）');ylabel('幅值');title(' 加入单频干扰后语音信号');axis tight; subplot(2,2,4);plot(f,Y);grid on;

xlabel('频率（单位：Hz）');ylabel('幅度谱');title('加入干扰后的语音信号幅度谱图');axis([-inf inf 0 250]);

fpd=4700;fsd=4850;fsu=5150;fpu=5300;Rp=1;As=60; %带阻滤波器设计指标

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(10.8*pi/dw)+1; %计算GAUSSWIN窗设计该滤波器时需要的阶数

n=0:M-1; %定义时间范围

gauss=gausswin(M);

hd_bs=ideal_lp(wcd,M)+ideal_lp(pi,M)-ideal_lp(wcu,M);% 调用自编函数计算理想带阻滤波器的脉冲响应

h_bs=gauss'.*hd_bs; % 用窗口法计算实际滤波器脉冲响应

[db,mag,pha,grd,w]=freqz_m(h_bs,1); %调用自编函数计算滤波器的频率特性

figure (2);

subplot(2,2,1);plot(w/pi,db,'LineWidth',2);axis([-inf,inf,-80,10]);

a=[0 3.15];b=[-60 -60];line(a, b,'Color','r','LineStyle','--', 'LineWidth',2);

a=[0 3.15];b=[-1 -1];line(a, b,'Color','r','LineStyle','--', 'LineWidth',2);

a=[wsd/pi wsd/pi];b=[-150 10];line(a, b,'Color','r','LineStyle','--', 'LineWidth',2);

a=[wsu/pi wsu/pi];b=[-150 10];line(a, b,'Color','r','LineStyle','--', 'LineWidth',2);

xlabel('w/pi');ylabel('dB');title('滤波器幅度响应图');axis([0,1,-70,10]);grid on;

subplot(2,2,2);plot(w/pi,mag); grid on;

xlabel('w/pi');ylabel('幅度mag');title('滤波器幅度响应图');axis([0,1,0,1.2]);

subplot(2,2,3);plot(w/pi,pha); grid on;

xlabel('w/pi');ylabel('相位pha');title('滤波器相位响应图');axis tight;

subplot(2,2,4);stem(n,h_bs, 'LineWidth',2);grid on;

xlabel('n');ylabel('h(n)');title('滤波器脉冲响应图');axis([-inf,inf,-0.1,1.1]);

y_fil=fftfilt(h_bs,y); %用设计好的滤波器对y进行滤波

Y_fil=abs(fft(y_fil));Y_fil=Y_fil(1:length(Y_fil)/2); %计算频谱取前一半

y_fil=fftfilt(h_bs,y); %用设计好的滤波器对y进行滤波

Y_fil=abs(fft(y_fil));Y_fil=Y_fil(1:length(Y_fil)/2); %计算频谱取前一半

figure (3);

subplot(3,2,1);plot(t,x);grid on;

xlabel('时间t');ylabel('幅度');title('原始语音信号时间x');axis tight;

subplot(3,2,2);plot(f,X);grid on;

xlabel('频率f');ylabel('幅度');title('原始语音信号幅度谱X');axis tight;

subplot(3,2,3);plot(t,y);grid on;

xlabel('时间t');ylabel('幅度');title('加干扰语音信号时间x1');axis tight;

subplot(3,2,4);plot(f,Y);grid on;

xlabel('频率f');ylabel('幅度');title('加干扰语音信号幅度谱X1');axis([-inf inf 0 250]); subplot(3,2,5);plot(t,y_fil);grid on;

xlabel('时间t');ylabel('幅度');title('滤波后语音信号时间y');axis tight;

subplot(3,2,6);plot(f,Y_fil);grid on;

xlabel('频率f');ylabel('幅度');title('滤波后语音信号幅度谱Y');axis tight;

sound(x,fs,bits);