语音端点检测（1）：双门限法（简单教学版）

为什么要有语音端点检测？或者换个角度说，静默检测、静音检测。

以下摘自百度。

语音活动检测(Voice Activity Detection,VAD)又称语音端点检测,语音边界检，是指在噪声环境中检测语音的存在与否,通常用于语音编码、语音增强等语音处理系统中,起到降低语音编码速率、节省通信带宽、减少移动设备能耗、提高识别率等作用。早先具有代表性的VAD方法有ITU-T的G.729 Annex B。

似乎很重要。

那怎么做？

最直接的，从能量的角度考虑。

如果大于某个阈值，就认为是语音段。当然，还要排除个别突发噪声的情况。

另外，声母的功率小，但也属于语音段啊。怎么办？考虑到它的ZCR比较大。

所以有个教学版的算法，简单，效果又挺好。

假定一段信号的前一段是纯净的噪声（噪声也不是没有用，还是有有用的时候）……

双门限的VAD算法，到处都有。百度一下就可以知道了。

以下贴出对应的matlab代码。

非常值得学习的是，3 sigma 原则。（sigma为标准差std）

测量环境噪声参数那一段时，用了个3 sigma的思想，对门限进行了选取。个人觉得挺巧妙的。

%% 语音静默检测
% 内容：时域特征+双阈值进行语音静默检测
% 作者：qcy

% 版本：v1.1
% 开始时间：2016年11月1日9:01:54
% 结束时间：2016年11月1日10:13:49
% 封装 get_st_zcr 函数
% 封装 get_st_energy 函数

% 版本：v1.0
% 参考文献：Rabiner L R, Schafer R W. Introduction to digital speech processing[J].
% Foundations & Trends in Signal Processing, 2007, 1(1):1-194.
% 开始时间：2016年10月24日21:10:58

clear;
close all;
clc

%% 导入文件
[x,fs] = audioread('bluesky1.wav');
% [x,fs] = audioread('静默检测.m4a');
x = x - mean(x);
x = x./max(abs(x));
T_duration = length(x)/fs;
t = 0:1/fs:(T_duration-1/fs);
% sound(x,fs);
figure(1);
subplot(411)
plot(t,x);
title('信号波形');

%% 计算短时过零率 short-term zero-crossing rate
% 1. 设置帧长、步长
wlen_time = 0.02; % [s]
step_time = 0.01; % [s]
wlen = round(wlen_time*fs);
nstep = round(step_time*fs);
nframes = fix((length(x)-wlen)/nstep)+1; % 帧数
frame_time = frame2time(nframes, wlen, nstep, fs); % 计算每帧对应的时间

Zr = get_st_zcr(x,fs,wlen_time,step_time);
figure(1);
subplot(412)
plot(frame_time,Zr);
title('短时过零率');

%% 计算短时能量 short-term energy
Er = get_st_energy( x,fs,wlen_time,step_time,'hamming','dB' );
figure(1);
subplot(413)
plot(frame_time,Er);
% plot(Er);
title('短时能量');

%%%%%%%%%%%%%%%%%%%%%     以下为静默检测    %%%%%%%%%%%%%%%%%%%%%
%% 1. 获取噪声特性，定义静默检测中所需要的一些参数
% (1) 测量环境噪声，获取一些噪声的参数
noise_frame_idx=floor((.1-wlen_time)/step_time)+1;  % 假定前100ms是环境噪声
eavg=mean(Er(1:noise_frame_idx)); % 计算环境噪声的 平均短时功率
esig=std(Er(1:noise_frame_idx)); % 计算环境噪声的 短时功率的标准差
zcavg=mean(Zr(1:noise_frame_idx)); % 计算环境噪声的 平均短时过零率
zcsig=std(Zr(1:noise_frame_idx));% 计算环境噪声的 平均过零率的标准差

% (2) 根据背景噪声设置阈值
IF=35;
IZCT=max([IF zcavg+3*zcsig]); % 均值往上加3倍标准差，作为短时过零率的阈值
ITU=-15; % constant in the range [-10 -20] dB (intensity threshold upper)
ITR=max([ITU-10 eavg+3*esig]); % 均值往上加3倍标准差，作为短时功率的阈值

% 2. 开始静默检测
B1 = 1; % 起点粗估计 量纲是帧序号
B2 = 1; % 起点精估计 量纲是帧序号
E1 = length(Er); % 终点粗估计 量纲同上
E2 = length(Er); % 终点精估计 量纲同上

