(完整版)数字音频处理

1、数字语音实验吕佩壕10024134一、实验要求1. 编程实现一句话语音的短时能量曲线，并比较窗长、窗口形状（以直角窗和和哈明窗为例）对短时平均能量的影响；2. 编程分析语音信号的短时谱特性，并比较窗长、窗口形状（以直角窗和和哈明窗为例）对语音短时谱的影响；3. 运用低通滤波器、中心削波和自相关技术估计一段男性和女性语音信号的基音周期，画出基音轨迹曲线，给出估计准确率。二、实验原理及实验结果1.窗口的选择通过对发声机理的认识，语音信号可以认为是短时平稳的。在550ms的范围内，语音频谱特性和一些物理特性参数基本保持不变。我们将每个短时的语音称为一个分析帧。一股帧长取1030ms我们采用一个长度有

2、限的窗函数来截取语音信号形成分析帧。通常会采用矩形窗和汉明窗。图1.1给出了这两种窗函数在窗长N=50时的时域波形。n211p1.8-0.9"*i1.6-0.8/J1.40.71.2-0.6-/-10.50.8-一w0.41n1f0.6-0.3r0.40.2/I0.2-0.101irh0窗hanming0406004060矩形窗20sample20sample图1.1矩形窗和hamming窗的时域波形矩形窗的定义：一个N点的矩形窗函数定义为如下：1,0nNw（n）0,Hamming的定义：一个N点的hamming®函数定义为如下:),0nNN10.540.46cos(2w(

3、n)=°,其他这两种窗函数都有低通特性，通过分析这两种窗的频率响应幅度特性可以发现（如图1.2）:矩形窗的主瓣宽度小（4*pi/N）,具有较高的频率分辨率，旁瓣峰值大（-13.3dB），会导致泄漏现象；汉明窗的主瓣宽8*pi/N,旁瓣峰值低（-42.7dB）,可以有效的克服泄漏现象，具有更平滑的低通特性。因此在语音频谱分析时常使用汉明窗，在计算短时能量和平均幅度时通常用矩形窗。表1.17Ei.Ura|yvn1VV1-ifX/ILfJVii1f-!'111.11yph1-500.10.20.30.40.50.60.70.80.91归一化频率(f/fs)-1000对比了这两种窗函

4、数的主瓣宽度和旁瓣峰值。0度幅归一化频率（f/fs）Hamming窗频率响应图1.2矩形窗和Hamming窗的频率响应2.短时能量由丁语音信号的能量随时间变化，活音和浊音之间的能量差别相当显著。因此对语音的短时能量进行分析，可以描述语音的这种特征变化情况。定义短时能量为：n22Enx(m)w(nm)x(m)w(nm)mmnN1,其中N为窗长特殊地，当采用矩形窗时，可简化为:Enx2(m)m图2.1和图2.2给出了不同矩形窗和hamming窗长，对所录的语音“我是吕佩壕”的短时能量函数：（1）矩形窗（从上至下依次为“我是吕佩壕”波形图，窗长分别为32,64,128,256,512的矩形窗的短时能

5、量函数）：图2.1矩形窗（2）hamming窗（从上至下依次为“我是吕佩壕”波形图，窗长分别为32,64,128,256,512的hamming®的短时能量函数）：图2.2hamming窗我们发现：在用短时能量反映语音信号的幅度变化时，不同的窗函数以及相应窗的长短均有影响。hamming®的效果比矩形窗略好。但是，窗的长短影响起决定性作用。窗过大（N很大），等效丁很窄的低通滤波器，不能反映幅度En的变化；窗过小（N很小），短时能量随时间急剧变化，不能得到平滑的能量函数。在11.025kHz左右的采样频率下，N选为100200比较合适。短时能量函数的应用：1）可用于区分活音段

6、与浊音段。En值大对应于浊音段，En值小对应于活音段。2）可用于区分浊音变为活音或活音变为浊音的时间（根据En值的变化趋势）。3）对高信噪比的语音信号，也可以用来区分有无语音（语音信号的开始点或终止点）。无信号（或仅有噪声能量）时，En值很小，有语音信号时，能量显著增大。Matlab程序：figure(3);a=wavread('C:audio.wav');subplot(6,1,1),plot(a);N=32;fori=2:6h=rectwin(2.A(i-2)*N);b=a.*a;En=conv2(h,b);%subplot(6,1,i),plot(En);i=i+1;if

