哈尔滨工程大学 语音信号处理实验报告

1、实 验 报 告实验课程名称： 语音信号处理实验 姓名： 班级： 20120811 学号： 实验序号实验名称实验过程实验结果实验成绩实验一语音信号的端点检测实验二语音信号的特征提取实验三语音信号的基频提取指导教师 张磊 实验教室 21B#293 实验时间 2015年4月12日 实验成绩 实验一 语音信号的端点检测一、实验目的1、掌握短时能量的求解方法2、掌握短时平均过零率的求解方法3、掌握利用短时平均过零率和短时能量等特征，对输入的语音信号进行端点检测。二、实验设备HP计算机、Matlab软件三、实验原理1、短时能量语音信号的短时能量分析给出了反应这些幅度变化的一个合适的描述方法。对于信号，短时

2、能量的定义如下：2、短时平均过零率短时平均过零率是指每帧内信号通过零值的次数。对于连续语音信号，可以考察其时域波形通过时间轴的情况。对于离散信号，实质上就是信号采样点符号变化的次数。过零率在一定程度上可以反映出频率的信息。短时平均过零率的公式为： 其中，是符号函数，即3、端点检测原理能够实现这些判决的依据在于，不同性质语音的各种短时参数具有不同的概率密度函数，以及相邻的若干帧语音应具有一致的语音特性，它们不会在S、U、V之间随机地跳来跳去。要正确判断每个输入语音的起点和终点，利用短时平均幅度参数和短时平均过零率可以做到这一点。首先，根据浊音情况下的短时能量参数的概率密度函数确定一个阈值参数，值

3、一般定的较高。当一帧输入信号的短时平均幅度参数超过时，就可以判定该帧语音信号不是无声，而有相当大的可能是浊音。根据可判定输入语音的前后两个点和。在和之间的部分肯定是语音段，但语音的精确起点、终点还要在之前和之后仔细查找，如图1所示。图1 端点检测原理示意图为此，再设定一个较低的阈值参数，由点向前找，当短时能量由大到小减至可以确定点。类似地，可以由点向后找，确定点。在和之间的仍能肯定是语音段。然后由向前和向后，利用短时平均过零率进行搜索。根据无声情况下的短时平均过零率，设置一个参数，如果由向前搜索时，短时平均过零率大于的3倍，则认为这些信号仍属于语音段，直到短时平均过零率下降到低于3倍的，这时的

4、点就是语音的精确的起点。对于终点也做类似的处理，可以确定终点。采用短时平均过零率的原因在于，点以前可能是一段清辅音，它的能量相当弱，依靠能量不可能将它们与无声段分开。而对于清辅音来说，它们的过零率明显高于无声段，因而能用这个参数将二者区分开来。存在噪声的情况下，正常的过零率计算会存在一定的误差，解决这个问题的一种方法是对上述过零率定义做一个简单的修改，即设立一个门限T，将过零率的含义修改为跨过正负门限的次数，如图2。图2 门限短时平均过零率于是，有： 这样计算的短时平均过零率就有一定的抗干扰能力。即使存在小的随机噪声，只要它不超过正负门限所构成的带，就不会产生虚假过零率。四、实验步骤及程序实验

5、程序：function point_check() s1=wavread('man.wav'); s=s1/max(abs(s1); s=filter(1 -0.9375,1,s);plot(s); len=length(s); round=160;repeat=80;inc=round-repeat;%b_len-2*repeat; ka=ceil(len-round)/(round-repeat)+1;s=s;zeros(round-repeat)*(ka-1)+round-len,1); len=length(s); w=zeros(round,ka); for i=1:

