《MATLAB语音信号分析与合成（第二版）》：第6章语音端点的检测（2）

文档中心

《MATLAB语音 信号分析与合成（第二版）》：第6章语音端点的检测（2）

前言
1. 数据与函数路径设置
2. MATLAB仿真一：语音信号谱距离法端点检测一
3. MATLAB仿真二：语音信号谱距离法端点检测二
4. MATLAB仿真三：语音信号谱距离法端点检测三
5. MATLAB仿真四：语音信号谱熵法端点检测一
6. MATLAB仿真五：语音信号谱熵法端点检测二
7. MATLAB仿真六：语音信号能零比法端点检测
8. MATLAB仿真七：语音信号能熵比法端点检测
9. MATLAB仿真八：语音信号小波变换法端点检测
10. MATLAB仿真九：语音信号EMD分解法端点检测
11. MATLAB仿真十：低信噪比时的端点检测一
12. MATLAB仿真十一：低信噪比时的端点检测二
13. MATLAB仿真十二：低信噪比时的端点检测三
小结

前言

《MATLAB语音信号分析与合成（第二版）》是中科院声学所的大佬宋知用老师数十年经验积累下的呕心之作，对于语音信号处理相关感兴趣的同学，日后希望在语音信号分析、处理与合成相关领域进行一定研究的话，可以以此进行入门。

语音信号处理是数字信号处理的一个重要分支。本书含有许多数字信号处理的方法和 MATLAB函数。全书共10章。第1-4章介绍语音信号处理的一些基本分析方法和手段,以及相应的MATLAB函数;第5-9章介绍语音信号预处理和特征的提取,包括消除趋势项和基本的减噪方法,以及端点检测、基音的提取和共振峰的提取,并利用语音信号处理的基本方法,给出了多种提取方法和相应的 MATLAB程序;第10章结合各种参数的检测介绍了语音信号的合成、语音信号的变速和变调处理,还介绍了时域基音同步叠加( TD PSOLA)的语音合成,并给出了相应的MATLAB程序。附录A中给出了调试复杂程序的方法和思路。本书可作为从事语音信号处理的本科高年级学生、研究生或科研工程技术人员的辅助读物,也可作为从事信号处理研究与应用的科研工程技术人员的参考用书。

我的研究生导师的主攻方向就是语音信号处理相关，虽然自己研究生期间的大论文方向是数字图像处理，但所谓语音图像不分家，自己在老师的研究生主讲课小波变换上虽然划水，但在后期导师的语音信号处理的课程设计和工程应用上自己在语音上还算入了一点小门道，在结课测试中拿到了小组第一，导师还特地发了三百大洋的伙食经费以资鼓励。

这次重新捡起语音识别，正好入手了宋老师的这本书，算是自己重新复习一遍吧，主要以介绍各章节中源码为主，这是本书的第六章的后12个仿真应用实例，话不多说，开始！

1. 数据与函数路径设置

书中经常会调用的一些函数(自编函数或取自其他应用工具箱中的函数)已集中在basic_tbx工具箱中，在运行本书的程序前请把该工具箱设置(用set path设置)在工作路径下；

当要运行EMD处理时，要把emd工具箱设置在工作路径下;

当要运行主体延伸基音检测时，要把Pitch_ztlib工具箱设置在工作路径下;

当要进行时域基音同步叠加语音合成时，要把psola_lib工具箱设置在工作路径下;

当要应用本书提供的语音数据时，最好把speech_signal设置在工作路径下。

本书的所有函数和程序都在MATLAB R2009a版本下调试通过。（我用的是MATLAB2015b，有些函数已经更新，所以我会进行修改，以便调试通过）

路径设置的方法如下：

打开MATLAB，点击“主页”，找到设置路径
在这里插入图片描述

将上述文件夹路径全部添加到MATLAB搜索路径中
在这里插入图片描述

添加完毕，保存，开始仿真。

2. MATLAB仿真一：语音信号谱距离法端点检测一

%% pr6_5_1 clear all; clc; close all;run Set_I   % 基本设置run PART_I  % 读入数据，分帧等准备Y=fft(y);   % FFT变换Y=abs(Y(1:fix(wlen/2)+1,:));     % 计算正频率幅值N=mean(Y(:,1:NIS),2);     % 计算前导无话段噪声区平均频谱NoiseCounter=0;for i=1:fn,     if i<=NIS      % 在前导无话段中设置为NF=1,SF=0 SpeechFlag=0; NoiseCounter=100; SF(i)=0; NF(i)=1;    else    % 检测每帧计算对数频谱距离 [NoiseFlag, SpeechFlag, NoiseCounter, Dist]=vad(Y(:,i),N,NoiseCounter,2.5,8);  SF(i)=SpeechFlag; NF(i)=NoiseFlag; D(i)=Dist;    endendsindex=find(SF==1);% 从SF中寻找出端点的参数完成端点检测voiceseg=findSegment(sindex);vosl=length(voiceseg);% 作图subplot 311; plot(time,x,'k'); title('纯语音波形');ylabel('幅值'); ylim([-1 1]);subplot 312; plot(time,signal,'k');title(['带噪语音 SNR=' num2str(SNR) '(dB)'])ylabel('幅值'); ylim([-1.2 1.2]);subplot 313; plot(frameTime,D,'k'); xlabel('时间/s'); ylabel('幅值'); title('对数频谱距离'); ylim([0 8]);for k=1 : vosl      % 标出语音端点    nx1=voiceseg(k).begin; nx2=voiceseg(k).end;    fprintf('%4d   %4d   %4d\n',k,nx1,nx2);    subplot 311    line([frameTime(nx1) frameTime(nx1)],[-1 1],'color','k','linestyle','-');    line([frameTime(nx2) frameTime(nx2)],[-1 1],'color','k','linestyle','--');    subplot 313    line([frameTime(nx1) frameTime(nx1)],[0 8],'color','k','linestyle','-');    line([frameTime(nx2) frameTime(nx2)],[0 8],'color','k','linestyle','--');end

