完整版数字音频处理.docx
《完整版数字音频处理.docx》由会员分享，可在线阅读，更多相关《完整版数字音频处理.docx（12页珍藏版）》请在冰豆网上搜索。
完整版数字音频处理
数字语音实验
吕佩壕10024134
一、实验要求
1.编程实现一句话语音的短时能量曲线，并比较窗长、窗口形状（以直角
窗和和哈明窗为例）对短时平均能量的影响；
2.编程分析语音信号的短时谱特性，并比较窗长、窗口形状（以直角窗和
和哈明窗为例）对语音短时谱的影响；
3.运用低通滤波器、中心削波和自相关技术估计一段男性和女性语音信号
的基音周期，画出基音轨迹曲线，给出估计准确率。
二、实验原理及实验结果
1.窗口的选择
通过对发声机理的认识，语音信号可以认为是短时平稳的。
在5~50ms的范
围内，语音频谱特性和一些物理特性参数基本保持不变。
我们将每个短时的语音称为一个分析帧。
一般帧长取10~30ms我们采用一个长度有限的窗函数来截取语音信号形成分析帧。
通常会采用矩形窗和汉明窗。
图1.1给出了这两种窗函数在窗长N=50时的时域波形。
图1.1矩形窗和hamming窗的时域波形矩形窗的定义：
一个N点的矩形窗函数定义为如下：
1,0nN
w（n）°,其他
"）,0
0.540.46cos（2
w（n）=0,其他
这两种窗函数都有低通特性，通过分析这两种窗的频率响应幅度特性可以发
现（如图1.2）:
矩形窗的主瓣宽度小（4*pi/N）,具有较高的频率分辨率，旁瓣峰值大（-13.3dB），会导致泄漏现象；汉明窗的主瓣宽8*pi/N，旁瓣峰值低
（-42.7dB），可以有效的克服泄漏现象，具有更平滑的低通特性。
因此在语音频谱分析时常使用汉明窗，在计算短时能量和平均幅度时通常用矩形窗。
表1.1
对比了这两种窗函数的主瓣宽度和旁瓣峰值。
图1.2矩形窗和Hamming窗的频率响应
2.短时能量
由于语音信号的能量随时间变化，清音和浊音之间的能量差别相当显著。
因此对语音的短时能量进行分析，可以描述语音的这种特征变化情况。
定义短时能量为：
n
22
En[x（m）w（nm）][x（m）w（nm）]
mmnN1
，其中N为窗长
特殊地，当采用矩形窗时，可简化为：
Enx2（m）
m
图2.1和图2.2给出了不同矩形窗和hamming窗长,对所录的语音“我是吕佩壕”的短时能量函数：
（1）矩形窗（从上至下依次为“我是吕佩壕”波形图，窗长分别为
32,64,128,256,512的矩形窗的短时能量函数）：
』一一屮・i一」亠一―LLi
D?
459IdFId
I||I
号二^亠一14-”一一丄t励般显仏
s2A韦4I?
14
xlf
*衬
图2.1矩形窗
（2）hamming窗（从上至下依次为“我是吕佩壕”波形图，窗长分别为
32,64,128,256,512的hamming窗的短时能量函数）：
图2.2hamming窗
我们发现：
在用短时能量反映语音信号的幅度变化时，不同的窗函数以及相
应窗的长短均有影响。
hamming窗的效果比矩形窗略好。
但是，窗的长短影响起决定性作用。
窗过大（N很大），等效于很窄的低通滤波器，不能反映幅度En
的变化；窗过小（N很小），短时能量随时间急剧变化，不能得到平滑的能量函数。
在11.025kHz左右的采样频率下，N选为100~200比较合适。
短时能量函数的应用:
1）可用于区分清音段与浊音段。
En值大对应于浊音段，En值小对应于清音段。
2）可用于区分浊音变为清音或清音变为浊音的时间（根据En值的变化趋势）。
3）对高信噪比的语音信号，也可以用来区分有无语音（语音信号的开始点或终止点）。
无信号（或仅有噪声能量）时，En值很小，有语音信号时，能量显著增大。
Matlab程序：
figure（3）;a=wavread（'C:
\audio.wav'）;subplot（6,1,1）,plot（a）;
N=32;
fori=2:
6
h=rectwin（24（i-2）*N）;
b=a.*a;
En=conv2（h,b）;%
subplot（6,1,i）,plot（En）;
i=i+1;
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
求短时能量函数En
figure（4）;
a=wavread（'C:
\audio.wav'）;subplot（6,1,1）,plot（a）;
N=32;
fori=2:
6
h=hamming（24（i-2）*N）;%
b=a.*a;
En=conv2（h,b）;%
subplot（6,1,i）,plot（En）;
i=i+1;
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
形成一个汉明窗，长度为2.A（i-2）*N
求短时能量函数En
3.短时谱
由于语音信号是短时平稳的随机信号，某一语音信号帧的短时傅立叶变换的定义为：
于是可以得到:
Xn（ejw）ejwnw（k）x（nk）ejwk
k
假定：
Xn（ejw）w（k）x（nk）ejwk
k
则可以得到：
Xn（ejw）ejwnXn（ejw）
同样，不同的窗口函数，将得到不同的傅立叶变换式的结果。
由上式可见，短时傅立叶变换有两个变量：
n和①，所以它既是时序n的离散函数，又是角频率①的连续函数。
与离散傅立叶变换逼近傅立叶变换一样，如令①=2nk/N，则
得离散的短时傅立叶变换如下：
Xn（ej2k/N）Xn（k）
x（m）w（nm）ej2km/N,（0kN1）
m
根据信号的时宽带宽积为一常数之一基本性质，可知W（ejw）主瓣宽度和窗口宽度成反比，N越大W（ejw）越窄。
尤其是N值大于语音音素长度时
w（ew）已不能反应语音音素的频谱了。
因此，应折衷选择窗的宽度n。
另外,窗的形状也对短时谱有影响，如矩形窗，虽然频率分辨率很高，但由于第一旁瓣的衰减很小，所以不适合用于频谱成分很宽的语音分析中，而汉明窗在频率范围中分辨率较高，而且旁瓣衰减大，具有频谱泄露少的优点，所以在求短时频谱时一般采用汉明窗。
图3.1到图3.6分别是不同窗长的汉明窗下的短时谱仿真图：
图3.1窗长N=4000
图3.1窗长N=4500
图3.1窗长N=7000
图3.1窗长N=10000
Matlab程序：
短时谱
figure⑴cleara=wavread（'C:
\audio.wav'）;
subplot（2,1,1）,plot（a）;title（'originalsignal'）;
grid
N=256;
h=hamming（N）;
form=1:
Nb（m）=a（m）*h（m）end
y=20*log（abs（fft（b）））
subplot（2,1,2）
plot（y）;title（'grid
短时谱'）;
4
理。
器,
.基于中心削波的基音检测
在基音检测的时候，为了改善基音检测器的性能我们都要进行与处理和后处在进行与处理的时候，具体的做法就是进行谱的平整处理。
谱平整从语音信号中排除共振峰结构，使每个谐波有相同的幅度。
主要方法有线性方法和非线性方法两种，线性方法是使用线性预测误差滤波非线性方法是使用中心削波技术。
下图为常用的三种中心削波函数：
x+CL,x<-CL
（a）
y=
:
clc（x）=
:
{x-C
；L,x>CL
0,-CL
0,|x|
（b）
y=
:
clp（x）=
={
x,|x|
>CL
-1,x
<-CL
（c）
y=
:
sgn（x）:
={1,x
>CL
0,|x|
其中CL为削波电平，由实际语音信号确定。
下图4.1为中心削波语音信号波形，图4.2为中心削波的自相关:
图4.2中心削波的自相关
由上图可知：
削波电平为样点中最大值的30%削波后的剩余信号只是位于原始基音周期上的几个脉冲；所得的自相关函数中引起混扰的外来峰相当少。
高的削波电平可以得到清楚的周期性指示；在整个语音段的持续时间内（如浊音语音的开始或终止处），信号幅度可能有相当大的变化，如果将中心削波电平置为语音段范围内最大幅度的高百分比（60~80%上,就会有更多的波形幅度低于削波电平而丢失，使基音周期估计出现问题。
无削波的自相关函数会带来很多基音估计错误，特别是对短基音周期；中心削波自相关函数消除了基音估计中的大多数误估；使用基音轨迹平滑可进一步减
少遗留的错误估计。
Matlab程序：
a=wavread（'C:
\audio.wav'）;%读取语音文件L=length（a）%测定语音长度
m=max（a）
fori=1:
L
a（i）=a（i）/m
（1）;%数据归一化
end
m=max（a）%找到最大正值n=min（a）%找到最小负值
ht=（m+n）/2;%保证幅度之余横坐标对称
fori=1:
L%数据中心下移，保持和横坐标轴对称
a（i）=a（i）-ht
（1）;
end
figure
（1）
subplot（2,1,1）
plot（a）
axis（[0,170000,-1,1]）;
title（'BeforeCenterClipping'）
xlabel（'TheSamplePoint'）ylabel（'Amplitude'）coeff=0.3;%中心削波函数系数取0.3th0=max（a）*coeff;%求中心削波函数门限fork=1:
L%中心削波
ifa（k）>=th0a（k）=a（k）-th0
（1）;
elseifa（k）<=th0a（k）=a（k）+th0
（1）;
elsea（k）=0;
endendm=max（a）;fori=1;%中心削波函数幅度归一化
a（i）=a（i）/m
（1）;
endsubplot（2,1,2）plot（a）axis（[0,17000,-2,2]）;title（'AfterCenterClipping'）xlabel（'TheSamplePoint'）ylabel（'Amplitude'）%没有经过中心削波的修正自相关
b=wavread（'C:
\audio.wav'）;
N=2048;%选择的窗长，加N=320的矩形窗
A=[];
fork=1:
2048%选择延时长度
sum=0;
form=1:
Nsum=sum+b（m+7500）*b（m+7500+k-1）;%计算自相关end
A（k）=sum;
end
fork=1:
2048
B（k）=A（k）/A
（1）;%归一化A（k）endfigure
（2）subplot（2,1,1）plot（B）axis（[0,2200
升级会员 