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

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

《大学本科语音信号处理实验讲义8学时.docx》由会员分享，可在线阅读，更多相关《大学本科语音信号处理实验讲义8学时.docx（26页珍藏版）》请在冰豆网上搜索。

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

语音信号处理实验讲义

时间：

2011-12

实验一语音信号生成模型分析……………………………………….3

实验二语音信号时域特征分析……………………………………….7

实验三语音信号频域特征分析………………………………………12

实验四语音信号的同态处理和倒谱分析……………………………16

实验一语音信号生成模型分析

一、实验目的

1、了解语音信号的生成机理，了解由声门产生的激励函数、由声道产生的调制函数和由嘴唇产生的辐射函数。

2、编程实现声门激励波函数波形及频谱，与理论值进行比较。

3、编程实现已知语音信号的语谱图，区分浊音信号和清音信号在语谱图上的差别。

二、实验原理

语音生成系统包含三部分：

由声门产生的激励函数

、由声道产生的调制函数

和由嘴唇产生的辐射函数

。

语音生成系统的传递函数由这三个函数级联而成，即

1、激励模型

发浊音时，由于声门不断开启和关闭，产生间隙的脉冲。

经仪器测试它类似于斜三角波的脉冲。

也就是说，这时的激励波是一个以基音周期为周期的斜三角脉冲串。

单个斜三角波的频谱表现出一个低通滤波器的特性。

可以把它表示成z变换的全极点形式

这里c是一个常数，T是脉冲持续时间。

周期的三角波脉冲还得跟单位脉冲串的z变换相乘：

这就是整个激励模型，

是一个幅值因子。

2、声道模型

当声波通过声道时，受到声腔共振的影响，在某些频率附近形成谐振。

反映在信号频谱图上，在谐振频率处其谱线包络产生峰值，把它称为共振峰。

一个二阶谐振器的传输函数可以写成

实践表明，用前3个共振峰代表一个元音足够了。

对于较复杂的辅音或鼻音共振峰要到5个以上。

多个

叠加可以得到声道的共振峰模型

3、辐射模型

从声道模型输出的是速度波，而语音信号是声压波。

二者倒比称为辐射阻抗，它表征了口唇的辐射效应，可用下式表示：

三、实验内容

1、设声门脉冲单个三角波的数学表达式为

试画出三角波波形图及其频谱。

取

。

参考程序：

%三角波及其频谱

n=linspace（0,25,125）;

g=zeros（1,length（n））;

i=0;

fori=0:

40

ifn（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）;

end

end

figure

（1）

subplot（1,2,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、给出语音段“数字信号处理”（speech.wav），画出它的语谱图。

clearall;

