《MATLAB语音信号分析与合成（第二版）》：第8章 基音周期的估算方法

• 前言
• 1. 数据与函数路径设置
• 2. MATLAB仿真一：基音周期提取的预处理
• 3. MATLAB仿真二：倒谱法的基音检测一
• 4. MATLAB仿真三：倒谱法的基音检测二
• 5. MATLAB仿真四：短时自相关的基音检测
• 6. MATLAB仿真五：短时平均幅度差的基音检测一
• 7. MATLAB仿真六：短时平均幅度差的基音检测二
• 8. MATLAB仿真七：线性预测的基音检测一
• 9. MATLAB仿真八：线性预测的基音检测二
• 10. MATLAB仿真九：主体-延伸法的基音检测
• 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仿真一：基音周期提取的预处理

%% 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.6f\n',b)fprintf('a=%5.6f   %5.6f   %5.6f   %5.6f   %5.6f   %5.6f   %5.6f\n',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. MATLAB仿真二：倒谱法的基音检测一

% % pr8_2_1  clc; close all; clear all;run Set_II;     % 参数设置run Part_II;    % 读入文件,分帧和端点检测lmin=fix(fs/500);      % 基音周期的最小值lmax=fix(fs/60);% 基音周期的最大值period=zeros(1,fn);    % 基音周期初始化for k=1:fn     if SF(k)==1 % 是否在有话帧中 y1=y(:,k).*hamming(wlen);    % 取来一帧数据加窗函数 xx=fft(y1);    % FFT a=2*log(abs(xx)+eps); % 取模值和对数 b=ifft(a);     % 求取倒谱  [R(k),Lc(k)]=max(b(lmin:lmax));     % 在lmin和lmax区间中寻找最大值 period(k)=Lc(k)+lmin-1;      % 给出基音周期    endend% 作图subplot 211, plot(time,x,'k');  title('语音信号')axis([0 max(time) -1 1]); ylabel('幅值');subplot 212; plot(frameTime,period,'k');xlim([0 max(time)]); title('基音周期'); grid; xlabel('时间/s'); ylabel('样点数');for k=1 : vosl  % 标出有话段    nx1=voiceseg(k).begin;    nx2=voiceseg(k).end;    nxl=voiceseg(k).duration;    fprintf('%4d   %4d   %4d   %4d\n',k,nx1,nx2,nxl);    subplot 211    line([frameTime(nx1) frameTime(nx1)],[-1 1],'color','k','linestyle','-');    line([frameTime(nx2) frameTime(nx2)],[-1 1],'color','k','linestyle','--');end

% Set_IIfiledir=[];      % 设置数据文件的路径filename='tone4.wav';   % 设置数据文件的名称fle=[filedir filename]  % 构成路径和文件名的字符串wlen=320; inc=80;% 分帧的帧长和帧移overlap=wlen-inc;% 帧之间的重叠部分T1=0.05;  % 设置基音端点检测的参数

% Part_II% [x,fs]=wavread(fle);   % 读入wav文件[x,fs]=audioread(fle);   % 读入wav文件x=x-mean(x);    % 消去直流分量x=x/max(abs(x));% 幅值归一化y  = enframe(x,wlen,inc)';    % 分帧fn  = size(y,2);% 取得帧数time = (0 : length(x)-1)/fs;  % 计算时间坐标frameTime = frame2time(fn, wlen, inc, fs);  % 计算各帧对应的时间坐标[voiceseg,vosl,SF,Ef]=pitch_vad1(y,fn,T1);   % 基音的端点检测

4. MATLAB仿真三：倒谱法的基音检测二

% % pr8_2_2 clc; close all; clear all;run Set_II;     % 参数设置run Part_II;    % 读入文件,分帧和端点检测lmin=fix(fs/500);      % 基音周期提取中最小值lmax=fix(fs/60);% 基音周期提取中最大值period=zeros(1,fn);    % 基音周期初始化for k=1:fn     if SF(k)==1 % 是否在有话帧中 y1=y(:,k).*hamming(wlen);    % 取来一帧数据加窗函数 xx=fft(y1);    % FFT a=2*log(abs(xx)+eps); % 取模值和对数 b=ifft(a);     % 求取倒谱  [R(k),Lc(k)]=max(b(lmin:lmax));     % 在lmin和lmax区间中寻找最大值 period(k)=Lc(k)+lmin-1;      % 给出基音周期    endendT0=zeros(1,fn); % 初始化T0和F0F0=zeros(1,fn);T0=pitfilterm1(period,voiceseg,vosl);% 对T0进行平滑处理求出基音周期T0Tindex=find(T0~=0);F0(Tindex)=fs./T0(Tindex);    % 求出基音频率F0% 作图subplot 311, plot(time,x,'k');  title('语音信号')axis([0 max(time) -1 1]); grid;  ylabel('幅值');subplot 312; line(frameTime,period,'color',[.6 .6 .6],'linewidth',3);xlim([0 max(time)]); title('基音周期'); hold on;ylim([0 150]); ylabel('样点数'); grid; for k=1 : vosl    nx1=voiceseg(k).begin;    nx2=voiceseg(k).end;    nxl=voiceseg(k).duration;    fprintf('%4d   %4d   %4d   %4d\n',k,nx1,nx2,nxl);    subplot 311    line([frameTime(nx1) frameTime(nx1)],[-1 1],'color','k','linestyle','-');    line([frameTime(nx2) frameTime(nx2)],[-1 1],'color','k','linestyle','--');endsubplot 312; plot(frameTime,T0,'k'); hold offlegend('平滑前','平滑后');subplot 313; plot(frameTime,F0,'k'); grid; ylim([0 450]);title('基音频率'); xlabel('时间/s'); ylabel('频率/Hz');

5. MATLAB仿真四：短时自相关的基音检测