%% (1) 起点和终点的粗估计 --> 能量阈值
% (1a) 粗估计：起始点
c = 1; % 游标
is_continue = 1; % 为了排除突发噪声
while is_continue
    while Er(c) < ITR % 从前往后扫，跳过所有低功率的帧
        if c>length(Er) % 合法性检查，不要越界
            break;
        end
        c = c+1;
    end
    is_continue = 0; % 扫描结束
    B1 = c; % 此时作为起点的粗估计
    
    % 为了防止是突发噪声，还要往后走M帧，看是否真正是语音
    M = 3;
    for k = c+1:c+M
        if k>length(Er) % 合法性检查，不要越界
            break;
        end
        if Er(k) < ITR % 后面一帧能量就下去了，说明前一帧多半是噪声
            % 所以必须有一个外层循环，让c跳过噪声的帧
            c = k+1;
            is_continue = 1;
        end
    end
end

% (1b) 粗估计：终止点
c = length(Er); % 游标
is_continue = 1;
while is_continue  % 为了排除突发噪声
    while Er(c) < ITR % 从后往前扫，跳过所有低功率的帧
        if c < 1 % 合法性检查，不要越界
            break;
        end
        c = c-1;
    end
    is_continue = 0; % 扫描结束
    E1 = c; % 此时作为起点的粗估计
    
    % 为了防止是突发噪声，还要往前走M帧，看是否真正是语音
    M = 3;
    for k = c-1:c-M
        if k < 1 % 合法性检查，不要越界
            break;
        end
        if Er(k) < ITR % 前面一帧能量就小于阈值了，说明刚刚后一帧多半是噪声
            % 所以必须有一个外层循环，让c跳过噪声的帧
            c = k-1;
            is_continue = 1;
        end
    end
end

%% (2) 起点和终点的精估计 --> 短时zcr阈值
% (2a) 精估计：起始点
B2 = B1;
% 累计从B1往左M_left帧，zcr大于zcr阈值的帧数 counter
% 如果这个counter确实大于某个数M_counter，就认为前一段也是语言
% 否则B2就认为B1
M_left = 20;
M_counter = 4;
nframes_zcr_higher_than_threshold_counter = 0;
possible_B2_idx = 0; % 起点精估计的可能的下标
for k = (B1-1):-1:(B1-M_left)
    if k < 1 % 合法性检查，不要越界
        break;
    end
    if Zr(k) > IZCT
        nframes_zcr_higher_than_threshold_counter = ...
            nframes_zcr_higher_than_threshold_counter+1;
        possible_B2_idx = k; 
    end
end
if nframes_zcr_higher_than_threshold_counter > M_counter
    B2 = possible_B2_idx;
end

% (2b) 精估计：终止点
E2 = E1;
% 累计从B1往左M_left帧，zcr大于zcr阈值的帧数 counter
% 如果这个counter确实大于某个数M_counter，就认为前一段也是语言
% 否则B2就认为B1
M_right = 20;
M_counter = 4;
nframes_zcr_higher_than_threshold_counter = 0;
possible_E2_idx = 0; % 终点精估计的可能的下标
for k = E1+1:E1+M_right
    if k > length(Zr) % 合法性检查，不要越界
        break;
    end
    if Zr(k) > IZCT
        nframes_zcr_higher_than_threshold_counter = ...
            nframes_zcr_higher_than_threshold_counter+1;
        possible_E2_idx = k; 
    end
end
if nframes_zcr_higher_than_threshold_counter > M_counter
    B2 = possible_E2_idx;
end

%% (3) 保守的起点和终点的再计算。
% 保守派：宁可多留一点噪声，也绝不漏掉一帧语音
% B2再往前走M_left2帧，E2再往后走M_right2帧
% 如果这些帧的能量还大于能量阈值
% 就认为这些还是语音段

% (3a) 起点
M_left2 = 20;M_right2 = 20;

for k = (B2-1) :-1: (B2-M_left2)
    if k < 1 % 合法性检查，不要越界
        break;
    end
    if Er(k) > ITR
        B2 = k;
    end
end

% (3b) 终点

for k = E2+1 : E2+M_right2
    if k > length(Er) % 合法性检查，不要越界
        break;
    end
    if Er(k) > ITR
        E2 = k;
    end
end

% 至此，B、E中分别保存的是起始、终止的帧序号
% 是否需要换算成秒
B_time = (B2 - 1) * step_time + wlen_time/2;
E_time = (E2 - 1) * step_time + wlen_time/2;

%%%%%%%%%%%%%    检测结束    %%%%%%%%%%%%%

