基于短时自相关函数的基因周期检测讲解

1、基音周期的检测什么是基音周期检测 人在发音时，根据声带是否震动可以将语音信号分为清音跟浊音两种。浊音又称有声语言，携带者语言中大部分的能量，浊音在时域上呈现出明显的周期性；而清音类似于白噪声，没有明显的周期性。发浊音时，气流通过声门使声带产生张弛震荡式振动，产生准周期的激励脉冲串。这种声带振动的频率称为基音频率，相应的周期就成为基音周期。 通常，基音频率与个人声带的长短、薄厚、韧性、劲度和发音习惯等有关系，在很大程度上反应了个人的特征。此外，基音频率还跟随着人的性别、年龄不同而有所不同。一般来说，男性说话者的基音频率较低，大部分在70200Hz的范围内，而女性说话者和小孩的基音频率相对较高，在

2、200450Hz之间。 目前，基音检测的算法有很多种，常用的检测算法有自相关法、平均幅度差函数法、并行处理法、倒谱法、简化逆滤波法等。 本次基音周期的检测，我们选用的是短时自相关函数法，包括四个模块。 第一个模块为基音的端点检测，主要为了区分浊音和清音。第二个模块为基音检测中的带通滤波器，主要为了减少共振峰的干扰。第三个模块为短时自相关函数法做基音检测，主要为了计算出基音周期。第四个模块为平滑处理，主要为了消除偏离值点。基音周期的检测（端点检测）1.设计数据文件名称:audioread功能：读入以.wav为扩展名的音频文件。调用格式：y,fs,nbits=audioread(wavFilena

3、me)说明：输出参数Y是数字化的音频；fs是信号采样频率；nbits是每个样点在编码时的二进制位数。2.分帧处理先定义一个帧数*帧长的数组，然后计算帧的偏移量来，按语音数据的位置把x放入y中。名称：enframe功能：把语音信号按帧长和帧移进行分帧调用格式：f=enframe(x,win,inc)3.端点检测名称：pitch_vad1功能：用能熵比法进行端点检测。说明：输入参数y是分帧后的数组，一般是一列表示一帧数据;fn是信号的总帧数；T1是一个阈值；输出参数voiceseg是语音端点的信息；vosl是几个有话段，当该帧为有话段时SF=1；否则为0，Ef是能熵比值。pitch_vad1for