% % pr8_3_1 clc; close all; clear all;run Set_II   % 参数设置run Part_II  % 读入文件，分帧和端点检测% 滤波器系数b=[0.012280   -0.039508   0.042177   0.000000   -0.042177   0.039508   -0.012280];a=[1.000000   -5.527146   12.854342   -16.110307   11.479789   -4.410179   0.713507];xx=filter(b,a,x);   % 带通数字滤波yy  = enframe(xx,wlen,inc)';      % 滤波后信号分帧lmin=fix(fs/500);   % 基音周期的最小值lmax=fix(fs/60);    % 基音周期的最大值period=zeros(1,fn); % 基音周期初始化period=ACF_corr(yy,fn,voiceseg,vosl,lmax,lmin); % 用自相关函数提取基音周期T0=pitfilterm1(period,voiceseg,vosl);    % 平滑处理% 作图subplot 211, plot(time,x,'k');  title('语音信号')axis([0 max(time) -1 1]); grid;  ylabel('幅值');subplot 212; plot(frameTime,T0,'k'); hold on;xlim([0 max(time)]); title('平滑后的基音周期'); grid; xlabel('时间/s'); ylabel('样点数');for k=1 : vosl    nx1=voiceseg(k).begin;    nx2=voiceseg(k).end;    nxl=voiceseg(k).duration;    fprintf('%4d   %4d   %4d   %4d\n',k,nx1,nx2,nxl);    subplot 211    line([frameTime(nx1) frameTime(nx1)],[-1 1],'color','k','linestyle','-');    line([frameTime(nx2) frameTime(nx2)],[-1 1],'color','k','linestyle','--');end

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));   % 在Pmin～Pmax范围内寻找最大值 period(k+ixb-1)=lmin+tloc-1;      % 给出对应最大值的延迟量    endend

6. MATLAB仿真五：短时平均幅度差的基音检测一

% % pr8_4_1 clc; close all; clear all;run Set_II;  % 参数设置run Part_II; % 读入文件,分帧和端点检测% 滤波器系数b=[0.012280   -0.039508   0.042177   0.000000   -0.042177   0.039508   -0.012280];a=[1.000000   -5.527146   12.854342   -16.110307   11.479789   -4.410179   0.713507];xx=filter(b,a,x);   % 带通数字滤波yy  = enframe(xx,wlen,inc)';      % 滤波后信号分帧lmin=floor(fs/500); % 基音周期的最小值lmax=floor(fs/60);  % 基音周期的最大值period=zeros(1,fn); % 基音周期初始化period=AMDF_mod(yy,fn,voiceseg,vosl,lmax,lmin); % 用AMDF_mod函数提取基音周期T0=pitfilterm1(period,voiceseg,vosl);    % 基音周期平滑处理% 作图subplot 211, plot(time,x,'k');  title('语音信号')axis([0 max(time) -1 1]); grid;  ylabel('幅值'); xlabel('时间/s');subplot 212; hold online(frameTime,period,'color',[.6 .6 .6],'linewidth',2); axis([0 max(time) 0 120]); title('基音周期'); grid; xlabel('时间/s'); ylabel('样点数');subplot 212; plot(frameTime,T0,'k'); hold offlegend('初估算值','平滑后值'); box on;

function period=AMDF_mod(y,fn,vseg,vsl,lmax,lmin)pn=size(y,2);if pn~=fn, y=y'; end % 把y转换为每列数据表示一帧语音信号period=zeros(1,fn);  % 初始化wlen=size(y,1);      % 取得帧长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);     % 取来一帧数据 for m = 1:wlen      R0(m) = sum(abs(u(m:wlen)-u(1:wlen-m+1))); % 计算平均幅度差函数 end  [Rmax,Nmax]=max(R0);% 求取AMDF中最大值和对应位置 for i = 1 : wlen    % 进行线性变换     R(i)=Rmax*(wlen-i)/(wlen-Nmax)-R0(i); end [Rmax,T]=max(R(lmin:lmax));% 求出最大值 T0=T+lmin-1;  period(k+ixb-1)=T0; % 给出了该帧的基音周期    endend

7. MATLAB仿真六：短时平均幅度差的基音检测二

% % pr8_4_2 clc; close all; clear all;run Set_II;run Part_II;lmin=floor(fs/500);  % 基音周期的最小值lmax=floor(fs/60);   % 基音周期的最大值period=zeros(1,fn);  % 基音周期初始化T0=zeros(1,fn);      % 初始化period=ACFAMDF_corr(y,fn,voiceseg,vosl,lmax,lmin);  % 提取基音周期T0=pitfilterm1(period,voiceseg,vosl);     % 基音周期平滑处理% 作图subplot 211, plot(time,x,'k');  title('语音信号')axis([0 max(time) -1 1]); grid;  ylabel('幅值'); xlabel('时间/s');subplot 212; hold online(frameTime,period,'color',[.6 .6 .6],'linewidth',2); xlim([0 max(time)]); title('基音周期'); grid; xlabel('时间/s'); ylabel('样点数');plot(frameTime,T0,'k'); hold offlegend('初估算值','平滑后值'); box on

function period=ACFAMDF_corr(y,fn,vseg,vsl,lmax,lmin)pn=size(y,2);if pn~=fn, y=y'; end % 把y转换为每列数据表示一帧语音信号period=zeros(1,fn);  % 初始化wlen=size(y,1);      % 取得帧长Acm=zeros(1,lmax);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);  % 取延迟量为正值的部分 for m = 1:wlen      R(m) = sum(abs(u(m:wlen)-u(1:wlen-m+1))); % 计算平均幅度差函数(AMDF) end  R=R(1:length(ru));  % 取与ru等长  Rindex=find(R~=0); Acm(Rindex)=ru(Rindex)'./R(Rindex);% 计算ACF/AMDF [tmax,tloc]=max(Acm(lmin:lmax));  % 在Pmin～Pmax范围内寻找最大值 period(k+ixb-1)=lmin+tloc-1;      % 给出对应最大值的延迟量     endend

8. MATLAB仿真七：线性预测的基音检测一

