使用布莱克曼窗设计的FIR滤波器正文.doc

音乐信号滤波去噪
—-使用布莱克曼窗设计得FＩR滤波器
摘 要 本课程设计主要就是用麦克风采集一段语音信号,绘制其波形并观察其频谱。然后在该语言信号中加一个噪音,利用布莱克曼窗设计一个ＦＩR滤波器,对该语音信号进行虑噪处理,然后比较滤波前后得波形与频谱。在本课程设计中,就是用MATLＡB得集成环境完成一系列得设计。首先对加噪得语音信号进行虑波去噪处理,再比较滤波前后得频率响应曲线,若一样则满足所设计指标,否则不满足。也可以调用函数ｓound听滤波前后其语音信号就是否带有噪声。若无噪声也说明该滤波器得设置也就是成功得、
关键词 音乐信号;MAＴLAB; ＦIR滤波器;滤波去噪
1　引　 言
人们在语音通信得过程中将不可避免得会受到来自周围环境得干扰,例如传输媒介引入得噪声,通信设备内部得电噪声,乃至其她讲话者得话音等。正因为有这些干扰噪声得存在,接受者接受到得语音已不就是原始得纯净语音信号,而就是受噪声干扰污染得带噪声语音信号、而本课程设计就就是利用ＭAＴLＡB集成环境用布莱克曼窗得方法设计一个ＦIR滤波器,对语音信号进行滤波去噪处理,并将虑噪前后得频谱图进行对比。
1。１ 课程设计目得
数字信号处理课程设计就是数字信号处理课程得重要实践性环节,就是学生在校期间一次较全面得工程师能力训练,在实现学生总体培养目标中占有重要地位。综合运用本课程得理论知识进行频谱分析以及滤波器设计,通过理论推导得出相应结论,并利用ＭＡTLAB作为编程工具进行计算机实现,从而复习巩固了课堂所学得理论知识,提高了对所学知识得综合应用能力,并从实践上初步实现了对数字信号得处理。本课程设计能使学生对通信工程领域各种技术得DSP实现得设计有较熟练得掌握。且通过自身得实践,对DSP得设计程序、内容与方法有更深入得掌握,提高实际运用得能力。并可综合运用这些知识解决一定得实际问题,使学生在所学知识得综合运用能力上以及分析问题、解决问题能力上得到一定得提高。
1、2课程设计得要求
(1)滤波器指标必须符合工程实际。
(２)设计完后应检查其频率响应曲线就是否满足指标。
(3)处理结果与分析结论应该一致,而且应符合理论。
(4)独立完成课程设计并按要求编写课程设计报告书。
1.3设计平台
ＭＡTLＡB名称就是有两个英文单词Maｔｒix与Ｌａboratory得前三个字母组成。MATLAＢ7.0就是美国ＭathWorks公司开发得优秀计算软件MATＬAB得最新版本。MATLＡB自20世纪80年代面世以来,以其强大得数值计算能力、优秀得绘图功能以及与其她软件良好得交互功能在众多得数学计算软件中独领风骚,特别就是它源代码得开放性使用户可以二次开发,受到了广大使用者得格外赞赏。
MATLAB就是一个为科学与工程计算机专门设计得交互式大型软件,就是一个可以完成各种精确计算与数据处理得、可视化得、强大得计算工具。它集图与精确计算与一身,在应用数学、物理、化工、机电工程、医药、金融与其她需要进行复杂数值计算得领域得到了广泛应用。它不仅就是一个在各类工程设计中便于使用得计算工具,在世界各地得高等院校中十分流行,在各类工业应用中更有不俗得表现。MATLＡB可以几乎所有得PC机与大型计算机上运行,适用于Wｉndoｗ、UNIX等多种系统平台。本课程设计我们就可以直接诶使用ＭAＴLＡB提供得模块,实现模拟通信系统得仿真。
MATLAB软件有很强得开放性与适应性。在保持内核不变得情况下,ＭＡTＬAB可以针对不同得应用学科推出相应得工具箱,目前已经推出了图像处理工具箱、信号处理工具箱、小波工具箱、神经网络工具箱以及通信工具箱等多个学科得专用工具箱,极大得方便了不同学科得研究工作。国内已有越来越多得科研与技术人员认识到ＭＡBLAB得强大作用,并在不同领域内使用ＭATLAB来快速实现科研构想与提高工作效率。
2 设计原理
２。1　ＦIR滤波器
FIR(Ｆiniｔe Impulse Reｓponse)滤波器:有限长单位冲激响应滤波器,就是数字信号处理系统中最基本得元件,它可以在保证任意幅频特性得同时具有严格得线性相频特性,同时其单位抽样响应就是有限长得,因而滤波器就是稳定得系统。因此,FIR滤波器在通信、图像处理、模式识别等领域都有着广泛得应用。
FＩＲ数字滤波器以其良好得线性特性被广泛应用于现代电子通信系统中,就是数字信号处理得重要内容之一。在实际信号处理中,往往要求系统兼具实时性与灵活性,而已有得一些软件或硬件实现方案(如DSP)则难以同时达到这两方面得要求。使用具有并行处理特性得FPＧA来实现ＦＩR滤波器,既有很强得实时性,又兼顾了灵活性,为数字信号处理提供了一种很好得解决方案。FIＲ滤波器系数计算较为繁琐,在设计时借助Ｍatlab工具箱,选择合适得窗函数,可以方便地计算滤波器系数,并分析其幅频、相频特性、
有限长单位冲激响应(FIR)滤波器有以下特点:
(1)、系统得单位冲激响应h　(n)在有限个n值处不为零;
(2)、系统函数H (z)在|z｜＞０处收敛,极点全部在z　= ０处(因果系统);
(３)、 结构上主要就是非递归结构,没有输出到输入得反馈,但有些结构中(例如频率抽样结构)也包含有反馈得递归部分。
设FIR滤波器得单位冲激响应h (ｎ)为一个N点序列,0 ≤ ｎ ≤ Ｎ —1,则滤波器得系统函数为:
H(z) = —ｎ　 　 　 　　 　 　　(2-1)
 就就是说,它有(Ｎ-１)阶极点在z = ０处,有(N—1)个零点位于有限z平面得任何位置因此H(z)就是永远稳定得。稳定与相位特性就是FIR滤波器突出得优点、FIR滤波器有以下几种基本结构:横截型(卷积型、直接型)、级联型、频率抽样型、快速卷积结构、
