MATLAB语音信号的处理

1、附件1：学 号： 课 程 设 计题 目基于MATLAB的基本运算和语音信号的处理 学 院信息工程学院专 业班 级 姓 名指导教师2011年11月29日课程设计任务书学生姓名： 专业班级： 指导教师： 工作单位： 题 目: 基于MATLAB的基本运算和语音信号的处理 初始条件：Matlab应用软件基本运算函数Windows自带的录音功能要求完成的主要任务: （包括课程设计工作量及其技术要求，以及说明书撰写等具体要求）1.必做题（1）极限的计算、微分的计算、积分的计算、级数的计算、求解代数方程、求解常微分方程；（2）矩阵的最大值、最小值、均值、方差、转置、逆、行列式、特征值的计算、矩阵的相乘、右除

2、、左除、幂运算（3）多项式加减乘除运算、多项式求导、求根和求值运算、多项式的部分分式展开、多项式的拟合、插值运算。2.选做题（1）完成语音信号的采集，利用windows自带的录音机或其他软件，录制一段语音，时间在1s以内，并对信号进行采样，画出采样信号的时域和频域波形。（2）要求用窗函数法和双线性变换法设计以下三种数字滤波器，其性能指标下a.低通滤波器性能指标 fb=1000hz ,fc=1200hz ,最大衰减As=1000db ,最小衰减fb=1000hz ；b.高通滤波器性能指标fs=4800hz,fb=5000hz, 最大衰减 As=100db ，最小衰减Ap=1db ； c.带通滤波

3、器性能指标 fb1=1200hz，fb2=3000hz ， fc1=1000hz, fc2=3200hz，最大衰减 As=100db ，最小衰减Ap=1db ； (3)用设计好滤波器的对语音信号进行滤波处理，画出处理后的时域和频域波形。(4)对滤波前后的语音信号频谱进行对比，并对设计结果进行分析。时间安排：第11周安排任务，分组。第12周设计仿真，撰写报告。第13周完成设计，提交报告，答辩。指导教师签名： 年 月 日系主任（或责任教师）签名： 年 月 日 目 录摘要4Abstract51 MATLAB的基本运算61.1 基本微积分运算61.1.1 极限值的计算61.1.2 微分的计算61.1.

4、3 积分的计算71.1.4级数的计算71.1.5求解代数方程81.1.6 求解常微分方程81.2 矩阵的基本运算91.2.1 矩阵最大值91.2.2 矩阵最小值91.2.3 矩阵的均值101.2.4 矩阵的方差101.2.5 矩阵的转置111.2.6 矩阵的逆矩阵111.2.7 矩阵的行列式121.2.8 矩阵特征值的计算121.2.9 矩阵的相乘131.2.10 矩阵的左除和右除141.2.11 矩阵的幂运算141.3 多项式的四则运算151.3.1 多项式加减乘除运算151.3.2 多项式求导、求根和求值运算161.3.3 多项式的部分分式展开181.3.4 多项式的拟合、插值运算182基

5、于MATLAB的语音信号处理212.1 设计原理212.1.1 采样212.1.2 时域信号的FFT分析方法222.1.3 数字滤波器222.2 设计方案与运行结果232.2.1语音信号的读取232.2.2 低通滤波器242.2.3 高通滤波器282.2.4 带通滤波器313 小结与体会36参考文献37摘要MATLAB是一种用于科学工程计算的高级语言，其基本的数据单位是矩阵，它的命令表达式与数学、工程中常用的习惯表达形式十分相似；MATLAB高效方便的矩阵和数组运算，使之在数字信号处理、建模、系统识别、自动控制、优化等领域的问题时，显得十分简洁、高交，具有其他高级语言不可比拟的优势；其典型应用

6、包括以下几个方面：数学计算、算法设计、建模与仿真、数据分析与可视化、科学及工程绘图和应用开发；语音信号处理是研究用数字信号处理技术和语音学知识对语音信号进行处理的新兴的学科，是目前发展最为迅速的信息科学研究领域的核心技术之一。通过语音传递信息是人类最重要、最有效、最常用和最方便的交换信息形式。本设计针对现在大部分语音处理软件内容繁多、操作不便等问题，采用MATLAB7.0综合运用GUI界面设计、各种函数调用等来实现语音信号的变频、变幅、傅里叶变换及滤波，程序界面简练，操作简便，具有一定的实际应用意义。关键字：Matlab；基本运算；滤波；语音信号处理；AbstractMATLAB is a u