% % pr8_5_1 clc; close all; clear all;run Set_II;     % 参数设置run Part_II;    % 读入文件,分帧和端点检测lmin=fix(fs/500);      % 基音周期的最小值lmax=fix(fs/60);% 基音周期的最大值period=zeros(1,fn);    % 基音周期初始化p=12;    % 设置线性预测阶数for k=1:fn     if SF(k)==1 % 是否在有话帧中 u=y(:,k).*hamming(wlen);     % 取来一帧数据加窗函数 ar = lpc(u,p); % 计算LPC系数 z = filter([0 -ar(2:end)],1,u);     % 一帧数据LPC逆滤波输出 E = u - z;     % 预测误差 xx=fft(E);     % FFT a=2*log(abs(xx)+eps); % 取模值和对数 b=ifft(a);     % 求取倒谱  [R(k),Lc(k)]=max(b(lmin:lmax));     % 在Pmin～Pmax区间寻找最大值 period(k)=Lc(k)+lmin-1;      % 给出基音周期    endendT1=pitfilterm1(period,voiceseg,vosl);% 基音周期平滑处理% 作图subplot 211, plot(time,x,'k');  title('语音信号')axis([0 max(time) -1 1]); grid;  ylabel('幅值'); xlabel('时间/s');subplot 212; hold online(frameTime,period,'color',[.6 .6 .6],'linewidth',2); axis([0 max(time) 0 150]); title('基音周期'); ylabel('样点数'); xlabel('时间/s'); grid; plot(frameTime,T1,'k'); hold offlegend('初估算值','平滑后值'); box on;

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;     % 赋值给yend

9. MATLAB仿真八：线性预测的基音检测二

% % pr8_5_2 clc; close all; clear all;run Set_II;   % 参数设置run Part_II;  % 读入文件,分帧和端点检测% 数字带通滤波器的设计Rp=1; Rs=50; fs2=fs/2;      % 通带波纹1dB,阻带衰减50dB Wp=[60 500]/fs2;     % 通带为60～500HzWs=[20 1000]/fs2;    % 阻带为20和1000Hz[n,Wn]=ellipord(Wp,Ws,Rp,Rs);      % 选用椭圆滤波器[b,a]=ellip(n,Rp,Rs,Wn);    % 求出滤波器系数x1=filter(b,a,x);    % 带通滤波x1=x1/max(abs(x1));  % 幅值归一化x2=resample(x1,1,4); % 按4:1降采样率lmin=fix(fs/500);    % 基音周期的最小值lmax=fix(fs/60);     % 基音周期的最大值period=zeros(1,fn);  % 基音周期初始化wind=hanning(wlen/4);% 窗函数y2=enframe(x2,wind,inc/4)'; % 再一次分帧p=4;   % LPC阶数为4for i=1 : vosl% 只对有话段数据处理    ixb=voiceseg(i).begin;  % 取一段有话段    ixe=voiceseg(i).end;    % 求取该有话段开始和结束位置及帧数    ixd=ixe-ixb+1;    for k=1 : ixd    % 对该段有话段数据处理 u=y2(:,k+ixb-1);    % 取来一帧数据 ar = lpc(u,p);      % 计算LPC系数 z = filter([0 -ar(2:end)],1,u);   % 一帧数据LPC逆滤波输出 E = u - z;   % 预测误差 ru1= xcorr(E, 'coeff');    % 计算归一化自相关函数 ru1 = ru1(wlen/4:end);     % 取延迟量为正值的部分 ru=resample(ru1,4,1);      % 按1:4升采样率 [tmax,tloc]=max(ru(lmin:lmax));   % 在Pmin～Pmax范围内寻找最大值 period(k+ixb-1)=lmin+tloc-1;      % 给出对应最大值的延迟量    endendT1=pitfilterm1(period,voiceseg,vosl);     % 基音周期平滑处理% 作图subplot 211, plot(time,x,'k');  title('语音信号')axis([0 max(time) -1 1]); grid;  ylabel('幅值'); xlabel('时间/s');subplot 212; hold online(frameTime,period,'color',[.6 .6 .6],'linewidth',2); xlim([0 max(time)]); title('基音周期'); grid; ylim([0 150]);  ylabel('样点数'); xlabel('时间/s'); plot(frameTime,T1,'k'); hold offlegend('初估算值','平滑后值'); box on

10. MATLAB仿真九：主体-延伸法的基音检测