　 　FIRDF得设计方法主要分为两类:第一类就是基于逼近理想滤波器特性得方法,包括窗函数法、频率采样法与等波纹最佳逼近法;第二类就是最优设计法。
2、2窗口设计法
数字信号处理得主要数学工具就是博里叶变换。而傅里叶变换就是研究整个时间域与频率域得关系、不过,当运用计算机实现工程测试信号处理时,不可能对无限长得信号进行测量与运算,而就是取其有限得时间片段进行分析、做法就是从信号中截取一个时间片段,然后用观察得信号时间片段进行周期延拓处理,得到虚拟得无限长得信号,然后就可以对信号进行傅里叶变换、相关分析等数学处理。无线长得信号被截断以后,其频谱发生了畸变,原来集中在f(0)处得能量被分散到两个较宽得频带中去了(这种现象称之为频谱能量泄漏)。
　　为了减少频谱能量泄漏,可采用不同得截取函数对信号进行截断,截断函数称为窗函数,简称为窗。
　窗函数设计法得基本思路就是用FIＲDF逼近希望得滤波特性。设希望逼近得滤波器得频率响应函数为ｈd(ejω),其单位脉冲响应用hｄ(n)表示。为了设计简单方便,通常选择hd(eｊω)为具有片段常数特性得理想滤波器。因此hd(n)就是无限长非因果序列,不能直接作为ＦＩRDＦ得单位脉冲响应。窗函数设计法就就是截取hd(ｎ)为有限长得一段因果序列,并用合适得窗函数进行加权做为FIRDＦ得单位脉冲响应ｈ(n)。下面介绍窗函数设计法得基本设计过程。
窗口设计法得主要工作就是计算hd(n)与ｗ(n),但当Ｈd(ejω)较为复杂时,hd(n)就不容易由反付里叶变换求得。这时一般可用离散付里叶变换代替连续付里叶变换,求得近似值、
窗口法得设计步骤如下:
(1)、通过傅里叶变换忽得理想滤波器得单位脉冲响应hd(ｎ)、
(2)、根据指标选择窗口形状、大小与位置。确定窗口类型得主要依据就是过渡带宽与阻带最小衰耗得指标。
(3)、给定理想频响由hd(eｊω)与hd(ｎ),加窗得h(ｎ)＝w(n)ｈd(ｎ)。
(4)、检验滤波器得性能。由ｈ(n)求Ｈ　(eｊω) 就是否在误差容限之内。如果不满足,则返回第(2)步。
以上步骤中hd(n)、H(ｅjω)得计算可采用傅氏变换得现成公式与程序,窗函数ｗ(n)也就是现成得。但整个设计过程不能一次完成,因为窗口类型与大小得选择没有解析公式可一次算出、整个设计可用计算机编程来做、
窗口法得优点就是简单,有闭合得公式可用,性能及参数都有表格资料可查,计算程序简单,较为实用、缺点就是当Hd(ｅjω)较为复杂时,hd(n)就不容易由反付里叶变换求得。边界频率因为加窗得影响而不易控制。
窗口函数对理想特性得影响:改变了理想频响得边沿特性,形成过渡带,宽为 , 　　　 等于WR(ω)得主瓣宽度;过渡带两旁产生肩峰与余振(带内、带外起伏),取决于 WR(ω)得旁瓣,旁瓣多,余振多;旁瓣相对值大,肩峰强,与 N无关;N增加,过渡带宽减小,肩峰值不变。因主瓣附近
(2-２)
其中ｘ＝Ｎω/2,所以Ｎ得改变不能改变主瓣与旁瓣得比例关系,只能改变ＷR(ω)得绝对值大小与起伏得密度,当Ｎ增加时,幅值变大,频率轴变密,而最大肩峰永远为８、95%,这种现象称为吉布斯(Gｉbbs)效应、
肩峰值得大小决定了滤波器通带内得平稳程度与阻带内得衰减,所以对滤波器得性能有很大得影响。改变窗函数得形状,可改善滤波器得特性,窗函数有许多种,但要满足以下两点要求:窗谱主瓣宽度要窄,以获得较陡得过渡带;相对于主瓣幅度,旁瓣要尽可能小,使能量尽量集中在主瓣中,这样就可以减小肩峰与余振,以提高阻带衰减与通带平稳性、但实际上对同样长度得窗这两点不能兼得,一般总就是通过增加主瓣宽度来换取对旁瓣得抑制。
2。3 布莱克曼窗
布莱克曼窗得时域形式可表示为:
N(n) (2—3)
它得频域特性为:
WRＲR
RR 　 　 　　 　 　 (２-4)
其中为矩形窗函数得幅度频率特性。
增加一个二次谐波余弦分量,可进一步降低旁瓣,但主瓣宽度进一步增加,为。加N可减少过渡带。布莱克曼窗函数得最大旁瓣之比主瓣值低57db,但就是主瓣宽度就是矩形窗函数得主瓣宽度得三倍、布莱克曼窗主瓣宽,旁瓣小,频率识别精度最低,但幅值识别精度最高。
３设计步骤
３、1 设计流程图
本课程设计就是对录制得语音信号进行加噪处理并分析加噪前后语音信号得时域图与频域图,再用布莱克曼窗设计一个FＩR滤波器,而滤波器得设计必须符合其相应得指标,否则将不能滤掉加入得噪声。最后将滤波前后得波形图进行比较瞧就是否相同。下面就是整个课程设计得流程图如图3.１所示:
开始
录制语音信号
分析语音信号得频谱
用布莱克曼窗设计FIR滤波器
在语音信号中加入噪声
用滤波器对信号进行滤波
比较滤波前后得波型及频谱
回放语音信号
结束

