MATLAB课程设计---基于MATLAB的语音信号处理.doc

武汉理工大学基于matlab的语音信号处理课程设计说明书课程设计任务书学生姓名：陈欢 专业班级： 通信0902班 指导教师：徐文君 工作单位： 信息工程学院 题 目: 基于matlab的语音信号处理设计任务与要求:1、matlab的基本运算（1） 极限的计算、微分的计算、积分的计算、级数的计算、求解代数方程、求解常微分方程；（2） 矩阵的最大值、最小值、均值、方差、转置、逆、行列式、特征值的计算、矩阵的相乘、右除、左除、幂运算；2、基于matlab的语音信号处理（1）完成语音信号的采集，利用windows自带的录音机或其他软件，录制一段午语音，时间在1s以内，并对信号进行采样，画出采样信号的时域和频域波、（2）要求用窗函数法和双线性变换法设计三种数字滤波器 （3）用设计好滤波器的对语音信号进行滤波处理，画出处理后的时域和频域波形。（4）对滤波前后的语音信号频谱进行对比，并对设计结果进行分析。参考文献：1、易克初，田斌，付强.语音信号处理.北京：国防工业出版社，2000年2、贾永红.数字图像处理.武汉大学出版社，2003年9月时间安排：第16周安排任务，分组。第17周设计仿真，撰写报告。第18周完成设计，提交报告，答辩。指导教师签名： 2011 年 月 日系主任（或责任教师）签名： 2011 年 月 日目录摘要1abstract21 课程设计的基本原理31.1 设计要求31.1.1 采用matlab基本运算如下31.1.2 基于matlab的语音信号处理31.2.matlab的基本运算41.2.1极限的计算41.2.2微分的计算41.2.3积分的计算51.2.4 级数的计算51.2.5 求解代数方程51.2.6 求解常微分方程61.2.7 矩阵的基本运算61.2.8 多项式的基本运算61.3 语音信号的录入与打开81.4 时域信号的fft分析81.5 数字滤波器设计原理81.5.1 用窗函数法设计fir滤波器91.5.2 用双线性变换法设计iir数字滤波器92 语音信号分析和处理过程102.1 语音信号的采集102.2 窗函数法设计112.2.1窗函数法低通滤波器112.2.2 窗函数法高通滤波器142.2.3 窗函数法带通滤波器152.3双极性变换法设计162.3.1 双极性变换法设计低通滤波器182.3.2 双极性变换法设计高通滤波器202.3.3 双极性变换法设计带通滤波器222.4 结果分析243 总结25参考文献26摘要matlab是矩阵实验室（matrix laboratory）的简称，和mathematica、maple并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。matlab可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。matlab的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用matlab来解算问题要比用c，fortran等语言完相同的事情简捷得多，并且mathwork也吸收了像maple等软件的优点,使matlab成为一个强大的数学软件。在新的版本中也加入了对c，fortran，c+ ，java的支持。可以直接调用,用户也可以将自己编写的实用程序导入到matlab函数库中方便自己以后调用，此外许多的matlab爱好者都编写了一些经典的程序，用户可以直接进行下载就可以用。关键字：matlab、矩阵运算、图像处理abstractmatlab is matrix lab (matrix laboratory) abbreviation, and mathematica, maple and called three mathematical software. it in math class technology applications in the numerical calculation is second to none. matlab matrix operations can be drawn function and data, and realize the algorithm, create the user interface, the link to the rest of the programming language procedure, etc. are mainly used in the engineering, control design, signal processing and communication, image processing and signal detection, financial modeling design and analysis, etc.matlab unit is the basic data of the matrix, its instruction expression and mathematics, engineering in the form of common is very similar, so matlab to the solution to problems than with c, fortran language such as the same things simple, and mathwork also absorbed like maple software such as the advantages of the matlab to become a powerful mathematical software. in the new version joining in the c, fortran, c + +, java support. can directly calls, users can take the practical program written into matlab function library to their future calls, moreover many of matlab lovers are writing some classic program, users can directly download can use.key word: matlab, matrix computation, image processing1 课程设计的基本原理1.1 设计要求1.1.1 采用matlab基本运算如下（3） 极限的计算、微分的计算、积分的计算、级数的计算、求解代数方程、求解常微分方程；（4） 矩阵的最大值、最小值、均值、方差、转置、逆、行列式、特征值的计算、矩阵的相乘、右除、左除、幂运算；（5） 多项式加减乘除运算、多项式求导、求根和求值运算、多项式的部分分式展开、多项式的拟合、插值运算。1.1.2 基于matlab的语音信号处理（1）.完成语音信号的采集，利用windows自带的录音机或其他软件，录制一段午语音，时间在1s以内，并对信号进行采样，画出采样信号的时域和频域波形。（2）.要求用窗函数法和双线性变换法设计以下三种数字滤波器，其性能指标如下a.低通滤波器性能指标 fb=1000hz ,fc=1200hz ,最大衰减as=1000db ,最小衰减fb=1000hz ；b.高通滤波器性能指标fs=4800hz,fb=5000hz, 最大衰减 as=100db ，最小衰减ap=1db ； c.带通滤波器性能指标 fb1=1200hz，fb2=3000hz ， fc1=1000hz, fc2=3200hz，最大衰减 as=100db ，最小衰减ap=1db ； (3)用设计好滤波器的对语音信号进行滤波处理，画出处理后的时域和频域波形。(4).对滤波前后的语音信号频谱进行对比，并对设计结果进行分析。1.2.matlab的基本运算1.2.1极限的计算matlab提供的命令函数limit()可以完成极限运算，其调用格式如下：limit(f,x,a,left)该命令对表达式f求极限，独立变量x从左边趋于a，函数中除f外的参数均可省略，left可换成right。举例如下：f=sym(1+a/x)x)f =(a/x + 1)xlimit(f,x,inf,left)ans =exp(a)1.2.2微分的计算matlab提供的函数diff()可以完成对给定函数求导函数的运算，其调用格式如下：diff(fun,x,n)其意义是求函数fun关于变量x的n阶导数，n为1时可省略。clear syms x y=log(x+2)/(1-x); dy=diff(y,x)dy =(1/(x - 1) - (x + 2)/(x - 1)2)*(x - 1)/(x + 2) dy3=diff(y,x,3)dy3 =(2*(1/(x - 1) - (x + 2)/(x - 1)2)*(x - 1)/(x + 2)3 - (2*(2/(x - 1)2 - (2*(x + 2)/(x - 1)3)/(x + 2) - (2*(1/(x - 1) - (x + 2)/(x - 1)2)/(x + 2)2 + (2*(2/(x - 1)2 - (2*(x + 2)/(x - 1)3)*(x - 1)/(x + 2)2 + (6/(x - 1)3 - (6*(x + 2)/(x - 1)4)*(x - 1)/(x + 2) pretty(dy3) %其功能是使它作用的表达式更符合数学上的书写习惯1.2.3积分的计算int(f)函数f对符号变量x或接近字母x的符号变量求不定积分； int(f,t)函数f对符号变量t 求不定积分； int(f,a,b)函数f 对符号变量x 或接近字母x的符号变量求从a到b的定积分； int(f,t,a,b)函数f 对符号变量t 求从a 到b 的定积分。 syms a x f=sin(a*x) g=int(f,0,pi) f =sin(a*x)g =(2*sin(pi*a)/2)2)/a1.2.4 级数的计算matlab中级数常用函数为symsum，格式如下：r = symsum(s)；r = symsum(s,v)；r = symsum(s,a,b)；r = symsum(s,v,a,b)函数表达的意义是表达式s关于变量v从a到b求和。taylor(f,n,a)函数f 对符号变量x （或最接近字母x 的符号变量）在a 点的n-1阶泰勒多项式（n 缺省时值为6，a 缺省值为0） taylor(sin(x) ans =x5/120 - x3/6 + x1.2.5 求解代数方程在matlab中，我们使用solve(f,t)命令解代数方程，所要做的事是把方程用的单引号引起来然后敲回车。其中solve的第二中调用方式，我们可以告诉它要解那个符号，语法如下：solve(equation,variable)g1=x+y+z-1;g2=x-y+z-2;g3=2*x-y-z-1;x,y,z=solve(g1,g2,g3)x =2/3 y =-1/2 z =5/61.2.6 求解常微分方程dsolve(eq,s1,s2,.,x),其中eq为方程；s1,s2,为初始条件缺省时给出含任意常数c1,c2,的通解；x为自变量，缺省时默认为t。y=dsolve(dy=-2*y+2*x2+2*x,y(0)=1,x)y =1/exp(2*x) + x21.2.7 矩阵的基本运算矩阵计算是线性代数中的核心内容，其对于整个数学系统的计算方面的意义是十分巨大的，集中它的基本运算包括最大值、最小值、均值、方差、转置、逆、行列式、特征值的计算、矩阵的相乘、右除、左除、幂运算等等，下面将具体介绍。矩阵的运算都是要以矩阵为基础的，本报告中决定选用一组矩阵来完成几乎全部可以完成的计算，那么首先就得生成矩阵了。矩阵的定义和分配可以 有多种方法。最简单的方法是有方括号包围的逐行给定元素。若定义一个标量，则方括号就不需要了。相同行中的元素是由一行或多个空格或一个逗号，分隔，列由分号; 或回车键分隔。没有结尾分号的每个命令在屏幕上显示出其结果。若结尾带分号，就执行计算，但计算结果并不显示。函数矩阵的最大值、最小值、均值、方差、逆、行列式、特征值的计算是通过函数max（）、min（）、mean（）、var（）、inv（）、det（）、eig（）来实现的；矩阵的相乘、右除、左除、幂运算可以直接通过相关运算符来实现.1.2.8 多项式的基本运算多项式的运算，主要包括多项式加减乘除、多项式求导、求根和求值运算、多项式的部分分式展开、多项式的拟合、插值运算。1.2.8.1 多项式的四则运算多项式求导、求根和求值运算是通过函数polyder（）、roots（）、polyval（）来实现的；1.2.8.2多项式的部分分式展开及拟合函数residue可以将多项式之比用部分分式展开，也可以将一个部分分式表示为多项式之比。其调用格式如下：r,p,k=residue(a,b)返回多项式之比a/b的部分分式展开，参照下面公式。a,b=residue(r,p,k)返回部分分式的多项式向量。通过函数residue（）、polyfit（）实现1.2.8.3 多项式的插值拟合插值函数通常是分段的，插值数据通过给定的数据点x，y。插值函数一般地可表示为yi=interpi（x,y,xi,method）其中i代表几维插值可取1、2,xi为插值范围内的任意点集的x坐标,yi是插值后对应数据点集的坐标,method为插值函数的类型选项，有linear为线性，也是缺省项，cubic和cubic spline为三次样条等三样。一维多项式插值：yi=interp1(x,y,xi,method);x必须是向量，y可是是向量也可以是矩阵。如果y是向量，则必须与x具有相同的长度，这时xi可以是标量、向量和任意维矩阵，yi与xi具有相同的大小；如果y是矩阵，则其大小必须是n,d1,d2,.,dk(n是向量x的长度)，函数对d1*d2*d3*dk组y值都进行插值。yi=interp1(y,xi);默认x为1:n，其中n是向量y的长度。yi=interp1(x,y,xi,method);输入变量method用于指定插值方法yi=interp1(x,y,xi,method,extrap);对超出插值范围的数据指定外推方法extrapyi=interp1(x,y,xi,method,extrapval);对超出差值范围的数值返回extrapval的值，一般为0或者nan(not a number)pp=interp1(x,y,method,pp);返回值pp为数据y的分段多项式形式。method指定产生多项式的方法。1.3 语音信号的录入与打开在matlab中，y,fs,bits=wavread(blip,n1 n2);用于读取语音，采样值放在向量y中，fs表示采样频率(hz)，bits表示采样位数。n1 n2表示读取从n1点到n2点的值（若只有一个n的点则表示读取前n点的采样值）。 sound(x,fs,bits); 用于对声音的回放。向量y则就代表了一个信号（也即一个复杂的“函数表达式”）也就是说可以像处理一个信号表达式一样处理这个声音信号。1.4 时域信号的fft分析在matlab的信号处理工具箱中函数fft和ifft用于快速傅立叶变换和逆变换。函数fft用于序列快速傅立叶变换，其调用格式为y=fft(x)，其中，x是序列，y是序列的fft，x可以为一向量或矩阵，若x为一向量，y是x的fft且和x相同长度；若x为一矩阵，则y是对矩阵的每一列向量进行fft。如果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进行处理。1.5 数字滤波器设计原理数字滤波是数字信号分析中最重要的组成部分之一，与模拟滤波相比，它具有精度和稳定性高、系统函数容易改变、灵活性强、便于大规模集成和可实现多维滤波等优点。在信号的过滤、检测和参数的估计等方面，经典数字滤波器是使用最广泛的一种线性系统。数字滤波器的作用是利用离散时间系统的特性对输入信号波形(或频谱)进行加工处理，或者说利用数字方法按预定的要求对信号进行变换。1.5.1 用窗函数法设计fir滤波器如果所希望的滤波器的理想频率响应函数为 hd(e j)，则其对应的单位脉冲响应为： （1）用窗函数w(n)将hd(n)截断，并进行加权处理，得到： （2）h(n)就作为实际设计的fir数字滤波器的单位脉冲响应序列，其频率响应函数h(e j)为： （3）如果要求线性相位特性，则h(n)还必须满足： （4）根据上式中的正、 负号和长度n的奇偶性又将线性相位fir滤波器分成四类。 要根据所设计的滤波特性正确选择其中一类。例如，要设计线性相位低通特性， 可选择h(n)=h(n-1-n)一类，而不能选h(n)=-h(n-1-n)一类。然后用freqz函数计算数字滤波器在频域中的特性，如果不满足要求，可根据具体情况调整窗函数类型或长度，知道满足要求。1.5.2 用双线性变换法设计iir数字滤波器双线性变换法从频域出发，实现s平面和z平面的一一映射关系，可消除频谱的混叠现象。2 语音信号分析和处理过程2.1 语音信号的采集 利用windows下得录音机（开始程序附件娱乐录音机，文件属性立即转换8khz，8位，单声道）录制一段自己的话音，或者采用windows自带的声音文件（默认为22050hz），时间控制在几秒左右，设置文件名为ding.wav，并保存在d:downloadmatlabbin下。然后在matlab软件平台下，利用函数wavread对语音信号进行分析。图2.1是基于pc机的语音信号采集过程，声卡可以完成语音波形的a/d转换，获得wave文件，为后续的处理储备原材料。调节录音机保存界面的“更改”选项，可以存储各种格式的wave文件。windows自带的录音机声音麦克风声卡滤波采样a/d转换 wav图2.1 基于pc机的语音信号采集过程采集到语音信号之后，需要对语音信号进行分析，如语音信号的时域分析、频谱分析、语谱图分析以及加噪滤波等处理。以下为一段程序，是语音信号在matlab中实现语音的读入打开，以及绘出了语音信号的波形频谱图。x,fs,bits=wavread(ding.wav,1024 6000); sound(x,fs,bits); x=fft(x,6000);magx=abs(x);angx=angle(x);subplot(221);plot(x);title(原始信号波形);subplot(222);plot(x); title(原始信号频谱);subplot(223);plot(magx);title(原始信号幅值);subplot(224);plot(angx);title(原始信号相位);程序运行可以听到声音，得到的结果如图2.2所示：图2.2 语音信号的读入与打开2.2 窗函数法设计2.2.1窗函数法低通滤波器窗函数设计低通滤波器程序如下：clc;close allz1,fs,bits=wavread(ding.wav)y1=z1(1:8192);y1=fft(y1);fp=1000;fc=1200;as=1000;ap=1;fs=8000;wc=2*pi*fc/fs; wp=2*pi*fp/fs;wdel=wc-wp;beta=0.112*(as-8.7);n=ceil(as-8)/2.285/wdel);wn= kaiser(n+1,beta); ws=(wp+wc)/2/pi;b=fir1(n,ws,wn);figure(1);freqz(b,1);x=fftfilt(b,z1);x=fft(x,8192);figure(2);subplot(2,2,1);plot(abs(y1);axis(0,1000,0,1.0);title(滤波前信号频谱);subplot(2,2,2);plot(abs(x);axis(0,1000,0,1.0);title(滤波后信号频谱);subplot(2,2,3);plot(z1);title(滤波前信号波形);subplot(2,2,4);plot(x);title(滤波后信号波形);sound(x,fs,bits);程序运行可以听到声音，得到的结果如图2.3所示图2.3 语音信号通过低通滤波器后的波形2.2.2 窗函数法高通滤波器窗函数设计高通滤波器程序设计如下：clc;close allz1,fs,bits=wavread(ding.wav)y1=z1(1:8192);y1=fft(y1);fp=4800;fc=5000;as=100;ap=1;fs=8000;wc=2*pi*fc/fs; wp=2*pi*fp/fs;wdel=wc-wp;beta=0.112*(as-8.7);n=ceil(as-8)/2.285/wdel);wn= kaiser(n,beta); ws=(wp+wc)/3/pi;b=fir1(n-1,ws,high,wn);figure(1);freqz(b,1);x=fftfilt(b,z1);x=fft(x,8192);figure(2);subplot(2,2,1);plot(abs(y1);axis(0,5000,0,1.0);title(滤波前信号频谱);subplot(2,2,2);plot(abs(x);axis(0,5000,0,1.0);title(滤波后信号频谱);subplot(2,2,3);plot(z1);title(滤波前信号波形);subplot(2,2,4);plot(x);title(滤波后信号波形);sound(x,fs,bits);程序运行可以听到声音，得到的结果如图2.4所示图2.4 语音信号通过高通滤波器后的波形2.2.3 窗函数法带通滤波器窗函数设计带通滤波器程序设计如下：clear;close allz1,fs,bits=wavread(ding.wav)y1=z1(1:8192);y1=fft(y1);fp1=1200 ;fp2=3000 ;fc1=1000 ;fc2=3200 ;as=100 ;ap=1 ;fs=8000 ;wp1=2*pi*fp1/fs; wc1=2*pi*fc1/fs; wp2=2*pi*fp2/fs; wc2=2*pi*fc2/fs;wdel=wp1-wc1;beta=0.112*(as-8.7);n=ceil(as-8)/2.285/wdel);ws =(wp1+wc1)/2/pi,(wp2+wc2)/2/pi;wn= kaiser(n+1,beta); b=fir1(n,ws,wn);figure(1); freqz(b,1)x=fftfilt(b,z1);x=fft(x,8192);figure(2);subplot(2,2,1);plot(abs(y1);axis(0,2000,0,1.0);title(滤波前信号频谱);subplot(2,2,2);plot(abs(x);axis(0,2000,0,0.03);title(滤波后信号频谱)subplot(2,2,3);plot(z1);title(滤波前信号波形);subplot(2,2,4);plot(x);title(滤波后信号波形);sound(x,fs,bits);程序运行可以听到声音，得到的结果如图2.5所示图2.5 语音信号通过带通滤波器后的波形2.3双极性变换法设计2.3.1 双极性变换法设计低通滤波器程序设计如下：clear;close allz1,fs,bits=wavread(ding.wav)y1=z1(1:8192);y1=fft(y1);fp=1000;fc=1200;as=100;ap=1; ;fs=8000;wc=2*fc/fs;wb=2*fp/fs;n,wp=cheb1ord(wc,wb,ap,as);b,a=cheby1(n,ap,wp);figure(1);freqz(b,a);x=filter(b,a,z1);x=fft(x,8192);figure(2);subplot(2,2,1);plot(abs(y1);axis(0,1000,0,1.0);title(滤波前信号频谱);subplot(2,2,2);plot(abs(x);axis(0,4000,0,0.03);title(滤波后信号频谱);subplot(2,2,3);plot(z1);title(滤波前信号波形);subplot(2,2,4);plot(x);title(滤波后信号波形);sound(x,fs,bits);程序运行可以听到声音，得到的结果如图2.6所示：图2.6 语音信号通过低通滤波器2.3.2 双极性变换法设计高通滤波器程序设计如下：clear;close allz1,fs,bits=wavread(ding.wav)y1=z1(1:8192);y1=fft(y1);fc=4800 ;fp=5000 ;as=100;ap=1; fs=8000;wc=fc/fs;wb=fp/fs;n,wp=cheb1ord(wc,wb,ap,as);b,a=cheby1(n,ap,wp,high);figure(1);freqz(b,a);x=filter(b,a,z1);x=fft(x,8192);figure(2);subplot(2,2,1);plot(abs(y1);axis(0,5000,0,1.0);title(滤波前信号频谱);subplot(2,2,2);plot(abs(x);axis(0,5000,0,0.03);title(滤波后信号频谱);subplot(2,2,3);plot(z1);title(滤波前信号波形);subplot(2,2,4);plot(x);title(滤波后信号波形);sound(x,fs,bits);程序运行可以听到声音，得到的结果如图2.7所示图2.7 语音信号通过高通滤波器2.3.3 双极性变换法设计带通滤波器程序设计如下clear;close allz1,fs,bits=wavread(ding.wav)y1=z1(1:8192);y1=fft(y1);fp1=1200 ;fp2=3000; fc1=1000 ;fc2=3200 ;as=100;ap=1; fs=8000; wc=2*fc1/fs,2* fc2/fs;wb=2*fp1/fs,2*fp2/fs; n,wp=cheb1ord(wc,wb,ap,as);b,a=cheby1(n,ap,wp,stop);figure(1);freqz(b,a);x=filter(b,a,z1);x=fft(x,8192);figure(2);subplot(2,2,1);plot(abs(y1);axis(0,1000,0,1.0);title(滤波前信号频谱);subplot(2,2,2);plot(abs(x);axis(0,4000,0,0.03);title(滤波后信号频谱);subplot(2,2,3);plot(z1);title(滤波前信号波形);subplot(2,2,4);plot(x);title(滤波后信号波形);sound(x,fs,bits);程序运行可以听到声音，得到的结果如图2.8所示图2.8 语音信号通过带通滤波器2.4 结果分析 从结构上看，iir滤波器必须采用递归结构来配置极点，并保证极点位置在单位圆内。由于有限字长效应，运算过程中将对系数进行舍入处理，引起极点的偏移。这种情况有时会造成稳定性问题，甚至产生寄生振荡。相反，fir滤波器只要采用非递归结构，不论在理论上还是在实际的有限精度运算中都不存在稳定性问题，因此造成的频率特性误差也较小。此外fir滤波器可以采用快速傅里叶变换算法，在相