%% pr8_6_1    clear all; clc; close allfiledir=[];   % 设置数据文件的路径filename='tone4.wav';% 设置数据文件的名称fle=[filedir filename]      % 构成路径和文件名的字符串% [xx,fs]=wavread(fle);% 读取文件[xx,fs]=audioread(fle);% 读取文件xx=xx-mean(xx);      % 消除直流分量xx=xx/max(abs(xx));  % 幅值归一化N=length(xx); % 信号长度time = (0 : N-1)/fs; % 设置时间刻度wlen=320;     % 帧长inc=80;% 帧移overlap=wlen-inc;    % 两帧重叠长度  lmin=floor(fs/500);  % 最高基音频率500Hz对应的基音周期lmax=floor(fs/60);   % 最高基音频率60Hz对应的基音周期 yy=enframe(xx,wlen,inc)';   % 第一次分帧fn=size(yy,2);% 取来总帧数frameTime = frame2time(fn, wlen, inc, fs);% 计算每一帧对应的时间Thr1=0.1;     % 设置端点检测阈值r2=0.5;% 设置元音主体检测的比例常数ThrC=[10 15]; % 设置相邻基音周期间的阈值% 用于基音检测的端点检测和元音主体检测[voiceseg,vosl,vseg,vsl,Thr2,Bth,SF,Ef]=pitch_vads(yy,fn,Thr1,r2,10,8);% 60～500Hz的带通滤波器系数b=[0.012280   -0.039508   0.042177   0.000000   -0.042177   0.039508   -0.012280];a=[1.000000   -5.527146   12.854342   -16.110307   11.479789   -4.410179   0.713507];x=filter(b,a,xx);    % 数字滤波x=x/max(abs(x));     % 幅值归一化y=enframe(x,wlen,inc)';     % 第二次分帧[Extseg,Dlong]=Extoam(voiceseg,vseg,vosl,vsl,Bth);   % 计算延伸区间和延伸长度T1=ACF_corrbpa(y,fn,vseg,vsl,lmax,lmin,ThrC(1));     % 对元音主体进行基音检测% 对语音主体的前后向过渡区延伸基音检测 T0=zeros(1,fn);      % 初始化F0=zeros(1,fn);for k=1 : vsl % 共有vsl个元音主体    ix1=vseg(k).begin;      % 第k个元音主体开始位置    ix2=vseg(k).end; % 第k个元音主体结束位置    in1=Extseg(k).begin;    % 第k个元音主体前向延伸开始位置    in2=Extseg(k).end;      % 第k个元音主体后向延伸结束位置    ixl1=Dlong(k,1); % 前向延伸长度    ixl2=Dlong(k,2); % 后向延伸长度    if ixl1>0 % 需要前向延伸基音检测 [Bt,voiceseg]=back_Ext_shtpm1(y,fn,voiceseg,Bth,ix1,... ixl1,T1,k,lmax,lmin,ThrC);    else      % 不需要前向延伸基音检测 Bt=[];    end    if ixl2>0 % 需要后向延伸基音检测 [Ft,voiceseg]=fore_Ext_shtpm1(y,fn,voiceseg,Bth,ix2,... ixl2,vsl,T1,k,lmax,lmin,ThrC);    else      % 不需要后向延伸基音检测 Ft=[];    end    T0(in1:in2)=[Bt T1(ix1:ix2) Ft];      % 第k个元音主体完成了前后向延伸检测 endtindex=find(T0>lmax);% 限制延伸后最大基音周期值不超越lmaxT0(tindex)=lmax;tindex=find(T0<lmin & T0~=0);      % 限制延伸后最小基音周期值不低于lminT0(tindex)=lmin;tindex=find(T0~=0);F0(tindex)=fs./T0(tindex);TT=pitfilterm1(T0,Extseg,vsl);      % T0平滑滤波FF=pitfilterm1(F0,Extseg,vsl);      % F0平滑滤波% STFT分析,绘制语谱图nfft=512;d=stftms(xx,wlen,nfft,inc);W2=1+nfft/2;n2=1:W2;freq=(n2-1)*fs/nfft;% 作图figure(1)pos = get(gcf,'Position');set(gcf,'Position',[pos(1), pos(2)-100,pos(3),pos(4)+100]);subplot 411; plot(time,xx,'k');title('原始信号波形'); axis([0 max(time) -1 1]); ylabel('幅值');for k=1 : vosl line([frameTime(voiceseg(k).begin) frameTime(voiceseg(k).begin)],... [-1 1],'color','k'); line([frameTime(voiceseg(k).end) frameTime(voiceseg(k).end)],... [-1 1],'color','k','linestyle','--');endsubplot 412; plot(frameTime,Ef,'k'); hold ontitle('能熵比'); axis([0 max(time) 0 max(Ef)]); ylabel('幅值');line([0 max(frameTime)],[Thr1 Thr1],'color','k','linestyle','--');plot(frameTime,Thr2,'k','linewidth',2);for k=1 : vsl line([frameTime(vseg(k).begin) frameTime(vseg(k).begin)],... [0 max(Ef)],'color','k','linestyle','-.'); line([frameTime(vseg(k).end) frameTime(vseg(k).end)],... [0 max(Ef)],'color','k','linestyle','-.');endtext(3.2,0.2,'Thr1');text(2.9,0.55,'Thr2');subplot 413; plot(frameTime,TT,'k'); title('基音周期'); grid; xlim([0 max(time)]); ylabel('样点数');subplot 414; plot(frameTime,FF,'k'); title('基音频率'); grid; xlim([0 max(time)]); xlabel('时间/s'); ylabel('频率/Hz')figure(2)pos = get(gcf,'Position');set(gcf,'Position',[pos(1), pos(2)-100,pos(3),(pos(4)-240)]);imagesc(frameTime,freq,abs(d(n2,:)));  axis xym = 64; LightYellow = [0.6 0.6 0.6];MidRed = [0 0 0]; Black = [0.5 0.7 1];Colors = [LightYellow; MidRed; Black];colormap(SpecColorMap(m,Colors));hold onplot(frameTime,FF,'w');ylim([0 1000]);title('语谱图上的基音频率曲线');xlabel('时间/s'); ylabel('频率/Hz')

function period=ACF_corrbpa(y,fn,vseg,vsl,lmax,lmin,ThrC)pn=size(y,2);if pn~=fn, y=y'; end % 把y转换为每列数据表示一帧语音信号wlen=size(y,1);      % 帧长period=zeros(1,fn);  % 初始化c1=ThrC(1);   % 取得阈值for i=1 : vsl % 只对有话段数据处理    ixb=vseg(i).begin;    ixe=vseg(i).end;    ixd=ixe-ixb+1;   % 求取一段有话段帧的帧数    Ptk=zeros(3,ixd);% 初始化    Vtk=zeros(3,ixd);    Ksl=zeros(1,3);    for k=1 : ixd     u=y(:,k+ixb-1);     % 取来一帧信号 ru=xcorr(u,'coeff');% 计算自相关函数 ru=ru(wlen:end);    % 取正延迟量部分 [Sv,Kv]=findmaxesm3(ru,lmax,lmin);% 获取三个极大值的数值和位置 lkv=length(Kv); Ptk(1:lkv,k)=Kv';   % 把每帧三个极大值位置存放在Ptk数组中    end    Kal=Ptk(1,:);    meanx=mean(Kal); % 计算Kal均值    thegma=std(Kal); % 计算Kal标准差    mt1=meanx+thegma;    mt2=meanx-thegma;    if thegma>5,  %判断基音候选组中有哪几个在第一次置信区域内 Ptemp=zeros(size(Ptk)); Ptemp(1,(Ptk(1,:)<mt1 & Ptk(1,:)>mt2))=1; Ptemp(2,(Ptk(2,:)<mt1 & Ptk(2,:)>mt2))=1; Ptemp(3,(Ptk(3,:)<mt1 & Ptk(3,:)>mt2))=1; % 如果第一组或(和)其他组都有值在第一次置信区内,取数值大的一个值赋于Pam Pam=zeros(1,ixd); for k=1 : ixd     if Ptemp(1,k)==1  Pam(k)=max(Ptk(:,k).*Ptemp(:,k));     end end mdex=find(Pam~=0);      % 在Pam非零的数值中 meanx=mean(Pam(mdex));  % 计算Pam均值 thegma=std(Pam(mdex));  % 计算Pam标准差 if thegma<0.5, thegma=0.5; end mt1=meanx+thegma; mt2=meanx-thegma;% 计算了第二次置信区 pindex=find(Pam<mt1 & Pam>mt2);% 寻找在置信区间的数据点 Pamtmp=zeros(1,ixd); Pamtmp(pindex)=Pam(pindex);    % 设置Pamtmp if length(pindex)~=ixd     bpseg=findSegment(pindex); % 计算置信区间内的数据分段信息     bpl=length(bpseg);  % 置信区间内的数据分成几段      bdb=bpseg(1).begin; % 置信区间内第一段的开始位置     if bdb~=1    % 如果置信区间内第一段开始位置不为1  Ptb=Pamtmp(bdb);% 置信区间内第一段开始位置的基音周期  Ptbp=Pamtmp(bdb+1);  pdb=ztcont11(Ptk,bdb,Ptb,Ptbp,c1);% 将调用ztcont11函数处理  Pam(1:bdb-1)=pdb;      % 把处理后的数据放入Pam中     end     if bpl>=2  for k=1 : bpl-1 % 如果有中间数据在置信区间外      pdb=bpseg(k).end;      pde=bpseg(k+1).begin;      Ptb=Pamtmp(pdb);      Pte=Pamtmp(pde);      pdm=ztcont21(Ptk,pdb,pde,Ptb,Pte,c1);  % 将调用ztcont21函数处理      Pam(pdb+1:pde-1)=pdm;     % 把处理后的数据放入Pam中  end     end     bde=bpseg(bpl).end;     Pte=Pamtmp(bde);     Pten=Pamtmp(bde-1);     if bde~=ixd  % 如果置信区间内最后一段结束位置不为ixd  pde=ztcont31(Ptk,bde,Pte,Pten,c1);% 将调用ztcont31函数处理  Pam(bde+1:ixd)=pde;    % 把处理后的数据放入Pam中     end end period(ixb:ixe)=Pam; else     period(ixb:ixe)=Kal;    endend