4、i=1:fnSp=abs(fft(y(:,i);%FFT取幅值(:表示所有y(:,1)表示第一列)Sp=Sp(1:wlen/2+1); %只取正频率部分Esum(i)=sum(Sp.*Sp);%计算能量值（能量放入Esum里）prob=Sp/(sum(Sp); %计算概率H(i)=-sum(prob.*log(prob+eps);%求谱熵值（eps表示很小的数，避免为零）（谱熵放入H里）end谱熵值很类似与过零率值，在有话段内的谱熵值小于噪声段的谱熵值，所以如同能零比一样，用能熵比更能突出噪声段的差别。基音周期的检测（基音检测中的带通滤波器）什么是带通滤波器 带通滤波器（band-passfi

5、lter）是一个允许特定频段的波通过同时屏蔽其他频段的滤波器。带通滤波器的作用 在使用相关法和AMDF法的基音检测之前常用到低通滤波器和带通滤波器，其主要目的是减少共振峰的干扰。我用了500Hz作为滤波器的上限频率，选择截止频率高不利于减少噪声和共振峰的影响。所以在本次的基音检测中的预滤波器选择的带宽为60500Hz，高频截止频率选择500Hz，是因为基频区间的高端就在这个区域中，低频截止频率选择60Hz是为了减少工频和低频噪声的干扰。什么是共振峰 共振峰是指在声音的频谱中能量相对集中的一些区域，共振峰不但是音质的决定因素，而且反映了声道（共振腔）的物理特征。 声音在经过共振腔时，受到腔体的滤

6、波作用，使得频域中不同频率的能量重新分配，一部分因为共振腔的共振作用得到强化，另一部分则受到衰减。由于能量分布不均匀，强的部分犹如山峰一般，故而称之为共振峰。为什么要减小共振峰 声道特性与基音周期估计有关，会对基音周期检测造成一定干扰。从语音信号中去除声道的影响，直接取出仅与声带振动有关的声源信息并非易事。声道共振峰可能强烈改变声门波形的结构，从而严重影响激励信号的谐波结构，给基因检测造成困难。我选用IIR滤波器中的椭圆滤波器，因为IIR滤波器的运算量比FIR少，当然IIR滤波器会带来延迟，也就是相位的变化，但语音信号是对相位不敏感的信号；又选用椭圆滤波器，因为它在经典滤波器设计中相同过渡带和

7、带宽条件下，需要的阶数比较小。椭圆滤波器（Ellipticfilter）是在通带和阻带等波纹的一种滤波器。椭圆滤波器在通带和阻带内都具有等波纹幅频响应特性。由于其极点位置与经典场论中的椭圆函数有关，所以由此取名为椭圆滤波器。对于给定的阶数和波纹要求，椭圆滤波器能得到较其它滤波器更窄的过渡带宽，可以获得对理想滤波器幅频响应的最好逼近，是一种性价比很高的滤波器，椭圆滤波器相比其他类型的滤波器，在阶数相同的条件下有着最小的通带和阻带波动。滤波器的要求为采样频率8000Hz，通带是60500Hz，通带波纹为1dB，阻带分别为30Hz和2000Hz，阻带衰减为40dB。 滤波器的阶数，就是指过滤谐波的次

8、数，一般来讲，同样的滤波器，其阶数越高，滤波效果就越好，但是，阶数越高，成本也就越高，因此，选择合适的阶数是非常重要的。所以本次实验选用椭圆6阶带通滤波器。实现代码%pr8_1_1clearall;clc;closeall;fs=8000;fs2=fs/2;%采样频率Wp=60500/fs2;%滤波器通带Ws=202000/fs2;%滤波器阻带Rp=1;Rs=40;%通带的波纹和阻带的衰减n,Wn=ellipord(Wp,Ws,Rp,Rs);%计算滤波器的阶数b,a=ellip(n,Rp,Rs,Wn);%计算滤波器的系数fprintf(b=%5.6f%5.6f%5.6f%5.6f%5.6f%5

9、.6f%5.6fn,b)fprintf(a=%5.6f%5.6f%5.6f%5.6f%5.6f%5.6f%5.6fn,a)代码理解 Ellipord函数的功能是求滤波器的最小阶数，其调用格式为n,Wn=ellipord(Wp,Ws,Rp,Rs,s) 其中各参量分别为：n-椭圆滤波器最小阶数；Wp-椭圆滤波器通带；Ws-椭圆滤波器阻带；Rp-通带波纹（dB）；Rs-阻带衰减(dB)； Ellip函数的功能是用来设计椭圆滤波器，其调用格式： b,a=ellip(n,Rp,Rs,Wp) 由程序计算出滤波器系数为： b=0.012280-0.0395080.0421770.000000-0.04217

10、70.039508 a=1.000000-5.52714612.854342-16.11030711.479789-4.4101790.713507基音周期的检测（短时自相关函数法一 算法原理 设语音信号的时间序列是x(n)，它的自相关函数定义为： R(k)=x(n)x(n+k)，K是时间的延迟量。 短时自相关函数有以下重要性质 如果信号x(n)具有周期性，周期是P，那么它的自相关函数R(k) 也具有周期性，而且周期与信号x(n)的周期性相同。 当k=0时，短时自相关函数具有最大值，即在延迟量为0，P，2P，时，周期信号的自相关函数也达到最大值。 短时自相关函数是偶函数，即R(k) =R(-k

11、)。短时自相关函数基因检测的主要原理大都是利用它的这些性质，通过比较原始信号和它的延迟后的信号之间的类似性质来确定基音周期。如果延迟量等于基音周期，那么两个信号具有最大类似性，或是直接找出短时自相关函数的两个最大值间的距离，即为基音周期的初估值。在用短时自相关函数检测基音时，常用归一化的自相关函数，表达式为r(k)=R(k)/R(0)以上性质中已指出，k=0，R(0)为最大值。所以r(k)的模值永远小于或等于1。已知基音频率范围为60Hz500Hz之间，采样率为fs时，则基音周期（样点值）范围在fs/500fs/60之间，然后在这个范围之间寻找归一化相关函数的最大值，对应的延迟量就是基音周期。

12、二代码编写用短时自相关函数法得基音周期的程序：functionperiod=ACF_corr(y,fn,vseg,vsl,lmax,lmin)pn=size(y,2);ifpn=fn,y=y;end%把y转换为每列数据表示一帧语音信号wlen=size(y,1);%取得帧长period=zeros(1,fn);%初始化fori=1:vsl%只对有话段数据处理ixb=vseg(i).begin;ixe=vseg(i).end;ixd=ixe-ixb+1;%求取一段有话段的帧数fork=1:ixd%对该段有话段数据处理u=y(:,k+ixb-1);%取来一帧数据ru=xcorr(u,coeff);

13、%计算归一化自相关函数ru=ru(wlen:end);%取延迟量为正值的部分tmax,tloc=max(ru(lmin:lmax);%在PminPmax范围内寻找最大值period(k+ixb-1)=lmin+tloc-1;%给出对应最大值的延迟量endend三实验结果四结果分析在端点检测和滤波后，再经过短时自相关函数法处理得到上面的实验结果图，可看出基音周期有很多的野点，需要进行平滑处理，也可以看出基音周期大约为40个样点数，即基音周期为40/fs,采样率为8000Hz，得周期为0.005s,基音频率为200Hz.基音周期的检测（平滑处理） 线性滤波算法 让基音轨迹序列通过一个“低通”FIR

14、滤波器。由于基音轨迹中的野点都是一些突发跳动点，它们具有很高的“频率分量”，采用“低通”滤波即可以将其去除。 线性平滑是用滑动窗口进行线性滤波处理，即 式中，w(m)，m=-L,-L+1,0,1,2,L为（2L+1）点平滑窗，满足平滑处理中值滤波算法：一种非线性滤波在被平滑点的左右各取L个样点。连同被平滑点共同构成一组信号采样值（共（2L+1）个样值），然后将这（2L+1）个样值按大小次序排成一对，取此队列中间者作为平滑器的输出。L值一般取为1或2，即中值平滑的“窗口”一般套住3或5个样值。中值平滑的优点是既可以有效地去除少量的野点，又不会破坏基音周期轨迹中两个平滑段之间的阶跃性变化。matl

15、ab函数y=medfilt1(x,n)x为输入序列；k为窗长，即套住的样点数，一般取3或5。y是中值滤波后的输出序列。设有一个一维序列f1,f2,fn,取窗口长度（点数）为m（m为奇数）,对其进行中值滤波，就是从输入序列中相继抽出m个数fi-v,fi-1,fi,fi+1,fi+v(其中fi为窗口中心值,v=(m-1)/2),再将这m个点按其数值大小顺序排序，取其序号的中心点的那个数作为滤波输出。数学公式表示为：Yi=Medfi-v,fi-1,fi,fi+1,fi+viNv=(m-1)/2；Yi称为序列fi-v,fi-1,fi,fi+1,fi+v的中值举例来说，输入：Y1-10:1,2,3,4,

16、5,6,7,8,9,10.取区间2k=4,所以k=2;执行中值滤波K=中值滤波(Y)由x-k+1=1，所以当k=2时，x=2滤波时：K1=Y1K2=(Y1、Y2、Y3、Y4)的中间值，即为2或3 T0=pitfilterm1(period,voiceseg,vosl); functiony=pitfilterm1(x,vseg,vsl) y=zeros(size(x);%初始化 fori=1:vsl%有段数据 ixb=vseg(i).begin;%该段的开始位置 ixe=vseg(i).end;%该段的结束位置 u0=x(ixb:ixe);%取来一段数据 y0=medfilt1(u0,5);%5点的中值滤波 v0=linsmoothm(y0,5);%线性平滑 y(ixb:ixe)=v0;%赋值给y end程序functiony=linsmoothm(x,n)ifnargin2n=3;endwin=han