大学本科语音信号处理实验讲义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
升级会员 