11. MATLAB仿真十：带噪语音的基音检测一

% % pr8_7_1   clc; close all; clear all;filedir=[]; % 设置数据文件的路径filename='tone4.wav';     % 设置数据文件的名称fle=[filedir filename]    % 构成路径和文件名的字符串% [xx,fs]=wavread(fle);     % 读取文件[xx,fs]=audioread(fle);     % 读取文件xx=xx-mean(xx);    % 消除直流分量xx=xx/max(abs(xx));% 幅值归一化SNR=5;      % 设置信噪比x=Gnoisegen(xx,SNR);      % 叠加噪声wlen=320;  inc=80; % 帧长和帧移overlap=wlen-inc;  % 两帧重叠长度 N=length(x);% 信号长度time=(0:N-1)/fs;   % 设置时间y=enframe(x,wlen,inc)';   % 第一次分帧fn=size(y,2);      % 取来帧长frameTime = frame2time(fn, wlen, inc, fs);  % 计算每一帧对应的时间T1=0.23;    % 设置T1[voiceseg,vsl,SF,Ef]=pitch_vad1(y,fn,T1); % 基音的端点检测% 60～500Hz的带通滤波器系数b=[0.012280   -0.039508   0.042177   0.000000   -0.042177   0.039508   -0.012280];a=[1.000000   -5.527146   12.854342   -16.110307   11.479789   -4.410179   0.713507];z=filter(b,a,x);   % 带通数字滤波yy  = enframe(z,wlen,inc)';      % 滤波后信号分帧lmin=floor(fs/500);% 基音周期的最小值lmax=floor(fs/60); % 基音周期的最大值period=zeros(1,fn);% 基音周期初始化F0=zeros(1,fn);    % 初始化 period=Wavelet_corrm1(yy,fn,voiceseg,vsl,lmax,lmin); % 小波-自相关函数基音检测tindex=find(period~=0);F0(tindex)=fs./period(tindex);   % 求出基音频率TT=pitfilterm1(period,voiceseg,vsl);    % 平滑处理FF=pitfilterm1(F0,voiceseg,vsl); % 平滑处理% 作图pos = get(gcf,'Position');set(gcf,'Position',[pos(1), pos(2)-150,pos(3),pos(4)+100]);subplot 511; plot(time,xx,'k');axis([0 max(time) -1 1]); title('原始语音'); ylabel('幅值');for k=1 : vsl    line([frameTime(voiceseg(k).begin) frameTime(voiceseg(k).begin)],...    [-1 1],'color','k','linestyle','-');    line([frameTime(voiceseg(k).end) frameTime(voiceseg(k).end)],...    [-1 1],'color','k','linestyle','--');endsubplot 512; plot(frameTime,Ef,'k')line([0 max(frameTime)],[T1 T1],'color','k','linestyle','-.');text(3.25, T1+0.05,'T1');axis([0 max(time) 0 1]); title('能熵比'); ylabel('幅值')for k=1 : vsl    line([frameTime(voiceseg(k).begin) frameTime(voiceseg(k).begin)],...    [-1 1],'color','k','linestyle','-');    line([frameTime(voiceseg(k).end) frameTime(voiceseg(k).end)],...    [-1 1],'color','k','linestyle','--');endsubplot 513; plot(time,x,'k');  title('带噪语音')axis([0 max(time) -1 1]); ylabel('幅值')subplot 514; plot(frameTime,TT,'k','linewidth',2);axis([0 max(time) 0 120]); title('基音周期'); grid; ylabel('样点数')subplot 515; plot(frameTime,FF,'k','linewidth',2);axis([0 max(time) 0 500]); title('基音频率'); grid; xlabel('时间/s');  ylabel('频率/Hz')

