《数字语音处理》实验指导书2015

1、数字语音处理 实验指导书 主编 王都生大连民族大学信息与通信工程学院2015年10月 目 录实验一 汉明窗加窗的时域波形以及频谱的显示.3实验二WAV格式文件的录制及其短时能量、短时平均幅度分析.5实验三 语音端点检测.7实验四 语音信号的倒谱的编程实现.11实验五 基于短时自相关函数的基音周期估计.15实验六 基于短时自相关函数的基音周期估计续.19实验七 声道冲激响应及激励信号的分离实现.23实验一 汉明窗加窗的时域波形以及频谱的显示一、实验环境 1Windows XP，Windows 98，Windows2000 都行。 2Matlab 二、实验目的 1 掌握加窗时域波形显示程序的编写。

2、 2了解清音和浊音的时域波形特征。 3了解清音和浊音的频谱特征。 三、实验原理 1窗函数 a 汉明窗函数 b 矩形窗函数 2. 给输入语音信号加窗 3.MATLAB程序设计 分为两步，对于浊音的时域波形显示及其频谱，对于清音的时域波形及其频谱的计算和显示。 其中，语音信号必须放在该程序所在的文件夹下。汉明窗的短时傅里叶变换y,Fs,bits=wavread(1.wav) %读出信号，采样率和采样位数。y=y(:,1) %这里假设声音是双声道，只取单声道作分析。y=y(4000:4400) %分别这里选择浊音信号和清音信号分别进行讨论。y=y.*hamming(length(y)figure s

3、ubplot(2,1,1)plot(y) sigLength=length(y)Y=zeros(Fs,1) Y=20*log10(abs(fft(y,sigLength)subplot(2,1,2) plot(Y) % 一点对应 1HZ xlabel(Frequency(Hz) ylabel(20lg|Xn(ejw)|)四、 实验步骤 1汉明窗加窗的时域波形及其短时频谱显示的源代码。 录音 wavrecord (N,Fs,CH) 读语音 wavread(“sound.wav”) 截取清浊音 加窗 y=y1.*hamming(length(y1) 2利用编写的源代码分别显示浊音和清音的时域波形及

4、频谱。求频谱 20lg|xn(e jw )|五、实验报告要求 总结语音信号的频谱特征，并指出加矩形窗和汉明窗的浊音信号、清音信号的频谱特征。实验二 WAV格式文件的录制及其短时能量、短时平均幅度分析一、实验环境 1、Windows XP，Windows 98，Windows 2000 都行。 2、Matlab 二、实验目的 1、掌握利用MATLAB 录制 WAV格式的语音文件编程。 2、掌握清音和浊音的短时能量和短时平均幅度计算的编程。 3、了解清音和浊音在不同的窗函数下的短时能量和短时平均幅度的特点。 4、了解在不同窗长的情况下，语音信号的短时能量及平均幅度的变化的规律。 三、实验原理 1、

5、短时能量定义为：其中，N为窗长。可见，短时能量为一帧样点值的加权平方和。若令：， 则 短时能量可以理解为：语音信号各个样点值平方，然后样点通过一个冲激响应为 h(n)的滤波器，输出为由短时能量构成的时间序列。因此，冲激响应 h(n)的选择或者说窗函数的选择直接影响着短时能量的计算。 短时平均幅度函数 Mn用来衡量语音幅度的变化，定义为：2、MATLAB 程序设计 a） 利用 MATLAB 进行录音。WAVRECORD 利用 WINDOWS 的音频输入设备进行录音。 格式：WAVRECORD(N,FS,CH) 其中 N 为产生语音序列的长度，FS 为选择的采样频率，CH 为输入声道数。 标准的采

6、样频率为 8000, 11025, 22050, and 44100 Hz. CH 是 1或者 2（单声道或双声道）. 采样序列的长度为 N CH. 如没特殊说明， FS=11025 Hz, and CH=1. b）短时能量及短时平均幅度的参考程序：短时能量的参考程序： a=wavread(10.wav) subplot(6,1,1) plot(a) N=32 for i=2:6 h=linspace(1,1, (i1)*N) %形成一个矩形窗，长度为N En=conv(h,a.*a) %求卷积得其短时能量函数 En subplot(6,1,i) plot(En) if(i=2) legend

7、(N=32) elseif(i=3) legend(N=64) elseif(i=4) legend(N=128) elseif(i=5) legend(N=256) elseif(i=6) legend(N=512) end end % 加 hamming窗的短时能量函数： % 把 h=linspace(1,1, (i1)*N) % 改为 h1=hamming(i1)*N)短时平均幅度的参考程序： % 加矩形窗的短时平均幅度： a=wavread(10.wav) subplot(6,1,1),plot(a) N=32for i=2:6h=linspace(1,1,(i1)*N)%形成一个矩形

8、窗，长度为 N En=conv(h,abs(a) %求卷积得其短时能量函数En subplot(6,1,i),plot(En) if(i=2) legend(N=32) elseif(i=3) legend(N=64) elseif(i=4) legend(N=128) elseif(i=5) legend(N=256) elseif(i=6) legend(N=512)end end 四、实验步骤 1、分别在矩形窗和汉明窗加窗时的时域波形的程序。 2、编写计算短时能量及短时平均幅度的程序。 五、实验报告要求 修改程序， 不用滤波器卷积的方法实现而用一帧样点值加权平方和实现。实验三 语音端点检

9、测一、实验环境 1、Windows XP，Windows 98，Windows 2000 都行。 2、Matlab 二、实验目的 1、掌握短时能量计算的编程及了解清、浊音波形特点。 2、掌握短时平均过零率计算的编程及了解清、浊音波形特点。 3、掌握语音端点检测的判断方法。三、实验原理 1、短时能量定义为：其中，N为窗长。可见，短时能量为一帧样点值的加权平方和。若令：， 则 短时能量可以理解为：语音信号各个样点值平方，然后样点通过一个冲激响应为 h(n)的滤波器，输出为由短时能量构成的时间序列。因此，冲激响应 h(n)的选择或者说窗函数的选择直接影响着短时能量的计算。2、短时平均过零率 在矩形窗

10、条件下，短时平均过零率定义为：短时平均过零率可以粗略估计语音的频谱特性。 由于高频率对应着高的过零率，低频率对应着低过零率，那么过零率与语音的清浊音特征就存在着 对应关系。一般经验结论是，清音和浊音的过零率分布大致为高斯分布，清音每 10ms 的短时平均过零次数的均值为 49 次， 浊音每 10ms 的短时平均过零次数的均值为 14 次。 3、语音端点的判断 在语音端点检测中，可以采用多种特征参数综合判决的方法，如使用短时平均能量和短时平均过零率来进行检测。 由于语音一般都存在能量较高浊音段,因此考察语音的短时能量的轮廓可以设定一个较高的门限 T1，使语音的起点和终点落在 T1 所确定的时间间

11、隔 AB 之外。然后根据背景噪声的平 均能量确定个门限较低的 T2，并从 A点往起点方向、从 B 点往终点方向搜索，分别找到与门限 T2 相交的两个点 C、D。这样我们就用双门限完成了第级粗判。第二级判决要利用短时平均过零率。同样根据背景噪声的 Zn 可以设定一个较低的门限 T3，从 C 点和 D 点分别向起点和终点方向搜索。可以找到 Zn与门限 T3 相交的两个点 E、F。这样就确定出了语音的端点 E、F。 4、MATLAB 程序设计参考程序： % 短时过零率： a=wavread(6.wav) n=length(a) N=320 subplot(3,1,1)；plot(a) h=linsp

12、ace(1,1,N) %形成一个矩形窗，长度为N En=conv(h, a.*a) %求卷积得其短时能量函数En subplot(3,1,2);plot(En) for i=1:n-1if a(i)=0 b(i)= 1 else b(i) =-1 end if a(i+1)=0 b(i+1)=1else b(i+1)=1 end w(i)=abs(b(i+1)b(i) end %求出每相邻两点符号的差值的绝对值 k=1 j=0while (k+N1)n Zm(k)=0for i=0:N1 Zm(k)=Zm(k)+w(k+i) endj=j+1 k=k+160 %每次移动半个窗 end for

13、w=1:j Q(w)=Zm(160*(w1)+1)/640 %短时平均过零率 end subplot(3,1,3),plot(Q) 四、实验步骤 1、短时能量计算的编程。 2、短时平均过零率计算的编程。 3、判断语音的端点（起始点及终点）。 五、实验报告要求 1 利用短时平均幅度代替 En计算； 2 修改程序是的 Zn与 En的样点值一一对应； 3 将 Zn与 En的横轴用时间表示； 4 实验报告上画出语音端点的判断过程。实验五 基于短时自相关函数的基音周期估计一、实验环境： 1Windows XP，Windows 98，Windows2000 都行。 2 matlab. 3、耳机，话筒。 二

14、、实验目的 1编写计算短时自相关函数的程序。 2从一段语音信号中求出该段语音信号的基音周期。 三、实验原理 1、基于短时自相关函数的基音周期估计原理 短时自相关函数在基音周期的整数倍位置存在较大的峰值。如果找出第 一个最大峰值的位置就可以估计出基音周期。但由于共振峰的影响，通常会 产生“倍频”或者“半频”的错误。一般可由两种方法解决： A带通滤波的方法 将输入信号通过一个频率范围为60,900Hz 的带通滤波器后，再进行 基音估计。因为最高基音频率为 450Hz,所以上截频设为 900Hz 可以保留语 音的一二次谐波。下截频为 60Hz是为了抑制 50Hz 电源干扰。 B中心削波法 它采用下面

15、的中心削波函数进行预处理：一般削波电平 T取本帧语音最大幅度的 6070。 将削波后的序列 y(n)用短时自相关函数估计基音周期，在基音周期位置的峰值更加尖锐，可 以有效减少倍频或半频错误。 2、软件实现 参考程序如下： function jiyinzhouqi(filename,shift) %短时自相关分析 %filename语音文件*.wav %zhouqi 基音周期 shift=10 signal,fs=wavread(FemaleVoice.wav)shift=round(fs*shift) %帧移 n1=fix(fs*0.97)+1 %分析起点970ms,帧长 30ms n2=fi

16、x(fs*1)+1ii=1 for ii=1:(length(signal)n1)/shift %分析次数 if n2aver/5) value(index)=0 %去除大野点的影响len=lenlength(index) for jj=1:3:len/3 %中值平滑，滑动窗口宽度 3，精度为中值 1/4（剔除野点） average=(value(jj)+value(jj+1)+value(jj+2)/3 for kk=1:3 if abs(value(jj1+kk)average)average/4 value(jj1+kk)=0 %将野点置零，同时数组长度减一 len=len1 enden

17、d end figure(2) stem(value) axis(0 length(value) 0 max(value) Tp=sum(value)/len/fs %求基音周期（Tp） function nmax=find_maxn(r) %寻找峰值最大的 n值及基音周期 %r,自相关序列 %maxn,为峰值最大的 n zer=find(r=0) %找第一个零点如果存在 jiaocha=0 %找第一近零点 ii=1 while (jiaocha0 & r(ii+1)0 & (ii+1)0 %检查是否存在零点 if zer(1)jiaocha %存在，则和 jiaocha 比较大小，用于祛除前

18、点的对基音周 期的查找带来的影响 jiaocha=zer(1) end end r(1:jiaocha)=0 %祛除影响 maxn=max(r) %找最大值 temp=find(r=maxn) %返回第一个最大值nmax=temp(1) 四、实验步骤 1应用 MATLAB 函数进行录音。 2用 matlab 编写实现代码。 3. 计算出所录语音的基音周期 五、实验报告要求 利用 matlab 实验 AMDF 基音检测法的实现。实验六 基于短时自相关函数的基音周期估计续一、实验环境： 1Windows XP，Windows 98，Windows2000 都行。 2 matlab. 3、耳机，话筒

19、。 二、实验目的 1编写计算短时自相关函数的程序。 2. 掌握中心削波法的短时自相关函数的基音周期估计程序编写。 3从一段语音信号中求出该段语音信号的基音周期。 三、实验原理 1、基于短时自相关函数的基音周期估计原理 短时自相关函数在基音周期的整数倍位置存在较大的峰值。如果找出第 一个最大峰值的位置就可以估计出基音周期。但由于共振峰的影响，通常会 产生“倍频”或者“半频”的错误。一般可由两种方法解决： A带通滤波的方法 将输入信号通过一个频率范围为60,900Hz 的带通滤波器后，再进行 基音估计。因为最高基音频率为 450Hz,所以上截频设为 900Hz 可以保留语 音的一二次谐波。下截频为

20、 60Hz是为了抑制 50Hz 电源干扰。 B中心削波法 它采用下面的中心削波函数进行预处理： 一般削波电平 T 取本帧语音最大幅度的 6070。将削波后的序列 y(n) 用短时自相关函数估计基音周期，在基音周期位置的峰值更加尖锐，可以有效减少倍频或半频错误。 2、软件实现 参考程序如下： function jiyinzhouqi2(filename,shift) %短时自相关分析 %filename语音文件*.wav %zhouqi 基音周期 %wld2shift=10 signal,fs=wavread(MaleVoice.wav) kk=0 while kk0.2 signal(kk)=

21、signal(kk)0.2 else if signal(kk)0.2&signal(kk)0.2 signal(kk)=0 end end end shift=round(fs*shift) %帧移 n1=fix(fs*0.97)+1 %分析起点970ms,帧长 30ms n2=fix(fs*1)+1ii=1 for ii=1:(length(signal)n1)/shift %分析次数 if n2aver/5) value(index)=0 %去除大野点的影响 len=lenlength(index) for jj=1:3:len/3 %中值平滑，滑动窗口宽度 3，精度为 中值 1/4（剔

22、除野点） average=(value(jj)+value(jj+1)+value(jj+2)/3 for kk=1:3 if abs(value(jj1+kk)average)average/4 value(jj1+kk)=0 %将野点置零，同时数组长度减一 len=len1 end end end %figure(4) %plot(average) figure(2) stem(value) axis(0 length(value) 0 max(value) Tp=sum(value)/len/fs %求基音周期（Tp）function nmax=find_maxn(r) %寻找峰值最大的

23、 n值及基音周期 %r,自相关序列 %maxn,为峰值最大的 n zer=find(r=0) %找第一个零点如果存在 jiaocha=0 %找第一近零点 ii=1 while (jiaocha0 & r(ii+1)0 & (ii+1)0 %检查是否存在零点 if zer(1)jiaocha %存在，则和 jiaocha 比较大小，用于祛除前点的对基音周期的查找带来的影响 jiaocha=zer(1) end end r(1:jiaocha)=0 %祛除影响 maxn=max(r) %找最大值 temp=find(r=maxn) %返回第一个最大值 nmax=temp(1) 四、实验步骤 1应用

24、 MATLAB 函数进行录音。 2用 matlab 编写实现代码。 3. 计算出所录语音的基音周期。 五、实验报告要求 编程实现三电平削波法的基音周期检测附录 MATLAB基本使用方法一、 MATLAB基本知识1、 MATLAB命令窗口运行MATLAB软件后首先看到的就是命令窗口，在命令窗口中，在MATLAB提示符下可键入MATLAB命令。例如输入一个33的矩阵：a=1 2 3; 4 5 6; 7 8 9按回车键后显示：a= 1 2 3 4 5 6 7 8 9但如果你忘不了C语言，输入的是：a=1 2 3; 4 5 6; 7 8 9;则按回车键后什么都不会显示。因为分号在MATLAB命令窗口中

25、是作为抑制显示符号，尽管变量a已经存在并被赋值，但你看不见它。抑制显示符号的一个重要作用是批处理运算：你可以先在文本文档中编写好由多条语句组成的一段程序，除最后一条外，所有其它语句都加上抑制显示符号，然后粘贴到命令窗口并回车，那么中间运算结果就不会显示，只看到最终结果。也许你对命令窗口的诸多内容感到眼花缭乱，那就试一下“Edit | Clear Command window”吧，它会像橡皮一样擦干净整张白纸，但是工作空间中已经存在的变量是MATLAB环境全局变量，它并不会消失，只是看不到而已，当你敲入变量名并回车后它又会显示出来。MATLAB命令窗口的以上特性使其被称为演算纸（我们俗称草稿纸，

26、而且是一张可以反复利用并且记性很好的草稿纸）。命令窗口除了作为演算纸外，它还是其它MATLAB功能的出发点。例如：可以通过“File”菜单新建或打开M文件、图形文件和图形用户界面（用于设计交互式程序的特殊图形文件）；可以通过“Help”菜单获取MATLAB帮助知识和演示程序。2、 M文件由MATLAB语言编写的程序文件称为M文件，扩展名为.m。M文件可以在命令窗口提示符下键入文件名来直接调用（不需要编辑、调试时），也可以通过命令窗口的文件菜单打开M文件编辑器（需要编辑、调试时），使用M文件的最大好处是它可以调试。从功能上看，M文件可分为两类：底稿文件和函数文件。（1）底稿文件底稿文件中的语句可

27、使用工作空间中的全部数据（包括命令窗口产生的数据），例如：有一包含以下MATLAB命令的底稿文件fibon.m：% an M file to calculate Fibonacci numbers（斐波纳契数列）f=1 1;I=1;While f(I)+f(I+1)1000 f(I+2)=f(I)+f(I+1); I=I+1;Endf其中“%”右边的语句为说明语句，它们只起到注释或帮助的作用。在MATLAB提示符下，如键入fibon，则MATLAB会自动执行这一文件中的每条命令，并产生执行结果：输入：fibon结果：f= 1 1 2 3 5 8 13 21 34 55 89 144 233 3

28、77 610 987注意，在底稿文件中的变量I和f都将保存在工作区中作为全局变量而存在，这一点与函数文件是不同的。（2）函数文件函数文件的第一行必须包含关键字“function”。函数文件与底稿文件的区别在于：函数文件可以传递参数，底稿文件不具备参数传递功能；在函数文件中定义及使用的变量都是局部变量，只在本函数的内有效，一旦退出该函数，则为无效变量，而底稿文件中定义或使用的变量都是全局变量，在退出文件后仍为有效变量。例如，函数文件mean.m 包含以下语句：function y=mean(x)% MEAN average or mean value% For vectors, MEAN(x)

29、return the mean value% For matrices, MEAN(x) is a row vector% containing the mean value of each columnm n=size(x);if m=1 m=n;endy=sum(x)/m;这个M文件定义了一个新函数mean，它的引用与其它MATLAB函数一样，其功能是计算向量或矩阵的平均值，例如：输入：z=1:99;m=mean(z)结果：m= 50附录 wavread函数用途读取微软WAVE波形声音文件（.wav）。句法y = wavread(filename)y,Fs,bits = wavread(f

30、ilename). = wavread(filename,N). = wavread(filename,N1 N2). = wavread(filename,size)描述Wavread函数支持多通道数据，最多可支持32位采样并支持读取24位和32位的.wav 文件。 y = wavread(filename) 载入由filename字符串指定的一个WAVE 文件，向y矢量返回采样数据。如果没有给出文件扩展名，函数自动附加.wav扩展名。样点幅度值介于-1,+1范围。y,Fs,bits = wavread(filename) 返回赫兹单位的采样率（Fs），以及用于文件数据编码的每样点比特数（b

31、its）。. = wavread(filename,N) 仅返回文件各通道数据的前N个样点。. = wavread(filename,N1 N2) 仅返回文件各通道数据的第N1到第N2样点。Great！siz = wavread(filename,size) 返回文件所含音频数据的大小而不是实际的音频数据，返回矢量siz的格式为样点数 通道数。参见auread，wavwrite，wavplay，wavrecord附录4 specgram函数用途时频分析（产生声谱图）。 句法B = specgram(a)B = specgram(a,nfft)B,f = specgram(a,nfft,fs)B

32、,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两者较小的一个，即：min(nffs,256)l fs = 2l window是长度为nff

33、t的周期性汉宁窗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