%% 画图
figure(1);
subplot(414)
plot(t,x);
title('静默检测结果');
hold on;
line([B_time B_time], [-1 1], 'Color', [1 0 0], 'LineWidth', 2); 
line([E_time E_time], [-1 1], 'Color', [1 0 0], 'LineWidth', 2);

%% 听效果
x_speech = x(round(B_time*fs):round(E_time*fs));
sound(x_speech,fs)

其中，计算短时能量、短时过零率的函数，已经在前文中讲了。

短时能量、短时过零率

本算法非常简单且有效。但不能检测出语音内的一段停顿。

检测结果如下图。

如果要改进，估计有点繁琐。

此外，本算法只能用于后处理（post processing）、离线处理（offline processing）。

要读完所有数据，才能检测语音的结尾。因为它的结尾首先是从右往左扫描的。

不过，简单嘛。可以接受。

补充：frame2time

之前的CSDN有代码片，现在新版系统好像找不到代码片了。。。之前的代码片引用好像也全部不见了！

-_-!! 什么东西嘛……

function frameTime=frame2time(frameNum,framelen,inc,fs)
% 分帧后计算每帧对应的时间
frameTime=(((1:frameNum)-1)*inc+framelen/2)/fs;

function zcr = get_st_zcr( x,fs,wlen_time,step_time,win_type )
%function zcr = get_st_zcr(x,fs,wlen_time,step_time,win_type )
%   获取短时过零率。
%   输入参数
%           x：语音信号 --> 单声道
%           fs：采样速率
%           wlen_time：窗口时间（s）
%           step_time：步进时间（s）
%           win_type：'hamming','hanning',...,默认'hamming'
%   返回参数
%           zcr：短时过零率（横坐标是帧序号）
%
% 作者：qcy
% 版本：v1.0
% 版本说明：计算短时过零率。
% 如果分帧时，不能整除，则抛弃最后一帧，不予以计算
% 时间：2016年10月31日21:08:22

if(min(size(x))>1) % 如果不是单声道
    % ...
end

wlen = round(wlen_time * fs);
nstep = round(step_time * fs);

if nargin < 5
    win = hamming(wlen);
elseif narmin == 5
    if strcmp(win_type, 'hamming')
        win = hamming(wlen);
    elseif strcmp(win_type, 'hanning')
        win = hanning(wlen);
    else
        win = hamming(wlen);
    end
else
    win = hamming(wlen);
end


nFrames = floor((length(x) - wlen)/nstep) + 1; % 总帧数
zcr = [];

for k = 1:nFrames
    idx = (k-1) * nstep + (1:wlen);
    x_sub = x(idx) .* win;
    x_sub1 = x_sub(1:end-1);
    x_sub2 = x_sub(2:end);
    zcr(k) = sum(abs(sign(x_sub1) - sign(x_sub2))) / 2 / length(x_sub1); 
end

end

function E = get_st_energy( x,fs,wlen_time,step_time,win_type,energy_unit )
%function zcr = get_st_energy( x,fs,wlen_time,step_time,win_type,energy_unit )
%   获取短时能量（没有除以帧长，所以不是计算的功率）。
%   输入参数
%           x：语音信号 --> 单声道
%           fs：采样速率
%           wlen_time：窗口时间（s）
%           step_time：步进时间（s）
%           win_type：'hamming','hanning',...,默认'hamming'
%           energy_unit：'dB'，以归一化的能量显示 （单位：dB）。否则是线性刻度。
%   返回参数
%           E：短时能量（横坐标是帧序号）
% 作者：qcy
% 版本：v1.0
% 版本说明：计算短时能量。
% 如果分帧时，不能整除，则抛弃最后一帧，不予以计算
% 时间：2016年10月31日21:21:23

wlen = round(wlen_time * fs);
nstep = round(step_time * fs);

if nargin < 5
    win = hamming(wlen);
elseif nargin == 5
    if strcmp(win_type, 'hamming')
        win = hamming(wlen);
    elseif strcmp(win_type, 'hanning')
        win = hanning(wlen);
    else
        win = hamming(wlen);
    end
else
    win = hamming(wlen);
end

nFrames = floor((length(x) - wlen)/nstep) + 1; % 总帧数
E = [];

for k = 1:nFrames
    idx = (k-1) * nstep + (1:wlen);
    x_sub = x(idx) .* win;
    E(k) = sum(x_sub.^2); 
end

% 是否需要化成dB
if nargin == 6
    if strcmp(energy_unit, 'dB') 
        E = 10*log10(E/max(E)+eps);
    end
end
end