function [voiceseg,vosl,SF,Ef]=pitch_vad1(y,fn,T1,miniL)if nargin<4, 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取幅值    Sp = Sp(1:wlen/2+1);    % 只取正频率部分    Esum(i) = sum(Sp.*Sp);      % 计算能量值    prob = Sp/(sum(Sp));    % 计算概率    H(i) = -sum(prob.*log(prob+eps));  % 求谱熵值endhindex=find(H<0.1);H(hindex)=max(H);Ef=sqrt(1 + abs(Esum./H));      % 计算能熵比Ef=Ef/max(Ef);    % 归一化zindex=find(Ef>=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    endendvosl=length(voiceseg);   % 有话段的段数 

function period=Wavelet_corrm1(y,fn,vseg,vsl,lmax,lmin)pn=size(y,2);if pn~=fn, y=y'; end % 把y转换为每列数据表示一帧语音信号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);     % 取来一帧数据 [ca1,cd1] = dwt(u,'db4');  % 用dwt做小波变换  a1 = upcoef('a',ca1,'db4',1);     % 用低频系数重构信号 ru=xcorr(a1, 'coeff');     % 计算归一化自相关函数 aL=length(a1); ru=ru(aL:end);      % 取延迟量为正值的部分 [tmax,tloc]=max(ru(lmin:lmax));   % 在lmin至lmax范围内寻找最大值  period(k+ixb-1)=lmin+tloc-1;      % 给出对应最大值的延迟量    endend

12. MATLAB仿真十一：带噪语音的基音检测二

%% pr8_7_2  clear all; clc; close all;filedir=[]; % 设置数据文件的路径filename='tone4.wav';     % 设置数据文件的名称fle=[filedir filename]    % 构成路径和文件名的字符串% [x,fs]=wavread(fle);      % 读取文件[x,fs]=audioread(fle);      % 读取文件x=x-mean(x);% 消除直流分量x=x/max(abs(x));   % 幅值归一化SNR=0;      % 设置信噪比signal=Gnoisegen(x,SNR);  % 叠加噪声snr1=SNR_singlech(x,signal)      % 计算初始信噪比值N=length(x);% 信号长度time = (0 : N-1)/fs;      % 设置时间刻度wlen=320; inc=80;  % 帧长和帧移overlap=wlen-inc;  % 两帧重叠长度  IS=0.15;    % 设置前导无话段长度NIS=fix((IS*fs-wlen)/inc +1);    % 求前导无话段帧数a=3; b=0.001;      % 设置参数a和boutput=simplesubspec(signal,wlen,inc,NIS,a,b); % 谱减运算snr2=SNR_singlech(x,output)      % 计算谱减后信噪比值y  = enframe(output,wlen,inc)';  % 分帧fn  = size(y,2);   % 取得帧数time = (0 : length(x)-1)/fs;     % 计算时间坐标frameTime = frame2time(fn, wlen, inc, fs);% 计算每一帧对应的时间T1=0.12;     % 设置基音端点检测的参数[voiceseg,vosl,SF,Ef]=pitch_vad1(y,fn,T1);   % 基音的端点检测% 60～500Hz的带通滤波器系数b=[0.012280   -0.039508   0.042177   0.000000   -0.042177   0.039508   -0.012280];a=[1.000000   -5.527146   12.854342   -16.110307   11.479789   -4.410179   0.713507];z=filter(b,a,output);     % 带通数字滤波yy  = enframe(z,wlen,inc)';      % 滤波后信号分帧lmin=floor(fs/500);% 基音周期的最小值lmax=floor(fs/60); % 基音周期的最大值period=zeros(1,fn);% 基音周期初始化period=ACF_corr(yy,fn,voiceseg,vosl,lmax,lmin);  % 用自相关函数提取基音周期tindex=find(period~=0);F0=zeros(1,fn);    % 初始化 F0(tindex)=fs./period(tindex);   % 求出基音频率TT=pitfilterm1(period,voiceseg,vosl);   % 对基音周期进行平滑滤波FF=pitfilterm1(F0,voiceseg,vosl);% 对基音频率进行平滑滤波% 作图pos = get(gcf,'Position');set(gcf,'Position',[pos(1), pos(2)-150,pos(3),pos(4)+100]);subplot 611; plot(time,x,'k'); ylabel('幅值');title('原始信号'); axis([0 max(time) -1 1]);for k=1 : vosl line([frameTime(voiceseg(k).begin) frameTime(voiceseg(k).begin)],... [-1 1],'color','k'); line([frameTime(voiceseg(k).end) frameTime(voiceseg(k).end)],... [-1 1],'color','k','linestyle','--');endsubplot 612; plot(time,signal,'k'); ylabel('幅值');title('带噪信号'); axis([0 max(time) -1 1]);subplot 613; plot(time,output,'k');title('消噪信号'); axis([0 max(time) -1 1]); ylabel('幅值');subplot 614; plot(frameTime,Ef,'k'); hold on; ylabel('幅值');title('能熵比'); axis([0 max(time) 0 max(Ef)]);line([0 max(frameTime)],[T1 T1],'color','k','linestyle','--');for k=1 : vosl line([frameTime(voiceseg(k).begin) frameTime(voiceseg(k).begin)],... [-1 1],'color','k'); line([frameTime(voiceseg(k).end) frameTime(voiceseg(k).end)],... [-1 1],'color','k','linestyle','--');endtext(3.2,0.2,'T1');subplot 615; plot(frameTime,TT,'k');  ylabel('样点数');title('基音周期'); grid; axis([0 max(time) 0 80]);subplot 616; plot(frameTime,FF,'k'); ylabel('频率/Hz')title('基音频率'); grid; axis([0 max(time) 0 450]); xlabel('时间/s'); 

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));   % 在Pmin～Pmax范围内寻找最大值 period(k+ixb-1)=lmin+tloc-1;      % 给出对应最大值的延迟量    endend

13. MATLAB仿真十二：带噪语音的基音检测三