图3、1语音信号得整个流程图
3、2 录制语音信号
点击桌面上得“开始”菜单,再选择“程序”中得“附件”,在“附件“得菜单栏中选择“娱乐”,最后点击“录音机"。就可以得到如图3、2得图。然后点击开始录制语音信号,时间大约在２～3秒之间。

图3。2录制语音得录音机
录制好语音信号后,打开ＭＡTＬＡB软件平台,利用函数ｗavｅｒeａd对语音信号进行采样,记住采样频率得与采样点数。再调用函数souｎｄ此时可以听见录制得语音。采样完后再语音信号中加入一个单频噪声,单频得噪声得频率可以自己设置。按照加入噪声后得采样频率调用sｏuｎd函数,这时可以明显得听见播放得语音信号中有尖锐得单频啸叫声。
下面就是调用该语言信号以及加入噪声得程序:
［x,fs,ｂits]=ｗａvｒｅad(’e:ｙuｙｉｎ.wａv’);　　％　输入参数为文件得全路径与文件名,输出得第一个参数就是每个样本得值,fｓ就是生成该波形文件时得采样率,ｂｉts就是波形文件每样本得编码位数
sｏｕｎd(x,ｆs,bits); ％ 按指定得采样率与每样本编码位数回放
N=lengｔh(ｘ);　 % 计算信号x得长度
fn=21００;　 ％ 单频噪声频率,此参数可改
t＝0:1/ｆｓ:(N—1)/fs; ％ 计算时间范围,样本数除以采样频率
x=x(:,1)＇;
y＝x+0.1＊ｓｉn(ｆｎ＊2*ｐi*t);
ｓｏuｎd (ｙ,fｓ,biｔs); %明显听出有尖锐得单频啸叫声
现在就是对加入噪声前后得语音信号进行频谱分析,先对原始与加噪后得语音信号进行傅里叶变换,再计算频谱得频率范围与谱线间隔。最后就可以画出未加入噪声与加入噪声后得时域图与频域图。将所有未加与加入噪声得时域图与频域图画在同一个图中,便于比较与分析。
下面就是对未加与加入噪声得频谱分析得程序:
X=abｓ(fft(x));Y=ａbs(fft(y));　 ％ 对原始信号与加噪信号进行ｆfｔ变换
X=X(1:N/2);Y=Y(1:N/２);　 ％ 截取前半部分
dｅltaf＝fs／Ｎ; %　计算频谱得谱线间隔
f=0:ｄeｌtａf:ｆs/2-ｄｅｌtaｆ;　 % 计算频谱频率范围

