《数字语音处理》实验指导书2005上课讲义.doc

Good is good, but better carries it.精益求精，善益求善。数字语音处理实验指导书2005-数字语音处理实验指导书学校：中国地质大学院系：机械与电子工程学院专业：电子信息工程编制：电信教研室版本：2005目录前言单元1MATLAB练习单元2综合分析实验附录1MATLAB软件简介附录2MATLAB基本使用方法附录3wavread函数附录4specgram函数参考文献单元1MATLAB练习实验目的l 了解信号谱分析特点，掌握各种窗函数功率谱的异同及其对短时信号谱分析的作用。l 明确短时傅立叶变换与声谱图的关系，理解其对语音信号的时频分析作用。l 学习音频文件的基本操作方法。l 学习用MATLAB实现基本的语音信号处理，提高自学和动手能力，培养学习兴趣。实验原理1、 离散频谱、幅度谱和功率谱离散频谱是信号序列的离散傅立叶变换，常用快速傅立叶变换（FFT）计算。离散频谱描述的是信号在频率采样点上的幅度和相位值，包括幅度-频率关系和相位频率关系，幅度谱就是前者。功率谱是幅度平方频率关系。由于语音的对数听觉效应，幅度谱和功率谱常用对数形式表达，单位为分贝。对数幅度谱=10lg幅度谱，对数功率谱=10lg幅度谱2=20lg幅度谱。由此可见，对数幅度谱与对数功率谱在谱形上是相似的，只不过谱的幅度方向相差一倍。对于窗序列的离散频谱、幅度谱和功率谱，结论也是类似的。2、 短时傅立叶变换和声谱图短时傅立叶变换的实质，是在时间和频率上同时对信号序列进行分解，使得信号被分解成若干个时段和频段上的分量（即分信号）。根据微分原理，若每个时段和频段都很小，对具体某个时段和频段上的分信号进行抽样，无论具体的抽样时点和频点如何选取，抽样值都不会有大的差异，因而一个分信号可以仅用一个该分信号的抽样值来表示。短时傅立叶变换的值X(n,k)实际上就是n时段和k频段上信号分量的幅度和相位抽样，这个抽样值大体代表了n时段和k频段上信号分量的幅度范围和相位范围。n和k的含义见下表。表1数字参数与模拟参数的对应关系数字时间模拟时间（s）数字频率模拟频率（HZ）nt=n/fsKf=k·fs/L注：L窗序列长度fs信号采样频率应用在语音处理方面，声谱图就是语音信号的短时傅立叶变换，其图形化显示往往只表示各个时段和频段上的分信号幅度抽样值。3、WAVE文件的基本操作WAVE是Windows操作系统支持的一种标准的数字波形音频文件格式。通常我们用麦克风和Windows录音机采集到计算中的声音就是WAVE格式的，WAVE文件后缀为.wav，多采用PCM格式或ADPCM格式，支持单声道或多声道数字音频数据采集和回放。与MP3等其它格式相比，WAVE格式尽管压缩比不好，却非常直观（除文件头外，实际上就是所有通道的采样数据按时间顺序交错排列），适合于数据量不大的音频录放（例如Windows提示音）。因此，对WAVE文件进行读取、保存、剪辑、播放等基本操作相对容易些。软件基础较好，有志于开发自己的音频应用的读者可到网上下载WAVE等文件的存储格式。实验器材l 多媒体计算机l 预装Windows98/2000/XP操作系统l 预装MATLAB6.5软件l 预装“数字语言处理实验”软件包l 自带移动存储实验内容1、输入并运行以下MATLAB代码：%（1）信号波形和对数功率谱示例%波形：100样点，抽样频率Fs=8kHz，信号频率f=2kHzx=sin(2*pi*(0:99)*2000/8000);%频谱：是波形的快速傅立叶变换X=fft(x);%对数功率谱：功率谱取10lg或幅度谱取20lgP=20*log10(abs(X);%绘制波形和对数功率谱（注意观察25和75处的尖峰）figure(1);plot(x);title('信号的波形');figure(2);plot(P);title('信号的对数功率谱');%（2）窗函数的对数幅度谱示例%避免对0取对数floor=0.1;%图1：矩形窗W_rec=10*log10(abs(fft(boxcar(512),1024)+floor);figure(1);plot(W_rec(1:end/2);title('矩形窗的对数幅度谱');%图2：海明窗W_hamming=10*log10(abs(fft(hamming(512),1024)+floor);figure(2);plot(W_hamming(1:end/2);title('海明窗的对数幅度谱');%图3：汉宁窗W_hanning=10*log10(abs(fft(hanning(512),1024)+floor);figure(3);plot(W_hanning(1:end/2);title('汉宁窗的对数幅度谱');2、编程绘制以下窗序列的对数幅度谱：（1）指数窗（2）矩形窗（3）海明窗（4）汉宁窗选N=512，用MATLAB编程绘制以上窗序列的幅度谱。3、编程读取WAVE文件并显示声谱图（取1024点汉宁窗）提示：4、 主要涉及的M函数有figure、plot、wavread、specgram等，理解如何使用它们。编程实现对WAVE文件的剪辑、保存和播放提示：主要涉及的M函数wavread、wavwrite、wavplay等，使用MATLAB帮助系统了解其用法。报告要求1、包括实验名称、实验目的、实验内容标题。2、在相应内容标题下：（1） 拷贝内容1的各个图形，标出谱零点。（2） 写出内容2的程序清单，说明几种窗函数谱的不同特点和对短时谱分析性能的影响。（3） 写出内容3的程序清单，拷贝声谱图图形，说明此声谱图形的物理意义。写出内容4的程序清单，说明你的程序是如何使用M函数实现WAVE文件的剪辑、保存和播放的。单元2综合分析实验实验目的1、 掌握利用短时自相关检测基音周期的原理，了解基音检测在清音段、谐音段、纯音段及其混杂信号段的适用情况。2、 掌握短时傅立叶变换的实现原理，理解并掌握数字频率与模拟频率间的对应关系，理解SFT对语音信号的时频分析作用。3、 初步了解线性预测的特点。实验原理1、 短时自相关基音检测设信号采样序列长度为M：x(m)|0mM，采样周期为fs矩形窗序列长度为N：w(m)|0mn，N<M则短时自相关为：其有效计算区间为：l 因此，对于矩形窗而言，短时自相关Rn(k)相当于取序列段x(m)|n-Nmn作自相关计算，最大值为Rn(0)。l 若所取序列段具有周期性，则Rn(k)将在k为序列段周期Tn时取得次大值Rnmax，因而取得Rnmax时的k值即为所取序列段周期：当Rn(k)=Rnmax时，Tn=kTs=k/fs。l 若Rnmax与Rn(0)相比相差太大，则表明所取序列段没有相关性，为清音帧。2、 短时傅立叶变换SFT的傅立叶解释X(m)w(n-m)w(m)x(m)0n-N窗选时间段nm注：SFT是窗选信号段的离散傅立叶变换，它反映窗选信号段的离散频谱；SFT系数Xn(ejk)是n-Nn时段信号在数字角频率k处的频谱抽样，含有该单频分量的幅度和相位信息，该单频分量用以反映所在频段的频谱特征。图1SFT的傅立叶解释SFT的滤波器解释W(n)X(n)Xn(ejk)e-jkn窗选信号频谱频率窗移频-f0Bff+B注：Xn(ejk)是数字信号X(n)被载波e-jkn调制移频后经过数字低通滤波器W(n)的输出。设k对应的模拟频率为f（HZ），低通滤波器通频带为0B（HZ）则只有频率范围在ff+B（HZ）的信号分量能通过低通滤波器。因此Xn(ejk)是频率范围在ff+B（HZ）內的信号频谱抽样。图2SFT的傅立叶解释（3）数字信号与模拟信号的参数对应关系对于矩形窗而言，数字参数与模拟参数存在以下对应关系：表1数字参数与模拟参数的对应关系数字时间模拟时间（s）数字频率模拟频率（HZ）数字角频率模拟角频率（rad/s）nt=n/fskf=k·fs/Lk=2·k/L=2·k·fs/L注：L窗序列长度fs信号采样频率（4）时频分辨率关系l 对于矩形窗，时间窗宽度=L/fs；频率窗宽度B=fs/L；时频窗面积常数K=·B=1。l Xn(ejk)是(n-L)/fsn/fs时段的信号在kB(k+1)B频段的频谱抽样，抽样频率为kB。l Xn(ejk)反映了在上述模拟时段和频段內的信号分量特征，若时间窗宽度越小，则时间分辨率越高，频率分辨率越低；反之，若时间窗宽度越大，则时间分辨率越低，频率分辨率越高。l SFT系数的集合Xn(ejk)是关于数字时间n和数字频率k的二维阵列，n、k的取值为：Xn(ejk)|n=N，3N/2，2N.M-1；k=0，1，2,N-1。其中M是语音信号的序列长度。3、线性预测（1）线性预测定义用短时窗选取一段信号序列x(m)|0mN，利用如下线性预测方程进行预测：其中N是窗序列长度，权系数ak|1kp是预测器参数（2）预测系数选取原则满足信号序列与其预测值序列间的均方预测误差最小准则。（3）窗长和预测阶数N过大，预测时段超出了语音信号的短时平稳区间，将使预测精度降低；而N过小，则不能由过少的信号样点得到其变化规律，同样将使预测精度降低。类似地预测阶数P的选取也应恰当。实验器材l 多媒体计算机l 预装Windows98/2000/XP操作系统l 预装VisualBasic6.0软件l 预装VisualC+6.0软件l 预装“2005数字语言处理实验”软件包l 自带移动存储实验内容1、 短时自相关基音检测（1） 打开“2005数字语音处理实验单元2：综合分析实验”文件夹，运行其中的VoiceProcess.exe实验程序，打开Sample.raw采样序列数据文件。（2） 点击“基音检测”菜单项，拷贝基音检测结果；根据基音检测结果并对比波形，指出哪些语段是浊音，哪些语段是清音，哪些语段两者都有，说明理由。2、 短时傅立叶变换（1） 点击“短时傅立叶变换”菜单项，观察SFT系数Xn(ejk)和各时段|Xn(ejk)|短时幅度谱。（2） 选取较短的窗长，在SFT谱值表中记录完整的一帧数据，同时拷贝下相对应的短时幅度谱图形。分析说明：l 短时谱的横、纵坐标分别代表什么。l 短时频谱与短时幅度谱的异同。l 根据所选时段、窗型和窗长，在SFT谱值表中任取一个短时谱谱值，确定其所分析的信号时段（秒）和频段（赫兹）。（3） 选取信号的单频正弦区间，改变窗长，观察并拷贝所显示的2-3个不同短时幅度谱图；注意窗长改变时，短时幅度谱对该单频波的频率分析区间宽度和幅度谱分布规律有何变化，从而理解频率分辨率概念，就此说明窗长改变对频率分辨率的影响。3、 线性预测（1） 点击“线性预测”菜单项，改变各观察时段，根据信号波形和预测波形的对比、以及预测信噪比统计数据，举例说明在信号段具有什么特征时预测准确性较高，何时较低。（2） 选取预测性能较好的一个时间段，改变窗长，观察并记录预测图形和预测信噪比数据，说明窗长的改变对预测准确性有什么影响，想想为什么？报告要求1、包括实验名称、实验目的、实验内容标题。2、在相应内容标题下，写出或拷贝实验内容所要求记录的数据和图形，完成实验内容所要求的分析说明。附录MATLAB软件简介MATLAB是什么？MATLAB是一套高性能的科学计算可视化软件，它是由操作界面、MATLAB语言、通用函数库、工具箱和帮助系统组成的一个科学计算集成实验环境。在这个环境下，对所要求解的问题，用户只需简单地列出数学表达式，其结果便以数值或图形方式显示出来。MATLAB有哪些主要特点？l 基于矩阵的运算单元。MATLAB的含义是矩阵实验室，现已发展成线性代数课程的标准工具，其基本计算单元是矩阵，即使是标量运算，也是作为1行1列的矩阵进行的。MATLAB针对矩阵运算做了很多优化，若使用循环语句将矩阵分解成标量进行运算，将在运算效率上比直接进行矩阵运算有显著降低。l 易扩展性。MATLAB函数库是由一系列函数（称为M文件）组成。除通用函数库外，MATLAB还包括了被称为ToolBox（工具箱）的特定学科问题求解工具，工具箱实际上是对MATLAB进行扩展应用的专用函数库。由于用户可自行建立M文件进而构成特定领域的工具箱，工具箱种类和内容都在不断扩充中，现已发展了信号处理、图象处理、控制系统、神经网络、财经等几十种工具箱，使得MATLAB具有广泛的应用领域。易学易用。与其它高级语言一样，MATLAB语言有其内定的规则，但其更接近于数学表示，因此其使用更为简便，避免了其它语言如C、FORTRAN中的许多限制，如变量、矩阵无须定义。基于功能齐全的多维向量运算功能、丰富实用的通用和专用函数库，一条MATLAB语句往往可实现相当于几十条甚至几百条C语句才能实现的复杂功能，使得用户无须深入了解算法的具体执行过程和编程技巧，通过简单的编程就可得到结果。演算纸和程序调试相结合，使得验证语句用法和定位错误变得简单。其可视化操作环境、丰富的图形函数和图形用户界面（GUI），便于实现科学计算过程的交互控制和结果的可视化。此外，MATLAB还提供了内容翔实的帮助系统。MATLAB能做什么？l 由于MATLAB具有如此简单而强大的科学演算和可视化功能，非常适合于多学科领域的算法设计和工程实践的仿真实验。在学习和研究方法上：在当今知识爆炸的时代，任何人不可能穷尽所有的知识；在工程实践中，人们往往是在确定问题定义和研究方案后，再对涉及的技术细节作深入细致的研究；因此，提出问题和尽可能了解各种可能解法的效果，可以避免研究工作的盲目性，是优先于掌握解法的实现细节的。MATLAB具有各个领域丰富实用的工具箱函数库，可以避免在复杂的编程技巧和具体的算法实现上过多地分散精力，迅速直观地观察到算法的运行结果，加深对算法核心物理意义的理解，十分适合于概念性学习和预研究工作。正因为如此，在美国和许多西方国家，MATLAB早已作为本科理工科专业的必修课程。附录MATLAB基本使用方法写在开始的话：一、 对于任何一种编程软件的使用，自学和实践是最好的老师：首先，找一本入门书了解基本知识和操作方法，然后更多地是举一反三，利用方便的帮助系统即用即学所需知识，并在实际编程和调试中验证你的理解。由于MATLAB的功能和函数十分丰富，要在一本实验指导书中完整地加以描述是不可能的，这里只介绍了一些最基本的和与实验内容相关的知识。MATLAB基本知识1、 MATLAB命令窗口运行MATLAB软件后首先看到的就是命令窗口，在命令窗口中，在MATLAB提示符下可键入MATLAB命令。例如输入一个3×3的矩阵：a=123;456;789按回车键后显示：a=123456789但如果你忘不了C语言，输入的是：a=123;456;789;则按回车键后什么都不会显示。因为分号在MATLAB命令窗口中是作为抑制显示符号，尽管变量a已经存在并被赋值，但你看不见它。抑制显示符号的一个重要作用是批处理运算：你可以先在文本文档中编写好由多条语句组成的一段程序，除最后一条外，所有其它语句都加上抑制显示符号，然后粘贴到命令窗口并回车，那么中间运算结果就不会显示，只看到最终结果。也许你对命令窗口的诸多内容感到眼花缭乱，那就试一下“Edit|ClearCommandwindow”吧，它会像橡皮一样擦干净整张白纸，但是工作空间中已经存在的变量是MATLAB环境全局变量，它并不会消失，只是看不到而已，当你敲入变量名并回车后它又会显示出来。MATLAB命令窗口的以上特性使其被称为演算纸（我们俗称草稿纸，而且是一张可以反复利用并且记性很好的草稿纸）。命令窗口除了作为演算纸外，它还是其它MATLAB功能的出发点。例如：可以通过“File”菜单新建或打开M文件、图形文件和图形用户界面（用于设计交互式程序的特殊图形文件）；可以通过“Help”菜单获取MATLAB帮助知识和演示程序。2、 M文件由MATLAB语言编写的程序文件称为M文件，扩展名为.m。M文件可以在命令窗口提示符下键入文件名来直接调用（不需要编辑、调试时），也可以通过命令窗口的文件菜单打开M文件编辑器（需要编辑、调试时），使用M文件的最大好处是它可以调试。从功能上看，M文件可分为两类：底稿文件和函数文件。（1）底稿文件底稿文件中的语句可使用工作空间中的全部数据（包括命令窗口产生的数据），例如：有一包含以下MATLAB命令的底稿文件fibon.m：%anMfiletocalculateFibonaccinumbers（斐波纳契数列）f=11;I=1;Whilef(I)+f(I+1)<1000f(I+2)=f(I)+f(I+1);I=I+1;Endf其中“%”右边的语句为说明语句，它们只起到注释或帮助的作用。在MATLAB提示符下，如键入fibon，则MATLAB会自动执行这一文件中的每条命令，并产生执行结果：输入：fibon结果：f=1123581321345589144233377610987注意，在底稿文件中的变量I和f都将保存在工作区中作为全局变量而存在，这一点与函数文件是不同的。（2）函数文件函数文件的第一行必须包含关键字“function”。函数文件与底稿文件的区别在于：函数文件可以传递参数，底稿文件不具备参数传递功能；在函数文件中定义及使用的变量都是局部变量，只在本函数的内有效，一旦退出该函数，则为无效变量，而底稿文件中定义或使用的变量都是全局变量，在退出文件后仍为有效变量。例如，函数文件mean.m包含以下语句：functiony=mean(x)%MEANaverageormeanvalue%Forvectors,MEAN(x)returnthemeanvalue%Formatrices,MEAN(x)isarowvector%containingthemeanvalueofeachcolumnmn=size(x);ifm=1m=n;endy=sum(x)/m;这个M文件定义了一个新函数mean，它的引用与其它MATLAB函数一样，其功能是计算向量或矩阵的平均值，例如：输入：z=1:99;m=mean(z)结果：m=50（3）M文件的调试MATLAB图形窗口3、剪切板的使用4、编程指南（1）命令窗口编程（2）M文件编程和调试（3）图形用户界面（4）控制流语句6、帮助系统二、 基本操作命令1、基本知识（1）简单矩阵的输入（2）矩阵元素（3）复数和复数矩阵（4）MATLAB语句和变量MATLAB语言是与大小写有关的语言，即变量A和a是完全不同的变量，应该注意所有的函数名均由小写字母构成，变量可由字母（大小写随意）和数字组成。（5）退出和保存工作空间（6）常数与算术运算符（7）函数（8）帮助命令2、矩阵运算（1） 矩阵转置（2） 矩阵加减（3） 矩阵乘法（4） 矩阵除法（5） 矩阵乘方（6） 矩阵超越函数1、 数组运算（1） 数组加减（2） 数组乘除（3） 数组乘方（4） 关系运算（5） 逻辑运算（6） 数学函数2、 向量和矩阵操作（1） 向量产生（2） 下标（3） 空矩阵（4） 构造大矩阵3、 绘图功能（1）建立图形对象（2）绘制简单的曲线（3）给图形添加标题6、文件操作附录wavread函数用途读取微软WAVE波形声音文件（.wav）。句法y=wavread('filename')y,Fs,bits=wavread('filename').=wavread('filename',N).=wavread('filename',N1N2).=wavread('filename','size')描述Wavread函数支持多通道数据，最多可支持32位采样并支持读取24位和32位的.wav文件。y=wavread('filename')载入由filename字符串指定的一个WAVE文件，向y矢量返回采样数据。如果没有给出文件扩展名，函数自动附加.wav扩展名。样点幅度值介于-1,+1范围。y,Fs,bits=wavread('filename')返回赫兹单位的采样率（Fs），以及用于文件数据编码的每样点比特数（bits）。.=wavread('filename',N)仅返回文件各通道数据的前N个样点。.=wavread('filename',N1N2)仅返回文件各通道数据的第N1到第N2样点。siz=wavread('filename','size')返回文件所含音频数据的大小而不是实际的音频数据，返回矢量siz的格式为样点数通道数。参见auread，wavwrite，wavplay，wavrecord附录4specgram函数用途时频分析（产生声谱图）。句法B=specgram(a)B=specgram(a,nfft)B,f=specgram(a,nfft,fs)B,f,t=specgram(a,nfft,fs)B=specgram(a,nfft,fs,window)B=specgram(a,nfft,fs,window,numoverlap)specgram(a)B=specgram(a,f,fs,window,numoverlap)描述specgram函数使用滑动窗计算信号的短时傅立叶变换。声谱图是该函数的幅度值。（1）B=specgram(a)计算矢量a所表示信号的短时傅立叶变换。该句法使用默认参数值：l nfft取矢量a长度和256两者较小的一个l fs=2l window是长度为nfft的周期性汉宁窗l numoverlap等于窗长的1/2其中：参数nfft指定specgram函数所用FFT长度，该值决定了计算短时傅立叶变换的各频点；参数fs是指定采样频率的一个标量；参数window指定了一个窗函数以及specgram函数用以分割矢量a的样点数；参数numoverlap是分割区域重叠的样点数。你从最后一个句法的输入参数列表中省略的任意参数将使用以上默认值。如果矢量a是实数，specgram函数仅在正频点计算短时傅立叶变换。如果nfft为偶数，specgram函数返回nfft/2+1行SFT值（包含0和奈奎斯特带宽）；如果nfft为奇数，specgram函数返回nfft/2行SFT值。返回矢量B中的列数为：k=fix(n-numoverlap)/(length(window)-numoverlap)如果矢量a是复数，specgram函数在正频点和负频点都计算短时傅立叶变换。在这种情况下，返回矢量B是一个nfft行的复数矩阵。从B中第1列的第1个样点开始，时间沿列号递增方向线性递增；频率从0开始，沿行号递增方向线性递增。（2）B=specgram(a,nfft)使用指定的FFT长度nfft进行其计算。（3）B,f=specgram(a,nfft,fs)返回的矢量f是该函数计算短时傅立叶变换的各个频点。参数fs对于输出B没有影响，它只不过是频率尺度倍乘系数。（4）B,f,t=specgram(a,nfft,fs)分别返回频率矢量f和时间矢量t。参数t是一个尺度化时间的列矢量，其长度等于矩阵B的列数。t(j)是第j个窗口截取矢量a的最早时刻，t(1)总是等于0。（5）B=specgram(a,nfft,fs,window)指定了窗函数和矢量a每个窗口分区的样点数。如果你提供了窗口尺度，specgram函数使用该长度的汉宁窗。窗长必须小于等于nfft。（6）B=specgram(a,nfft,fs,window,numoverlap)以numoverlap个样点重叠a的各分区。你可以利用空矩阵为任意输入参数指定默认值。例如：B=specgram(x,10000)，除了具有采样频率10000Hz而不是默认的2Hz以外，其它都等同于B=specgram(x)。（7）specgram(.)没有输出参数，使用函数imagesc(t,f,20*log10(abs(b),axisxy,colormap(jet)在当前图形窗口显示对数尺度声谱图。模式axisxy在坐标系的左下角显示信号第1部分的低频成分。Specgram函数根据真实的时间和频率，利用来刻度坐标轴。（8）B=specgram(a,f,fs,window,numoverlap)利用高斯z变换（多于20个平坦分布的频点）或者10取其1多相滤波器组，在矢量f所指定的频点处计算声谱图。参数f是以赫兹为单位的频点矢量，它至少应含有2个元素。示例显示数字化语音信号的声谱图：loadmtlbspecgram(mtlb,512,Fs,kaiser(500,5),475)title('Spectrogram')注意：你可以在信号处理工具箱的Demo中观看和操作类似的声谱图，可从MATLAB基础（查看菜单的LaunchPad菜单项）访问该Demo。算法specgram函数按如下过程计算给定信号的声谱图：1、 它应用窗口参数所指定的窗口，将信号分割成重叠的分区。2、 用nfft长度的FFT计算各分区的短视时傅立叶变换，以产生对信号短期频率成分的一个估计；这些变换组成了矩阵B的各个列。数值(length(window)-numoverlap)指定经过多少个样点specgram函数改变一次窗口。对于实数输入，specgram函数将声谱图截短到开始的nfft/2+1个点（若nfft为偶数）或(nfft+1)/2个点（若nfft为奇数）。诊断当使用了错误的参数时回显示相应的诊断消息：Requireswindow'slengthtobenogreaterthantheFFTlength.RequiresNOVERLAPtobestrictlylessthanthewindowlength.RequirespositiveintegervaluesforNFFTandNOVERLAP.Requiresvectorinput.参见cohere,csd,pwelch,tfe参考文献1楼顺天等编著.MATLAB程序设计语言.西安电子科技大学出版社,1998.2楼顺天等编著.基于MATLAB的系统分析与设计信号处理.西安电子科技大学出版社,1998.3Oppenheim,A.V.,andR.W.Schafer,Discrete-TimeSignalProcessing,Prentice-Hall,EnglewoodCliffs,NJ,1989,pp.713-718.4Rabiner,L.R,andR.W.Schafer,DigitalProcessingofSpeechSignals,Prentice-Hall,EnglewoodCliffs,NJ,1978.-