function [NoiseFlag, SpeechFlag, NoiseCounter, Dist]=vad(signal,noise,NoiseCounter,NoiseMargin,Hangover)%[NOISEFLAG, SPEECHFLAG, NOISECOUNTER, DIST]=vad(SIGNAL,NOISE,NOISECOUNTER,NOISEMARGIN,HANGOVER)%Spectral Distance Voice Activity Detector%SIGNAL is the the current frames magnitude spectrum which is to labeld as%noise or speech, NOISE is noise magnitude spectrum template (estimation),%NOISECOUNTER is the number of imediate previous noise frames, NOISEMARGIN%(default 3)is the spectral distance threshold. HANGOVER ( default 8 )is%the number of noise segments after which the SPEECHFLAG is reset (goes to%zero). NOISEFLAG is set to one if the the segment is labeld as noise%NOISECOUNTER returns the number of previous noise segments, this value is%reset (to zero) whenever a speech segment is detected. DIST is the%spectral distance. %Saeed Vaseghi%edited by Esfandiar Zavarehei%Sep-04if nargin<4    NoiseMargin=3;endif nargin<5    Hangover=8;endif nargin<3    NoiseCounter=0;end    FreqResol=length(signal);SpectralDist= 20*(log10(signal)-log10(noise));SpectralDist(find(SpectralDist<0))=0;Dist=mean(SpectralDist); if (Dist < NoiseMargin)     NoiseFlag=1;     NoiseCounter=NoiseCounter+1;else    NoiseFlag=0;    NoiseCounter=0;end% Detect noise only periods and attenuate the signal     if (NoiseCounter > Hangover)     SpeechFlag=0;    else     SpeechFlag=1; end

%% Set_I.mIS=0.25;    % 设置前导无话段长度wlen=200;   % 设置帧长为25msinc=80;     % 求帧移filedir=[]; % 设置文件路径filename='bluesky1.wav';  % 设置文件名称fle=[filedir filename]    % 构成文件路径和名称SNR=10;     % 设置信噪比

%% PART_I% [xx,fs]=wavread(fle);     % 读入数据[xx,fs]=audioread(fle);     % 读入数据xx=xx-mean(xx);    % 消除直流分量x=xx/max(abs(xx)); % 幅值归一化N=length(x);% 取信号长度time=(0:N-1)/fs;   % 设置时间signal=Gnoisegen(x,SNR);  % 叠加噪声wnd=hamming(wlen); % 设置窗函数overlap=wlen-inc;  % 求重叠区长度NIS=fix((IS*fs-wlen)/inc +1);    % 求前导无话段帧数y=enframe(signal,wnd,inc)';      % 分帧fn=size(y,2);      % 求帧数frameTime=frame2time(fn, wlen, inc, fs);% 计算各帧对应的时间

function soundSegment=findSegment(express)if express(1)==0    voicedIndex=find(express);% 寻找express中为1的位置else    voicedIndex=express;endsoundSegment = [];k = 1;soundSegment(k).begin = voicedIndex(1);     % 设置第一组有话段的起始位置for i=1:length(voicedIndex)-1,if voicedIndex(i+1)-voicedIndex(i)>1,   % 本组有话段结束soundSegment(k).end = voicedIndex(i);      % 设置本组有话段的结束位置soundSegment(k+1).begin = voicedIndex(i+1);% 设置下一组有话段的起始位置  k = k+1;endendsoundSegment(k).end = voicedIndex(end);     % 最后一组有话段的结束位置% 计算每组有话段的长度for i=1 :k    soundSegment(i).duration=soundSegment(i).end-soundSegment(i).begin+1;end

在这里插入图片描述

3. MATLAB仿真二：语音信号谱距离法端点检测二

%% pr6_5_2 clear all; clc; close all;run Set_I     % 基本设置run PART_I    % 读入数据，分帧等准备for i=1 : fn    u=y(:,i); % 取来一帧数据    U(:,i)=rceps(u); % 求取倒谱endC0=mean(U(:,1:5),2); % 计算出前5帧倒谱系数的平均值作为背景噪声倒谱系数的估算值for i=6 : fn  % 从第6帧开始计算每帧倒谱系数与背景噪声倒谱系数的距离    Cn=U(:,i);   Dst0=(Cn(1)-C0(1)).^2;    Dstm=0;    for k=2 :12 Dstm=Dstm+(Cn(k)-C0(k)).^2;    end    Dcep(i)=4.3429*sqrt(Dst0+2*Dstm);   % 倒谱距离endDcep(1:5)=Dcep(6);Dstm=multimidfilter(Dcep,10);    % 平滑处理dth=max(Dstm(1:(NIS)));   % 阈值计算T1=1*dth;T2=1.5*dth;[voiceseg,vsl,SF,NF]=vad_param1D(Dstm,T1,T2);% 倒谱距离双门限的端点检测% 作图subplot 311; plot(time,x,'k');title('纯语音波形');ylabel('幅值'); axis([0 max(time) -1 1]);subplot 312; plot(time,signal,'k');title('加噪语音波形(信噪比10dB)');ylabel('幅值'); axis([0 max(time) -1 1]);subplot 313; plot(frameTime,Dstm,'k');title('短时倒谱距离值'); axis([0 max(time) 0 1.2*max(Dstm)]);xlabel('时间/s'); ylabel('幅值'); line([0,frameTime(fn)], [T1 T1], 'color','k','LineStyle','--');line([0,frameTime(fn)], [T2 T2], 'color','k','LineStyle','-');for k=1 : vsl      % 标出语音端点    nx1=voiceseg(k).begin; nx2=voiceseg(k).end;    fprintf('%4d   %4d   %4d\n',k,nx1,nx2);    subplot 311;     line([frameTime(nx1) frameTime(nx1)],[-1 1],'color','k','LineStyle','-');    line([frameTime(nx2) frameTime(nx2)],[-1 1],'color','k','LineStyle','--');    subplot 313;     line([frameTime(nx1) frameTime(nx1)],[0 1.2*max(Dstm)],'color','k','LineStyle','-');    line([frameTime(nx2) frameTime(nx2)],[0 1.2*max(Dstm)],'color','k','LineStyle','--');end