6、ka for k=1:round w(k,i)=s(k+(i-1)*(round-repeat),1); end end ST=0.01;F=zeros(ka,round); for i=1:ka for k=0:round-1 for j=1:round-k F(ka,k+1)=F(ka,k+1)+abs(w(j,i)-w(j+k,i); end if(k<round-1) if(w(k+1,i)>0&&w(k+2,i)<0) z=sign(w(k+1,i)-ST)-sign(w(k+2,i)+ST); elseif(w(k+1,i)<0&&a

7、mp;w(k+2,i)>0) z=sign(w(k+2,i)-ST)-sign(w(k+1,i)+ST); end if(z=2) f(i,1)=f(i,1)+1; end end end end amp=sum(abs(w).2,1);实验结果截图：五、实验结果与分析 这次实验利用MATLAB对信号进行分析和处理，利用短时过零率和短时能量，对语音信号的端点进行检测。我之前虽然接触过MATLAB软件，但从未进行过有关语音信号的操作，在实验过程中欠缺独立性，代码是参考网上的，在理解的基础上借鉴。学习是无止境的，现在的动手能力不强也是由于之前的松懈造成的，以后一定要多动手，毕竟能力是立身之本

8、。实验二 语音信号的特征提取一、实验目的1、掌握语音信号的Mel倒谱特征（MFCC）的求解方法2、掌握语音信号的线性预测原理以及LPC特征的求解方法二、实验设备HP计算机、Matlab软件三、实验原理1、MFCC特征原理及求解方法在语音识别和说话人识别中，常用的语音特征是基于Mel频率的、倒谱系数（Mel Frequency Cepstrum Coefficient，简称MFCC）。由于MFCC参数是将人耳的听觉感知特性和语音的产生机制相结合，因此目前大多数语音识别系统中广泛使用这种特征。研究者根据心理学实验得到了类似于耳蜗作用的一组滤波器组，这就是Mel频率滤波器组。Mel频率可以用公式表达

9、如下：在实际应用中，MFCC倒谱系数计算过程如下：（1） 将信号进行分帧，预加重和加汉明窗处理，然后进行短时傅立叶变换得到其频谱；（2） 求出频谱平方，即能量谱，并用M个Mel带通滤波器进行滤波，由于每一个频带中分量的作用在人耳中是叠加的，因此将每个滤波频带内的能量进行叠加，这时第k个滤波器输出功率谱；（3） 将每个滤波器的输出取对数，得到相应频带的对数功率谱；并进行反离散余弦变换，得到个MFCC系数，如下式所示。一般取1216个左右； （4） 这种直接得到的MFCC特征作为静态特征，将这种静态特征做一阶和二阶差分，得到相应的动态特征。2、线性预测原理及LPC特征求解方法根据参数模型功率谱估计

10、的思想，可以将语音信号看作是由一个输入序列激励一个全极点的系统（模型）而产生的输出，如图3所示。系统的传递函数为：其中为常数，为实数，p为模型的阶数。显而易见，这种模型是以系数和增益G为模型参数的全极点模型，即AR模型。 用系数可以定义一个阶线性预测器这个阶预测器从时域角度可以理解为，用信号的前个样本来预测当前的样本得到预测值，因为预测器是用AR模型的系数来构造的，而AR模型是在最小均方意义上对数据的拟合，所以预测器必然是一个最佳预测器，即此时预测器的预测误差短时总能量最小。语音信号的线性预测分析就是根据这一性质，从语音信号出发，依据最小均方误差准则，估计出一组线性预测器的系数，它就是我们所要

11、求的信号AR模型的系数。被称为线性预测系数或LPC系数。预测器的预测误差为：由上式可知，是输入为，且具有如下形式传递函数的滤波器的输出因此称为预测误差滤波器。比较式(4-1)和式(4-5)可知，即预测误差滤波器是系统的逆滤波器。为了在最小均方误差意义上计算一组最佳预测系数，定义短时预测均方误差为由于语音信号的时变特性，线性预测分析应该在短时的语音段上进行，即按帧进行。因此上式的求和通常也是在一帧语音的范围内进行。为了使达到最小，必须满足。则有：这样可以得到以为变量的线性方程组, 要构造信号的AR模型，还应估算增益因子。将系统传递函数转化成差分方程的形式：则预测误差和残差能量分别计算如下：激励信

12、号无法准确计算，但根据前文所述的语音产生模型，在浊音情况下，激励可以看作是准周期的脉冲串；在清音时，可以看作为高斯白噪声。因此输入信号总能量可以认为近似为1，这样估算增益因子 。这种计算只是一种近似的方法，。为了使模型的假定能够很好地符合语音产生的模型，需要考虑如下两个因素。（1）模型阶数的选择阶数应与共振峰的个数相吻合，通常一对极点对应一个共振峰，因此当共振峰数为5时，应取=10。考虑到一些情况，一般按如下的方式计算模型阶数，其中为共振峰的个数。（2）通过预加重进行高频提升由于声门脉冲形状和口唇辐射的影响，语音信号的频谱在总趋势上会产生高频衰落的现象，大约每倍程下降6dB。要抵消这一影响，通

13、常在LPC分析之前采用一个非常简单的一阶FIR滤波器进行预加重以进行高频提升，其传递函数为，其中为预加重系数，对10kHz采样的语音，一般取=0.95。大量的实践证明：LPC参数是反映语音信号特征的良好参数。附：莱文逊杜宾递推算法该算法过程如下：（1）计算自相关系数；（2）；（3）；（4）开始按如下公式进行递推运算： ， （5） 。若则算法结束退出，否则返回第（4）步。 四、实验步骤及程序MFCC特征提取程序：close all clear clc x=wavread('N:secondman.wav'); bank=melbankm(24,256,8000,0,0.5,

14、9;m'); bank=full(bank); bank=bank/max(bank(:); for k=1:12 n=0:23; dctcoef(k,:)=cos(2*n+1)*k*pi/(2*24); end w = 1 + 6 * sin(pi * 1:12 ./ 12); w = w/max(w);% 预加重滤波器xx=double(x); xx=filter(1 -0.9375,1,xx); % 语音信号分帧xx=enframe(xx,256,80); % 计算每帧的MFCC参数for i=1:size(xx,1) y = xx(i,:); s = y' .* ham

15、ming(256); t = abs(fft(s); t = t.2; c1=dctcoef * log(bank * t(1:129); c2 = c1.*w'm(i,:)=c2' end figure imagesc(m);xlabel('帧数');ylabel('幅度');title('MFCC'); 实验结果截图：LPC特征提取程序：clear clcI = wavread('N:secondman.wav');%读入原始语音%subplot(3,1,1), plot(I); title('原始语音

16、波形') %对指定帧位置进行加窗处理Q = I' N = 256; % 窗长Hamm = hamming(N); % 加窗frame = 60;%需要处理的帧位置M = Q(frame - 1) * (N / 2) + 1):(frame - 1) * (N / 2) + N); Frame = M .* Hamm'%加窗后的语音帧B,F,T = specgram(I,N,N/2,N); m,n = size(B); for i = 1:m FTframe1(i) = B(i,frame); end P =50;%input('请输入预测器阶数 = ')

17、; ai = lpc(Frame,P); % 计算lpc系数LP = filter(0 -ai(2:end),1,Frame); % 建立语音帧的正则方程FFTlp = fft(LP); E = Frame - LP; % 预测误差figuresubplot(2,1,1),plot(1:N,Frame,1:N,LP,'-r');grid; title('原始语音和预测语音波形') subplot(2,1,2),plot(E);grid; title('预测误差'); figurefLength(1 : 2 * N) = M,zeros(1,N);

18、 Xm = fft(fLength,2 * N); X = Xm .* conj(Xm); Y = fft(X , 2 * N); Rk = Y(1 : N); PART = sum(ai(2 : P + 1) .* Rk(1 : P); G = sqrt(sum(Frame.2) - PART); A = (FTframe1 - FFTlp(1 : length(F') ./ FTframe1 subplot(2,1,1),plot(F',20*log(abs(FTframe1),F',(20*log(abs(1 ./ A),'-r');grid;xl

19、abel('频率/dB');ylabel('幅度'); title('短时谱'); subplot(2,1,2),plot(F',(20*log(abs(G ./ A);grid; xlabel('频率/dB');ylabel('幅度'); title('LPC谱'); figure%求出预测误差的倒谱pitch = fftshift(rceps(E); M_pitch = fftshift(rceps(Frame); subplot(2,1,1),plot(M_pitch);grid; x

20、label('语音帧');ylabel('/dB'); title('原始语音帧倒谱'); subplot(2,1,2),plot(pitch);grid; xlabel('语音帧');ylabel('/dB'); title('预测误差倒谱'); figure %画出语谱图ai1 = lpc(I,P); % 计算原始语音lpc系数LP1 = filter(0 -ai(2:end),1,I); % 建立原始语音的正则方程subplot(2,1,1); specgram(I,N,N/2,N); tit

21、le('原始语音语谱图'); subplot(2,1,2); specgram(LP1,N,N/2,N); title('预测语音语谱图');实验结果截图： 五、实验结果与分析这次实验利用MATLAB对信号进行分析和处理，主要就是进行语音线性预测，语音线性预测的基本思想是：一个语音信号的抽样值可以用过去若干个取样值的线性组合来逼近。通过使实际语音抽样值与线性预测抽样值的均方误差达到最小，可以确定唯一的一组线性预测系数。做了LPC特征提取的实验之后，促进了我对MFCC特征提取进行探究和验证，此次试验让我对相关知识加深了理解。实验三 语音信号的基频提取一、实验目的1

22、、掌握语音信号基频的概念，加深对基频刻画声调特征作用的理解。2、掌握语音信号基频特征的典型求解方法二、实验设备HP计算机、Matlab软件三、实验原理1、基频概念的理解基音检测的方法大致可分为三类： 波形估计法。直接由语音波形来估计基音周期，分析出波形上的周期峰值。包括并行处理法、数据减少法等； 相关处理法。这种方法在语音信号处理中广泛使用，这是因为相关处理法抗波形的相位失真能力强，另外它在硬件处理上结构简单。包括波形自相关法、平均振幅差分函数法（AMDF）、简化逆滤波法（SIFT）等； 变换法。将语音信号变换到频域或倒谱域来估计基音周期，利用同态分析方法将声道的影响消除，得到属于激励部分的信

23、息，进一步求取基音周期，比如倒谱法。虽然倒谱分析算法比较复杂，但基音估计效果较好。这里我们只考虑自相关方法的基频提取。2、自相关方法的基频提取浊音信号的自相关函数在基音周期的整数倍位置上出现峰值，而清音的自相关函数没有明显的峰值出现，因此检测自相关函数是否有峰值就可以判断是清音或浊音，峰峰值之间对应的就是基音周期。为了提高自相关方法检测基音周期的准确性，需要进行一些前期的预处理。1）预处理语音信号的低幅值部分包含大量的共振峰信息，而高幅值部分包含较多的基音信息。因此，任何削减或者抑制语音低幅度部分的非线性处理都会使自相关方法的性能得到改善。中心削波即是一种非线性处理，用以消除语音信号的低幅度部

24、分，其削波特性如图4所示，数学表达形式为其中削波电平由语音信号的峰值幅度来确定，它等于语音段最大幅度的一个固定百分数，一般取最大信号幅度的60%70%。这个门限的选择是重要的，一般在不损失基音信息的情况下应尽可能选得高些，以达到较好的效果。经过中心削波后只保留了超过削波电平的部分，其结果是削去了许多和声道响应有关的波动。中心削波后的语音再计算自相关函数，这样在基音周期位置呈现大而尖的峰值，而其余的次要峰值幅度都很小。据报道使用这种方法，对电话带宽的语音在信噪比低至18dB的情况下获得了良好的性能。图4 中心削波函数计算自相关函数的运算量是很大的，其原因是传统的计算机进行乘法运算非常费时。尽管近

25、年来随着数字信号处理器的广泛使用，实时地计算自相关函数已经不是问题，但在基音检测中仍然有一些减少短时自相关运算的有效方法。如可对中心削波函数进行修正，采用三电平中心削波的方法，如图5所示。其削波函数为即削波器的输出在时为1，时为，除此以外均为零。虽然这一处理会增加刚刚超过削波电平峰的重要性，但大多数次要的峰被滤除掉了，而只保留了明显的周期性峰。图5 三电平削波函数此外，还可以用一个通带为900Hz的线性相位低通滤波器滤除高次谐波分量。这样处理后的信号，基本上只含有第一共振峰以下的基波和谐波分量。实验表明，用这种方法做预处理，对改善自相关和平均幅度差函数法的基音检测都有明显的效果。2）基于自相关

26、函数的基音检测短时自相关函数在基音周期的各个整数倍点上有很大的峰值，只要找到第一最大峰值点的位置，并计算它与原点的间隔，便能估计出基音周期。但实际上并不是这么简单，第一个最大峰值点的位置有时不能与基音周期相吻合。产生这种情况的原因有以下两个方面。第一方面是与窗的长度有关。一般认为窗长至少应大于两个基音周期，才可能有较好的效果。第二方面与声道特性的影响有关。有的情况下，即使窗长已经选得足够长，第一个最大峰值点与基音周期仍不一致，这就是声道共振峰特性的干扰。经过上述带通滤波的预处理，就可以消除大部分的共振峰的影响。但是，如果希望减少自相关计算中的乘法运算，可以把上述中心削波后的信号的自相关用两个信

27、号的互相关代替，其中一个信号是，令一个信号是对进行三电平量化产生的结果。显然，只有-1，0，+1三种可能的取值，因而这里的互相关计算只需做加减法，而这个互相关序列的周期性与的自相关序列近似相同。下面结合L.R.Rabiner一篇论文中介绍的具体例子来叙述关于自相关函数的基音检测方法。假设信号的采样率为10kHz，窗序列采用300点的矩形窗，帧叠200点。这时对每一帧进行基音周期估计的步骤如下：（1） 用900Hz低通滤波器对一帧语音信号进行滤波，并去掉开头的20个输出值不用，得到；（2） 分别求的前部100个样点和后部100个样点的最大幅度，并取其中较小的一个，乘以因子0.68作为门限电平；（

28、3） 对分别进行中心削波得到和三电平量化得到；（4） 求这两个信号的互相关值。其中，此处的取值范围20150相应于基音频率范围60Hz500Hz，相应于短时能量；得到互相关值后，可以得到中的最大值，如果，则认为本帧为清音，令其基音周期值为0，否则基音周期即为使为最大值时位置的值，即。3）后处理无论采用哪一种算法求得的基音周期轨迹与真实的基音周期轨迹不可能完全吻合。大部分段落是吻合的，而在一些局部段落或区域中有一个或几个基音周期的估计值偏离了正常的轨迹这时称为基音轨迹产生了若干“野点”。为了去除这些野点，可以采用平滑算法，其中最常用的是中值平滑算法和线性平滑算法。在中值滤波平滑算法中，被平滑点的

29、左右各取个样点，连同被平滑点共同构成一组个信号样点值。将这些样点值按大小次序排成一队，取此队列中间者作为平滑器的输出。值一般取为1或2，即中值平滑的“窗口”一般套住3或5个样值。中值平滑的优点是既可以有效地去除少量野点，又不会破坏基音周期轨迹中的两个平滑段之间的阶跃性变化。线性平滑是用滑动窗进行线性滤波处理，即：其中为点平滑窗，满足：例如三点窗的值可取为。线性平滑在纠正输入信号中不平滑处样点的同时，也使附近的样点值做了修改，所以窗长不易过长。四、实验步骤及程序实验程序：x=wavread('man.wav');%读取声音文件 figure; stem(x,'.'

30、); %显示声音信号的波形 n=160; %取20ms的声音片段，即160个样点 for m=1:length(x)/n; %对每一帧求短时自相关函数 for k=1:n; Rm(k)=0; for i=(k+1):n; Rm(k)=Rm(k)+x(i+(m-1)*n)*x(i-k+(m-1)*n); end end p=Rm(10:n); %防止误判，去掉前边10个数值较大的点 Rmax,N(m)=max(p); %读取第一个自相关函数的最大点 end %补回前边去掉的10个点 N=N+10; T=N/8; %算出对应的周期 figure;stem(T,'.');axis(0

31、 length(T) 0 10); xlabel('帧数(n)');ylabel('周期(ms)');title('各帧基音周期'); T1= medfilt1(T,5); %去除野点 figure;stem(T1,'.');axis(0 length(T1) 0 10); xlabel('帧数(n)');ylabel('周期(ms)');title('各帧基音周期');a=wavread('1.wav');L=length(a);m=max(a);n=min(a);ht=(m+n)/2;for i=1:L; a(i)=a(i)-ht;endfor i=1:L; a(i)=a(i)/m;endN=320;A=;for k=1:320; sum=0; fo