语音信号的分析及滤波

课程设计说明书PAGE71引言1.1数字滤波器的背景及意义数字滤波器(digitalfilter)是由数字乘法器、加法器和延时单元组成的一种装置。其功能是对输入离散信号的数字代码进行运算处理，以达到改变信号频谱的目的。由于电子计算机技术和大规模集成电路的发展，数字滤波器已可用计算机软件实现，也可用大规模集成数字硬件实时实现。数字滤波器广泛用于数字信号处理中，如电视、VCD、音响等。数字滤波器是一个离散时间系统（按预定的算法，将输入离散时间信号转换为所要求的输出离散时间信号的特定功能装置）。应用数字滤波器处理模拟信号时，首先须对输入模拟信号进行限带、抽样和模数转换。数字滤波器输入信号的抽样率应大于被处理信号带宽的两倍，其频率响应具有以抽样频率为间隔的周期重复特性，且以折叠频率即1／2抽样频率点呈镜像对称。为得到模拟信号，数字滤波器处理的输出数字信号须经数模转换、平滑。数字滤波器具有高精度、高可靠性、可程控改变特性或复用、便于集成等优点。数字滤波器在语言信号处理、图像信号处理、医学生物信号处理以及其他应用领域都得到了广泛应用。数字滤波器有低通、高通、带通、带阻和全通等类型。它可以是时不变的或时变的、因果的或非因果的、线性的或非线性的。应用最广的是线性、时不变数字滤波器。数字滤波器是指完成信号滤波处理功能的，用有限精度算法实现的离散时间线性非时变系统，其输入是一组数字量，其输出是经过变换的另一组数字量。因此，它本身即可以是用数字硬件装配成的一台完成给定运算的专用数字计算机，也可以是将所需运算编成程序，让通用计算机来执行。数字滤波器具有稳定性高、精度高、灵活性大等优点。随着数字技术的发展，用数字技术实现滤波器的功能越来越受到人们的注意和广泛的应用。1.2设计要求1.完成语音信号的采集，利用windows自带的录音机或其他软件，录制一段语音，时间在1s内；2.进行语音信号的频谱分析；3.进行数字滤波器的设计，滤波器的性能指标可以根据实际情况作调整，要求用窗函数法和双线性变换法设计以下三种数字滤波器：(1)低通滤波器:，，最大衰减，最小衰减；(2)高通滤波器:，，最大衰减，最小衰减；(3)带通滤波器:，，，，最大衰减，最小衰减；4.对语音信号进行滤波处理；5.对滤波前后的语音信号频谱进行对比，并对设计结果进行独立思考和分析；6.在基本要求的基础上，学生可以根据个人对该课程设计的理解，添加一些新的内容，如设计系统的人机对话界面。1.3设计目的掌握双线性变换法及窗函数设计IIRFIR数字滤波器的具体设计方法及其原理，观察双线性变换及脉冲响应不变法设计的滤波器的频域特性，了解双线性变换法及窗函数法的特点。2数字滤波器的设计2.1数字滤波器的基本概念2.1.1数字滤波器的分类数字滤波器从实现的网络结构或者从单位脉冲响应分类，可以分成无限脉冲响应(IIR)滤波器和有限脉冲响应(FIR)滤波器。它们的系统函数分别为：2.1.2数字滤波器的技术要求我们通常用的数字滤波器一般属于选频滤波器。假设数字滤波器的传输函数用下式表示：通带内和阻带内允许的衰减一般用dB数表示，通带内允许的最大衰减用Ap表示，阻带内允许的最小衰减用As表示，Ap和As分别定义为：如将归一化为1后得：2.2利用窗函数法设计FIR数字滤波器窗函数设计法是一种通过截短和计权的方法使无限长非因果序列成为有限长脉冲响应序列的设计方法。通常在设计滤波器之前，应该先根据具体的工程应用确定滤波器的技术指标。在大多数实际应用中，数字滤波器常常被用来实现选频操作，所以指标的形式一般为在频域中以分贝值给出的相对幅度响应和相位响应。用窗函数法设计FIR滤波器的步骤如下：(1)根据过渡带宽及阻带衰减要求，选择窗函数的类型并估计窗口长度N（或阶数M=N-1）。窗函数类型可根据最小阻带衰减AS独立选择，因为窗口长度N对最小阻带衰减AS没有影响。在确定窗函数类型以后，可根据过渡带宽小于给定指标确定所拟用的窗函数的窗口长度N。设待求滤波器的过渡带宽为△ω，它与窗口长度N近似成反比。窗函数类型确定后，其计算公式也确定了，不过这些公式是近似的，得出的窗口长度还要在计算中逐步修正。原则是在保证阻带衰减满足要求的情况下，尽量选择较小的N。在N和窗函数类型确定后，即可调用MATLAB中的窗函数求出窗函数。(2)根据待求滤波器的理想频率响应求出理想单位脉冲响应。如果给出待求滤波器的频率响应为，则理想的单位脉冲响应可以用傅里叶反变换式得到，在一般情况下，是不能用封闭公式表示的，需要采用数值方法表示。从ω=0到ω=2π采样N点，采用离散傅里叶反变换(IDFT)即可求出。(3)计算滤波器的单位脉冲响应h(n)。它是理想单位脉冲响应和窗函数的乘积，即，在MATLAB中用点乘命令表示为。(4)验算技术指标是否满足要求。为了计算数字滤波器在频域中的特性，可调用freqz子程序，如果不满足要求，可根据具体情况，调整窗函数类型或长度，直到满足要求为止。使用窗函数法设计时要满足以下两个条件：(1)窗谱主瓣尽可能地窄，以获得较陡的过渡带；(2)尽量减少窗谱的最大旁瓣的相对幅度，也就是使能量尽量集中于主瓣，减小峰肩和纹波，进而增加阻带的衰减。在实际工程中常用的窗函数有五种，即矩形窗、三角窗、汉宁窗、海明窗和凯泽窗。2.3利用双线性变换法设计IIR数字滤波器2.3.1变换原理双线性变换法是使数字滤波器的频率响应与模拟滤波器的频率响应相似的一种变换方法。为了克服脉冲响应不变法的多值映射这一缺点，我们首先把整个s平面压缩变换到某一中介的平面的一横带里（宽度为，即从到），然后再通过上面讨论过的标准变换关系将此横带变换到整个z平面上去，这样就使s平面与z平面是一一对应的关系，消除了多值变换性，也就消除了频谱混叠现象，基本原理如图所示。将s平面整个平面压缩到平面的到，可采用以下的变换关系其中C为常数；这样变为，变为，可将上式写成图双线性变换法的映射关系图双线性变换法的映射关系令，则可得再将平面通过以下标准变化关系映射到z平面后可得2.3.2变换常数C的选择为了使模拟滤波器与数字滤波器在低频处有较确切的对应关系，即在低频处有，当较小时有由上式可知因而得则Z和S的表达式可重新写成：即2.3.3模拟滤波器的数字化由于双线性变换法中，s与z之间有简单的代数关系，故可由模拟系统函数通过代数置换直接得到数字滤波器的系统函数。即可见数字滤波器的极点数等于模拟滤波器的极点数。频率响应也可用直接置换得到这一公式可用于将滤波器的数字域指标，转换为模拟域指标。3基于MATLAB的数字滤波器的设计3.1MATLAB软件介绍Matlab是由matrix和laboratory两个英文单词的前三个字母组合而成，是“MatrixLaboratory”的缩写，意为“矩阵实验室”，它是MathWorks公司于1984年推出的一套高性能的数值计算可视化软件，集数值分析、矩阵运算、信号处理和图形显示于一体，被誉为“巨人肩上的工具”，是当今美国很流行的科学计算软件。信息技术、计算机技术发展到今天，科学计算在各个领域得到了广泛的应用．在许多诸如控制论、时间序列分析、系统仿真、图像信号处理等方面产生了大量的矩阵及其相应的计算问题。自己去编写大量的繁复的计算程序，不仅会消耗大量的时间和精力，减缓工作进程，而且往往质量不高．美国Mathwork软件公司推出的Matlab软件就是为了给人们提供一个方便的数值计算平台而设计的。Matlab是一个交互式的系统，它的基本运算单元是不需指定维数的矩阵，按照IEEE的数值计算标准（能正确处理无穷数Inf(Infinity)、无定NaN(not-a-number)及其运算）进行计算。系统提供了大量的矩阵及其它运算函数，可以方便地进行一些很复杂的计算，而且运算效率极高。Matlab命令和数学中的符号、公式非常接近，可读性强，容易掌握，还可利用它所提供的编程语言进行编程完成特定的工作。除基本部分外，Matlab还根据各专门领域中的特殊需要提供了许多可选的工具箱，如应用于自动控制领域的ControlSystem工具箱和神经网络中NeuralNetwork工具箱等。MATLAB主要特点：(1)语言简洁紧凑，使用方便灵活，库函数极其丰富。MATLAB程序书写形式自由，利用起丰富的库函数避开繁杂的子程序编程任务，压缩了一切不必要的编程工作。由于库函数都由本领域的专家编写，用户不必担心函数的可靠性。可以说，用MATLAB进行科技开发是站在专家的肩膀上。(2)运算符丰富。由于MATLAB是用C语言编写的，MATLAB提供了和C语言几乎一样多的运算符，灵活使用MATLAB的运算符将使程序变得极为简短。(3)MATLAB既具有结构化的控制语句（如for循环，while循环，break语句和if语句），又有面向对象编程的特性。(4)程序限制不严格，程序设计自由度大。例如，在MATLAB里，用户无需对矩阵预定义就可使用。(5)程序的可移植性很好，基本上不做修改就可以在各种型号的计算机和操作系统上运行。(6)MATLAB的图形功能强大。在FORTRAN和C语言里，绘图都很不容易，但在MATLAB里，数据的可视化非常简单。MATLAB还具有较强的编辑图形界面的能力。(7)MATLAB的缺点是，它和其他高级程序相比，程序的执行速度较慢。由于MATLAB的程序不用编译等预处理，也不生成可执行文件，程序为解释执行，所以速度较慢。(8)功能强大的工具箱是MATLAB的另一特色。MATLAB包含两个部分：核心部分和各种可选的工具箱。核心部分中有数百个核心内部函数。其工具箱又分为两类：功能性工具箱和学科性工具箱。功能性工具箱主要用来扩充其符号计算功能，图示建模仿真功能，文字处理功能以及与硬件实时交互功能。功能性工具箱用于多种学科。而学科性工具箱是专业性比较强的，control,toolbox,signlproceessingtoolbox,commumnicationtoolbox等。这些工具箱都是由该领域内学术水平很高的专家编写的，所以用户无需编写自己学科范围内的基础程序，而直接进行高、精、尖的研究。(9)源程序的开放性。开放性也许是MATLAB最受人们欢迎的特点。除内部函数以外，所有MATLAB的核心文件和工具箱文件都是可读可改的源文件，用户可通过对源文件的修改以及加入自己的文件构成新的工具箱。3.2MATLAB常用函数简介(1)wavread函数：[y,fs,bits]=wavread('Blip')用于读取语音，采样值放在向量y中，fs表示采样频率(Hz)，bits表示采样位数。sound(x,fs,bits);用于对声音的回放。向量y则就代表了一个信号（也即一个复杂的“函数表达式”）也就是说可以像处理一个信号表达式一样处理这个声音信号。(2)FFT函数：FFT用于序列快速傅立叶变换。其调用格式为：y=fft(x)。其中，x是序列，y是序列的FFT，x可以为一向量或矩阵：若x为一向量，y是x的FFT，且和x相同长度。若x为一矩阵，则y是对矩阵的每一列向量进行FFT。(3)Fir1函数：fir1函数是基于窗函数的FIR滤波器设计——标准频率响应形状。其调用格式为b=fir1(N,we,’ftype’,window)。ftype和window可以默认。b=fir1(N,wc)可得到截止频率为wc且满足线性相位条件的N阶FIR低通滤波器，window默认选用hamming窗。其单位脉冲响应h(n)为：h(n)=b(n+1)，n=0,1,2,…,N。当wc=[wc1,wc2]时，得到的是通带为wc1<w<wc2的带通滤波器。b=fir1(N,wc,’ftype’)，可设计高通和带阻滤波器。当ftype=high时，设计高通FIR滤波器；当ftype=stop时，设计带阻FIR滤波器。(4)sound(x,fs,bits)：该函数用于播放生硬，向量y就代表了一个信号（也即一个复杂的“函数表达式”），也就是说可以像处理一个信号表达式一样处理这个声音信号。3.3FIR滤波器设计3.3.1用kaiser窗设计FIR滤波器的整体流图如图1所示：作出kaiser作出kaiser相应滤波器的幅频响应曲线确定模拟滤波器的技术指标技术指标从模拟到数字的转换计算中心频率，过渡带宽，设定旁瓣高度alfa利用fir1函数确定滤波器的实际单位脉冲响应h(n)利用freqz函数确定幅频响应开始结束图1kaiser窗滤波器设计流图3.3.2FIR低通滤波器设计其基本程序为：fb=1000;fc=1200;%确定技术指标alfa=60;Fs=22050;%alfa=50wb=2*pi*fb/Fs;%求归一化频率wc=2*pi*fc/Fs;wo=wc-wb;%确定过渡带beta=0.1102*(alfa-0.87);%求凯泽窗函数参数N=ceil((alfa-8)/2.285/wo);%由过渡带决定N的值ww=(wb+wc)/2/pi;b=fir1(N,ww,kaiser(N+1,beta));[Hw,w]=freqz(b,1);显示出其滤波器的特性分析如图2所示：图2滤波器频谱分析读入信号对信号傅里叶变换进行频谱分析，得到滤波前后的波形与频谱图如图3所示：图图3滤波前后波形图与频谱图比较分析滤波前后信号的变化，波形图变窄，频谱变化表明:信号的低频率段被保留，高频率段被滤除。分析滤波前后的声音变化，滤波后声音明显变低而沉闷，这说明滤波器设计基本符合指标要求。3.3.3FIR带通滤波器设计其基本程序为：wp1=2*pi*fc1/Fs;ws1=2*pi*fb1/Fs;wp2=2*pi*fb2/Fs;ws2=2*pi*fc2/Fs;%求归一化频率wo1=wp1-ws1;%求过渡带wo2=ws2-wp2;wo=max(wo1,wo2);beta=0.1102*(alfa-0.87);%求凯泽窗的参数N=ceil((alfa-8)/2.285/wo);%有过渡带决定Nwn1=(fb1+fc1)/Fs;wn2=(fb2+fc2)/Fs;wn=[wn1,wn2];b=fir1(N,wn,'band',kaiser(N+1,beta));[Hw,w]=freqz(b,1,512,Fs);Kaiser窗带通滤波器的特性分析如图4所示：图图4滤波器性能分析读入信号对信号傅里叶变换进行频谱分析，得到滤波前后的波形与频谱图如图5所示：图图5滤波前后波形图与频谱图比较分析滤波前后信号的变化，波形图变窄。频谱变化表明:信号的低频率段和高频率段被滤除，中间频率段被保留；分析滤波前后的声音变化，滤波后声音一定程度上变的尖锐，与高通滤波器滤波后的声音相比较低。这说明滤波器设计基本符合指标要求。3.3.4FIR高通滤波器设计其基本程序为：wb=2*pi*fb/Fs;%求归一化频率ws=2*pi*fs/Fs;wo=wb-ws;%过渡带beta=0.1102*(alfa-0.87);%求凯泽窗函数参数N=ceil((alfa-8)/2.285/wo);%有过渡带决定Nwc=(wb+ws)/2/pi;b=fir1(N,wc,'high',kaiser(N+1,beta));[Hw,w]=freqz(b,1);显示出其滤波器的特性分析如图6所示：图图6滤波器性能分析读入信号对信号傅里叶变换进行频谱分析，得到滤波前后的波形与频谱图如图7所示：图图7滤波前后波形图与频谱图比较分析滤波前后信号的变化，波形图变窄，频谱变化表明:信号的高频率段被保留，低频率段被滤除。分析滤波前后的声音变化，滤波后声音明显变高而尖锐，这说明滤波器设计基本符合指标要求。3.4IIR滤波器设计3.4.1用双线性变换法设计IIR滤波器的整体流图如图8所示：模拟低通滤波器的系统函数模拟低通滤波器的系统函数计算归一化滤波器的零点，极点和增益确定滤波器的阶数N和wc技术指标从数字到模拟的转换归一化，求得数字边缘频率令wc=wap，反归一化双线性变换法：实现H(s)到H(z)的转化作出(z)的幅频特性曲线确定模拟滤波器的技术指标H(z)的幅频特性结束开始图图8用双线性变换法设计IIR滤波器的流图3.4.2IIR低通滤波器设计其基本程序为：wb=2*pi*fb/Fs;%归一化，求得数字边缘频率ws=2*pi*fc/Fs;wap=2/Ts*tan(wb/2);%将数字指标转换成模拟指标was=2/Ts*tan(ws/2);[N,wc]=buttord(wap,was,Ap,As,'s');%选择滤波器的最小阶数[z,p,k]=buttap(N);%创建butterworth模拟滤波器[Bap,Aap]=zp2tf(z,p,k);[b,a]=lp2lp(Bap,Aap,wc);[bz,az]=bilinear(b,a,Fs);%用双线性变换法实现模拟滤波器到数字滤波器的转换[H,w]=freqz(bz,az);显示出其滤波器的特性分析如图9所示：图图9滤波器性能分析读入信号对信号傅里叶变换进行频谱分析，与滤波前进行比较。并在滤波前后用sound函数播放语音信号，观察滤波前后的变化。结果如图10所示：图图10滤波前后波形图与频谱图比较分析滤波前后信号的变化，波形图变窄，频谱变化表明:信号的低频率段被保留，高频率段被滤除。分析滤波前后的声音变化，滤波后声音明显变低而沉闷，这说明滤波器设计基本符合指标要求。3.4.3其基本程序为：wp1=2*pi*fp1/Fs;wp2=2*pi*fp2/Fs;ws1=2*pi*fs1/Fs;ws2=2*pi*fs2/Fs;%归一化，求得数字边缘频率Ts=1/Fs;wap1=2/Ts*tan(wp1/2);wap2=2/Ts*tan(wp2/2);%将数字指标换成模拟指标was1=2/Ts*tan(ws1/2);was2=2/Ts*tan(ws2/2);wap=[wap1,wap2];was=[was1,was2];W1=(wap1+was1)/2;W2=(wap2+wap2)/2;Bw=W2-W1;%带宽wo=sqrt(W1*W2);[N,wc]=buttord(wap,was,Ap,As,'s');%选择滤波器的最小阶数[z,p,k]=buttap(N);%创建巴特沃兹模拟滤波器[Bap,Aap]=zp2tf(z,p,k);[b,a]=lp2bp(Bap,Aap,wo,Bw);[bz,az]=bilinear(b,a,Fs);%用双线性变换法实现模拟滤波器到数字滤波器的转换[H,w]=freqz(bz,az);显示其滤波器的特性分析如图11所示：图图11滤波器的性能分析读入信号对信号傅里叶变换进行频谱分析，与滤波前进行比较。并在滤波前后用sound函数播放语音信号，观察滤波前后的变化。结果如图12所示：图图12滤波前后波形图与频谱图比较分析滤波前后信号的变化，波形图变窄。频谱变化表明:信号的低频率段和高频率段被滤除，中间频率段被保留；分析滤波前后的声音变化，滤波后声音一定程度上变的尖锐，与高通滤波器滤波后的声音相比较低。这说明滤波器设计基本符合指标。3.4.4IIR高通滤波器设计其基本程序为：wb=2*pi*fb/Fs;%归一化，的边缘数字频率ws=2*pi*fs/Fs;wap=2/Ts*tan(wb/2);%数字指标转换成模拟指标was=2/Ts*tan(ws/2);[N,wc]=buttord(wap,was,Ap,As,'s');%选择滤波器的最小阶数[z,p,k]=buttap(N);%创建巴特沃兹模拟滤波器[Bap,Aap]=zp2tf(z,p,k);[b,a]=lp2hp(Bap,Aap,wc);[bz,az]=bilinear(b,a,Fs);%用双线性变换法实现模拟滤波器到数字滤波器的转换[H,w]=freqz(bz,az);显示其滤波器的特性分析如图13所示：图图13波器的性能分析读入信号对信号傅里叶变换进行频谱分析，与滤波前进行比较。并在滤波前后用sound函数播放语音信号，观察滤波前后的变化。结果如图14所示：图图14滤波前后波形图与频谱图比较分析滤波前后信号的变化，波形图变窄，频谱变化表明:信号的高频率段被保留，低频率段被滤除。分析滤波前后的声音变化，滤波后声音明显变高而尖锐，这说明滤波器设计基本符合指标要求。4图形用户界面GUI4.1GUI简介图形用户界面GUI是一种提供人机交互的工具和方法。GUI是包含图形对象，如窗口，图标,菜单和文本等图文并茂的用户界面。一个优秀的GUI能非常直观的使用户知道如何操作MATLAB界面，了解设计者的开发意图。用MATLAB图形用户界面开发环境设计GUI的一般步骤是：(1)进行界面设计；(2)设置控件属性；(3)进行M语言编译；4.2GUIDE及其组成部分GUIDE提供一个界面设计工具集实现图形用户界面的创建工作，这些工具包括：界面设计编辑器：添加并排列图形窗口的组件对象；属性检查器：检查并设置组件的属性值；对象浏览器：观察此次的MATLAB运行过程中图形对象的句柄继承关系表；菜单编辑器：创建窗口菜单和文本菜单。以上这些工具通过界面设计编辑器的相应菜单进行调用，使用GUIDE命令来启动界面设计编辑器，界面设计编辑器的外观如图15所示：图图15GUI编辑器的外观4.3滤波器用GUI界面显示的设计滤波器共有六个，一个滤波器出五个图形，故设置六个pushbutton键，五个axes图形窗。每一个部分都回调相应的滤波器，分别用来观察滤波器的幅频特性和信号在滤波前后的时域波形和幅频特性。设计好并运行的GUI界面如图16所示：图16设计并运行的GUI界面图16设计并运行的GUI界面5总结通过这次课设，我认识到了MATLAB功能非常的强大，使得我们在使用的时候用户直接调用这些库函数并赋予实际参数就能解决实际问题，具有极高的变成效率。我也熟悉了MATLAB的工作环境，可以很熟练的对MATLAB进行常规的操作，快速进行程序编辑和仿真。本文次课设通过一个设计实例，利用MATLAB实现IIR数字滤波器设计与滤波的三种方法，从仿真结果可以看出它们均可以达到技术指标要求，而且方法简单、快捷，大大减轻了工作量。滤波器的设计工作完成后，可以借助于MATLAB的export操作导出所设计滤波器的系统函数H(z)。由于MATLAB具有强大的接口功能，仿真后的结果可以很方便的移植到DSP、CPLD。在实际应用中,只需按要求修改滤波器参数,并对程序作较少的改动,即可实现不同截止频率的FIR滤波器,实用性较强。完成任务后才发现这次训练的内容并不算难。首先要做的工作就是把设计条件写入到MATLAB中，其次就是对条件进行分析，主要是得出一些参数。最后就是对程序经行仿真，这里采用的是FDATool,Simulink仿真。以上的内容全都要通过程序的编写来完成，不过主要的操作是对函数的调用，这就要我们熟悉MATLAB的一些函数的基本功能以及它们的基本语句。当然，经过了课程设计，我也发现了自己的很多不足。但是通过自己的动手动脑，既增加了知识，又给了我专业知识以及专业技能上的提升，我也会更加努力，认真学习，争取在以后的课程中做得更好！在课设的过程中我们都遇到了很多的难题，请教老师都给了细心地指导，帮助我们解决了很多的问题。老师还教我们意识到自己动手对自己的成长和发展过程中的重要性，还敦促我们多上网查找问题的解决方法，使我们的在网上获取信息的能力得到了增长，再此特别感谢董老师和郭老师对我们的认真指导和教诲，希望老师们生活幸福并向老师保证，我一定继续学习自己不会的知识，不让老师失望。参考文献[1]范寿康.DSP技术与DSP芯片.北京:电子工业出版社，2008[2]王艳芬.数字信号处理原理及实现.清华大学出版社，2008[3]王金龙.DSP设计与实验教程.机械工业出版社，2007[4]赵健,李毅.数字信号处理学习辅导.北京:清华大学出版社出版,2006[5]丁玉美.数字信号处理学习指导与题解.北京:电子工业出版社出版,2007[6]丛玉良.数字信号处理原理及其MATLAB实现.北京:电子工业出版社,2005[7]陈生谭,郭宝龙等.信号与系统.西安:西安电子科技大学出版社,2001附录%FIR低通滤波器fb=1000;fc=1200;%确定技术指标alfa=60;Fs=22050;%alfa>50wb=2*pi*fb/Fs;%求归一化频率wc=2*pi*fc/Fs;wo=wc-wb;%确定过渡带beta=0.1102*(alfa-0.87);%求凯泽窗函数参数N=ceil((alfa-8)/2.285/wo);%由过渡带决定N的值ww=(wb+wc)/2/pi;b=fir1(N,ww,kaiser(N+1,beta));[Hw,w]=freqz(b,1);figure(1);plot(w*Fs/2/pi,20*log10(abs(Hw)));%绘制频率响应曲线grid;title('FIR低通滤波器');ylabel('幅度/db');xlabel('频率/Hz');axis([0,12000,-180,50]);x1=wavread('3.wav');%利用函数wavread对语音信号进行采样sound(x1,22050);%回放语音信号x=x1(:,1);y=filter(b,1,x);%语音信号经滤波器滤波sound(y,22050);%播放滤波后的信号figure(2);subplot(2,2,1);plot(x);%输出原始语音信号的脉冲响应title('输入信号时域波形');xlabel('f/Hz');ylabel('幅度');subplot(2,2,2);plot(y(:,1));%输出原始语音信号经滤波后的的脉冲响应title('经低通滤波器后的时域波形');xlabel('f/Hz');ylabel('幅度');X=fft(x,1024);%对原始信号进行傅立叶变换并进行1024点采样Y=fft(y,1024);%对滤波后的信号进行傅立叶变换并进行1024点采样f=Fs*[0:1023]/1024;subplot(2,2,3);plot(f,abs(X));%输出原始语音信号的频谱图title('原始语音信号的频谱图');xlabel('频率/Hz');ylabel('幅度/db');axis([0,8000,0,10]);subplot(2,2,4);plot(f,abs(Y));%输出原始语音信号经滤波器后的频谱图title('经低通滤波器后的频谱图');xlabel('频率/Hz');ylabel('幅度/db');axis([0,8000,0,10]);%FIR带通滤波器fb1=1200;fc1=1000;fb2=3000;fc2=3200;%确定技术指标alfa=120;Fs=22050;wp1=2*pi*fc1/Fs;%求归一化频率ws1=2*pi*fb1/Fs;wp2=2*pi*fb2/Fs;ws2=2*pi*fc2/Fs;wo1=wp1-ws1;%求过渡带wo2=ws2-wp2;wo=max(wo1,wo2);beta=0.1102*(alfa-0.87);%求凯泽窗的参数N=ceil((alfa-8)/2.285/wo);%有过渡带决定Nwn1=(fb1+fc1)/Fs;wn2=(fb2+fc2)/Fs;wn=[wn1,wn2];b=fir1(N,wn,'band',kaiser(N+1,beta));[Hw,w]=freqz(b,1,512,Fs);figure(3);plot(w,20*log10(abs(Hw)));%绘制频率响应曲线grid;title('FIR带通滤波器');ylabel('幅度/db');xlabel('频率/Hz');axis([0,12000,-250,10]);x2=wavread('3.wav');%利用函数对原始信号进行采样sound(x2,22050);%回放原始语音信号x=x2(:,1);y=filter(b,1,x);%语音信号经滤波器滤波sound(y,22050);%播放滤波后的信号figure(4);subplot(2,2,1);plot(x(:,1));%输出原始语音信号的脉冲响应title('输入信号时域波形');xlabel('f/Hz');ylabel('幅度');subplot(2,2,2);plot(y(:,1));%输出原始语音信号经滤波器后的脉冲响应title('经带通滤波器后的时域波形');xlabel('f/Hz');ylabel('幅度');X=fft(x,1024);%对原始语音信号进行傅立叶变换并进行1024点采样Y=fft(y,1024);%对滤波后的信号进行傅立叶变换并进行1024点采样f=Fs*[0:1023]/1024;subplot(2,2,3);plot(f,abs(X));%输出原始信号的频谱图title('原始语音信号的频谱图');xlabel('频率/Hz');ylabel('幅度/db');axis([0,8000,0,10]);subplot(2,2,4);plot(f,abs(Y));%输出原始语音信号的经滤波后的频谱图title('经带通滤波器后的频谱图');xlabel('频率/Hz');ylabel('幅度/db');axis([0,8000,0,10]);%FIR高通滤波器fb=5000;%确定技术指标fs=4800;alfa=60;Fs=22050;wb=2*pi*fb/Fs;%求归一化频率ws=2*pi*fs/Fs;wo=wb-ws;%过渡带beta=0.1102*(alfa-0.87);%求凯泽窗函数参数N=ceil((alfa-8)/2.285/wo);%有过渡带决定Nwc=(wb+ws)/2/pi;b=fir1(N,wc,'high',kaiser(N+1,beta));[Hw,w]=freqz(b,1);figure(5);plot(w*Fs/2/pi,20*log10(abs(Hw)));grid;title('FIR高通滤波器的幅频特性');ylabel('幅度/db');xlabel('频率/Hz');axis([0,10000,-200,10]);x1=wavread('3.wav');%对原始信号进行采样sound(x1,22050);%回放原始信号x=x1(:,1);y=filter(b,1,x);%语音信号经滤波器滤波sound(y,22050);%播放滤波后的信号figure(6);subplot(2,2,1);plot(x);%输出原始语音信号的脉冲响应title('输入信号时域波形');xlabel('f/Hz');ylabel('幅度');subplot(2,2,2);plot(y);%输出原始语音信号经滤波后的脉冲响应title('高通滤波后的时域波形');xlabel('f/Hz');ylabel('幅度');X=fft(x,1024);%对原始的语音信号进行傅立叶变换并进行1024点采样Y=fft(y,1024);%对滤波后的信号进行傅立叶变换并进行1024点采样f=Fs*[0:1023]/1024;subplot(2,2,3);plot(f,abs(X));%输出原始语音信号的频谱图title('原始语音信号的频谱图');xlabel('频率/Hz');ylabel('幅度/db');axis([0,10000,0,10]);subplot(2,2,4);plot(f,abs(Y));%输出原始语音信号经滤波器后的频谱图title('经高通滤波器后的频谱图');xlabel('频率/Hz');ylabel('幅度/db');axis([0,10000,0,10]);%iir低通滤波器fb=1000;fc=1200;Fs=22050;%确定技术指标Ap=1;As=15;Ts=1/Fs;wb=2*pi*fb/Fs;%归一化，求得数字边缘频率ws=2*pi*fc/Fs;wap=2/Ts*tan(wb/2);%将数字指标转换成模拟指标was=2/Ts*tan(ws/2);[N,wc]=buttord(wap,was,Ap,As,'s');%选择滤波器的最小阶数[z,p,k]=buttap(N);%创建butterworth模拟滤波器[Bap,Aap]=zp2tf(z,p,k);[b,a]=lp2lp(Bap,Aap,wc);[bz,az]=bilinear(b,a,Fs);%用双线性变换法实现模拟滤波器到数字滤波器的转换[H,w]=freqz(bz,az);figure(7);plot(w*Fs/2/pi,20*log10(abs(H)));%绘制频率响应曲线grid;axis([0,12000,-400,100]);xlabel('频率/Hz');ylabel('幅度/db')title('巴特沃兹低通滤波器');x1=wavread('3.wav');%把语音信号加入matlab仿真平台中sound(x1,22050);x=x1(:,1);figure(8);f1=filter(bz,az,x);%语音信号经滤波器滤波sound(f1,22050);%播放滤波后的信号subplot(2,2,1);plot(x);%画出滤波前的时域图title('滤波前时域波形');xlabel('f/Hz');ylabel('幅度');subplot(2,2,2);plot(f1);%画出滤波后的时域图title('滤波后的时域波形');xlabel('f/Hz');ylabel('幅度');y=fft(f1,1024);%对滤波后的信号进行傅立叶变换并进行1024点的采样f=Fs*[0:511]/1024;%？？？？y1=fft(x1,1024);%对滤波前的信号进行傅立叶变换并进行1024点的采样subplot(2,2,3);plot(f,abs(y1(1:512)));%滤波前的频谱图title('原始语音信号滤波前的频谱图');xlabel('频率/Hz');ylabel('幅度/db');axis([0,8000,0,10]);subplot(2,2,4);F1=plot(f,abs(y(1:512)));%经巴特沃兹滤波器后的频谱图title('经巴特沃兹低通滤波器滤波后的频谱图');xlabel('频率/Hz');ylabel('幅度/db');axis([0,8000,0,10]);%iir带通滤波器fp1=1200;fp2=3000;Fs=22050;%确定技术指标fs1=1000;fs2=3200;Ap=1;As=10;wp1=2*pi*fp1/Fs;%归一化，求得数字边缘频率wp2=2*pi*fp2/Fs;ws1=2*pi*fs1/Fs;ws2=2*pi*fs2/Fs;Ts=1/Fs;wap1=2/Ts*tan(wp1/2);%将数字指标换成模拟指标wap2=2/Ts*tan(wp2/2);was1=2/Ts*tan(ws1/2);was2=2/Ts*tan(ws2/2);wap=[wap1,wap2];was=[was1,was2];W1=(wap1+was1)/2;W2=(wap2+wap2)/2;Bw=W2-W1;%带宽wo=sqrt(W1*W2);[N,wc]=buttord(wap,was,Ap,As,'s');%选择滤波器的最小阶数[z,p,k]=buttap(N);%创建巴特沃兹模拟滤波器[Bap,Aap]=zp2tf(z,p,k);[b,a]=lp2bp(Bap,Aap,wo,Bw);[bz,az]=bilinear(b,a,Fs);%用双线性变换法实现模拟滤波器到数字滤波器的转换[H,w]=freqz(bz,az);figure(9);plot(w*Fs/2/pi,20*log10(abs(H)));%绘制频率响应曲线grid;axis([0,8000,-300,50]);xlabel('频率/Hz');yla