function [voiceseg,vsl,SF,NF]=vad_param1D(dst1,T1,T2)fn=size(dst1,2);  % 取得帧数maxsilence = 8;   % 初始化  minlen  = 5;    status  = 0;count   = 0;silence = 0;%开始端点检测xn=1;for n=2:fn   switch status   case {0,1}      % 0 = 静音, 1 = 可能开始      if dst1(n) > T2     % 确信进入语音段  x1(xn) = max(n-count(xn)-1,1);  status  = 2;  silence(xn) = 0;  count(xn)   = count(xn) + 1;      elseif dst1(n) > T1 % 可能处于语音段%      zcr(n) < zcr2  status = 1;  count(xn)  = count(xn) + 1;      else  % 静音状态  status  = 0;  count(xn)   = 0;  x1(xn)=0;  x2(xn)=0;      end   case 2,  % 2 = 语音段      if dst1(n) > T1     % 保持在语音段  count(xn) = count(xn) + 1;  silence(xn) = 0;      else  % 语音将结束  silence(xn) = silence(xn)+1;  if silence(xn) < maxsilence    % 静音还不够长，尚未结束     count(xn)  = count(xn) + 1;  elseif count(xn) < minlen      % 语音长度太短，认为是噪声     status  = 0;     silence(xn) = 0;     count(xn)   = 0;  else      % 语音结束     status  = 3;     x2(xn)=x1(xn)+count(xn);  end      end   case 3,  % 语音结束，为下一个语音准备 status  = 0;    xn=xn+1;  count(xn)   = 0; silence(xn)=0; x1(xn)=0; x2(xn)=0;   endend   el=length(x1);if x1(el)==0, el=el-1; end% 获得x1的实际长度if el==0, return; endif x2(el)==0% 如果x2最后一个值为0，对它设置为fn    fprintf('Error: Not find endding point!\n');    x2(el)=fn;endSF=zeros(1,fn);    % 按x1和x2，对SF和NF赋值NF=ones(1,fn);for i=1 : el    SF(x1(i):x2(i))=1;    NF(x1(i):x2(i))=0;endspeechIndex=find(SF==1);  % 计算voicesegvoiceseg=findSegment(speechIndex);vsl=length(voiceseg);

在这里插入图片描述

4. MATLAB仿真三：语音信号谱距离法端点检测三

%% pr6_5_3 clear all; clc; close all;run Set_I   % 基本设置run PART_I  % 读入数据，分帧等准备ccc=mfcc_m(signal,fs,16,wlen,inc);      % 计算MFCCfn1=size(ccc,1);   % 取帧数1frameTime1=frameTime(3:fn-2);    % MFCC对应的时间刻度Ccep=ccc(:,1:16);  % 取得MFCC系数C0=mean(Ccep(1:5,:),1);   % 计算噪声平均MFCC倒谱系数的估计值for i=6 : fn1    Cn=Ccep(i,:);  % 取一帧MFCC倒谱系数    Dstu=0;    for k=1 : 16   % 从第6帧开始计算每帧MFCC倒谱系数与 Dstu=Dstu+(Cn(k)-C0(k))^2;      % 噪声MFCC倒谱系数的距离    end    Dcep(i)=sqrt(Dstu);endDcep(1:5)=Dcep(6);Dstm=multimidfilter(Dcep,2);     % 平滑处理dth=max(Dstm(1:NIS-2));   % 阈值计算T1=1.2*dth;T2=1.5*dth;[voiceseg,vsl,SF,NF]=vad_param1D(Dstm,T1,T2);% MFCC倒谱距离双门限的端点检测% 作图subplot 311; plot(time,x,'k');title('纯语音波形');ylabel('幅值'); axis([0 max(time) -1 1]);subplot 312; plot(time,signal,'k');title('加噪语音波形(信噪比10dB)');ylabel('幅值'); axis([0 max(time) -1 1]);subplot 313; plot(frameTime1,Dstm,'k');title('短时MFCC倒谱距离值'); axis([0 max(time) 0 1.2*max(Dstm)]);xlabel('时间/s'); ylabel('幅值'); line([0,frameTime(fn)], [T1 T1], 'color','k','LineStyle','--');line([0,frameTime(fn)], [T2 T2], 'color','k','LineStyle','-');for k=1 : vsl      % 标出语音端点    nx1=voiceseg(k).begin; nx2=voiceseg(k).end;    fprintf('%4d   %4d   %4d\n',k,nx1,nx2);    subplot 311;     line([frameTime(nx1) frameTime(nx1)],[-1 1],'color','k','LineStyle','-');    line([frameTime(nx2) frameTime(nx2)],[-1 1],'color','k','LineStyle','--');    subplot 313;     line([frameTime(nx1) frameTime(nx1)],[0 1.2*max(Dstm)],'color','k','LineStyle','-');    line([frameTime(nx2) frameTime(nx2)],[0 1.2*max(Dstm)],'color','k','LineStyle','--');end

在这里插入图片描述

5. MATLAB仿真四：语音信号谱熵法端点检测一

%% pr6_6_1 clear all; clc; close allrun Set_I   % 基本设置run PART_I  % 读入数据，分帧等准备for i=1:fn    Sp = abs(fft(y(:,i)));% FFT变换取幅值    Sp = Sp(1:wlen/2+1);     % 只取正频率部分    Ep=Sp.*Sp;     % 求出能量    prob = Ep/(sum(Ep)); % 计算每条谱线的概率密度    H(i) = -sum(prob.*log(prob+eps));  % 计算谱熵endEnm=multimidfilter(H,10); % 平滑处理Me=min(Enm);% 计算阈值 eth=mean(Enm(1:NIS));     Det=eth-Me;T1=0.98*Det+Me;T2=0.93*Det+Me;[voiceseg,vsl,SF,NF]=vad_param1D_revr(Enm,T1,T2);% 用双门限法反向检测端点% 作图subplot 311; plot(time,x,'k'); hold ontitle('语音波形');ylabel('幅值'); axis([0 max(time) -1 1]);subplot 312; plot(time,signal,'k'); hold ontitle(['带噪语音波形 SNR=' num2str(SNR) 'dB'] );ylabel('幅值'); axis([0 max(time) -1 1]);top=Det*1.1+Me;botton=Me-0.1*Det;subplot 313; plot(frameTime,Enm,'k');  axis([0 max(time) botton top]); line([0 fn],[T1 T1],'color','k','LineStyle','--');line([0 fn],[T2 T2],'color','k','LineStyle','-');title('短时谱熵'); xlabel('时间/s'); ylabel('谱熵值');for k=1 : vsl      % 标出语音端点    nx1=voiceseg(k).begin; nx2=voiceseg(k).end;    fprintf('%4d   %4d   %4d\n',k,nx1,nx2);    subplot 311    line([frameTime(nx1) frameTime(nx1)],[-1 1],'color','k','LineStyle','-');    line([frameTime(nx2) frameTime(nx2)],[-1 1],'color','k','LineStyle','--');    subplot 313    line([frameTime(nx1) frameTime(nx1)],[botton top],'color','k','LineStyle','-');    line([frameTime(nx2) frameTime(nx2)],[botton top],'color','k','LineStyle','--');end

