基于短时自相关函数法的基音周期检测

项目总结报告 基音周期的检测1. 项目整体框架1.1目标了解语音基音周期估计方法，掌握自相关法估计基音周期的原理。1.2主要内容 本次基音周期的估算，我们选用的是短时自相关函数法，包括四个模块。第一个模块为基音的端点检测，主要为了区分浊音和清音。第二个模块为基音检测中的带通滤波器，主要为了减少共振峰的干扰。第三个模块为短时自相关函数法做基音检测，主要为了计算出基音周期。第四个模块为平滑处理，主要为了消除偏离值点。2. 模块一（端点检测）2.1主要负责工作利用能熵比法进行语音端点检测，区分语音帧的起点以及终点。2.2具体实现方法 2.2.1实验步骤1）取一段语音“tone4.wav”，该语音内容是“妈妈，好吗，上马，骂人”，语音长度为3.5秒，采样率Fs=8000. 进行简单的去除直流分量，然后幅值归一化，时域波形如图1所示。2）设置好分帧参数，帧长wlen=320,帧移inc=80，调用函数y =enframe(x,wlen,inc);对语音信号x分帧处理。最后帧数Fn=337。3）设置端点检测门限值T1=0.05，使用能熵比法进行端点检测。对分帧后的语音y 每一帧进行FFT运算，然后计算每一帧的能熵比值。从而计算出语音y中的语音端点。结果如图2所示。2.2.2能熵比法 设语音信号时域波形为X(n),加窗分帧处理后得到的第i帧语音信号为Xi(m),则FFT后表示为Xi(k),其中下标i表示为第i帧，而k表示为第k条谱线。该语音帧在频域中的短时能量为Ei=k=0n2Xik*Xi(k)式中，N为FFT的长度，只去正频率部分。而对于某一谱线k的能量谱为Yik=Xik*Xi(k),则每个频率分量的归一化谱概率密度函数定义为Pik=Yi(k)l=0N/2Yi(l)=Yi(k)Ei k=0,1,N-1该语音帧的短时谱熵定义为Hi=-k=0N/2Pi(k)logPi(k)其中，只取正频率部分的谱熵，对应的能熵比表示为EEFi=1+|Ei/Hi|2.2.3代码编写在主程序中，执行的是voiceseg,vosl,SF,Ef=pitch_vad1(y,fn,T1);而调用了以下函数function voiceseg,vosl,SF,Ef=pitch_vad1(y,fn,T1,miniL)if nargin4, miniL=10; endif size(y,2)=fn, y=y; end % 把y转换为每列数据表示一帧语音信号wlen=size(y,1); % 取得帧长for i=1:fn Sp = 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里）endhindex=find(H=T1); % 寻找Ef中大于T1的部分zseg=findSegment(zindex); % 给出端点检测各段的信息zsl=length(zseg); % 给出段数j=0;SF=zeros(1,fn);for k=1 : zsl % 在大于T1中剔除小于miniL的部分 if zseg(k).duration=miniL j=j+1; in1=zseg(k).begin; in2=zseg(k).end; voiceseg(j).begin=in1; voiceseg(j).end=in2; voiceseg(j).duration=zseg(k).duration; SF(in1:in2)=1; % 设置SF（有话段的SF帧为1） endendvosl=length(voiceseg); % 有话段的段数 2.2.4实验结果图1 “妈妈，好吗，上马，骂人”语音的时域波形图2 图中实线代表语音起始点，虚线代表语音终点2.2.5结果分析基音周期是语音信号的重要参数之一，它描述了语音激励源的一个重要特征。语音的头、尾部并不能具有声带振动那样的周期性，也就是检测不到相应的基音周期。如果要进行基音周期的计算，首要的必然是检测语音端点，将有话段从整段语音中分离出来，才能开始对语音的基音周期进行计算.利用能熵比法检测语音端点后，该语音段可分为四段，每段的起点以及终点如图二所示.2.3总结语音的预处理是非常关键的一环，做好这个实验的关键在于理解matlab语句的作用，录入语音信号，设定好参数，进行分帧处理，并通过能熵值的大小标比较，确定好语音信号的端点。通过这次的检测，我了解了Matlab函数库非常的多，理解起来需要查阅很多公式以及资料，非常的耗时。做实验务必讲究认真专心，一步一个脚印的去做好每一个环节就可以了，与队友的配合以及沟通也是重要的一环，及时了解队友的情况，做好衔接。3. 模块二（基音检测中的带通滤波器） 3.1主要负责工作用带通滤波器减少共振峰的干扰。 3.2具体实现方法 3.2.1算法原理 在使用相关法和AMDF法的基音检测之前常用到低通滤波器和带通滤波器，其主要目的是减少共振峰的干扰。在文献中作者选用了500Hz作为滤波器的上限频率，并指出选择截止频率高不利于减少噪声和共振峰的影响。所以在本书的基音检测中的预滤波器选择的带宽为60500Hz，高频截止频率选择500Hz，是因为基频区间的高端就在这个区域中，低频截止频率选择60Hz是为了减少工频和低频噪声的干扰。 我们选用IIR滤波器中的椭圆滤波器，因为IIR滤波器的运算量比FIR少，当然IIR滤波器会带来延迟，也就是相位的变化，但语音信号是对相位不敏感的信号；又选用椭圆滤波器，因为它在经典滤波器设计中相同过渡带和带宽条件下，需要的阶数比较小。在阶数相同的条件下,椭圆滤波器相比于其他类型的滤波器，能获得更窄的过渡带宽和较小的阻带波动, 就这点而言, 椭圆滤波器是最优的。滤波器的要求为采样频率8000Hz，通带是60500Hz，通带波纹为1dB，阻带分别为30Hz和2000Hz，阻带衰减为40dB。 3.2.2代码编写% pr8_1_1 clear all; clc; close all; fs=8000; fs2=fs/2; % 采样频率Wp=60 500/fs2; % 滤波器通带Ws=20 2000/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.6f %5.6fn,b)fprintf(a=%5.6f %5.6f %5.6f %5.6f %5.6f %5.6f %5.6fn,a) db, mag, pha, grd,w=freqz_m(b,a); % 求取频率响应曲线plot(w/pi*fs/2,db,k); % 作图grid; ylim(-90 10);xlabel(频率/Hz); ylabel(幅值/dB); title(椭圆6阶带通滤波器频率响应曲线);3.2.3实验结果由程序计算出滤波器系数为：b=0.012280 -0.039508 0.042177 0.000000 -0.042177 0.039508 -0.012280a=1.000000 -5.527146 12.854342 -16.110307 11.479789 -4.410179 0.713507实验结果图：图3 滤波前和后的语音信号3.2.4结果分析带通滤波器使原始信号变得更加平整，滤去了杂音，为基音检测减少干扰，使检测更加顺利进行。 3.3总结通过这次课程实践，我深刻体会到自己在课程理论方面的不足，MATLAB的应用也相当生疏。在设计过程中，看不懂的东西很多，但是学会了通过上网搜查资料和自主从书本寻找知识。实践中遇到许多问题，但是和同学分工合作与讨论让我学会了很多平常学不到的知识。这次实践让我明白知识就像大海一样宽广，我们要像海绵吸水一样吸取有用的知识。4模块三（短时自相关函数法做基音检测） 4.1主要负责工作 在前面两个步骤之后用自相关函数法去实现基音周期的检测 4.2具体实现方法 4.2.1算法设语音信号的时间序列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)。短时自相关函数基因检测的主要原理大都是利用它的这些性质，通过比较原始信号和它的延迟后的信号之间的类似性质来确定基音周期。如果延迟量等于基音周期，那么两个信号具有最大类似性，或是直接找出短时自相关函数的两个最大值间的距离，即为基音周期的初估值。在用短时自相关函数检测基音是，常用归一化的短时自相关函数，表达式为 r(k)= R(k)/ R(0)以上性质中已指出，k=0，R(0)为最大值。所以r(k)的模值永远小于或等于1。已知基音频率范围为60Hz500Hz之间，采样率为fs时，则基音周期（样点值）范围在fs/500fs/60之间，然后在这个范围之间寻找归一化相关函数的最大值，对应的延迟量就是基音周期。4.2.2代码编写步骤 根据前面步骤所得分帧后的语音数组，对它进行处理，把它转换为每一列数据表示一阵语音信号，行数为总的帧长。 对每一帧语音信号的基音周期进行初始化，让其为零保证清音段无基音周期。 利用前面端点检测后所得的有话段和有话段段数，对每一段的有话段的每一帧数据计算归一化自相关函数。 最后对每一帧所得的归一化自相关函数取正延迟量的部分，然后在基音频率范围内得最大值，从而可得最大值对应的样点数，利用函数循环计算，最后可以得出每段有话段内每一帧的基音周期。4.2.3代码编写用短时自相关函数法的程序：function period=ACF_corr(y,fn,vseg,vsl,lmax,lmin)pn=size(y,2);if pn=fn, y=y; end % 把y转换为每列数据表示一帧语音信号wlen=size(y,1); % 取得帧长period=zeros(1,fn); % 初始化for i=1 : vsl % 只对有话段数据处理 ixb=vseg(i).begin; ixe=vseg(i).end; ixd=ixe-ixb+1; % 求取一段有话段的帧数 for k=1 : ixd % 对该段有话段数据处理 u=y(:,k+ixb-1); % 取来一帧数据 ru= xcorr(u, coeff); % 计算归一化自相关函数 ru = ru(wlen:end); % 取延迟量为正值的部分 tmax,tloc=max(ru(lmin:lmax); % 在基音频率范围内寻找最大值 period(k+ixb-1)=lmin+tloc-1; % 给出对应最大值的延迟量 endend4.2.4实验结果 图4 端点检测的语音信号 图5 有话段的基音周期4.2.5结果分析在端点检测和滤波后，经过短时自相关函数法处理后，由上面的实验结果图可看出基音周期有很多的偏离点即野点，需要进行平滑处理，也可以看出基音周期大约为40个样点数，即基音周期为40/fs,采样率为8000Hz，得周期为0.005s,基音频率为200Hz.4.3总结通过这次设计，我深刻体会到在信号处理方面我的严重不足，MATLAB的应用也相当生疏。在设计过程中，很多东西都看不懂，对原理也一塌糊涂。搜的一些材料也没有发挥他们的作用。有些公式原理看懂了可程序又不行了，看不出每一步的作用及实现的功能。设计中遇到的问题通过同学的交流也解决了一部分，总的来说在这次设计中我还是学到了一些知识，即有独立完成设计的喜悦又有对自己不足的担忧。在校学习的机会也不多了，我会在余下的时间里努力学习，改正自己的缺点，弥补自己。5.模块四（平滑处理）5.1主要负责工作主要负责第四个模块，把中值滤波算法和线性平滑算法相结合来去除偏离值点，进行后处理。5.2具体实现方法5.2.1算法中值滤波算法：在被平滑点的左右各取L个样点。连同被平滑点共同构成一组信号采样值（共（2L+1）个样值），然后将这（2L+1）个样值按大小次序排成一对，取此队列中间者作为平滑器的输出。L值一般取为1或2，即中值平滑的“窗口”一般套住3或5个样值。中值平滑的优点是既可以有效地去除少量的野点，又不会破坏基音周期轨迹中两个平滑段之间的阶跃性变化。matlab函数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+v iN v=(m-1)/2 ；Yi称为序列fi-v,fi-1,fi,fi+1,fi+v的中值线性滤波算法让基音轨迹序列通过一个“低通”FIR滤波器。由于基音轨迹中的野点都是一些突发跳动点，它们具有很高的“频率分量”，采用“低通”滤波即可以将其去除。线性平滑是用滑动窗口进行线性滤波处理，即yn=m=-LLxn-mw(m)式中，w(m)，m=-L,-L+1,0,1,2,L为（2L+1）点平滑窗，满足m=-LLwm=15.2.2 实验步骤 在模块三得到的基音周期（图6）上进行平滑处理； 先单独进行中值滤波处理，得到图7； 再单独进行线性平滑处理，得到图8； 对模块三得到的原始基音周期（图6）进行中值滤波处理和线性平滑结合处理，得到图9； 对图6到图9进行对比并分析。5.2.3 代码编写在主程序中，执行的是T0= pitfilterm1(period,voiceseg,vosl);而调用了以下函数function y=pitfilterm1(x,vseg,vsl)y=zeros(size(x); % 初始化for i=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; % 赋值给yendfunction y = linsmoothm(x,n)if nargin 2 n=3;endwin=hanning(n); % 用hanning窗win=win/sum(win); % 归一化 x=x(:); % 把x转换为行序列len=length(x);y= zeros(len,1); % 初始化y% 对x序列前后补n个数,以保证输出序列与x相同长if mod(n, 2) =0 l=n/2; x = ones(1,1)* x(1) x ones(1,l)* x(len);else l=(n-1)/2; x = ones(1,1)* x(1) x ones(1,l+1)* x(len);end% 线性平滑处理for k=1:len y(k) = win* x(k:k+ n- 1);end当只进行中值滤波平滑处理时，调用的函数改为：function y=pitfilterm1(x,vseg,vsl)y=zeros(size(x); % 初始化for i=1 : vsl % 有段数据 ixb=vseg(i).begin; % 该段的开始位置 ixe=vseg(i).end; % 该段的结束位置 u0=x(ixb:ixe); % 取来一段数据 y0=medfilt1(u0,5); % 5点的中值滤波 y(ixb:ixe)=y0; % 赋值给yend当只进行线性平滑处理时，调用的函数改为：function y=pitfilterm1(x,vseg,vsl)y=zeros(size(x); % 初始化for i=1 : vsl % 有段数据 ixb=vseg(i).begin; % 该段的开始位置 ixe=vseg(i).end; % 该段的结束位置 u0=x(ixb:ixe); % 取来一段数据 v0=linsmoothm(u0,5); % 线性平滑 y(ixb:ixe)=v0; % 赋值给yendfunction y = linsmoothm(x,n)if nargin 2 n=3;endwin=hanning(n); % 用hanning窗win=win/sum(win); % 归一化 x=x(:); % 把x转换为行序列len=length(x);y= zeros(len,1); % 初始化y%