7、sed in scientific engineering calculation languages, and its basic data unit is matrix, it commands expression and mathematics, engineering in the habit of commonly used expression is very similar to form; MATLAB effective convenient matrix and array of operation, make in digital signal processing,

8、modeling, system identification, automatic control and optimization, etc, is very concise, high pay, with other senior language incomparable advantages; Its typical applications include the following aspects: mathematical calculation, algorithm design and modeling and simulation, data analysis and v

9、isualization, science and engineering drawing and application development;The speech signal processing is using digital signal processing technology and knowledge of phonetics the speech signal processing new discipline, is the most rapid development of the information science research in the field

10、of core technical one. Through the speech information is the most important and the most effective human, the most common and most convenient exchange information form.This design for now most speech processing various content, operating software inconvenience of, use MATLAB7.0 comprehensive use of

11、GUI interface design, all sorts of function calls, to realize the speech signal frequency conversion, derricking and Fourier transform and filter, and the program interface concise, simple operation, and has certain actual application significance.1 MATLAB的基本运算1.1 基本微积分运算1.1.1 极限值的计算MATLAB中极限函数limit

12、格式如下：limit(F,x,a) ；limit(F,a)；limit(F)；limit(F,x,a,'right')；limit(F,x,a,'left')其中F表示的是函数式，x表示求极限的变量，a表示的是变量取的值，left或right表示是取左极限还是右极限。例：求limx2，y0 x3+cosy+100示例程序如下：syms x y limit(limit(x3+cos(y)+100,x,2),y,0)计算结果： 1.1.2 微分的计算MATLAB中微分函数diff格式如下：Y = diff(X)；Y = diff(X,n)；Y = diff(X,n,

13、dim)其中X表示待微分的变量，n表示n次微分，第三式表示沿着定维dim的n阶微分。例：求y=12cos(3x-4)+x3的微分;示例程序如下：syms x y dyy=12*cos(3*x-4)+x3;dy=diff(y)计算结果：1.1.3 积分的计算MATLAB中积分常用函数为int，格式如下：int(f,x,a,b)其中f表示待积分的函数，x表示积分变量，而a，b则分别表示积分起始终止点。例：求下列积分-115x3+1dx示例程序如下：syms xint(5*x3+1,-1,1)计算结果：1.1.4级数的计算MATLAB中级数常用函数为symsum，格式如下：r = symsum(s)