function [voiceseg,vsl,SF,NF]=vad_param1D_revr(dst1,T1,T2)fn=length(dst1);  % 取得帧数maxsilence = 8;   % 初始化  minlen  = 5;    status  = 0;count   = 0;silence = 0;%开始端点检测xn=1;for n=2:fn   switch status   case {0,1}      % 0 = 静音, 1 = 可能开始      if dst1(n) < T2     % 确信进入语音段  x1(xn) = max(n-count(xn)-1,1);  status  = 2;  silence(xn) = 0;  count(xn)   = count(xn) + 1;      elseif dst1(n) < T1 % 可能处于语音段  status = 1;  count(xn)  = count(xn) + 1;      else  % 静音状态  status  = 0;  count(xn)   = 0;  x1(xn)=0;  x2(xn)=0;      end   case 2,  % 2 = 语音段      if dst1(n) < T1     % 保持在语音段  count(xn) = count(xn) + 1;  silence(xn) = 0;      else  % 语音将结束  silence(xn) = silence(xn)+1;  if silence(xn) < maxsilence    % 静音还不够长，尚未结束     count(xn)  = count(xn) + 1;  elseif count(xn) < minlen      % 语音长度太短，认为是噪声     status  = 0;     silence(xn) = 0;     count(xn)   = 0;  else      % 语音结束     status  = 3;     x2(xn)=x1(xn)+count(xn);  end      end   case 3,  % 语音结束，为下一个语音准备 status  = 0;    xn=xn+1;  count(xn)   = 0; silence(xn)=0; x1(xn)=0; x2(xn)=0;   endend   el=length(x1);if x1(el)==0, el=el-1; end% 获得x1的实际长度if el==0, return; endif x2(el)==0% 如果x2最后一个值为0，对它设置为fn    fprintf('Error: Not find endding point!\n');    x2(el)=fn;endSF=zeros(1,fn);    % 按x1和x2，对SF和NF赋值NF=ones(1,fn);for i=1 : el    SF(x1(i):x2(i))=1;    NF(x1(i):x2(i))=0;endspeechIndex=find(SF==1);  % 计算voicesegvoiceseg=findSegment(speechIndex);vsl=length(voiceseg);

在这里插入图片描述

6. MATLAB仿真五：语音信号谱熵法端点检测二

%% pr6_6_2 clear all; clc; close allrun Set_I   % 基本设置run PART_I  % 读入数据，分帧等准备df=fs/wlen; % 求出FFT后频率分辨率fx1=fix(250/df)+1; fx2=fix(3500/df)+1;  % 找出250Hz和3500Hz的位置km=floor(wlen/8);  % 计算出子带个数K=0.5;      % 常数Kfor i=1:fn    A=abs(fft(y(:,i)));   % 取来一帧数据FFT后取幅值    E=zeros(wlen/2+1,1);  E(fx1+1:fx2-1)=A(fx1+1:fx2-1);      % 只取250～3500Hz之间的分量    E=E.*E; % 计算能量    P1=E/sum(E);   % 幅值归一化    index=find(P1>=0.9);  % 寻找是否有分量的概率大于0.9    if ~isempty(index), E(index)=0; end % 若有,该分量置0    for m=1:km     % 计算子带能量 Eb(m)=sum(E(4*m-3:4*m));    end    prob=(Eb+K)/sum(Eb+K);% 计算子带概率    Hb(i) = -sum(prob.*log(prob+eps));  % 计算子带谱熵end   Enm=multimidfilter(Hb,10);% 平滑处理Me=min(Enm);% 设置阈值eth=mean(Enm(1:NIS));Det=eth-Me;T1=0.99*Det+Me;T2=0.96*Det+Me;[voiceseg,vsl,SF,NF]=vad_param1D_revr(Enm,T1,T2);% 用双门限法反向端点检测  % 作图top=Det*1.1+Me;botton=Me-0.1*Det;subplot 311; plot(time,x,'k'); hold ontitle('语音波形');ylabel('幅值'); axis([0 max(time) -1 1]);subplot 312; plot(time,signal,'k'); hold ontitle(['带噪语音波形 SNR=' num2str(SNR) 'dB'] );ylabel('幅值'); axis([0 max(time) -1 1]);subplot 313; plot(frameTime,Enm,'k');  axis([0 max(time) botton top]); line([0 fn],[T1 T1],'color','k','LineStyle','--');line([0 fn],[T2 T2],'color','k','LineStyle','-');title('短时改进子带谱熵'); xlabel('时间/s'); ylabel('谱熵值');for k=1 : vsl      % 标出语音端点    nx1=voiceseg(k).begin; nx2=voiceseg(k).end;    fprintf('%4d   %4d   %4d\n',k,nx1,nx2);    subplot 311    line([frameTime(nx1) frameTime(nx1)],[-1 1],'color','k','LineStyle','-');    line([frameTime(nx2) frameTime(nx2)],[-1 1],'color','k','LineStyle','--');    subplot 313    line([frameTime(nx1) frameTime(nx1)],[botton top],'color','k','LineStyle','-');    line([frameTime(nx2) frameTime(nx2)],[botton top],'color','k','LineStyle','--');end

在这里插入图片描述

7. MATLAB仿真六：语音信号能零比法端点检测