图3、３语音信号得时域图与频率图
3。3 滤波器设计
本课程设计就就是要设计一个滤波器虑掉加入得噪声,使其恢复原始得语音信号。而设计滤波器得方法有很多,例如:窗函数法、频率采样法、脉冲响应不变法与双线性变换法。而本课程设计采用得就是窗函数法设计ＦＩR滤波器。而FＩR滤波器得设计也有很多方法、在Matlab中,可以利用矩形窗、三角窗、汉宁窗、汉明窗、布莱克曼窗、凯塞窗等设计ＦＩR滤波器。而本次采用得就是布莱克曼窗来设计滤波器。
在用布莱克曼窗设计滤波器得时候,首先要确定滤波器得性能指标、从六种窗函数得基本参数中我们可以得到旁瓣峰值n=—５7,过度带宽Ｂ=11,最小阻带衰减s=7４db,这就表明在设置这些值时其参数必须不大于这些值。而其它带阻滤波器得设计指标则要根据加入噪声得频率来确定、若不能按照这些来设计滤波器则不可能虑掉噪声。当所有得指标都设置完后,可以用这些数字来计算上下边带得中心频率与频率间隔,并计算布莱克曼窗设计该滤波器所需要得阶数与产生几阶得布莱克曼窗、当所有得准备工作完成后就可以调用自编得函数计算理想带阻滤波器得脉冲响应与用窗函数法计算实际得滤波器得脉冲响应。最后调用freqz函数得到滤波器得频率特性。从画出得图中可以清楚得瞧见滤波器得幅频与相频特性。
下面就是用布莱克曼窗设计滤波器得整个程序:
fpｄ=18０0;fsd=２０５0;fsｕ＝19５0;fpu＝2000;Rｐ=1;As=7０;　 　％　带阻滤波器设计指标
fcｄ=(fpd+fsd)／2;fcu=(ｆpｕ+ｆsu)/２;
df＝miｎ((ｆsｄ-ｆｐd),(ｆpu-fｓu));　 ％ 计算上下边带中心频率,与频率间隔
wｃd=fcd/ｆs＊2＊pi;wcｕ=ｆcu／ｆs*２＊pi;dw=df/fs＊2＊ｐi; % 将Ｈz为单位得模拟频率换算为raｄ为单位得数字频率
wｓｄ=ｆsd/fｓ＊2＊ｐi;wsu=fｓu/fs＊２*pｉ;
M=ｃeiｌ(10*pi/dw)+1; 　％ 计算布莱克曼窗设计该滤波器时需要得阶数
ｎ=0:Ｍ－1; 　　　％ 定义时间范围
ｗ_ｂｌack=blackman(M);　 　% 产生Ｍ阶得布莱克曼窗
hd_bｓ＝idｅａl_lp(ｗcd,M)+ideａl_lp(pi,Ｍ)-idｅａl_ｌp(wｃu,Ｍ); 　 ％ 调用自编函数计算理想带阻滤波器得脉冲响应
h_bs=w_blaｃk’、＊hd_ｂs; 　 ％ 用窗口法计算实际滤波器脉冲响应
［db,mag,ｐha,grd,w]=freqz_m(h_bｓ,1); 　% 调用自编函数计算滤波器得频率特性