[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）;

forb=0:

h:

（s-n）

u=win'.*x（（b+1）:

（b+n））;

t=fft（u）;

d（:

c）=t（1:

（1+n/2））;

c=c+1;

end

tt=[0:

h:

（s-n）]/sr;

ff=[0:

（n-2）]*sr/n;

imagesc（tt,ff/1000,20*log10（abs（d）））;

colormap（gray）;

axisxy

xlabel（'时间/s'）

ylabel（'频率/kHz'）

图1-2语谱图

四、思考题

1、声门激励脉冲信号是高频衰减的还是高频增强的？

2、画语谱图时为什么要给语音信号加汉明窗？

若加矩形窗会有什么区别？

3、在语谱图上观察，浊音信号的和清音信号的频谱有什么区别？

实验二语音信号时域特征分析

一、实验目的

1、了解自相关函数及自相关函数在语音信号处理中的应用。

2、编写程序分析语音信号的短时自相关特征，计算语音信号的基音周期。

3、编写修正短时自相关函数的程序，并与未修正的函数进行比较。

二、实验原理

自相关函数用于衡量信号自身时间波形的相似性。

由前面的讨论可知，清音和浊音的发声机理不同，因而在波形上也存在着较大的差异。

浊音的时间波形呈现出一定的周期性，波形之间相似性较好；清音的时间波形呈现出随机噪声的特性，杂乱无章，样点间的相似性较差，这样，可以用短时自相关函数来测定语音的相似特性。

时域离散确定信号的自相关函数定义为：

对于语音信号来说，采用短时分析方法，可以定义短时自相关函数为

因为

，所以

定义

，则上式可以写成

如果长基音周期用窄的窗，将得不到预期的基音周期；但是如果短的基音周期用长的窗，自相关函数将对多个基因周期做平均计算，从而模糊语音的短时特性，这是不希望的。

为了解决这个问题，可以采用修正的短时自相关函数，选择的窗长不一定要等于自相关函数的最大自变量取值。

这种方法可以采用较窄的窗，同时避免了短时自相关函数随k增加而衰减的不足。

三、实验内容

1、根据给出的浊音信号，分别画出浊音信号的时域波形、加矩形窗和加汉明窗后计算短时自相关归一化后的结果。

语音的抽样频率为8kHz，窗长为320。

参考程序：

[x,fs,nbits]=wavread（'speech_dsp.wav'）;

s1=x（2500:

2819）;

N=320;

A=[];

fork=1:

320

sum=0;

form=1:

N-k+1

sum=sum+s1（m）*s1（m+k-1）;

end

A（k）=sum

end

fork=1:

320

A1（k）=A（k）/A

（1）;

end

f=zeros（1,320）;

n=1,j=1;

whilej<=320

f（1,j）=s1（n）*[0.54-0.46*cos（2*pi*n/319）];

j=j+1;

n=n+1;

end

B=[];

fork=1:

320

sum=0;

form=1:

N-k+1

sum=sum+s1（m）*s1（m+k-1）;

end

B（k）=sum

end

fork=1:

320

B1（k）=B（k）/B

（1）;

end

%画图

s2=s1/max（s1）;

figure

（1）

subplot（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（'自相关函数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=[];

fork=1:

320

sum=0;

form=1:

N

sum=sum+b1（m）*b1（m+k-1）;

end

A（k）=sum;

end

fork=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=[];

fork=1:

160

sum=0;

form=1:

N

sum=sum+b2（m）*b2（m+k-1）;

end

B（k）=sum;

end

fork=1:

160

B1（k）=B（k）/B

（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=[];

fork=1:

70

sum=0;

form=1:

N

sum=sum+b3（m）*b3（m+k-1）;

end

C（k）=sum;

end

fork=1:

70

C1（k）=C（k）/C

（1）;

end

%画图

figure

（1）

subplot（3,1,3）

plot（C1）;

xlabel（'延时k'）

ylabel（'R（k）'）

legend（'N=70'）

axis（[0,320,-0.5,1]）

图2-2修正的自相关函数（参加自相关运算的点数N取不同值）

四、思考题

1、自相关函数的作用是什么？

互相关函数的作用是什么？

2、浊音信号分别加矩形窗和汉明窗时自相关函数有什么不同？

3、清音信号的自相关函数和浊音信号的有什么区别？

实验三语音信号频域特征分析

一、实验目的

1、了解语音信号进行短时傅里叶分析的基础。

2、画出浊音信号加不同窗时的时域波形和频谱，并进行比较。

3、画出清音信号加不同窗时的时域波形和频谱，并与浊音信号进行比较。

二、实验原理

语音信号本质上是非平稳信号，其非平稳特性是由发声器官的物理运动产生的。

发声器官的运动由于存在惯性，所以可以假设语音信号在10~30ms这样短的时间段内是平稳的，这是短时傅里叶分析的基础。

某一帧的短时傅里叶变换的定义式如下：

式中，

是窗函数，不同的窗函数，可得到不同的傅里叶变换的结果。

在式中，短时傅里叶变换有两个变量，即离散时间n及连续频率w，若令

，则可得到离散的短时傅里叶变换如下：

对于

窗来说，它除了具有选出

序列中被分析部分的作用外，它的形状对时变傅里叶变换的特性也有重要作用。

采用矩形窗时，基音谐波的各个峰都比较尖锐，而整个频谱图显得比较破碎，这是因为矩形窗的主瓣较窄，具有较高的频率分辨率，但它也具有较高的旁瓣，因而使基音的相邻谐波之间的干扰比较严重。

相邻谐波之间的这种“泄露”的现象，抵消了矩形窗主瓣窄的优点，因此，在语音短时频谱分析中极少采用矩形窗。

当加汉明窗时，得到的短时频谱要平滑得多，因而在语音分析中汉明窗用得比较普遍。

三、实验内容

1、根据给出的浊音语音，画出它分别加矩形窗和汉明窗的时域波形和短时频谱。

%函数文件fra.m，对语音信号进行帧长为frame，重叠为shift的分帧

functionoutput=fra（frame,shift,x）

len=length（x）;

frame_num=floor（（len-frame）/shift）+1;

y=zeros（frame_num,frame）;

fori=1:

frame_num

y（i,:

）=x（（shift*（i-1）+1）:

（shift*（i-1）+frame））';

end

output=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]）