% % pr6_7_1 clear all; clc; close all;run Set_I   % 基本设置run PART_I  % 读入数据，分帧等准备aparam=2; bparam=1;% 设置参数etemp=sum(y.^2);   % 计算能量etemp1=log10(1+etemp/aparam);    % 计算能量的对数值zcr=zc2(y,fn);     % 求过零点值Ecr=etemp1./(zcr+bparam); % 计算能零比Ecrm=multimidfilter(Ecr,2);      % 平滑处理dth=mean(Ecrm(1:(NIS)));  % 阈值计算T1=1.2*dth;T2=2*dth;[voiceseg,vsl,SF,NF]=vad_param1D(Ecrm,T1,T2);% 能零比法的双门限端点检测% 作图subplot 311; plot(time,x,'k');title('纯语音波形');ylabel('幅值'); axis([0 max(time) -1 1]);subplot 312; plot(time,signal,'k');title('加噪语音波形(信噪比10dB)');ylabel('幅值'); axis([0 max(time) -1 1]);subplot 313; plot(frameTime,Ecrm,'k');title('短时能零比值'); grid; ylim([0 1.2*max(Ecrm)]);xlabel('时间/s'); ylabel('能零比值'); line([0,frameTime(fn)], [T1 T1], 'color','k','LineStyle','--');line([0,frameTime(fn)], [T2 T2], 'color','k','LineStyle','-');for k=1 : vsl      % 标出语音端点    nx1=voiceseg(k).begin; nx2=voiceseg(k).end;    fprintf('%4d   %4d   %4d\n',k,nx1,nx2);    subplot 311;     line([frameTime(nx1) frameTime(nx1)],[-1 1],'color','k','LineStyle','-');    line([frameTime(nx2) frameTime(nx2)],[-1 1],'color','k','LineStyle','--');    subplot 313;     line([frameTime(nx1) frameTime(nx1)],[0 1.2*max(Ecrm)],'color','k','LineStyle','-');    line([frameTime(nx2) frameTime(nx2)],[0 1.2*max(Ecrm)],'color','k','LineStyle','--');end

在这里插入图片描述

8. MATLAB仿真七：语音信号能熵比法端点检测

%% pr6_7_2 clear all; clc; close allrun Set_I % 基本设置run PART_I% 读入数据，分帧等准备aparam=2; % 设置参数for i=1:fn    Sp = abs(fft(y(:,i)));     % FFT变换取幅值    Sp = Sp(1:wlen/2+1);   % 只取正频率部分    Esum(i) = log10(1+sum(Sp.*Sp)/aparam);   % 计算对数能量值    prob = Sp/(sum(Sp));      % 计算概率    H(i) = -sum(prob.*log(prob+eps)); % 求谱熵值    Ef(i) = sqrt(1 + abs(Esum(i)/H(i)));     % 计算能熵比end   Enm=multimidfilter(Ef,10);     % 平滑滤波 Me=max(Enm);     % Enm最大值eth=mean(Enm(1:NIS));   % 初始均值ethDet=Me-eth;      % 求出值后设置阈值T1=0.05*Det+eth;T2=0.1*Det+eth;[voiceseg,vsl,SF,NF]=vad_param1D(Enm,T1,T2); % 用能熵比法的双门限端点检测% 作图top=Det*1.1+eth;botton=eth-0.1*Det;subplot 311; plot(time,x,'k'); hold ontitle('语音波形');ylabel('幅值'); axis([0 max(time) -1 1]);subplot 312; plot(time,signal,'k'); hold ontitle(['带噪语音波形 SNR=' num2str(SNR) 'dB'] );ylabel('幅值'); axis([0 max(time) -1 1]);subplot 313; plot(frameTime,Enm,'k');  axis([0 max(time) botton top]); line([0,frameTime(fn)], [T1 T1], 'color','k','LineStyle','--');line([0,frameTime(fn)], [T2 T2], 'color','k','LineStyle','-');title('短时能熵比'); xlabel('时间/s');for k=1 : vsl      % 标出语音端点    nx1=voiceseg(k).begin; nx2=voiceseg(k).end;    fprintf('%4d   %4d   %4d\n',k,nx1,nx2);    subplot 311    line([frameTime(nx1) frameTime(nx1)],[-1 1],'color','k','linestyle','-');    line([frameTime(nx2) frameTime(nx2)],[-1 1],'color','k','linestyle','--');    subplot 313    line([frameTime(nx1) frameTime(nx1)],[botton top],'color','k','linestyle','-');    line([frameTime(nx2) frameTime(nx2)],[botton top],'color','k','linestyle','--');end

在这里插入图片描述

9. MATLAB仿真八：语音信号小波变换法端点检测

%% pr6_8_1 clear all; clc; close all;run Set_I   % 基本设置run PART_I  % 读入数据，分帧等准备% 小波分解后参数start=[ 1  8  15  22  29  37  47  60  79  110  165];send=[ 7  14  21  28  36  46  59  78  109  164  267];duration=[ 7  7  7  7  8  10  13  19  31  55  103];for i=1 : fn    u=y(:,i);      % 取一帧    [c,l]=wavedec(u,10,'db4');   % 用母小波db4进行10层分解    for k=1 : 10 E(11-k)=mean(abs(c(start(k+1):send(k+1))));% 计算每层的平均幅值    end    M1=max(E(1:5)); M2=max(E(6:10));    % 按式(6-8-2)求M1和M2    MD(i)=M1*M2;   % 按式(6-8-3)计算MDendMDm=multimidfilter(MD,10);% 平滑处理MDmth=mean(MDm(1:NIS));   % 计算阈值T1=2*MDmth;T2=3*MDmth;[voiceseg,vsl,SF,NF]=vad_param1D(MDm,T1,T2);% 用小波分解系数平均幅值积法的双门限端点检测% 作图subplot 311; plot(time,x,'k');title('纯语音波形');ylabel('幅值'); axis([0 max(time) -1 1]);subplot 312; plot(time,signal,'k');title(['加噪语音波形 信噪比=' num2str(SNR) 'dB']);ylabel('幅值'); axis([0 max(time) -1 1]);subplot 313; plot(frameTime,MDm,'k');title('小波分解短时系数平均幅值积'); grid; ylim([0 1.2*max(MDm)]); xlabel('时间/s'); ylabel('幅值');line([0,frameTime(fn)], [T1 T1], 'color','k','LineStyle','--');line([0,frameTime(fn)], [T2 T2], 'color','k','LineStyle','-');for k=1 : vsl    nx1=voiceseg(k).begin; nx2=voiceseg(k).end;    fprintf('%4d   %4d   %4d\n',k,nx1,nx2);    subplot 311;     line([frameTime(nx1) frameTime(nx1)],[-1 1],'color','k','LineStyle','-');    line([frameTime(nx2) frameTime(nx2)],[-1 1],'color','k','LineStyle','--');    subplot 313;     line([frameTime(nx1) frameTime(nx1)],[0 1.2*max(MDm)],'color','k','LineStyle','-');    line([frameTime(nx2) frameTime(nx2)],[0 1.2*max(MDm)],'color','k','LineStyle','--');end

在这里插入图片描述

10. MATLAB仿真九：语音信号EMD分解法端点检测

