matlab 端点检测 能零比法_基于短时能量与过零率的端点检测的matlab分析

基于语音端点检测的方法有很多，从历史的发展来看。

首先是基于短时能量和短视过零率的端点检测=〉各变换域=〉人工神经网络=〉基于倒谱距离的检测算法=〉基于谱熵的方法=〉几何门限的方法=〉sigma函数=〉近些年的从分形技术和混沌理论引入的端点检测。

作为最早的短时能量与过零率的检测方法，比较简单，当然也有很大的缺陷：在SNR比较低的情况下，准确度大打折扣，严重的话会失效。以目前而言，先研究这个简单的再说。

1，实现原理。

1.1基于两个公式，即短时能量和短时过零率。

1.2理论基础：语音信号一般可分为无声段、清音段和浊音段。无声段是背景噪声段,平均能量最低;浊音段为声带振动发出对应的语音信号段,平均能量最高;清音段是空气在口腔中的摩擦、冲击或爆破而发出的语音信号段,平均能量居于前两者之间。清音段和无声段的波形特点有明显的不同,无声段信号变化较为缓慢,而清音段信号在幅度上变化剧烈,穿越零电平次数也多。经验表明,通常清音段过零率最大。端点检测就是首先判断/有声0还是/无声0,如果有声,则还要判断是/清音0还是/浊音0。为正确地实现端点检测,一般综合利用短时能量和过零率两个特征,采用/双门限检测法0。

1.3基本思路：根据信号设置三个阈值：能量阈值，TL，TH；过零率阈值ZCR,当某帧信号大于TL或者大于ZCR时，认为信号的开始、起点，当大于TH时，则认为正式的语音信号，如果保持一段时间，则确认这信号即所需信号。

2，matlab实现

2.1 matlab语句分析。

function [x1,x2] = vad(x)

%幅度归一化到[-1,1]

x = double(x);

x = x / max(abs(x));

%常数设置

FrameLen = 240;%帧长为240点

FrameInc = 80;%帧移为80点

amp1 = 10;%初始短时能量高门限

amp2 = 2;%初始短时能量低门限

zcr1 = 10;%初始短时过零率高门限

zcr2 = 5;%初始短时过零率低门限

maxsilence = 8;  % 8*10ms  = 80ms

%语音段中允许的最大静音长度，如果语音段中的静音帧数未超过此值，则认为语音还没结束；如果超过了

%该值，则对语音段长度count进行判断，若count

%状态0；若count>minlen，则认为语音段结束；

minlen  = 15;    % 15*10ms = 150ms

%语音段的最短长度，若语音段长度小于此值，则认为其为一段噪音

status  = 0;     %初始状态为静音状态

count   = 0;     %初始语音段长度为0

silence = 0;     %初始静音段长度为0

%计算过零率

tmp1  = enframe(x(1:end-1), FrameLen, FrameInc);%分帧，所得矩阵为fix((x(1:end-1)-FrameLen+FrameInc)/FrameInc)*FrameLen

tmp2  = enframe(x(2:end)  , FrameLen, FrameInc);%分帧，所得矩阵为fix((x(2:end)-FrameLen+FrameInc)/FrameInc)*FrameLen

signs = (tmp1.*tmp2)<0;%tmp1.*tmp2所得矩阵小于等于零的赋值为1，大于零的赋值为0

diffs = (tmp1 -tmp2)>0.02;%tmp1-tmp2所得矩阵小于0.02的赋值为0，大于等于0.02的赋值为1

zcr   = sum(signs.*diffs, 2);

%计算短时能量

%amp = sum(abs(enframe(filter([1 -0.9375], 1, x), FrameLen, FrameInc)), 2);

amp = sum(abs(enframe(x, FrameLen, FrameInc)), 2);

%调整能量门限

amp1 = min(amp1, max(amp)/4);

amp2 = min(amp2, max(amp)/8);

%开始端点检测

x1 = 0;

x2 = 0;

for n=1:length(zcr) %length(zcr)得到的是整个信号的帧数

goto = 0;

switch status

case {0,1}                   % 0 = 静音, 1 = 可能开始

if amp(n) > amp1          % 确信进入语音段