14、；r = symsum(s,v)；r = symsum(s,a,b)；r = symsum(s,v,a,b)函数表达的意义是表达式s关于变量v从a到b求和。例：      计算级数 11n2 的值,分析：对于级数我们最先接触的就是常数项无穷级数，因为技术结果一般较复杂，所以用simple求解最简形式,还有Inf为无穷大。示例程序如下：syms xsimple(symsum(1/x2,1,Inf)计算结果：1.1.5求解代数方程代数方程分为很多种，有简单有复杂，方法也有很多。一般多项式方程的根可为实数，也可为复数，可用MATLAB符号工具箱中的s

15、olve( )函数， MATLAB中求解代数方程常用函数solve格式如下：x,=solve(eqn1,eqn2,eqnn，x,)其中eqn表示的是式子，x等表示的是变量。例：求解方程3x2+17x-6=0 示例程序如下：syms xx=solve('3*x2+17*x-6=0','x')计算结果：1.1.6 求解常微分方程未知的函数以及它的某些阶的导数连同自变量都由一已知方程联系在一起的方程称为微分方程。如果未知函数是一元函数，称为常微分方程。MATLAB中主要用dsolve求符号解析解。dsolve的基本格式是：s=dsolve(方程1, 方程2,初始条件1

16、，初始条件2 ,自变量)例：求解微分方程y+3y=0,y(0)=2求解析解。示例程序如下:syms x yr=dsolve('Dy+3*y=0','y(0)=2','x')计算结果：1.2 矩阵的基本运算1.2.1 矩阵最大值MATLAB中max函数可以表示求每一列的最大值，那么经过分析可以知道，先求出每一列的最大值然后求出这些最大值里面的最大值，下面以A矩阵为例。示例程序如下：A=2 -9 1;3 6 5;7 2 5y=max(A)X=max(y)运行结果如下： 1.2.2 矩阵最小值与求最大值类似，MATLAB中min函数可以表示求每一列的最

17、小值，那么经过分析可以知道，先求出每一列的最小值然后求出这些最小值里面的最小值，下面以A矩阵为例A=2 -9 1;3 6 5;7 2 5y=min(A)X=min(y)运行结果如下：1.2.3 矩阵的均值MATLAB中求解矩阵均值的函数是mean，它的具体用法如下：mean(A,1)表示对列取平均，mean(A,2）表示对行取平均，mean(A）则默认为mean(A,1)。下面以矩阵A分别举例，程序示例如下：A=2 -9 1;3 6 5;7 2 5y=mean(A)x=mean(A,2)运行结果：1.2.4 矩阵的方差MATLAB中求解矩阵方差的函数是var，它的常用格式是V = var(X)

18、，如果X是一个矩阵，var(X)返回一个包含矩阵X每一列方差的行向量。下面还是以矩阵A来示例，程序如下：A=2 -9 1;3 6 5;7 2 5d=var(var(A)计算结果：1.2.5 矩阵的转置矩阵的一个重要的运算是转置，如果A是一个实数矩阵，那么它被转置时，第1行变成第1列，第2行变成第2列，依此类推，一个m×n矩阵变为一个n×m矩阵。如果矩阵是方阵，那么这个矩阵在主对角线反映出来。MATLAB中求转置的函数是A，以A为例，编程如下：A=2 -9 1;3 6 5;7 2 5a=A'计算结果：1.2.6 矩阵的逆矩阵实际中求矩阵的逆跟行列式均要求矩阵是方阵，M

19、ATLAB中求逆的函数是inv，格式为Y = inv(X)。下面仍以矩阵A 为例来编程示例，如下：A=2 -9 1;3 6 5;7 2 5y=inv(A)运行结果如下：1.2.7 矩阵的行列式在MATLAB中det用于求行列式。以矩阵A 为例，编程如下A=2 -9 1;3 6 5;7 2 5y=det(A)运行结果如下：1.2.8 矩阵特征值的计算矩阵的特征值的求解，就是找到方程组的解：其中是一个标量，x是一个长度为n的列向量。标量是A的特征值，x是相对应的特征向量。对于实数矩阵A来说，特征值和特征向量可能是复数。一个 n×n的矩阵有n个特征值，表示为。求矩阵的特征值和特征向量可用e

20、ig函数。Eig(A)求包含矩阵A的特征值的向量。V,D =eig(A)产生一个矩阵A的特征值在对角线上的对角矩阵D和矩阵V，它的列是相应的特征向量，满足AV=VD，下面以矩阵A为例来演示。编程如下：A=2 -9 1;3 6 5;7 2 5V,D=eig(A)计算结果：1.2.9 矩阵的相乘 假定有两个矩阵A和B，若A为m×n矩阵，B为n×p矩阵，则C=A×B为m×p矩阵。元素是A的第i行和B的第j列的点积。对于方阵，也定义了积B A，但其结果通常与A B不同。MATLAB中求矩阵的乘积直接用符号*即可，下面以A、B矩阵为例来分别演示A*B与B*A区别。

21、示例程序如下：A=2 -9 1;3 6 5;7 2 5B=1 2 5;9 3 5;6 8 2Y=A*BX=B*A计算结果：1.2.10 矩阵的左除和右除在MATLAB中，有两个矩阵除法的符号，左除“ ”和右除“/”。如果A是一个非奇异方阵，那么A B和B / A对应A的逆与B的左乘和右乘，即分别等价于命令 inv(A)*B和B*inv(A)。可是，MATLAB执行它们时是不同的，且在MATLAB中求解一个系统用左除比用逆和乘法所需的运算次数要少。令X=AB, Y=B/A ， 下面仍然以A、B为例来演示。示例程序如下：A=2 -9 1;3 6 5;7 2 5B=1 2 5;9 3 5;6 8 2

22、X=ABY=B/A计算结果：1.2.11 矩阵的幂运算对于二维方阵，A的p次乘方可以用Ap实现。如果p是一个正整数，那么这个幂可以由许多矩阵乘法运算定义。对于 p= 0，得到与A维数相同的同一个矩阵；当 p< 0时，如果A- 1存在，可定义A p，它是与inv(A)(-p)相同。A0=A3，A1=A3，A2=A-3Ap0为3个A矩阵相乘，Ap1中的元素为A矩阵中相应元素的立方，矩阵Ap2为矩阵A的逆矩阵的乘积，A3为A0的逆矩阵。以矩阵A为例，分别编程实例如下：A=2 -9 1;3 6 5;7 2 5B=1 2 5;9 3 5;6 8 2A0=A3A1=A.3A2=A-3A3=A0-1计

23、算结果：1.3 多项式的四则运算1.3.1 多项式加减乘除运算多项式的四则运算就是包括加减乘除，其中加减运算可以直接用+、-来运算，它们的运算规则中注意要满足向量的长度相同，而乘除就得用函数了，其中乘法的计算函数是conv，它本来是卷积的意思，同时它也符合多项式函数的运算规则，除法运算是相乘的逆运算，但会有余子式。下面以A,B为例来进行四则运算。示例编程如下：a=5,1,8,18;b=4,1,32,2;x=a+by=a-bz=conv(a,b)q,r=deconv(a,b)计算结果：1.3.2 多项式求导、求根和求值运算多项式的求导、求根、求值运算是多项式运算的又一大板块，其中多项式求导数的函

24、数是polyder，调用格式是e=polyder（c），其中c表示的是待求导的函数式，然后求根运算的函数是roots或poly，其中roots是根据函数求多项式的根，它的调用格式是h=roots（c），c代表待求根的函数式，而poly函数是根据根求函数，格式是i=poly（h），表示根据h求函数i，然后求值运算的函数是polyval，将多项式的自变量赋予值z，则调用格式是j=polyval（f，z），表示当变量是1时，函数f的结果。下面编程演示。示例编程如下：a=2 9 16 5x=polyder(a) %求函数a的导数y=roots(a) %求函数a的根z=poly(x) %根据所求的根求函

25、数w=polyval(a,4) %当变量值为4时，函数a的值 计算结果：1.3.3 多项式的部分分式展开函数residue可以将多项式之比用部分分式展开，也可以将一个部分分式表示为多项式之比。其调用格式如下：r,p,k=residue(a,b)返回多项式之比a/b的部分分式展开，参照下面公式。a,b=residue(r,p,k)返回部分分式的多项式向量。示例程序如下：a=2 -4 10 5;b=6 1 12 9;r,p,k=residue(a,b) %返回多项式之比a/b的部分分式展开c,d=residue(r,p,k) %f返回部分分式的多项式向量计算结果：1.3.4 多项式的拟合、插值运算

26、在MATLAB中指令polyfit用于求数据的拟合多项式。可以通过polyval指令给求得的多项式赋值，数值化后就可以通过数值画出图形。这种方法可以用于对一些数据的预测。调用函数polyfit常用格式为p = polyfit(x,y,n)，x为变量，y为函数，n为阶数。示例程序如下：x=linspace(0,4*pi,100); %定义向量，从0到4派分为100份y=5*sin(x+2); %定义函数t=polyfit(x,y,7); %7次拟合y1=polyval(t,x); %根据拟合的结果求多项式的值subplot(2,1,1), plot(x,y,'ro',x,y1,&

27、#39;g-') %作图s=polyfit(x,y,4); %4次拟合y2=polyval(s,x); %根据拟合的结果求多项式的值subplot(2,1,2),plot(x,y,'ro',x,y2,'g-') %作图计算结果： 由图可以知道高次拟合的效果要好一些插值函数通常是分段的，插值数据通过给定的数据点x，y。插值函数一般地可表示为yi=interpi（x,y,xi,method）其中i代表几维插值可取1、2,xi为插值范围内的任意点集的x坐标,yi是插值后对应数据点集的坐标,method为插值函数的类型选项，有linear为线性，也是缺省项，cu

28、bic和cubic spline为三次样条等三样。例:已知某矩形温箱中3×5个测试点上的温度，求全箱的温度分布。给定：wid=1：5；dep=1：3；tem=60 61 66 62 64;79 62 74 78 66;70 62 69 77 71要求计算沿宽度和深度细分网格：di=1：0.2：3；wi=1：0.2：5；交点温度。示例程序如下：wid=1:5;dep=1:3;tem=60 61 66 62 64;79 62 74 78 66;70 62 69 77 71;di=1:0.2:3;wi=1:0.2:5;tc=interp2(wid,dep,tem,wi,di',&#

29、39;cubic'); %求温度,di必须变成列向量mesh(wi,di,tc) %画三维图计算结果:2基于MATLAB的语音信号处理2.1 设计原理语音信号是基于时间轴上的一维数字信号，它带着各种信息。在语音编码、语音合成、语音识别和语音增强等语音处理中无一例外需要提取语音中的各种信息。语音分析的目的在于方便有效地提取并表示语音信号所携带的信息。本课题就从时域和频域两个角度对信号进行分析，并通过分析频谱来设计出合适的滤波器。2.1.1 采样在进行模拟/数字信号的转换过程中，当采样频率fs.max大于信号中，最高频率fmax的2倍时，即：fs.max>=2fmax,则采样之后的数

30、字信号完整地保留了原始信号中的信息，一般实际应用中保证采样频率为信号最高频率的510倍；采样定理又称奈奎斯特定理。 1924年奈奎斯特(Nyquist)就推导出在理想低通信道的最高大码元传输速率的公式:理想低通信道的最高大码元传输速率=2W*log2 N (其中W是理想低通信道的带宽,N是电平强度)采样频率是指计算机每秒钟采集多少个声音样本，是描述声音文件的音质、音调，衡量声卡、声音文件的质量标准。采样频率越高，即采样的间隔时间越短，则在单位时间内计算机得到的声音样本数据就越多，对声音波形的表示也越精确。采样频率与声音频率之间有一定的关系，根据奎斯特理论，只有采样频率高于声音信号最高频率的两倍

31、时，才能把数字信号表示的声音还原成为原来的声音。这就是说采样频率是衡量声卡采集、记录和还原声音文件的质量标准。采样位数即采样值或取样值，用来衡量声音波动变化的参数，是指声卡在采集和播放声音文件时所使用数字声音信号的二进制位数。采样频率是指录音设备在一秒钟内对声音信号的采样次数，采样频率越高声音的还原就越真实越自然。采样位数和采样率对于音频接口来说是最为重要的两个指标，也是选择音频接口的两个重要标准。无论采样频率如何，理论上来说采样的位数决定了音频数据最大的力度范围。每增加一个采样位数相当于力度范围增加了6dB。采样位数越多则捕捉到的信号越精确。对于采样率来说你可以想象它类似于一个照相机，44.

32、1kHz意味着音频流进入计算机时计算机每秒会对其拍照达441000次。显然采样率越高，计算机摄取的图片越多，对于原始音频的还原也越加精确。2.1.2 时域信号的FFT分析方法FFT即为快速傅氏变换，是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。在MATLAB的信号处理工具箱中函数FFT和IFFT用于快速傅立叶变换和逆变换。函数FFT用于序列快速傅立叶变换，其调用格式为y=fft(x)，其中，x是序列，y是序列的FFT，x可以为一向量或矩阵，若x为一向量，y是x的FFT且和x相同长度；若x为一矩阵，则y是对矩阵的每一列向量进行FFT。

33、如果x长度是2的幂次方，函数fft执行高速基2FFT算法，否则fft执行一种混合基的离散傅立叶变换算法，计算速度较慢。函数FFT的另一种调用格式为y=fft(x,N)，式中，x，y意义同前，N为正整数。函数执行N点的FFT，若x为向量且长度小于N，则函数将x补零至长度N；若向量x的长度大于N，则函数截短x使之长度为N；若x 为矩阵，按相同方法对x进行处理。2.1.3 数字滤波器FIR:Finite Impulse response，有限冲击响应IIR:Infinite Impulse response，无限冲击响应数字滤波器的作用是利用离散时间系统的特性对输入信号波形(或频谱)进行加工处理，或

34、者说利用数字方法按预定的要求对信号进行变换。数字滤波器可以理解为是一个计算程序或算法，将代表输入信号的数字时间序列转化为代表输出信号的数字时间序列，并在转化过程中，使信号按预定的形式变化。数字滤波器有多种分类，根据数字滤波器冲激响应的时域特征，可将数字滤波器分为两种，即无限长冲激响应（IIR）滤波器和有限长冲激响应（FIR）滤波器。从性能上来说，IIR滤波器传输函数的极点可位于单位圆内的任何地方，因此可用较低的阶数获得高的选择性，所用的存贮单元少，所以经济而效率高。但是这个高效率是以相位的非线性为代价的。选择性越好，则相位非线性越严重。相反，FIR滤波器却可以得到严格的线性相位，然而由于FIR

35、滤波器传输函数的极点固定在原点，所以只能用较高的阶数达到高的选择性；对于同样的滤波器设计指标，FIR滤波器所要求的阶数可以比IIR滤波器高510倍，结果，成本较高，信号延时也较大；如果按相同的选择性和相同的线性要求来说，则IIR滤波器就必须加全通网络进行相位较正，同样要大增加滤波器的节数和复杂性。整体来看，IIR滤波器达到同样效果阶数少，延迟小，但是有稳定性问题，非线性相位；FIR滤波器没有稳定性问题，线性相位，但阶数多，延迟大。2.2 设计方案与运行结果2.2.1语音信号的读取下面的一段程序是语音信号在MATLAB中的最简单表现，它实现了语音的读入打开，以及绘出了语音信号的波形时域图和频谱图

36、。x,fs,bits=wavread('look.wav');sound(x,fs,bits);y=fft(x,1024);f=fs*(0:511)/1024;t=(0:length(x)-1)/fs;figure(1);plot(t,x);figure(2);plot(f,20*log10(abs(y(1:512);图 1 时域波形图 2 频域波形2.2.2 低通滤波器1 .窗函数法y,fs,bit=wavread('look.wav');sound(y,fs,bit);f=fs*(0:511)/1024;t=(0:length(y)-1)/fs;fr=120

37、0;wr=2*pi*fr;Window=boxcar(8);b=fir1(7,fr/fs*2,Window);h,w=freqz(b,1);y1=filter(b,1,y);subplot(221);plot(t,y);title('原始语音信号波形');xlabel('t/s');x1=fft(y,1024);subplot(222);plot(f,20*log10(abs(x1(1:512);title('原始语音信号频谱');xlabel('f/Hz');ylabel('dB');subplot(223);p

38、lot(t,y1);title('滤波后语音信号波形');xlabel('t/s');y2=fft(y1,1024);subplot(224);plot(f,20*log10(abs(y2(1:512);title('滤波后语音信号频谱');xlabel('f/Hz');ylabel('dB');sound(y1,fs,bit);滤波处理后的时域和频域波形：图 3 窗函数低通滤波结果比较上图的原语音频谱图和滤波处理后的语音频谱图，频率较高的语音信号被明显地衰减，低频的语音信号保留了下来，起到了低通滤波的效果。2 双

39、线性变换法y,fs,bit=wavread('look.wav');sound(y,fs,bit);x1=fft(y,1024);fp=1000;fr=1200;f=fs*(0:511)/1024;t=(0:length(y)-1)/fs;subplot(221);plot(t,y);title('原始语音信号波形');xlabel('t/s');subplot(222);plot(f,20*log10(abs(x1(1:512);title('原始语音信号频谱');xlabel('f/Hz');ylabel(&#

40、39;dB');wp=(fp/fs)*2*pi;ws=(fr/fs)*2*pi;OmegaP=2*fs*tan(wp/2);OmegaS=2*fs*tan(ws/2);ap=1;as=15;N,Wn=buttord(OmegaP,OmegaS,ap,as,'s'); %选择滤波器的最小阶数b,a=butter(N,Wn,'s');bz,az=bilinear(b,a,fs);h,w=freqz(bz,az);y1=filter(bz,az,y);subplot(223);plot(t,y1);title('滤波后语音信号波形');xlab

41、el('t/s');y2=fft(y1,1024);subplot(224);plot(f,20*log10(abs(y2(1:512);title('滤波后语音信号频谱');xlabel('f/Hz');ylabel('dB');sound(y1,fs,bit);滤波处理后的时域和频域波形：图 4 双极性变换法低通滤波结果比较上图的原语音频谱图和滤波处理后的语音频谱图，频率较高的语音信号被明显地衰减，低频的语音信号保留了下来，起到了低通滤波的效果。2.2.3 高通滤波器1 窗函数法y,fs,bit=wavread('lo

42、ok.wav');sound(y,fs,bit);f=fs*(0:511)/1024;t=(0:length(y)-1)/fs;subplot(221);plot(t,y);title('原始语音信号波形');xlabel('t/s');x1=fft(y,1024);subplot(222);plot(f,20*log10(abs(x1(1:512);title('原始语音信号频谱');xlabel('f/Hz');ylabel('dB');fr=5000;wr=2*pi*fr;Window=boxcar(

43、8);b=fir1(7,fr/(fs/2),'high');h,w=freqz(b,1);y1=fftfilt(b,y);subplot(223);plot(t,y1);title('滤波后语音信号波形');xlabel('t/s');y2=fft(y1,1024);subplot(224);plot(f,20*log10(abs(y2(1:512);title('滤波后语音信号频谱');xlabel('f/Hz');ylabel('dB');sound(y1,fs,bit);滤波处理后的时域和频域

44、波形：图 5 窗函数法高通滤波结果比较上图的原语音频谱图和滤波处理后的语音频谱图，频率较低的语音信号被明显地衰减，高频的语音信号保留了下来，起到了高通滤波的效果。2 双线性变换法y,fs,bit=wavread('look.wav');sound(y,fs,bit);f=fs*(0:511)/1024;t=(0:length(y)-1)/fs;subplot(221);plot(t,y);title('原始语音信号波形');xlabel('t/s');x1=fft(y,1024);subplot(222);plot(f,20*log10(abs(

45、x1(1:512);title('原始语音信号频谱');xlabel('f/Hz');ylabel('dB');fp=5000;fr=4800;wp=(fp/fs)*2*pi;ws=(fr/fs)*2*pi;OmegaP=2*fs*tan(wp/2);OmegaS=2*fs*tan(ws/2);ap=1;as=20;N,Wn=buttord(OmegaP,OmegaS,ap,as,'s'); %选择滤波器的最小阶数b,a=butter(N,Wn,'high','s');bz,az=bilinear(

46、b,a,fs);h,w=freqz(bz,az);y1=filter(bz,az,y);subplot(223);plot(t,y1);title('滤波后语音信号波形');xlabel('t/s');y2=fft(y1,1024);subplot(224);plot(f,20*log10(abs(y2(1:512);title('滤波后语音信号频谱');xlabel('f/Hz');ylabel('dB');sound(y1,fs,bit); 滤波处理后的时域和频域波形：比较上图的原语音频谱图和滤波处理后的语音频

47、谱图，频率较低的语音信号被明显地衰减，高频的语音信号保留了下来，起到了高通滤波的效果。2.2.4 带通滤波器1 窗函数法y,fs,bit=wavread('look.wav');sound(y,fs,bit);f=fs*(0:511)/1024;t=(0:length(y)-1)/fs;subplot(221);plot(t,y);title('原始语音信号波形');xlabel('t/s');x1=fft(y,1024);subplot(222);plot(f,20*log10(abs(x1(1:512);title('原始语音信号频谱

48、');xlabel('f/Hz');ylabel('dB');fr=1200;fp=3000;wr=2*pi*fr;Window=blackman(16);b=fir1(15,fr/(fs/2),fp/(fs/2),Window);h,w=freqz(b,1);y1=filter(b,1,y);subplot(223);plot(t,y1);title('滤波后语音信号波形');xlabel('t/s');y2=fft(y1,1024);subplot(224);plot(f,20*log10(abs(y2(1:512);title('滤波后语音信号频谱');xlabel('f/Hz');ylabel('