%% pr8_7_3   clear all; clc; close all;filedir=[];   % 设置语音文件路径filename='tone4.wav';% 设置文件名fle=[filedir filename]      % 构成路径和文件名的字符串SNR=0; % 设置信噪比IS=0.15;      % 设置前导无话段长度wlen=240;     % 设置帧长为25msinc=80;% 设置帧移10ms% [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;     % 设置时间lmin=floor(fs/500);  % 最高基音频率500Hz对应的基音周期lmax=floor(fs/60);   % 最高基音频率60Hz对应的基音周期 signal=Gnoisegen(x,SNR);    % 叠加噪声overlap=wlen-inc;    % 求重叠区长度NIS=fix((IS*fs-wlen)/inc +1);      % 求前导无话段帧数snr1=SNR_singlech(x,signal) % 计算初始信噪比值a=3; b=0.01;  % 设置a和boutput=simplesubspec(signal,wlen,inc,NIS,a,b);% 谱减法减噪snr2=SNR_singlech(x,output) % 计算谱减后信噪比值yy  = enframe(output,wlen,inc)';   % 滤波后信号分帧time = (0 : length(x)-1)/fs;% 计算时间坐标fn=size(yy,2);frameTime=frame2time(fn,wlen,inc,fs);Thr1=0.12;    % 设置端点检测阈值r2=0.5;% 设置元音主体检测的比例常数ThrC=[10 15]; % 设置相邻基音周期间的阈值% 用于基音检测的端点检测和元音主体检测[voiceseg,vosl,vseg,vsl,Thr2,Bth,SF,Ef]=pitch_vads(yy,fn,Thr1,r2,10,8);% 60～500Hz的带通滤波器系数b=[0.012280   -0.039508   0.042177   0.000000   -0.042177   0.039508   -0.012280];a=[1.000000   -5.527146   12.854342   -16.110307   11.479789   -4.410179   0.713507];x=filter(b,a,xx);    % 数字滤波x=x/max(abs(x));     % 幅值归一化y=enframe(x,wlen,inc)';     % 第二次分帧lmax=floor(fs/60);   % 基音周期的最小值lmin=floor(fs/500);  % 基音周期的最大值[Extseg,Dlong]=Extoam(voiceseg,vseg,vosl,vsl,Bth);% 计算延伸区间和延伸长度T1=ACF_corrbpa(y,fn,vseg,vsl,lmax,lmin,ThrC(1));  % 对元音主体进行基音检测% 对语音主体的前后向过渡区延伸基音检测 T0=zeros(1,fn);      % 初始化F0=zeros(1,fn);for k=1 : vsl % 共有vsl个元音主体    ix1=vseg(k).begin;      % 第k个元音主体开始位置    ix2=vseg(k).end; % 第k个元音主体结束位置    in1=Extseg(k).begin;    % 第k个元音主体前向延伸开始位置    in2=Extseg(k).end;      % 第k个元音主体后向延伸结束位置    ixl1=Dlong(k,1); % 前向延伸长度    ixl2=Dlong(k,2); % 后向延伸长度    if ixl1>0 % 需要前向延伸基音检测 [Bt,voiceseg]=back_Ext_shtpm1(y,fn,voiceseg,Bth,ix1,... ixl1,T1,k,lmax,lmin,ThrC);    else      % 不需要前向延伸基音检测 Bt=[];    end    if ixl2>0 % 需要后向延伸基音检测 [Ft,voiceseg]=fore_Ext_shtpm1(y,fn,voiceseg,Bth,ix2,... ixl2,vsl,T1,k,lmax,lmin,ThrC);    else      % 不需要后向延伸基音检测 Ft=[];    end    T0(in1:in2)=[Bt T1(ix1:ix2) Ft];      % 第k个元音主体完成了前后向延伸检测 endtindex=find(T0>lmax);% 限止延伸后最大基音周期值不超越lmaxT0(tindex)=lmax;tindex=find(T0<lmin & T0~=0);      % 限止延伸后最小基音周期值不低于lminT0(tindex)=lmin;tindex=find(T0~=0);F0(tindex)=fs./T0(tindex);   % 求出基音频率TT=pitfilterm1(T0,Extseg,vsl);      % T0平滑滤波FF=pitfilterm1(F0,Extseg,vsl);      % F0平滑滤波% STFT分析,绘制语谱图nfft=512;d=stftms(xx,wlen,nfft,inc);W2=1+nfft/2;n2=1:W2;freq=(n2-1)*fs/nfft;% 作图pos = get(gcf,'Position');set(gcf,'Position',[pos(1), pos(2)-150,pos(3),pos(4)+150]);subplot 611; plot(time,xx,'k'); ylabel('幅值');title('原始信号'); axis([0 max(time) -1 1]);for k=1 : vosl line([frameTime(voiceseg(k).begin) frameTime(voiceseg(k).begin)],... [-1 1],'color','k'); line([frameTime(voiceseg(k).end) frameTime(voiceseg(k).end)],... [-1 1],'color','k','linestyle','--');endsubplot 612; plot(time,signal,'k'); ylabel('幅值');title('带噪信号'); axis([0 max(time) -1 1]);subplot 613; plot(time,output,'k'); ylabel('幅值');title('消噪信号'); axis([0 max(time) -1 1]);subplot 614; plot(frameTime,Ef,'k'); hold on; ylabel('幅值');title('能熵比'); axis([0 max(time) 0 max(Ef)]);line([0 max(frameTime)],[Thr1 Thr1],'color','k','linestyle','--');plot(frameTime,Thr2,'k','linewidth',2);for k=1 : vsl line([frameTime(vseg(k).begin) frameTime(vseg(k).begin)],... [0 max(Ef)],'color','k','linestyle','-.'); line([frameTime(vseg(k).end) frameTime(vseg(k).end)],... [0 max(Ef)],'color','k','linestyle','-.');endtext(3.2,0.2,'Thr1');text(2.9,0.55,'Thr2');subplot 615; plot(frameTime,TT,'k'); ylabel('样点数'); title('基音周期'); grid; axis([0 max(time) 0 80]);subplot 616; plot(frameTime,FF,'k'); ylabel('频率/Hz'); title('基音频率'); grid; axis([0 max(time) 0 350]); xlabel('时间/s'); figure(2)pos = get(gcf,'Position');set(gcf,'Position',[pos(1), pos(2)-100,pos(3),(pos(4)-240)]);imagesc(frameTime,freq,abs(d(n2,:)));  axis xym = 64; LightYellow = [0.6 0.6 0.6];MidRed = [0 0 0]; Black = [0.5 0.7 1];Colors = [LightYellow; MidRed; Black];colormap(SpecColorMap(m,Colors));hold onplot(frameTime,FF,'w');ylim([0 1000]);title('语谱图上的基音频率曲线');xlabel('时间/s'); ylabel('频率/Hz')

function period=ACF_corrbpa(y,fn,vseg,vsl,lmax,lmin,ThrC)pn=size(y,2);if pn~=fn, y=y'; end % 把y转换为每列数据表示一帧语音信号wlen=size(y,1);      % 帧长period=zeros(1,fn);  % 初始化c1=ThrC(1);   % 取得阈值for i=1 : vsl % 只对有话段数据处理    ixb=vseg(i).begin;    ixe=vseg(i).end;    ixd=ixe-ixb+1;   % 求取一段有话段帧的帧数    Ptk=zeros(3,ixd);% 初始化    Vtk=zeros(3,ixd);    Ksl=zeros(1,3);    for k=1 : ixd     u=y(:,k+ixb-1);     % 取来一帧信号 ru=xcorr(u,'coeff');% 计算自相关函数 ru=ru(wlen:end);    % 取正延迟量部分 [Sv,Kv]=findmaxesm3(ru,lmax,lmin);% 获取三个极大值的数值和位置 lkv=length(Kv); Ptk(1:lkv,k)=Kv';   % 把每帧三个极大值位置存放在Ptk数组中    end    Kal=Ptk(1,:);    meanx=mean(Kal); % 计算Kal均值    thegma=std(Kal); % 计算Kal标准差    mt1=meanx+thegma;    mt2=meanx-thegma;    if thegma>5,  %判断基音候选组中有哪几个在第一次置信区域内 Ptemp=zeros(size(Ptk)); Ptemp(1,(Ptk(1,:)<mt1 & Ptk(1,:)>mt2))=1; Ptemp(2,(Ptk(2,:)<mt1 & Ptk(2,:)>mt2))=1; Ptemp(3,(Ptk(3,:)<mt1 & Ptk(3,:)>mt2))=1; % 如果第一组或(和)其他组都有值在第一次置信区内,取数值大的一个值赋于Pam Pam=zeros(1,ixd); for k=1 : ixd     if Ptemp(1,k)==1  Pam(k)=max(Ptk(:,k).*Ptemp(:,k));     end end mdex=find(Pam~=0);      % 在Pam非零的数值中 meanx=mean(Pam(mdex));  % 计算Pam均值 thegma=std(Pam(mdex));  % 计算Pam标准差 if thegma<0.5, thegma=0.5; end mt1=meanx+thegma; mt2=meanx-thegma;% 计算了第二次置信区 pindex=find(Pam<mt1 & Pam>mt2);% 寻找在置信区间的数据点 Pamtmp=zeros(1,ixd); Pamtmp(pindex)=Pam(pindex);    % 设置Pamtmp if length(pindex)~=ixd     bpseg=findSegment(pindex); % 计算置信区间内的数据分段信息     bpl=length(bpseg);  % 置信区间内的数据分成几段      bdb=bpseg(1).begin; % 置信区间内第一段的开始位置     if bdb~=1    % 如果置信区间内第一段开始位置不为1  Ptb=Pamtmp(bdb);% 置信区间内第一段开始位置的基音周期  Ptbp=Pamtmp(bdb+1);  pdb=ztcont11(Ptk,bdb,Ptb,Ptbp,c1);% 将调用ztcont11函数处理  Pam(1:bdb-1)=pdb;      % 把处理后的数据放入Pam中     end     if bpl>=2  for k=1 : bpl-1 % 如果有中间数据在置信区间外      pdb=bpseg(k).end;      pde=bpseg(k+1).begin;      Ptb=Pamtmp(pdb);      Pte=Pamtmp(pde);      pdm=ztcont21(Ptk,pdb,pde,Ptb,Pte,c1);  % 将调用ztcont21函数处理      Pam(pdb+1:pde-1)=pdm;     % 把处理后的数据放入Pam中  end     end     bde=bpseg(bpl).end;     Pte=Pamtmp(bde);     Pten=Pamtmp(bde-1);     if bde~=ixd  % 如果置信区间内最后一段结束位置不为ixd  pde=ztcont31(Ptk,bde,Pte,Pten,c1);% 将调用ztcont31函数处理  Pam(bde+1:ixd)=pde;    % 把处理后的数据放入Pam中     end end period(ixb:ixe)=Pam; else     period(ixb:ixe)=Kal;    endend

小结

基音周期是语音信号最重要的参数之一，它描述了语音激励源的一个重要特征。因为汉语是一种有调语言，基音的变化模式称为声调，它携带非常重要的具有辨义作用的信息，有区别意义的功能。所以，基音的提取和估算对汉语分析更是一个十分重要的课题。

但由于人的声道的易变性以及声道特征是因人而异的，而基音周期的范围又很宽，且同一个人在不同情感条件下发音的基音周期也不同，加之基音周期还受到单词发音音调的影响，因而基音周期的精确检测实际上是一件比较困难的事情。基音提取的主要困难反映在以下几个方面：

（1）声门激励信号并不是一个完全周期的序列，在某个音的头、尾部并不具有声带振动那样的周期性，有些清音和浊音的过渡帧是很难被确切地判断为是周期性还是非周期性的。

（2）基音频率处在100~200 Hz 的情况占多数，浊音信号往往可能包含有三四十次谐波分戴，而其基波分量往往不是最强的分量，这造成了检测基音时会把谐波当做基波。

（3）语音的第一共振峰通常在300~1 000 Hz 范围内，形成2~8次谐波分量常常比基波分
量还强，这也是造成误判的因素。

（4） 基音周期变化范围大，从老年男性的50 Hz到儿童和女性的500 Hz，接近三个多倍频程，给基音检测带来了一定的困难。由千这些困难，所以迄今为止尚未找到一个完善的方法可以对千各类人群（包括男、女、老、幼及不同语种）、各类应用领域和各种环境条件都能获得满意的检测方法。

尽管基音检测有许多困难，但因为它的重要性，基音的检测提取一直是一个重要的研究课题，为此人们提出了各种各样的基音检测算法。本章先介绍纯语音信号的基音检测方法：自相关函数(ACF)法、平均幅度差函数CAMDF)法、倒谱法、线性预测法、小波法等，主体－延伸法的基音检测模式，最后介绍了带噪语音的基音检测。

对本章内容感兴趣或者想充分学习了解的，建议去研习书中第八章节的内容。后期会对其中一些知识点在自己理解的基础上进行讨论补充，欢迎大家一起学习交流。

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

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

网赚站