发布时间 : 星期一 文章(620512681) 自相关基频提取算法的MATLAB实现 - 图文更新完毕开始阅读8ab3e380e518964bce847c1a
第31卷总第80期 西北民族大学学报(自然科学版) V01．31．No．4
1 0年1 2 0 2月 Journal of Nonhw铭t University for Nationalities(Natural Science) Dec，2010
自相关基频提取算法的MATLAB实现 马
效敏1，郑文思2，陈 琪2 (西北民族大学科研处，甘肃兰州730030；2．中国民族信息技术研究院，甘肃兰州730030)
孽膊-聘妒掣摩护-庸、护岬—b序妒峰庸护峥高、。，砖、—气—昏尊舟峥舟、-，南
[摘要]在对自相关基音检测算法进行分析的基础上，对自相关基频提取算法进行了探讨．考虑到检测准确度和 检测速率两方面的因素，结合带通滤波及三电平削波模块，设计了基于Madab的估计方法，然后通过对一段具体的语音 进行处理，得到了比较准确的浊音语音的基音周期．实验证明。该方法简单有效．
【关键词]Matlab；基音检测；自相关函数；三电平削波
[中图分类号】TP391 [文献标识码]A [文章编号]1009—2102{2010104—0054—05
O引言 基频是指浊音发声过程中声带振动的频率．它是语音信号中的重要参数之一，在语音识别、语音合
成和语音编码中有广泛的应用．由于语音信号变化十分复杂，声门激励信号并不是一个完全周期的序 列、不同发声类型声带振动的复杂性和不同信号在反映声带振动时的差异有时会严重影响到激励信号 的基频检测．同时，基音周期变化范围大，从老年男性的80 Hz到儿童女性的500 Hz。接近三倍频程．这 些困难使得基频检测一直没有得到很好的解决，所以虽然已经提出各种各样的基音检测算法，如自相关 函数(ACF)法、峰值提取算法(PPA)、平均幅度差函数(AMDF)法等及基于它们的改进算法，但没有一 种算法，在各种环境条件下对各类人群都能获得满意的检测结果．
此外，声道特性与基音周期估计有关，会对基频提取造成一定干扰．从语音信号中去除声道的影 响，直接取出仅与声带振动有关的声源信息并非易事．声道共振峰可能强烈改变声门波形的结构，从而 严重影响激励信号的谐波结构?，给基音检测造成困难．因此，减少声遭共振峰的影响对基频提取尤为 重要．减少共振峰影响通常可用的方法：①采用带通滤波．②采用中心削波．
本文在采用滤波和削波的基础上，对用自相关函数提取语音基频进行了探讨．对于一段具体的语 音信号先通过带通滤波去掉语音信号中的高次共振峰，然后从语音信号中选取一段样本并将样本中获 得的语音样值点分为若干个语音帧；再经过三电平削波滤除大多数次要的峰，只保留明显的周期性的 峰．再用自相关函数求出自相关序列，并把自相关序列的周期转化为频率．最后通过中值平滑进行后处 理，平滑掉半频点、倍点和随机错误点，较准确地提取出浊音语音信号的基音频率．
1
基频提取算法原理 自相关函数的性质：若原信号具有周期性，那么它的自相关函数也具有周期性，并且周期与原信号
的周期相同．清音信号没有周期性，它的自相关函数R(”会随着k的增大呈衰减趋势．浊音信号具有 准周期性，它的R(k)在基音周期整数倍上有很大的峰值，通常取第一最大峰值点为基音周期点．自相 关函数法基音检测正是利用这一性质对语音信号进行基音检测的，同时结合滤波和三电平削波及
[收稿日期]2010—10—20 [作者简介]马效敏(1967一)，男，甘肃临夏人，讲师，主要从事计算机应用技术方面的研究
_。——54--——
万方数据
后处理使基频提取更准确、高效． 1．1带通滤 波
语音基音的频率范围是50 Hz～450 Hz，为了保留语音的一二次谐波以保证正确的处理结果，一般 先使语音信号通过一个频率范围为60 Hz～900 Hz的带通滤波器，下截止频率为60 Hz可以抑制电 源 的干扰，上截止频率为900 Hz不但可以保留基音的一二次谐波，还可以滤掉高次谐波【2|．在语音信 号 中，基音的频率一般在100 Hz～200 Hz之间，而第一共振峰在300 Hz～1 000 Hz之间，通过带通滤 波 器，可以去掉语音信号中的高次共振峰．实验表明，滤波处理后的信号只含有第一共振峰以下的基波 和 谐波分量，明显改善了检测效果．
1．2采样分帧 这里的“采样”是指从语音信号中选取一段样本，一般取样点数为帧长的整数倍．每秒钟 的采样样
本数叫做采样频率．分帧主要完成将取样模块中获得的语音样值点分为若干个语音帧． 语音是不平稳的 时变信号，在时间足够短的情况下，可以近似认为是平稳的，短时分析将语音流分
为一段一段来处理，每一段就被称为一“帧”．分帧时需对语音信号进行加窗操作，即用一个有限长度的 窗 序列截取一段语音信号来进行分析．该窗函数可以按时间方向滑动，以便分析任一时刻附近的信号． 常 见的窗函数有：方窗、Hamming窗及Harm窗．如果把窗函数理解成为某个滤波器的单位冲激响应， 由于窗函数一般是中间大两头小的光滑函数，因此该滤波器具有低通特性．窗口长度的选择非常重要， 窗长过短会使分析窗内没有包含足够的数据点来进行周期判断，且短时能量变化剧烈窗长过长，短时能 量是一段长时间的平均，不但不能反映语音信号基频的细节变化部分，而且使得计算量增大．窗口长度 至少要大于基音周期的两倍．
1．3 三电平削波
自相关函数提基频，运算量很大，其原因是乘法运算非常费时．然而从提取基频的角度看，短时自 相关函数所包含的信息许多是多余的，用自相关函数提取基音频率时，关心的是时间，也就是自相关函 数峰值出现的位置，而与峰值本身大小没有什么关系．这样就可以采用三电平削波法，其输入输出关系 式为：
卜Y(玎)=C 7[z(n)]=<0，l z(行“行)l≤CL}?CL ，
1
【一1，z(n)<一CLJ
即削波器的输出在z(竹)>C／．时为1，z(n)<一CL时为一1，除此以外都为0．这样可以增加刚刚超 过 电平的峰的重要性，滤除大多数次要的峰，只保留了明显的具有周期性的峰，避免了对清音段不必要 的 分析处理，这对后期处理有很大的好处．
削波电平由语音信号的峰值幅度来确定，它等于语音段最大幅度的一个固定百分数．这个门限的 选择是重要的，由于在整个语音段的持续时间内，信号幅度可能有较大变化，因此，采用较高的削波电平 可能会把一些有用的波形削去，然而削波电平太低，又达不到应用的效果．所以在不损失基音信息的情 况下应尽可能选得高些，以达到较好的效果．
1．4 自相关函数
经三电平削波后：自相关函数的计算极为简单．用Y(n)表示削波器的输出，则短时自相关函数如 下所示：
N一1一 K
R(k)=芝：[Y(咒+仇)硼7(优)][Y(，l+优+k)W 7(m+k)]． ，A=0
由自相关函数的定义可以看出其所具有的一些性质： 1)自相关函数是偶函数，满足R(k)=R(一k)； 2)k=0时函数取得最大值；
3)如果原序列是周期为丁的周期信号，那么自相关函数也是周期为T的周期函数，即
一55—
万方数据
R(愚)=R(T+k)?．
1．5
基频计算 因为浊音语音的自相关函数具有一定的周期性，在相隔一定的取样后，自相关函数达到 最大值．浊
音语音的周期可用自相关函数中的第一个峰值的位置来估算．得到语音的周期后，根据周期与频率的 倒数关系就可容易得到语音频率．
1．6 后处理
没有任何一种预测方法能够准确计算出所有的基音频率．往往大部分的点是比较准确的，但基音 频率轨迹中总是会有一个或几个基音频率的估计值偏离了正常的轨迹(通常是偏离到正常值的两倍或 1／2)，称这种偏离点为基音轨迹的“野点”．此时为了去除这些野点，可以采用各种平滑算法，其中最常 用 的是中值平滑算法、线性平滑算法及组合平滑算法．此处采用的是中值平滑处理．
中值平滑处理的基本原理是：设z(ft)为输入信号，Y(咒)为中值滤波器的输出，采用一个滑动窗， 输 出值就是将窗外的中心移到咒处时窗口输入样点的中值，及扎点的左右各取L个样点．连同被平滑点 中心共同构成一组信号采样点(共2L+1个样值)，然后将这(2+1)个样点按大小序列排成一排，此序 列中中间者作为平滑器的输出．L值一般为1或2，即中值平滑处的“窗口”一般套住3个或5个样点．中 值平滑的优点是既可以有效地除去少量的野点，又不会破坏基音周期轨迹中两个平滑段之间的阶跃性 变化[4I．
2算法实现与实验结果 Matlab作为mathworks公司的高性能的数值计算和可视化软件，集成数值分 析、矩阵计算、信号运
算、信号处理和图形显示于一体，构成了一个方便且界面友好的用户环境．以下，我们在voicebox工具 箱的帮助下完成对一段具体语音的基频提取．
2．1带通滤波
图1为语音信号波形图，图2为滤波前后波形图对照．很明显，滤波后谐振峰的影响被削弱了， 语 音波形得到了改善，而且信号能量比较集中．
图1语音信号波形图 图2滤波前后波形图对照
2．2三电平削波 三电平消波 Matlab程序如下： for JJ=1
iN
if(data(jj)>Amax)
data(jj)=1；
一56 一
万方数据
elseif(data(jj)<一Amax)
data(jj)=一l；
elseif(data(jj)(=Amax&&data(jj))=一Amax)
data(jj)=0；
end end
其中Amax取最大幅值的0．65，data存放的是加窗之后的数据，长度为一帧．
其中一帧语音经Matlab编程削波处理后的语音波形如图3所示． 图3削波前后波形图对照 图4基颏轨迹
2．3 自相关求基频 经过三电平滤波后，自相关函数的计算由繁复的乘法运算变为简单的逻辑组合，其 Matlab实现如 下：
for k=1：N一1 ％求自相关序列
for jj=1：N—k
R(k)=R(k)+data(jj)*data(jj+k)；
end end 其中N为
每一帧的长度． 由信号周期转换为频率的Matlab程序为：value(ii)=Fs／find-maxn(R)． 其中Fs为采样频率，读取语音时获得，find—maxn(R)为寻找最大峰值位置的函数，其具体实现为： zer=find(r==0)；％找第一个零点如果存在
clozer；％找第一近零点 ii=1： while(dozer<=O)
if(r(ii)>0&&r(ii+1)<0&&(ii+1) dozer=ii； end ii=ii十 l： if ii==length(r)％没有找到符合要求的点 dozer21； end end if length(zer)>0％检查是否存在零点 一57— 万方数据