数字音频处理.doc

数字语音实验
吕佩壕 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.2 hamming窗
我们发现：在用短时能量反映语音信号的幅度变化时，不同的窗函数以及相应窗的长短均有影响。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;
for i=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;
for i=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，则得离散的短时傅立叶变换如下：

根据信号的时宽带宽积为一常数之一基本性质，可知W（ejw）主瓣宽度和窗口宽度成反比，N越大 W（ejw）越窄。尤其是N值大于语音音素长度时W（ejw）已不能反应语音音素的频谱了。因此，应折衷选择窗的宽度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('original signal');
grid
N=256;
h=hamming(N);
for m=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=clcx=x+CL,x≤-CLx-CL,x≥CL0,-CL高的削波电平可以得到清楚的周期性指示；在整个语音段的持续时间内（如浊音语音的开始或终止处），信号幅度可能有相当大的变化，如果将中心削波电平置为语音段范围内最大幅度的高百分比（60~80%）上，就会有更多的波形幅度低于削波电平而丢失，使基音周期估计出现问题。
无削波的自相关函数会带来很多基音估计错误，特别是对短基音周期；中心削波自相关函数消除了基音估计中的大多数误估；使用基音轨迹平滑可进一步减少遗留的错误估计。

Matlab程序：
a=wavread('C:\audio.wav'); %读取语音文件
L=length(a) %测定语音长度
m=max(a)
for i=1:L
a(i)=a(i)/m(1); %数据归一化
end
m=max(a) %找到最大正值
n=min(a) %找到最小负值
ht=(m+n)/2; %保证幅度之余横坐标对称
for i=1:L %数据中心下移，保持和横坐标轴对称
a(i)=a(i)-ht(1);
end
figure(1)
subplot(2,1,1)
plot(a)
axis([0,170000,-1,1]);
title('Before Center Clipping')
xlabel('The Sample Point')
ylabel('Amplitude')
coeff0.3; %中心削波函数系数取0.3
th0=max(a)*coeff; %求中心削波函数门限
for k=1:L %中心削波
if a(k)>=th0
a(k)=a(k)-th0(1);
elseif a(k)<=th0
a(k)=a(k)+th0(1);
else a(k)=0;
end
end
m=max(a);
for i=1; %中心削波函数幅度归一化
a(i)=a(i)/m(1);
end
subplot(2,1,2)
plot(a)
axis([0,17000,-2,2]);
title('After Center Clipping')
xlabel('The Sample Point')
ylabel('Amplitude') %没有经过中心削波的修正自相关

b=wavread('C:\audio.wav');
N=2048; %选择的窗长，加N=320的矩形窗
A=[];
for k=1:2048 %选择延时长度
sum=0;
for m=1:N
sum=sum+b(m+7500)*b(m+7500+k-1); %计算自相关
end
A(k)=sum;
end
for k=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('Modified Autocorrelation Before Center Clipping')
xlabel('Delay')
ylabel('Amplitude')
N=2048; %选择的窗长，加N=320的矩形窗
A=[];
for k=1:2048 %选择延时长度
sum=0;
for m=1:N
sum=sum+a(m+7500)*a(m+7500+k-1); %计算自相关
end
A(k)=sum;
end
for k=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('Modified Autocorrelation After Center Clipping')
xlabel('Delay')
ylabel('Amplitude')

三、心得体会
在随着上课的进度做这样几个实验，确确实实让我对数字音频这门课程的学习有更深的理解，尤其是在做过仿真以后，能帮助我更深的理解平时上课的时候不能理解的内容，而且通过在试验中对于各个参数的不停修正，让我对一些公式的具体应用的理解更加深刻，并且在自己独立解决这些问题的时候，没有别人的帮助确实能够从各个方面锻炼自己的一些能力。虽然没能做完所有的实验比较可惜，但是从中确实获得了不少知识，并且解决问题的过程让我十分享受并且受益匪浅。

吕佩壕
10024134