完整word版数字音频处理.docx
《完整word版数字音频处理.docx》由会员分享，可在线阅读，更多相关《完整word版数字音频处理.docx（13页珍藏版）》请在冰豆网上搜索。
完整word版数字音频处理
数字语音实验
吕佩壕10024134
一、实验要求
1.编程实现一句话语音的短时能量曲线，并比较窗长、窗口形状（以直角窗和和哈明窗为例）对短时平均能量的影响；
2.编程分析语音信号的短时谱特性，并比较窗长、窗口形状（以直角窗和和哈明窗为例）对语音短时谱的影响；
3.运用低通滤波器、中心削波和自相关技术估计一段男性和女性语音信号的基音周期，画出基音轨迹曲线，给出估计准确率。
二、实验原理及实验结果
1.窗口的选择
通过对发声机理的认识，语音信号可以认为是短时平稳的。
在5~50ms的范围内，语音频谱特性和一些物理特性参数基本保持不变。
我们将每个短时的语音称为一个分析帧。
一般帧长取10~30ms。
我们采用一个长度有限的窗函数来截取语音信号形成分析帧。
通常会采用矩形窗和汉明窗。
图1.1给出了这两种窗函数在窗长N=50时的时域波形。
图1.1矩形窗和hamming窗的时域波形
矩形窗的定义：
一个N点的矩形窗函数定义为如下：
Hamming窗的定义：
一个N点的hamming窗函数定义为如下：
这两种窗函数都有低通特性，通过分析这两种窗的频率响应幅度特性可以发现（如图1.2）：
矩形窗的主瓣宽度小（4*pi/N），具有较高的频率分辨率，旁瓣峰值大（-13.3dB），会导致泄漏现象；汉明窗的主瓣宽8*pi/N，旁瓣峰值低（-42.7dB），可以有效的克服泄漏现象，具有更平滑的低通特性。
因此在语音频谱分析时常使用汉明窗，在计算短时能量和平均幅度时通常用矩形窗。
表1.1对比了这两种窗函数的主瓣宽度和旁瓣峰值。
图1.2矩形窗和Hamming窗的频率响应
2．短时能量
由于语音信号的能量随时间变化，清音和浊音之间的能量差别相当显著。
因此对语音的短时能量进行分析，可以描述语音的这种特征变化情况。
定义短时能量为：
，其中N为窗长
特殊地，当采用矩形窗时，可简化为：
图2.1和图2.2给出了不同矩形窗和hamming窗长,对所录的语音“我是吕佩壕”的短时能量函数：
（1）矩形窗（从上至下依次为“我是吕佩壕”波形图，窗长分别为32,64,128,256,512的矩形窗的短时能量函数）：
图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（2.^（i-2）*N）;
b=a.*a;
En=conv2（h,b）;%求短时能量函数En
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
figure（4）;
a=wavread（'C:
\audio.wav'）;
subplot（6,1,1）,plot（a）;
N=32;
fori=2:
6
h=hamming（2.^（i-2）*N）;%形成一个汉明窗，长度为2.^（i-2）*N
b=a.*a;
En=conv2（h,b）;%求短时能量函数En
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
3.短时谱
由于语音信号是短时平稳的随机信号，某一语音信号帧的短时傅立叶变换的定义为：
其中w（n-m）是实窗口函数序列，n表示某一语音信号帧。
令n-m=k'，则得到:
于是可以得到：
假定：
则可以得到：
同样，不同的窗口函数，将得到不同的傅立叶变换式的结果。
由上式可见，短时傅立叶变换有两个变量：
n和ω，所以它既是时序n的离散函数，又是角频率ω的连续函数。
与离散傅立叶变换逼近傅立叶变换一样，如令ω=2πk/N，则得离散的短时傅立叶变换如下：
根据信号的时宽带宽积为一常数之一基本性质，可知主瓣宽度和窗口宽度成反比，N越大越窄。
尤其是N值大于语音音素长度时已不能反应语音音素的频谱了。
因此，应折衷选择窗的宽度N。
另外，窗的形状也对短时谱有影响，如矩形窗，虽然频率分辨率很高，但由于第一旁瓣的衰减很小，所以不适合用于频谱成分很宽的语音分析中，而汉明窗在频率范围中分辨率较高，而且旁瓣衰减大，具有频谱泄露少的优点，所以在求短时频谱时一般采用汉明窗。
图3.1到图3.6分别是不同窗长的汉明窗下的短时谱仿真图：
图3.1窗长N=4000
图3.1窗长N=4100
图3.1窗长N=4500
图3.1窗长N=5000
图3.1窗长N=7000
图3.1窗长N=10000
Matlab程序：
短时谱
figure
（1）
clear
a=wavread（'C:
\audio.wav'）;
subplot（2,1,1）,
plot（a）;title（'originalsignal'）;
grid
N=256;
h=hamming（N）;
form=1:
N
b（m）=a（m）*h（m）
end
y=20*log（abs（fft（b）））
subplot（2,1,2）
plot（y）;title（'短时谱'）;
grid
4.基于中心削波的基音检测
在基音检测的时候，为了改善基音检测器的性能我们都要进行与处理和后处理。
在进行与处理的时候，具体的做法就是进行谱的平整处理。
谱平整从语音信号中排除共振峰结构，使每个谐波有相同的幅度。
主要方法有线性方法和非线性方法两种，线性方法是使用线性预测误差滤波器，非线性方法是使用中心削波技术。
下图为常用的三种中心削波函数：
（a）y=
（b）y=
（c）y=
其中CL为削波电平，由实际语音信号确定。
下图4.1为中心削波语音信号波形，图4.2为中心削波的自相关：
图4.1中心削波语音信号波形
图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.3
th0=max（a）*coeff;%求中心削波函数门限
fork=1:
L%中心削波
ifa（k）>=th0
a（k）=a（k）-th0
（1）;
elseifa（k）<=th0
a（k）=a（k）+th0
（1）;
elsea（k）=0;
end
end
m=max（a）;
fori=1;%中心削波函数幅度归一化
a（i）=a（i）/m
（1）;
end
subplot（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:
N
sum=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）
end
figure
（2）
subplot（2,1,1）
plot（B）
axis（[0,2200,-1,1]）;
title（'ModifiedAutocorrelationBeforeCenterClipping'）
xlabel（'Delay'）
ylabel（'Amplitude'）
N=2048;%选择的窗长，加N=320的矩形窗
A=[];
fork=1:
2048%选择延时长度
sum=0;
form=1:
N
sum=sum+a（m+7500）*a（m+7500+k-1）;%计算自相关
end
A（k）=sum;
end
fork=1:
2048
C（k）=A（k）/A
（1）%归一化A（k）
end
figure
（2）
subplot（2,1,2）
plot（C）
axis（[0,2200,-1,1.0]）;
title（'ModifiedAutocorrelationAfterCenterClipping'）
xlabel（'Delay'）
ylabel（'Amplitud
升级会员 