%矩形窗傅立叶变换

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）;

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'）

ylabel（'对数幅度/dB'）

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（[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（length（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'）

ylabel（'对数幅度/dB'）

title（'清音加汉明窗时语音谱'）

axis（[0,4000,-80,15]）

图3-1清音信号加不同窗时的频谱

四、思考题

1、浊音信号分别加矩形窗和汉明窗时的频谱有什么区别？

2、清音信号的频谱和浊音信号相比有什么不同?

3、对语音信号进行短时傅里叶分析的基础是什么？

实验四语音信号的同态处理和倒谱分析

一、实验目的

1、掌握同态处理的原理，倒谱与复倒谱的含义。

2、利用同态处理对语音信号进行基音检测。

3、利用同态处理对语音信号进行共振峰检测。

二、实验原理

1、基音检测

语音的倒谱是将语音的短时谱取对数后再进行IDFT得到的，所以浊音信号的周期性激励反映在倒谱上是同样周期的冲激。

借此，可从倒谱波形中估计出基音周期。

一般把倒谱波形中第二个冲激，认为是对应激励源的基频。

下面给出一种倒谱法求基音周期的框图。

先计算倒谱，然后在预期的基音周期附近寻找峰值。

如果倒谱的峰值超出了预先规定的门限，则输入语音段定为浊音，而峰的位置就是基音周期的良好估值。

如果没有超出门限的峰值，则输入语音段定为清音。

如果计算的是一个时变的倒谱，则可估计出激励源模型及基音周期随时间的变化。

一般每隔10~20ms计算一次倒谱，这是因为在一般语音中激励参数是缓慢变化的。

图4.1一种倒谱法求基音周期的实现框图

2、共振峰检测

倒谱将基音谐波和声道的频谱包络分离开来。

倒谱的低时部分可以分析声道、声门和辐射信息，而高频部分可用来分析激励源信息。

对倒谱进行低时窗选，通过语音倒谱分析系统的最后一级，进行DFT后的输出即为平滑后的对数模函数，这个平滑得对数谱显示了特定输入语音段的谐振结构，即谱的峰值基本上对应于共振峰频率，对平滑过的对数谱中的峰值进行定位，即可估计共振峰。

原理框图如图4.2所示。

图4.2共振峰检测框图

三、实验内容

1、利用同态处理计算给出语音信号的基音周期和基频。

[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）;

ylen=frameSize;

frameSize=length（frame）;

frame2=frame.*hamming（length（frame））;

rwy=rceps（frame2）;

cepstrum=rwy（1:

ylen/2）;

fori=1:

ylen/2

cepstrum1（i）=rwy（ylen/2+1-i）;

end

fori=（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_cepind]=max（cn）;

ifmx_cep>0.08

a=fs/（LF+ind）;

else

a=0;

end

pitch=a

%画图

figure

（1）;

subplot（3,1,1）;

plot（time1,y）;

title（'语音波形'）;

axistight

ylim=get（gca,'ylim'）;

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（'幅度'）;

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]=max（cn）;

%共振峰检测核心代码

%找到最大的突起的位置

NN=ind+LF;

ham=hamming（NN）;

cep=cepstrum（1:

NN）;

ceps=cep.*ham;

%formant1=20*log（abs（fft（ceps）））;

formant1=abs（fft（ceps））;

formant（1:

2）=formant1（1:

2）;

fort=3:

NN

%dosomemedianfiltering

z=formant1（t-2:

t）;

md=median（z）;

formant2（t）=md;

end

fort=1:

NN-1

ift<=2

formant（t）=formant1（t）;

else

formant（t）=formant2（t-1）*0.25+formant2（t）*0.5+formant2（t+1）*0.25;

end

end

subplot（3,1,1）

plot（cepstrum）;

title（'倒谱'）;

xlabel（'样点数'）;

ylabel（'幅度'）;

axis（[0,ylen/2,-0