%% pr6_8_2 clear all; clc; close all;run Set_I   % 基本设置run PART_I  % 读入数据，分帧等准备imf=emd(signal);   % EMD分解M=size(imf,1);     % 取得分解后IMF的阶数u=zeros(1,N);h=waitbar(0,'Running...');% 设置运行程序进度条图,初始化set(h,'name','端点检测 - 0%');    % 设置本图的名称“端点检测”for k=3 : M % 丢弃前2阶IMF重构语音信号    u=u+imf(k,:);endz=enframe(u,wnd,inc)';    % 重构语音信号的分帧for k=1 : fn    v=z(:,k);      % 取来一帧    imf=emd(v);    % EMD分解    L=size(imf,1); % 取得分解后IMF的阶数L    Etg=zeros(1,wlen);    for i=1 : L    % 计算每阶IMF的平均Teager能量 Etg=Etg+steager(imf(i,:));    end    Tg(k,:)=Etg;    Tgf(k)=mean(Etg);     % 计算本帧的平均Teager能量    waitbar(k/fn,h)% 显示运行的百分比,用红条表示% 显示本图的名称"端点检测",并显示运行的百分比数,用数字表示    set(h,'name',['端点检测 - ' sprintf('%2.1f',k/fn*100) '%'])endclose(h)    % 关闭程序进度条Zcr=zc2(y,fn);     % 计算过零率Tgfm=multimidfilter(Tgf,10);     % 平滑处理Zcrm=multimidfilter(Zcr,10);     % 平滑处理Mtg=max(Tgfm);     % 计算阈值Tmth=mean(Tgfm(1:NIS));Zcrth=mean(Zcrm(1:NIS));T1=1.5*Tmth;T2=3*Tmth;T3=0.9*Zcrth;T4=0.8*Zcrth;% 双参数双门限的端点检测[voiceseg,vsl,SF,NF]=vad_param2D_revr(Tgfm,Zcrm,T1,T2,T3,T4);% 作图pos = get(gcf,'Position');set(gcf,'Position',[pos(1), pos(2)-200,pos(3),pos(4)+150]) subplot 511; plot(time,x,'k');title('纯语音波形');ylabel('幅值'); axis([0 max(time) -1 1]);subplot 512; plot(time,signal,'k');title(['加噪语音波形 信噪比=' num2str(SNR) 'dB']);ylabel('幅值'); axis([0 max(time) -1 1]);subplot 513; plot(time,u,'k'); ylim([-1 1]);title('EMD后重构的语音波形');ylabel('幅值'); axis([0 max(time) -1 1]);subplot 514; plot(frameTime,Tgfm,'k'); title('短时EMD分解后Teager能量平均值'); ylim([0 1.2*Mtg]);ylabel('T能量平均值'); line([0 max(time)], [T1 T1],'color','k','lineStyle','--');line([0 max(time)], [T2 T2], 'color','k','lineStyle','-');subplot 515; plot(frameTime,Zcrm,'k');title('短时过零率值'); ylim([0 1.2*max(Zcrm)]);xlabel('时间/s'); ylabel('过零率值'); line([0 max(time)], [T3 T3], 'color','k','lineStyle','--');line([0 max(time)], [T4 T4], 'color','k','lineStyle','-');for k=1 : vsl    nx1=voiceseg(k).begin; nx2=voiceseg(k).end;    fprintf('%4d   %4d   %4d\n',k,nx1,nx2);    figure(1); subplot 511;     line([frameTime(nx1) frameTime(nx1)],[-1 1],'color','k','lineStyle','-');    line([frameTime(nx2) frameTime(nx2)],[-1 1],'color','k','lineStyle','--');    subplot 514;     line([frameTime(nx1) frameTime(nx1)],[0 1.2*Mtg],'color','k','lineStyle','-');    line([frameTime(nx2) frameTime(nx2)],[0 1.2*Mtg],'color','k','lineStyle','--');    subplot 515;    line([frameTime(nx1) frameTime(nx1)],[0 1.2*max(Zcrm)],'color','k','lineStyle','-');    line([frameTime(nx2) frameTime(nx2)],[0 1.2*max(Zcrm)],'color','k','lineStyle','--');end

function [voiceseg,vsl,SF,NF]=vad_param2D_revr(dst1,dst2,T1,T2,T3,T4)fn=length(dst1);  % 取帧数maxsilence = 8;   % 初始化  minlen  = 5;    status  = 0;count   = 0;silence = 0;%开始端点检测xn=1;for n=1:fn   switch status   case {0,1}      % 0 = 静音, 1 = 可能开始      if dst1(n) > T2  |  ...    % 确信进入语音段  ( nargin==6 & dst2(n) < T4 )   x1(xn) = max(n-count(xn)-1,1);  status  = 2;  silence(xn) = 0;  count(xn)   = count(xn) + 1;      elseif dst1(n) > T1 | ... % 可能处于语音段      dst2(n) < T3  status = 1;  count(xn)  = count(xn) + 1;      else  % 静音状态  status  = 0;  count(xn)   = 0;  x1(xn)=0;  x2(xn)=0;      end   case 2,  % 2 = 语音段      if dst1(n) > T1 | ...     % 保持在语音段  dst2(n) <  T3   count(xn) = count(xn) + 1;  silence(xn) = 0;      else  % 语音将结束  silence(xn) = silence(xn)+1;  if silence(xn) < maxsilence    % 静音还不够长，尚未结束     count(xn)  = count(xn) + 1;  elseif count(xn) < minlen      % 语音长度太短，认为是噪声     status  = 0;     silence(xn) = 0;     count(xn)   = 0;  else      % 语音结束     status  = 3;     x2(xn)=x1(xn)+count(xn);  end      end   case 3,  % 语音结束，为下一个语音准备 status  = 0;    xn=xn+1;  count(xn)   = 0; silence(xn)=0; x1(xn)=0; x2(xn)=0;   endend   el=length(x1);if x1(el)==0, el=el-1; end% 获得x1的实际长度if x2(el)==0% 如果x2最后一个值为0，对它设置为fn    fprintf('Error: Not find endding point!\n');    x2(el)=fn;endSF=zeros(1,fn);    % 按x1和x2，对SF和NF赋值NF=ones(1,fn);for i=1 : el    SF(x1(i):x2(i))=1;    NF(x1(i):x2(i))=0;endspeechIndex=find(SF==1);  % 计算voicesegvoiceseg=findSegment(speechIndex);vsl=length(voiceseg);

