大学本科语音信号处理实验讲义8学时

1、语音信号处理实验讲义时间：2011-12 目录实验一 语音信号生成模型分析.3实验二 语音信号时域特征分析.7实验三 语音信号频域特征分析12实验四 语音信号的同态处理和倒谱分析16实验一 语音信号生成模型分析一、 实验目的1、了解语音信号的生成机理，了解由声门产生的激励函数、由声道产生的调制函数和由嘴唇产生的辐射函数。2、编程实现声门激励波函数波形及频谱，与理论值进行比较。3、编程实现已知语音信号的语谱图，区分浊音信号和清音信号在语谱图上的差别。二、实验原理语音生成系统包含三部分：由声门产生的激励函数、由声道产生的调制函数和由嘴唇产生的辐射函数。语音生成系统的传递函数由这三个函数级联而成，即

2、 1、激励模型发浊音时，由于声门不断开启和关闭，产生间隙的脉冲。经仪器测试它类似于斜三角波的脉冲。也就是说，这时的激励波是一个以基音周期为周期的斜三角脉冲串。单个斜三角波的频谱表现出一个低通滤波器的特性。可以把它表示成z变换的全极点形式这里c是一个常数，T是脉冲持续时间。周期的三角波脉冲还得跟单位脉冲串的z变换相乘：这就是整个激励模型，是一个幅值因子。2、声道模型当声波通过声道时，受到声腔共振的影响，在某些频率附近形成谐振。反映在信号频谱图上，在谐振频率处其谱线包络产生峰值，把它称为共振峰。一个二阶谐振器的传输函数可以写成实践表明，用前3个共振峰代表一个元音足够了。对于较复杂的辅音或鼻音共振峰

3、要到5个以上。多个叠加可以得到声道的共振峰模型3、辐射模型从声道模型输出的是速度波，而语音信号是声压波。二者倒比称为辐射阻抗，它表征了口唇的辐射效应，可用下式表示：三、实验内容1、设声门脉冲单个三角波的数学表达式为试画出三角波波形图及其频谱。取。参考程序：%三角波及其频谱n=linspace(0,25,125);g=zeros(1,length(n);i=0;for i=0:40 if n(i+1)<=5 g(i+1)=0.5*(1-cos(n(i+1)*pi/5); else g(i+1)=cos(n(i+1)-5)*pi/8); endendfigure(1)subplot(1,2,

4、1)plot(n,g)xlabel('时间/ms')ylabel('幅度')axis(0,25,-0.4,1.2)r=fft(g,1024);r1=abs(r);yuanlai=20*log10(r1);signal(1:512)=yuanlai(1:512);pinlv=(0:1:511)*8000/1024;subplot(1,2,2)plot(pinlv,signal);xlabel('频率/Hz')ylabel('幅度/dB')axis(0,620,0,30) 图1-1 三角波及其频谱2、给出语音段“数字信号处理”（spe

5、ech.wav），画出它的语谱图。clear all;x,sr=wavread('speech_dsp.wav');s=length(x);w=round(44*sr/1000);n=w;shift=w/2;h=w-shift;%win=hanning(n)'win=hamming(n)'c=1;ncols=1+fix(s-n)/h);d=zeros(1+n/2),ncols);for b=0:h:(s-n) u=win'.*x(b+1):(b+n); t=fft(u); d(:,c)=t(1:(1+n/2); c=c+1;endtt=0:h:(s-n)

