%% 语音静默检测
% 内容：时域特征+双阈值进行语音静默检测
% 作者：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:58clear;
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_continuewhile Er(c) < ITR % 从前往后扫，跳过所有低功率的帧if c>length(Er) % 合法性检查，不要越界break;endc = c+1;endis_continue = 0; % 扫描结束B1 = c; % 此时作为起点的粗估计% 为了防止是突发噪声，还要往后走M帧，看是否真正是语音M = 3;for k = c+1:c+Mif k>length(Er) % 合法性检查，不要越界break;endif Er(k) < ITR % 后面一帧能量就下去了，说明前一帧多半是噪声% 所以必须有一个外层循环，让c跳过噪声的帧c = k+1;is_continue = 1;endend
end% (1b) 粗估计：终止点
c = length(Er); % 游标
is_continue = 1;
while is_continue  % 为了排除突发噪声while Er(c) < ITR % 从后往前扫，跳过所有低功率的帧if c < 1 % 合法性检查，不要越界break;endc = c-1;endis_continue = 0; % 扫描结束E1 = c; % 此时作为起点的粗估计% 为了防止是突发噪声，还要往前走M帧，看是否真正是语音M = 3;for k = c-1:c-Mif k < 1 % 合法性检查，不要越界break;endif Er(k) < ITR % 前面一帧能量就小于阈值了，说明刚刚后一帧多半是噪声% 所以必须有一个外层循环，让c跳过噪声的帧c = k-1;is_continue = 1;endend
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;endif Zr(k) > IZCTnframes_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_counterB2 = 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_rightif k > length(Zr) % 合法性检查，不要越界break;endif Zr(k) > IZCTnframes_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_counterB2 = 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;endif Er(k) > ITRB2 = k;end
end% (3b) 终点for k = E2+1 : E2+M_right2if k > length(Er) % 合法性检查，不要越界break;endif Er(k) > ITRE2 = 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)