x1 = max(n-count-1,1);

status  = 2;

silence = 0;

count   = count + 1;

elseif amp(n) > amp2 | ... % 可能处于语音段

zcr(n) > zcr2

status = 1;

count  = count + 1;

else                       % 静音状态

status  = 0;

count   = 0;

end

case 2,                       % 2 = 语音段

if amp(n) > amp2 | ...     % 保持在语音段

zcr(n) > zcr2

count = count + 1;

else                       % 语音将结束

silence = silence+1;

if silence 

count  = count + 1;

elseif count 

status  = 0;

silence = 0;

count   = 0;

else                    % 语音结束

status  = 3;

end

end

case 3,

break;

end

end

count = count-silence/2;

x2 = x1 + count -1;

subplot(311)    %subplot(3,1,1)表示将图排成3行1列，最后的一个1表示下面要画第1幅图

plot(x)

axis([1 length(x) -1 1])    %函数中的四个参数分别表示xmin,xmax,ymin,ymax，即轴的范围

ylabel('Speech');

line([x1*FrameInc x1*FrameInc], [-1 1], 'Color', 'red');

%这里作用为用直线画出语音段的起点和终点，看起来更直观。第一个[]中的两个参数为线起止点的横坐标，

%第二个[]中的两个参数为线起止点的纵坐标。最后两个参数设置了线的颜色。

line([x2*FrameInc x2*FrameInc], [-1 1], 'Color', 'red');

subplot(312)

plot(amp);

axis([1 length(amp) 0 max(amp)])

ylabel('Energy');

line([x1 x1], [min(amp),max(amp)], 'Color', 'red');

line([x2 x2], [min(amp),max(amp)], 'Color', 'red');

subplot(313)

plot(zcr);

axis([1 length(zcr) 0 max(zcr)])

ylabel('ZCR');

line([x1 x1], [min(zcr),max(zcr)], 'Color', 'red');

line([x2 x2], [min(zcr),max(zcr)], 'Color', 'red');

2.2由语句提出的程序流程。

由上边的程序可以看出程序流程为：

3，mat lab程序中的部分解释说明

3.1流程图的说明

amp2,amp1为能量的两个阈值，前者为小的，后者为大的，zcr2为过零率的阈值小值，当>amp2 or >zcr2，count开始加1，在此期间如果有不满足该条件的话，count立即为0，回到0状态。如果>amp1时，count加1，然后进入2状态。在2状态里边，当>amp2

or >zcr2，count也1，如果不满足条件，则Silence+1，如果Silence即在2状态期间处于静音状态满足结束时的静音条件，则判断所有计数的信号即count的值是否满足最小的语音信号长度值，如果满足，则找到结束点，否则认为是噪声，重新开始。如果Silence即在2状态期间处于静音状态不足结束时的静音条件，则count继续加1。

3.2起点和终点的判断

即判断X1,X2。根据程序x1=max(n-count-1,1);n为找到>amp1时，此时的贞的序列值。而count为在这之前的>amp2

or >zcr2，的贞的个数。一般突发信号从无慢慢到有，如果N=6，COUNT=2，则起点从第三帧开始。很好理解。

对于X2，有count = count-silence/2;x2 = x1 + count - 1;总count值为>amp2 or

>zcr2开始，到判断超过最小静音为止。而当进入silence加的时候，已经为静音阶段了。当silence=6时，便结束了，然后X2的计算方法，基本了解，但silence/2感觉可以不要也行，后边的-1也感觉可以不要。

3.3 enframe函数的说明

前边算能量，过零率都是基于帧来计算的，而enframe函数是用来把信号进行分帧的。

其代码如下：

function f=enframe(x,win,inc)%定义函数。

nx=length(x(:)); %x(:)的作用是把x给弄成一个向量，x为一行，则变成一列，如果为矩阵，则按每一列的顺序排成一列。得出的nx为序列的数据个数。

nwin=length(win);

if (nwin == 1)

len = win; %如果win中就一个数，则len就=该数，此例中为256个点。即每帧长

else

len = nwin; %如果有多个数，则len=个数。

end

if (nargin < 3) % nargin返回的是函数输入的个数，如果中间有变量，返回的是负值。