7、(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');endendfigure(4);a=wavread('C:audio.wav');subplot(6,1,1),plot(a);N=32;fori=2:6h=hamming(2.A(i-2)*N);%b=a.*a;En=conv2(h,b);%sub

8、plot(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求短时能量函数En形成一个汉明窗，长度为2.A（i-2）*N求短时能量函数Enend3.短时谱由于语音信号是短时平稳的随机信号，某一语音信号帧的短时傅立叶变换的定义为：Xn(ejw)x(m)w(nm)e

9、jwmm其中w(n-m)是实窗口函数序歹U,n表示某一语音信号帧。令n-m=k',则得到：Xn(ejw)w(k')x(nk')ejw(nk)k'于是可以得到：Xn(ejw)ejwnw(k)x(nk)ejwkk假定：Xn(ejw)w(k)x(nk)ejwkk则可以得到：Xn(ejw)ejwnXn(ejw)同样，不同的窗口函数，将得到不同的傅立叶变换式的结果。由上式可见，短时傅立叶变换有两个变量：n和3,所以它既是时序n的离散函数，乂是角频率3的连续函数。与离散傅立叶变换逼近傅立叶变换一样，如令3=2兀k/N,则得离散的短时傅立叶变换如下：片伦八小)Xn(k)x(m

10、)w(nm)ej2km/N,(0kN1)m根据信号的时宽带宽积为一常数之一基本性质，可知W(ejw)主瓣宽度和窗口宽度成反比，N越大W(ejw)越窄。尤其是N值大于语音音素长度时w(aw)已不能反应语音音素的频谱了。因此，应折衷选择窗的宽度n。另外,窗的形状也对短时谱有影响，如矩形窗，虽然频率分辨率很高，但由于第一旁瓣的衰减很小，所以不适合用于频谱成分很宽的语音分析中，而汉明窗在频率范围中分辨率较高，而且旁瓣衰减大，具有频谱泄露少的优点，所以在求短时频谱时一般采用汉明窗。图3.1到图3.6分别是不同窗长的汉明窗下的短时谱仿真图：II巧，地2制矿广图3.1窗长N=4100图3.1窗长N=4500

11、图3.1窗长N=7000图3.1窗长N=10000Matlab程序：短时谱figure(1)cleara=wavread('C:audio.wav');subplot(2,1,1),plot(a);title('originalsignal');gridN=256;h=hamming(N);form=1:Nb(m)=a(m)*h(m)endy=20*log(abs(fft(b)subplot(2,1,2)plot(y);title('grid短时谱');4理。器,.基丁中心削波的基音检测在基音检测的时候，为了改善基音检测器的性能我们都要进行与处理

12、和后处在进行与处理的时候，具体的做法就是进行谱的平整处理。谱平整从语音信号中排除共振峰结构，使每个谐波有相同的幅度。主要方法有线性方法和非线性方法两种，线性方法是使用线性预测误差滤波非线性方法是使用中心削波技术。下图为常用的三种中心削波函数：3)中心削波函数I1EehzireenterClipri-g2U0Q4W1QGQOO0(1QQ10-Q0Q*00。UOQQ160CWThePoinr(访三电平中心削波函数x+CL,x<-CL(a) y=clc(x)=(x-CL,x>CL0,-CL<x<CL0,|x|<CL(b) y=clp(x)=(x,|x|>CL-1,

13、x<-CL(c) y=sgn(x)=(1,x>CL0,|x|<CL其中CL为削波电平，由实际语音信号确定。下图4.1为中心削波语音信号波形,图4.2为中心削波的自相关:图4.1中心削波语音信号波形211040CbOUHOd11M012C014001曲12IflOU2(10022tflMDdrfiedAutocafFellsiiflriEafcreILenierClippingD«iayMudilivdAuiocDmtaLiotiAilwC*ni*(白ipfNngo51no1;即Ir1o图4.2中心削波的自相关由上图可知：削波电平为样点中最大值的30%削波后的剩余信号

14、只是位丁原始基音周期上的几个脉冲；所得的自相关函数中引起混扰的外来峰相当少。高的削波电平可以得到活楚的周期性指示；在整个语音段的持续时间内（如浊音语音的开始或终止处）,信号幅度可能有相当大的变化，如果将中心削波电平置为语音段范围内最大幅度的高白分比（6080%上，就会有更多的波形幅度低丁削波电平而丢失，使基音周期估计出现问题。无削波的自相关函数会带来很多基音估计错误，特别是对短基音周期；中心削波自相关函数消除了基音估计中的大多数误估；使用基音轨迹平滑可进一步减少遗留的错误估计。Matlab程序：a=wavread('C:audio.wav');%读取语音文件L=length(a

15、)%测定语音长度m=max(a)fori=1:La(i)=a(i)/m(1);%数据归一化endm=max(a)%找到最大正值n=min(a)%找到最小负值ht=(m+n)/2;%保证幅度之余横坐标对称fori=1:L%数据中心下移，保持和横坐标轴对称a(i)=a(i)-ht(1);endfigure(1)subplot(2,1,1)plot(a)axis(0,170000,-1,1);title('BeforeCenterClipping')xlabel('TheSamplePoint')ylabel('Amplitude')coeff=0.3

16、;%中心削波函数系数取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(&#

17、39;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);%计算自相关endA(k)=sum;endfork=1:2048B(k)=A(k)/A(1);%归一化A(k)endfigure(2)subplot(2,1,1)plot(B)axis(0,2200,-1,1);title('ModifiedAutocorrelationBef

18、oreCenterClipping')xlabel('Delay')ylabel('Amplitude')N=2048;%选择的窗长，加N=320的矩形窗A=；fork=1:2048%选择延时长度sum=0;form=1:Nsum=sum+a(m+7500)*a(m+7500+k-1);%计算自相关endA(k)=sum;endfork=1:2048C(k)=A(k)/A(1)%归一化A(k)endfigure(2)subplot(2,1,2)plot(C)axis(0,2200,-1,1.0);title('ModifiedAutocorrelationAfterCenterClippin