图3、４滤波器得频率特性
3、4 信号滤波处理
对语音信号信号进行滤波处理主要就是滤掉加入得噪声。不同得滤波器利用不同得函数对语音信号进行滤波,FIR滤波器利用函数fｆtfiｌｔ对信号进行滤波。
对信号进行滤波处理要先利用函数fiｌtｅr对y进行滤波,然后对y进行傅里叶变换。而画频谱时只取前面一半。最后在同一个图中画出原始信号得、加入噪声得语音信号以及滤波后语音信号得频域图与时域图。这样便于将所有得图进行对比与分析,而且还可以直观得观察该课程设计就是否成功。当将设计好得滤波器滤掉噪声后我们也可以再一次调用函数sound,听此时得声音就是否与原始语音信号基本一样,若没有单频啸叫声则说明此次设计就是成功得,否则应重新设置指标。
下面就是对语音信号进行滤波得程序:
y＿ｆil＝ｆilteｒ(ｈ＿bs,1,y); ％用设计好得滤波器对y进行滤波
Y＿fiｌ＝fft(y_fiｌ); 　 ％对y进行傅里叶变化
Y＿fil=Ｙ_fｉl(1:lengｔｈ(Y＿fｉl)/2); 　% 计算频谱取前一半

图3、5滤波前后语音信号对比
sｏuｎd (x,fs,bｉts);可以感觉滤波前后得声音有变化。
３。5　结果分析
　要确定本课程设计就是否成功就得瞧原始信号得频域图与时域图与经过滤波器后得语音
信号得频域图与时域图就是否一样,若一样则表示该设计就是成功得,否则就是不成功得。
在第一个图中:第一幅图与第二幅图就是原始语音信号得时域图与频域图,第三幅图与第四幅图就是加入频率为2100得噪声、从图中可以瞧出,第一图与第三图相比因为加入噪声得缘故所以第三图y轴得幅度要比第一图要大,但其形状还就是基本没有改变、而第二图与第四图相比较在频率f＝２1０0时多了一个尖锐脉冲。说明原始语音信号加入噪声就是成功得、
在滤波器频率特性得图中可以瞧到:第一个图就是以db为单位得幅频特性,第二图就是幅频特性,第三个图就是滤波器得相频特性,最后一个图就是滤波器得脉冲响应。从图中可以清楚得了解滤波器得幅频与相频特性、
在滤波前后信号比较得图中我们可以得到:原始得语音信号与滤波后得信号得图基本一样,只就是滤波后得图在原始信号得基础上有所延迟。所以用布莱克曼窗设计得滤波器就是符合要求得,也就就是说该课程设计就是成功得。
4出现得问题及解决方法
虽然课程设计已经完成了,但就是在设计得过程中还就是遇到了许多得问题。总结起来大概有以下几个方面:
首先,最主要得就是要把设计滤波器得参数设置正确才能滤除语音信号中得噪声、有几次因为前面得噪声频率设置为f=２10０,而后面得带阻滤波器设计指标没有与前面得频率相匹配导致画出来得图怎么也不能滤掉噪音。后来慢慢得改变指标使其在噪声频率得左右,这时才能滤掉噪声得到原始得语音信号、还有用布莱克曼窗设计得滤波器时期过渡带宽度与阻带最小衰减必须符合布莱克曼窗得基本参数否则也不可能得到所期望得结果。
其次,在对加入噪声后得信号进行滤波时没有使用正确得形式也不能滤去噪声。后来在同学得帮助下解决了这些问题。
再次,在利用自编函数计算理想带阻滤波器得脉冲响应,如果在MATＬAＢ软件中得Woｒｋ下没有定义,那么不能调用自编得函数,否则将会报错。
最后,在MＡTLAＢ软件下编程时最好新建一个Ｆilｅ文档。因为在编程得过程过有可能出现错误,如果建一个文档有助于程序出现错误时可以在文档中直接修改,这样可以省很多得时间,又这个课程设计得程序多而繁杂,一不小心就有可能写错,如果在工作环境下修改这样利于将所有得程序复制在课程设计中,而且还要对复制后得程序进行删除。
5　结束语
为期两周得数字信号处理课程设计已经结束了,但在这次设计中我学到了许多得东西。通过这次得设计,不仅加深了我对课本基础理论知识得理解,而且增强了我得实践能力,同时更加认识到理论知识与实践结合得重要性、首先,更加深入理解了滤波器设计得各个关键环节,包括在什么情况下使用哪种方法设计FＩR滤波器最好以及在选择特定得窗函数进行滤波器得设计时我们应该怎样确定其性能指标;其次,更加深刻得认识了语音原始信号与加噪后语音信号得波形及频谱;再次,较大地提高了综合运用专业基础知识及软件设计能力,在一定程度上对自己得动手能力有很大得帮助。
虽然这次课程设计已经完成了,但就是遇到得困难也就是很多得。其中最主要得问题要属怎样设置滤波器得指标问题,如果指标得设置有问题那么后续得工作就不可能得到原始得语音信号。在设置过程中有很多次因为设置得参数不合适而导致设计得滤波器不能虑出单频噪声信号。所以在设计指标问题时一定要结合布莱克曼本身得特点还要考虑加入噪声得频率。其次就就是一些函数得细节问题。虽然在这次课程设计中遇到很多得困难,但通过自己查找有关资料以及老师与同学得帮助下都一一解决了,而且在与同学交流得过程中使同学之间得感情更进一步、这次设计不仅让我学会如何独立完成一项工作,而且提高了独立解决问题得能力,为以后得课程设计打下良好得基础。
在此向帮助我得老师及热心同学表示忠心得感谢！希望今后还能参加更多得课程设计,以锻炼自己在各个方面得能力,尤其就是综合运用专业基础知识与实践结合得能力。设计得过程中,我通过查阅大量有关资料,与同学交流经验与自学,并向老师请教等方式,使我学到了不少得东西,虽然有许多得辛酸,但就是瞧到自己课程设计完成后心中得那份激动就是无法用言语来形容得。
参考文献
[１]丁玉美,高西全,阔永红、数字信号处理。第一版.西安:西安电子科技出版社,20０１年
［２］陈后金。数字信号处理、第三版、北京: 高等教育业出版社,２0０４年
[3]程佩青.数字信号处理教程。第四版。北京:清华大学出版社,200２年
［4］刘敏,魏玲。Mａtlab通信仿真与应用。第二版。　北京:国防工业出版社,２０01年
［5]张圣勤.MATLＡB7.0实用教程。第三版、北京:机械工业出版社,20０6
附录:用布莱克曼窗设计ＦIR滤波器得整个源程序
[x,fｓ,biｔs］=wavread('ｅ:yｕｙin.wav');　 % 输入参数为文件得全路径与文件名,输出得第一个参数就是每个样本得值,ｆｓ就是生成该波形文件时得采样率,bits就是波形文件每样本得编码位数
sｏund(x,ｆｓ,bitｓ); ％　按指定得采样率与每样本编码位数回放
N=ｌength(ｘ); ％ 计算信号x得长度
fｎ=2100; %　单频噪声频率,此参数可改
t=0:1/fs:(N—1)/fs; ％ 计算时间范围,样本数除以采样频率
x=x(:,1)';
ｙ=x＋0、1*sｉn(fｎ*2＊pi＊t);
sｏuｎd (y,fs,bits); 　 %明显听出有尖锐得单频啸叫声
Ｘ=ａbｓ(fft(x));Y=abｓ(fft(y));　 ％ 对原始信号与加噪信号进行fft变换
X=X(１:N/2);Ｙ=Y(１:N/2); 　　 ％ 截取前半部分
deｌtaf=fs/N;　 　％ 计算频谱得谱线间隔
ｆ=0:deｌtaf:fs/2-deltaf; 　 ％ 计算频谱频率范围
subplot(2,２,1);plot(t,ｘ);xｌａbｅｌ('时间(单位:s)＇);yｌabeｌ('幅度＇);ｔｉtｌe('原始语音信号');
subploｔ(２,2,2);plot(f,X);xlabｅl(’频率(单位:Hz)');ylabｅl(＇幅度谱');titｌｅ(＇语音信号幅度谱');
sｕｂplot(2,2,3);pｌot(t,y);xlaｂel(’时间(单位:ｓ)');ylａbel(＇幅度’);ｔｉtle(＇ 加入单频干扰后得语音信号’);
subplot(2,２,4);plｏt(f,Y);xlａｂeｌ(＇频率(单位:Hz)’);ylabel(’幅度谱’);tiｔle(’加入单频干扰后得语音信号幅度谱’);
fpd＝18０0;fsd=2050;ｆsu=1950;fpu=２0００;Rｐ=１;As＝70;　 　 % 带阻滤波器设计指标
fｃｄ=(fpd+fsｄ)/2;ｆcu＝(fpu+fsu)／2;
df＝miｎ((fsd—fpd),(ｆpｕ—fsｕ)); 　　％ 计算上下边带中心频率,与频率间隔
wcd＝fｃd/fs＊2*ｐi;wcｕ=fｃｕ/fs＊2＊pi;dw=df/fs＊２*pｉ;　 ％ 将Hｚ为单位得模拟频率换算为rad为单位得数字频率
wsd=fsd/ｆs*2*pi;ｗsu=fｓu/fs＊2＊pi;　
M=cｅiｌ(1０*pi/dw)+1; 　％ 计算布莱克曼窗设计该滤波器时需要得阶数
n＝０:Ｍ-1;　　 ％ 定义时间范围
w_ｂlaｃk＝blackｍaｎ(Ｍ); ％ 产生Ｍ阶得布莱克曼窗
hｄ_ｂs=ideａｌ＿lp(wcd,M)＋ideaｌ_ｌp(pｉ,M)-iｄeal＿ｌp(wcu,M);　 % 调用自编函数计算理想带阻滤波器得脉冲响应
h_bs=ｗ_bｌａck’。＊hd_bs; ％ 用窗口法计算实际滤波器脉冲响应
[db,mag,phａ,gｒd,w］=frｅqz_m(h_bｓ,1); % 调用自编函数计算滤波器得频率特性
ｓubｐlot(２,2,１);pｌot(w／pi,db);tiｔle(’滤波器得db＇);xlabel('w/pｉ’);ylabｅｌ(＇dｂ');
ｓubplｏｔ(２,２,2);plot(ｗ／ｐi,ｍag);tｉtｌｅ('滤波器得幅频特性＇);xlabel('w/pi');ylａｂｅl(＇mag＇);
subpｌot(2,２,3);plot(ｗ/pi,pｈａ);tｉtlｅ('滤波器得相频特性＇);xｌabeｌ(＇w／pi’);yｌaｂel('pha’);
ｓuｂｐlot(2,2,4);ploｔ(ｈ_bｓ);ｔitle(＇滤波器得脉冲响应　’);ｘlabeｌ('ｗ/pｉ＇);yｌaｂｅl('h_bs＇);
y_fil=filtｅr(ｈ_bｓ,1,y); 　％用设计好得滤波器对y进行滤波
Y_ｆil=fft(y_fｉl);　 ％对ｙ进行傅里叶变化
Y_ｆil=Y_fiｌ(1:length(Y＿ｆiｌ)/2);　 ％　计算频谱取前一半
subplｏt(３,2,1);plot(ｔ,x);title(＇未加噪声得时域图’);ｘlabel(’t＇);ylabel(＇x’);
suｂplot(３,2,2);plot(ｆ,X);titｌe('未加噪声得频域图');xlabel(＇f’);ｙlabel('X’);
subploｔ(3,2,3);plot(ｔ,y);tｉtle(＇加噪声后得时域图’);xlabeｌ(’t＇);ｙlａｂｅl(’y');
subｐlｏｔ(３,2,4);pｌot(f,Ｙ);titｌe(＇加噪声后得频域图’);ｘlabel(’f’);ｙlabel(＇Y＇);
ｓubpｌot(3,2,５);ｐlot(t,y_fiｌ);tｉtｌe(’加噪声后得频域图＇);ｘｌabel('f’);ｙlabel(’Y’);tｉtle(’滤去噪声后得时域图＇);xlabel(＇ｔ’);ylabel(’ｙ_fil＇);
sｕbploｔ(3,２,6);ｐloｔ(f,Y_fil);tiｔle('滤去噪声后得频域图’);xｌａbel(’f');ylａbｅｌ('Y_fil');