inc = len; %也就是说，如果函数enframe的输入只有两个的话，系统就自动赋inc

end

nf = fix((nx-len+inc)/inc);%这个比较关键，nf为分帧的组数，结合下边的可以分析这里

各参数的意义，len为帧长，inc为未覆盖的数据，nx为整个数

据量。假设数据为1：30，len为10，未覆盖为5，则nf=5，5

组，第一组为1，2，……10，第二组为6，7，……15，依次列

推，便可知其缘由，即(nx-len)/inc + 1；

f=zeros(nf,len); %构成以组数为行，帧长为列的矩阵。

indf= inc*(0:(nf-1)).'; %indf为一列nf个数据，即0到nf-1的inc倍，即分好了每幀起点。

inds = (1:len); %构成了长度为len的一行。

f(:) = x(indf(:,ones(1,len))+inds(ones(nf,1),:));

%上一条语句为整个算法的核心部分了，indf(:,ones(1,len))把indf的

第一列扩展了nf*len的矩阵。同理inds(ones(nf,1),:)把inds第一行扩

展为nf*len的矩阵，相加得到

[ 1 2 3……len

inc+1 inc+2 inc+3……inc+len

2*inc+1………………2*inc+len

。

。…………………………………]然后就按照这个矩阵从x中把数据

给选出来，达到分帧的目的。

if (nwin > 1) %nwin大于1的情况就不说了。

w = win(:)';

f = f .* w(ones(nf,1),:);

end

3.4过零率的计算

其语句如下：

tmp1 = enframe(x(1:end-1), FrameLen,inc);

tmp2 = enframe(x(2:end) , FrameLen,inc);

signs = (tmp1.*tmp2)<0;%对于tmp1.*tmp2算出来的矩阵，矩阵中<0的数都为1，其他

为0，后边也是一样的。

diffs = (tmp1 - tmp2)>0.02;

zcr = sum(signs.*diffs,2);

假设数据量为1，2……21，帧长为10，inc为5，则

tmp1为[1 2 3 4 5 6 7 8 9 10

6 7 8 9 10 11 12 13 14 15

11 12 13 14 15 16 17 18 19 20]

tmp2为[2 3 4 5 6 7 8 9 10 11

7 8 9 10 11 12 13 14 15 16

12 13 14 15 16 17 18 19 20 21]

在这里注意一个问题：即数组的乘法与矩阵的乘法是不一样的。数组乘法：A.*A,矩阵乘法：A*A’。前者有‘.’号，算出来的结果是在矩阵A中每一个数据平方，而后者成为另一个数组，行与列相乘然后相加作为一个值。

在这里tmp1.*tmp2为数组相乘，第一个数乘以第二个数，第二个数乘以第三个数，依次，从而判断两者的符号，<0的为1。

然后进行相减，第一个减第二个数，第二个减第三个数……，>0.02,为什么>0.02了？首先得到的signs是真正的过零率，但得限制能量，因为对于噪音的话，也会在0点附近上下摆动，但噪声能量显然是没有语音大的，根据实际情况，所以选择>0.02，其次感觉diffs

= abs((tmp1 - tmp2))>0.02;%?就是加个绝对值，因为对于负值-正值也会满足条件，结果验证也不影响。

3.5能量的计算

语句为：

amp = sum((abs(enframe( x, FrameLen, inc))).^2, 2);

通过对enframe函数的分析，就比较容易了，enframe对x分帧后，绝对值然后平方，最后是sum(x，2)2代表是各列相加最后得到的是一列数据，即各帧的平方和。从中可以看出矩阵处理数据的方便性，一个矩阵就把各帧的结果给弄出来了。

最后运行结果为图：

注：实验过程中发现能量出计算好像有问题，应为下面两句(filter()函数可不要，所以采用下面的，上面注释掉了)，特记，

%amp = sum((abs(enframe(filter([1 -0.9375], 1, x), FrameLen, inc))).^2, 2);

amp = sum((abs(enframe( x, FrameLen, inc))).^2, 2);

转自：http://home.eeworld.com.cn/my/space-uid-436810-blogid-82338.html