6、/sr;ff=0:(n-2)*sr/n;imagesc(tt,ff/1000,20*log10(abs(d);colormap(gray);axis xyxlabel('时间/s')ylabel('频率/kHz') 图1-2 语谱图四、思考题1、声门激励脉冲信号是高频衰减的还是高频增强的？2、画语谱图时为什么要给语音信号加汉明窗？若加矩形窗会有什么区别？3、在语谱图上观察，浊音信号的和清音信号的频谱有什么区别？实验二 语音信号时域特征分析一、实验目的1、了解自相关函数及自相关函数在语音信号处理中的应用。2、编写程序分析语音信号的短时自相关特征，计算语音信号的基音

7、周期。3、编写修正短时自相关函数的程序，并与未修正的函数进行比较。二、实验原理自相关函数用于衡量信号自身时间波形的相似性。由前面的讨论可知，清音和浊音的发声机理不同，因而在波形上也存在着较大的差异。浊音的时间波形呈现出一定的周期性，波形之间相似性较好；清音的时间波形呈现出随机噪声的特性，杂乱无章，样点间的相似性较差，这样，可以用短时自相关函数来测定语音的相似特性。时域离散确定信号的自相关函数定义为：对于语音信号来说，采用短时分析方法，可以定义短时自相关函数为因为，所以定义，则上式可以写成 如果长基音周期用窄的窗，将得不到预期的基音周期；但是如果短的基音周期用长的窗，自相关函数将对多个基因周期做

8、平均计算，从而模糊语音的短时特性，这是不希望的。为了解决这个问题，可以采用修正的短时自相关函数，选择的窗长不一定要等于自相关函数的最大自变量取值。这种方法可以采用较窄的窗，同时避免了短时自相关函数随k增加而衰减的不足。三、实验内容1、根据给出的浊音信号，分别画出浊音信号的时域波形、加矩形窗和加汉明窗后计算短时自相关归一化后的结果。语音的抽样频率为8kHz，窗长为320。参考程序：x,fs,nbits=wavread('speech_dsp.wav');s1=x(2500:2819);N=320;A=;for k=1:320 sum=0; for m=1:N-k+1 sum=su

9、m+s1(m)*s1(m+k-1); end A(k)=sumendfor k=1:320 A1(k)=A(k)/A(1);endf=zeros(1,320);n=1,j=1;while j<=320 f(1,j)=s1(n)*0.54-0.46*cos(2*pi*n/319); j=j+1; n=n+1;endB=;for k=1:320 sum=0; for m=1:N-k+1 sum=sum+s1(m)*s1(m+k-1); end B(k)=sumendfor k=1:320 B1(k)=B(k)/B(1);end%画图s2=s1/max(s1);figure(1)subplot

10、(3,1,1)plot(s2)title('一帧语音信号');xlabel('样点数');ylabel('幅度');axis(0,320,-1,1);subplot(3,1,2)plot(A1)title('加矩形窗的自相关函数')xlabel('延时k')ylabel('自相关函数R(k)')axis(0,320,-1,1);subplot(3,1,3)plot(B1)title('加汉明窗的自相关函数')xlabel('延时k')ylabel('自相关函数

11、R(k)')axis(0,320,-1,1); 图2-1 浊音信号加不同窗时的自相关函数2、仍选取上题中的语音信号，改变窗长和截取语音段的长度，计算修正的短时自相关函数。取值分别为：(1)N=320,M=640;(2)N=160,M=320;(3)N=70,M=140;x,fs,nbits=wavread('speech_dsp.wav');s1=x(2500:3139);b=s1;%窗长640，自相关运算取320个点。b1=b(1:640);N=320;A=;for k=1:320 sum=0; for m=1:N sum=sum+b1(m)*b1(m+k-1); e

12、nd A(k)=sum;endfor k=1:320 A1(k)=A(k)/A(1);end%画图figure(1)subplot(3,1,1)plot(A1);xlabel('延时k')ylabel('R(k)')legend('N=320')axis(0,320,-0.5,1)%窗长320，自相关运算取160个点。b2=b(1:320);N=160;A=;for k=1:160 sum=0; for m=1:N sum=sum+b2(m)*b2(m+k-1); end B(k)=sum;endfor k=1:160 B1(k)=B(k)/B(

13、1);end%画图figure(1)subplot(3,1,2)plot(B1);xlabel('延时k')ylabel('R(k)')legend('N=160')axis(0,320,-0.5,1)%窗长140，自相关运算取70个点。b3=b(1:140);N=70;A=;for k=1:70 sum=0; for m=1:N sum=sum+b3(m)*b3(m+k-1); end C(k)=sum;endfor k=1:70 C1(k)=C(k)/C(1);end%画图figure(1)subplot(3,1,3)plot(C1);xla

14、bel('延时k')ylabel('R(k)')legend('N=70')axis(0,320,-0.5,1) 图2-2 修正的自相关函数（参加自相关运算的点数N取不同值）四、思考题1、自相关函数的作用是什么？互相关函数的作用是什么？2、浊音信号分别加矩形窗和汉明窗时自相关函数有什么不同？3、清音信号的自相关函数和浊音信号的有什么区别？实验三 语音信号频域特征分析一、实验目的1、了解语音信号进行短时傅里叶分析的基础。2、画出浊音信号加不同窗时的时域波形和频谱，并进行比较。3、画出清音信号加不同窗时的时域波形和频谱，并与浊音信号进行比较。二、实验

15、原理语音信号本质上是非平稳信号，其非平稳特性是由发声器官的物理运动产生的。发声器官的运动由于存在惯性，所以可以假设语音信号在1030ms这样短的时间段内是平稳的，这是短时傅里叶分析的基础。某一帧的短时傅里叶变换的定义式如下：式中，是窗函数，不同的窗函数，可得到不同的傅里叶变换的结果。在式中，短时傅里叶变换有两个变量，即离散时间n及连续频率w，若令，则可得到离散的短时傅里叶变换如下： 对于窗来说，它除了具有选出序列中被分析部分的作用外，它的形状对时变傅里叶变换的特性也有重要作用。采用矩形窗时，基音谐波的各个峰都比较尖锐，而整个频谱图显得比较破碎，这是因为矩形窗的主瓣较窄，具有较高的频率分辨率，但

16、它也具有较高的旁瓣，因而使基音的相邻谐波之间的干扰比较严重。相邻谐波之间的这种“泄露”的现象，抵消了矩形窗主瓣窄的优点，因此，在语音短时频谱分析中极少采用矩形窗。当加汉明窗时，得到的短时频谱要平滑得多，因而在语音分析中汉明窗用得比较普遍。三、实验内容1、根据给出的浊音语音，画出它分别加矩形窗和汉明窗的时域波形和短时频谱。%函数文件fra.m，对语音信号进行帧长为frame，重叠为shift的分帧function output=fra(frame,shift,x)len=length(x);frame_num=floor(len-frame)/shift)+1;y=zeros(frame_num

17、,frame);for i=1:frame_num y(i,:)=x(shift*(i-1)+1):(shift*(i-1)+frame)'endoutput=y;%浊音的波形和短时频谱图(窗长256)x,fs,nbits=wavread('speech_dsp.wav');e=fra(256,128,x);ee=e(20,:);subplot(2,2,1)ee1=ee/max(ee);plot(ee1)xlabel('样点数')ylabel('幅度')title('原始语音')axis(0,256,-1.5,1.5)%矩

18、形窗傅立叶变换r=fft(ee,1024);r1=abs(r);r1=r1/max(r1);yuanlai=20*log10(r1);signal(1:512)=yuanlai(1:512);pinlv=(0:1:511)*8000/1024;subplot(2,2,2)plot(pinlv,signal);xlabel('频率/Hz')ylabel('对数幅度/dB')title('加矩形窗时语音谱')axis(0,4000,-80,15)%加汉明窗f=ee'.*hamming(length(ee);f1=f/max(f);subplo

19、t(2,2,3)plot(f1)axis(0,256,-1.5,1.5)xlabel('样点数')ylabel('幅度')title('窗选语音')%加汉明窗傅立叶变换r=fft(f,1024);r1=abs(r);r1=r1/max(r1);yuanlai=20*log10(r1);signal(1:512)=yuanlai(1:512);pinlv=(0:1:511)*8000/1024;subplot(2,2,4)plot(pinlv,signal);xlabel('频率/Hz')ylabel('对数幅度/dB

20、9;)title('加汉明窗时语音谱')axis(0,4000,-80,15) 图3-1 浊音信号加不同窗时的频谱2、根据给出的清音语音，画出它分别加矩形窗和汉明窗的时域波形和短时频谱。%清音的波形和短时频谱图(窗长256)x,fs,nbits=wavread('speech_dsp.wav');e=fra(256,128,x);ee=e(10,:);subplot(2,2,1)ee1=ee/max(ee);plot(ee1)xlabel('样点数')ylabel('幅度')title('清音原始语音')axis(

21、0,256,-1.5,1.5)%清音加矩形窗傅立叶变换r=fft(ee,1024);r1=abs(r);r1=r1/max(r1);yuanlai=20*log10(r1);signal(1:512)=yuanlai(1:512);pinlv=(0:1:511)*8000/1024;subplot(2,2,2)plot(pinlv,signal);xlabel('频率/Hz')ylabel('对数幅度/dB')title('清音加矩形窗时语音谱')axis(0,4000,-80,15)%清音加汉明窗f=ee'.*hamming(lengt

22、h(ee);f1=f/max(f);subplot(2,2,3)plot(f1)axis(0,256,-1.5,1.5)xlabel('样点数')ylabel('幅度')title('清音窗选语音')%清音加汉明窗傅立叶变换r=fft(f,1024);r1=abs(r);r1=r1/max(r1);yuanlai=20*log10(r1);signal(1:512)=yuanlai(1:512);pinlv=(0:1:511)*8000/1024;subplot(2,2,4)plot(pinlv,signal);xlabel('频率/Hz

23、')ylabel('对数幅度/dB')title('清音加汉明窗时语音谱')axis(0,4000,-80,15)图3-1 清音信号加不同窗时的频谱四、思考题1、浊音信号分别加矩形窗和汉明窗时的频谱有什么区别？2、清音信号的频谱和浊音信号相比有什么不同?3、对语音信号进行短时傅里叶分析的基础是什么？实验四 语音信号的同态处理和倒谱分析一、实验目的1、掌握同态处理的原理，倒谱与复倒谱的含义。2、利用同态处理对语音信号进行基音检测。3、利用同态处理对语音信号进行共振峰检测。二、实验原理1、基音检测语音的倒谱是将语音的短时谱取对数后再进行IDFT得到的，所以浊

24、音信号的周期性激励反映在倒谱上是同样周期的冲激。借此，可从倒谱波形中估计出基音周期。一般把倒谱波形中第二个冲激，认为是对应激励源的基频。下面给出一种倒谱法求基音周期的框图。先计算倒谱，然后在预期的基音周期附近寻找峰值。如果倒谱的峰值超出了预先规定的门限，则输入语音段定为浊音，而峰的位置就是基音周期的良好估值。如果没有超出门限的峰值，则输入语音段定为清音。如果计算的是一个时变的倒谱，则可估计出激励源模型及基音周期随时间的变化。一般每隔1020ms计算一次倒谱，这是因为在一般语音中激励参数是缓慢变化的。DFTIDFT峰值检测 基音检测 s(n)w(n) 图4.1 一种倒谱法求基音周期的实现框图2、

25、共振峰检测倒谱将基音谐波和声道的频谱包络分离开来。倒谱的低时部分可以分析声道、声门和辐射信息，而高频部分可用来分析激励源信息。对倒谱进行低时窗选，通过语音倒谱分析系统的最后一级，进行DFT后的输出即为平滑后的对数模函数，这个平滑得对数谱显示了特定输入语音段的谐振结构，即谱的峰值基本上对应于共振峰频率，对平滑过的对数谱中的峰值进行定位，即可估计共振峰。原理框图如图4.2所示。分帧加窗|FFT|logIFFT加窗|FFT|平滑处理共振峰log图4.2 共振峰检测框图三、实验内容1、利用同态处理计算给出语音信号的基音周期和基频。y,fs,nbits=wavread('speech_dsp.w

26、av');time1=1:length(y);time=(1:length(y)/fs;frameSize=floor(30*fs/1000);startIndex=round(2080);endIndex=startIndex+frameSize-1;frame=y(startIndex:endIndex);ylen=frameSize;frameSize=length(frame);frame2=frame.*hamming(length(frame);rwy=rceps(frame2);cepstrum=rwy(1:ylen/2);for i=1:ylen/2 cepstrum1

27、(i)=rwy(ylen/2+1-i);endfor i=(ylen/2+1):ylen cepstrum1(i)=rwy(i+1-ylen/2);end%基音检测LF=floor(fs/500);HF=floor(fs/70);cn=cepstrum(LF:HF);mx_cep ind=max(cn);if mx_cep>0.08 a=fs/(LF+ind);else a=0;endpitch=a%画图figure(1);subplot(3,1,1);plot(time1,y);title('语音波形');axis tightylim=get(gca,'ylim

28、');line(time1(startIndex),time1(startIndex),ylim,'color','r');line(time1(endIndex),time1(endIndex),ylim,'color','r');xlabel('样点数');ylabel('幅度');subplot(3,1,2);plot(frame);axis(0,240,-0.5,0.5)title('一帧语音');xlabel('样点数');ylabel('幅

29、度');subplot(3,1,3);time2=-119:1:-1,0:1:120;plot(time2,cepstrum1);axis(-120,120,-0.5,0.5)title('一帧语音的倒谱');xlabel('样点数');ylabel('幅度'); 图4-1 对浊音信号进行基音检测2、利用同态处理计算给出语音信号的共振峰。y,fs,nbits=wavread('speech_dsp.wav');time1=1:length(y);time=(1:length(y)/fs;frameSize=floor(30*fs/1000);startIndex=round(2080);endIndex=startIndex+frameSize-1;frame=y(startIndex:endIndex);%frameSize=length(frame);frame2=frame.*hamming(length(frame);rwy=rceps(frame2);ylen=length(rwy);cepstrum=rwy(1:ylen/2);%基音检测LF=floor(fs/500);HF=floor(fs/70);cn=cepstrum(LF:HF);mx_cep,ind=ma