在这里插入图片描述

11. MATLAB仿真十：低信噪比时的端点检测一

%% pr6_9_1 clear all; clc; close all;filedir=[]; % 设置文件路径filename='bluesky1.wav';  % 设置文件名称fle=[filedir filename]    % 构成文件路径和名称% [xx,fs]=wavread(fle);     % 读入数据[xx,fs]=audioread(fle);     % 读入数据x=xx/max(abs(xx)); % 幅值归一化N=length(x);time=(0:N-1)/fs;   % 设置时间IS=0.3;     % 设置前导无话段长度wlen=200;   % 设置帧长为25msinc=80;     % 求帧移SNR=20;     % 设置信噪比wnd=hamming(wlen); % 设置窗函数overlap=wlen-inc;  % 求重叠区长度NIS=fix((IS*fs-wlen)/inc +1);    % 求前导无话段帧数noisefile='destroyerops.wav';    % 指定噪声的文件名[signal,noise] = add_noisefile(x,noisefile,SNR,fs);% 叠加噪声y=enframe(signal,wnd,inc)';      % 分帧fn=size(y,2);      % 求帧数frameTime=frame2time(fn, wlen, inc, fs);Y=fft(y);   % FFTY=abs(Y(1:fix(wlen/2)+1,:));     % 计算正频率幅值N=mean(Y(:,1:NIS),2);     % 计算前导无话段噪声区平均频谱NoiseCounter=0;NoiseLength=9;for i=1:fn,     if i<=NIS      % 在前导无话段中设置为NF=1,SF=0 SpeechFlag=0; NoiseCounter=100; SF(i)=0; NF(i)=1; TNoise(:,i)=N;    else    % 检测每帧计算对数频谱距离 [NoiseFlag, SpeechFlag, NoiseCounter, Dist]=...     vad(Y(:,i),N,NoiseCounter,2.5,8);  SF(i)=SpeechFlag; NF(i)=NoiseFlag; D(i)=Dist; if SpeechFlag==0  % 如果是噪声段对噪声谱进行更新     N=(NoiseLength*N+Y(:,i))/(NoiseLength+1); %      end TNoise(:,i)=N;    end    SN(i)=sum(TNoise(:,i));      % 计算噪声谱幅值之和end% 作图pos = get(gcf,'Position');set(gcf,'Position',[pos(1), pos(2)-100,pos(3),pos(4)+100]) subplot 411; plot(time,x,'k',frameTime,SF,'k--'); title('纯语音波形');ylabel('幅值'); ylim([-1.2 1.2]);subplot 412; plot(time,signal,'k');title(['带噪语音 SNR=' num2str(SNR) '(dB)'])ylabel('幅值'); ylim([-1.5 1.5]);subplot 413; plot(frameTime,D,'k'); ylabel('幅值'); grid;title('对数频谱距离'); ylim([0 25]);subplot 414; plot(frameTime,SN,'k','linewidth',2); xlabel('时间/s'); ylabel('幅值'); grid;title('噪声幅值和');  ylim([7 10]);

function [y,noise] = add_noisefile(s,filepath_name,SNR,fs)s=s(:); % 把信号转换成列数据s=s-mean(s);   % 消除直流分量% [wavin,fs1,nbits]=wavread(filepath_name);   %读入噪声文件的数据[wavin,fs1]=audioread(filepath_name);   %读入噪声文件的数据wavin=wavin(:);% 把噪声数据转换成列数据if fs1~=fs     % 纯语音信号的采样频率与噪声的采样频率不相等    wavin1=resample(wavin,fs,fs1);  % 对噪声重采样，使噪声采样频率与纯语音信号的采样频率相同else    wavin1=wavin;endwavin1=wavin1-mean(wavin1);  % 消除直流分量ns=length(s);  % 求出s的长度noise=wavin1(1:ns);   % 把噪声长度截断为与s等长noise=noise-mean(noise);     % 噪声去除直流分量signal_power = 1/ns*sum(s.*s);      % 求出信号的能量noise_power=1/ns*sum(noise.*noise); % 求出噪声的能量noise_variance = signal_power / ( 10^(SNR/10) );   % 求出噪声设定的方差值noise=sqrt(noise_variance/noise_power)*noise;      % 调整噪声幅值y=s+noise;     % 构成带噪语音

在这里插入图片描述

12. MATLAB仿真十一：低信噪比时的端点检测二

%% pr6_9_2   clear all; clc; close all;run Set_I   % 基本设置SNR=0;      % 重新设置信噪比SNRrun PART_I  % 读入数据,分帧等准备snr1=SNR_singlech(x,signal)      % 计算初始信噪比值a=3; b=0.01;output=simplesubspec(signal,wlen,inc,NIS,a,b);  % 计算谱减snr2=SNR_singlech(x,output)      % 计算谱减后信噪比值y=enframe(output,wlen,inc)';     % 谱减后输出序列分帧nl2=wlen/2+1;Y=fft(y);   % FFT转成频域Y_abs=abs(Y(1:nl2,:));    % 取正频率域幅值M=floor(nl2/4);    % 计算子带数for k=1 : fn    for i=1 : M    % 每个子带由4条谱线相加 j=(i-1)*4+1; SY(i,k)=Y_abs(j,k)+Y_abs(j+1,k)+Y_abs(j+2,k)+Y_abs(j+3,k);    end    Dvar(k)=var(SY(:,k)); % 计算每帧子带分离的频带方差endDvarm=multimidfilter(Dvar,10);   % 平滑处理dth=max(Dvarm(1:(NIS)));  % 阈值计算T1=1.5*dth;T2=3*dth;[voiceseg,vsl,SF,NF]=vad_param1D(Dvarm,T1,T2);% 频域方差双门限的端点检测% 作图pos = get(gcf,'Position');set(gcf,'Position',[pos(1), pos(2)-150,pos(3),pos(4)+100]) subplot 411; plot(time,x,'k'); axis tight;title('纯语音波形'); ylabel('幅值'); xlabel('(a)');subplot 412; plot(time,signal,'k'); axis tight; xlabel('(b)');title(['带噪语音 信噪比=' num2str(SNR) 'dB']); ylabel('幅值')subplot 413; plot(time,output,'k'); xlabel('(c)');title(['谱减后语音波形 SNR=' num2str(round(snr2*100)/100) 'dB'] ); ylabel('幅值')subplot 414; plot(frameTime,Dvarm,'k');title('谱减短时均匀子带频带方差值'); xlabel(['时间/s' 10 '(d)']);  ylabel('方差值');line([0,frameTime(fn)], [T1 T1], 'color','k','LineStyle','--');line([0,frameTime(fn)], [T2 T2], 'color','k','LineStyle','-');ylim([0 max(Dvar)]);for k=1 : vsl    nx1=voiceseg(k).begin; nx2=voiceseg(k).end;    fprintf('%4d   %4d   %4d\n',k,nx1,nx2);    subplot 411;     line([frameTime(nx1) frameTime(nx1)],[-1 1],'color','k','LineStyle','-');    line([frameTime(nx2) frameTime(nx2)],[-1 1],'color','k','LineStyle','--');    subplot 414;     line([frameTime(nx1) frameTime(nx1)],[0 max(Dvar)],'color','k','LineStyle','-');    line([frameTime(nx2) frameTime(nx2)],[0 max(Dvar)],'color','k','LineStyle','--');end

function y=multimidfilter(x,m)a=x;for k=1 : m    b=medfilt1(a, 5);     a=b;endy=b;

在这里插入图片描述

13. MATLAB仿真十二：低信噪比时的端点检测三

%% pr6_9_3 clear all; clc; close allrun Set_I% 基本设置SNR=0;   % 重新设置信噪比SNRrun PART_I      % 读入数据,分帧等准备snr1=SNR_singlech(x,signal)   % 计算加噪后的信噪比alpha=2.8; beta=0.001; c=1;   % 设置参数alpha,beta和c% 调用多窗谱减函数Mtmpsd_ss,实现减噪处理output=Mtmpsd_ssb(signal,wlen,inc,NIS,alpha,beta,c);  snr2=SNR_singlech(x,output)   % 计算减噪后的信噪比 y=enframe(output,wlen,inc)';  % 对减噪后的信号分帧aparam=2;% 设置参数for i=1:fn      % 计算各帧能熵比    Sp = abs(fft(y(:,i)));    % FFT变换取幅值    Sp = Sp(1:wlen/2+1);  % 只取正频率部分    Esum(i) = log10(1+sum(Sp.*Sp)/aparam);  % 计算对数能量值    prob = Sp/(sum(Sp));     % 计算概率    H(i) = -sum(prob.*log(prob+eps));% 求谱熵值    Ef(i) = sqrt(1 + abs(Esum(i)/H(i)));    % 计算能熵比end   Enm=multimidfilter(Ef,10);    % 平滑滤波 Me=max(Enm);    % 取Enm的最大值eth=mean(Enm(1:NIS));  % 求均值ethDet=Me-eth;     % 设置阈值T1=0.05*Det+eth;T2=0.1*Det+eth;[voiceseg,vsl,SF,NF]=vad_param1D(Enm,T1,T2);% 用双门限法端点检测% 作图top=Det*1.1+eth;botton=eth-0.1*Det;pos = get(gcf,'Position');set(gcf,'Position',[pos(1), pos(2)-200,pos(3),pos(4)+100])subplot 411; plot(time,x,'k');title('语音波形');ylabel('幅值'); axis([0 max(time) -1 1]);subplot 412; plot(time,signal,'k');title(['带噪语音波形 SNR=' num2str(SNR) 'dB'] );ylabel('幅值'); axis([0 max(time) -1 1]);subplot 413; plot(time,output,'k');title(['谱减后语音波形 SNR=' num2str(round(snr2*100)/100) 'dB'] );ylabel('幅值'); axis([0 max(time) -1 1]);subplot 414; plot(frameTime,Enm,'k');  axis([0 max(time) botton top]); title('短时能熵比'); xlabel('时间/s'); line([0 fn],[T1 T1],'color','k','linestyle','--');line([0 fn],[T2 T2],'color','k','linestyle','-');for k=1 : vsl    nx1=voiceseg(k).begin; nx2=voiceseg(k).end;    fprintf('%4d   %4d   %4d\n',k,nx1,nx2);    subplot 411    line([frameTime(nx1) frameTime(nx1)],[-1 1],'color','k','linestyle','-');    line([frameTime(nx2) frameTime(nx2)],[-1 1],'color','k','linestyle','--');    subplot 414    line([frameTime(nx1) frameTime(nx1)],[botton top],'color','k','linestyle','-');    line([frameTime(nx2) frameTime(nx2)],[botton top],'color','k','linestyle','--');end

在这里插入图片描述

小结

语音信号的端点检测即检测出有效语音段，排除非有效语音的干扰，从而实现更好的语音识别效果。本章从最基础的双门限检测，到双门限检测法的改进方法，相关法、方差法、谱距离法、谱熵法、能零比方法、能熵比方法、小波变换方法、EMD分解方法，最后对低信噪比时的语音信号也进行了一定的讨论与研究。

语音信号的端点检测属于语音信号进入下一步语音识别的基础，对本章内容感兴趣或者想充分学习了解的，建议去研习书中第六章节的内容。后期会对其中一些知识点在自己理解的基础上进行讨论补充，欢迎大家一起学习交流。

关于宋老师：宋知用——默默传授MATLAB与信号处理知识的老人家

本系列文章列表如下：
《MATLAB语音信号分析与合成（第二版）》：第2章语音信号的时域、频域特性和短时分析技术
《MATLAB语音信号分析与合成（第二版）》：第3章语音信号在其他变换域中的分析技术和特性
《MATLAB语音信号分析与合成（第二版）》：第4章语音信号的线性预测分析
《MATLAB语音信号分析与合成（第二版）》：第5章带噪语音和预处理
《MATLAB语音信号分析与合成（第二版）》：第6章语音端点的检测（1）
《MATLAB语音信号分析与合成（第二版）》：第6章语音端点的检测（2）
《MATLAB语音信号分析与合成（第二版）》：第7章语音信号的减噪
《MATLAB语音信号分析与合成（第二版）》：第8章基音周期的估算方法
《MATLAB语音信号分析与合成（第二版）》：第9章共振峰的估算方法
《MATLAB语音信号分析与合成（第二版）》：第10章语音信号的合成算法

唱吧

《MATLAB语音信号分析与合成（第二版）》：第6章语音端点的检测（2）