lpana2.m

来自「这是一个用于语音信号处理的工具箱」· M 代码 · 共 189 行

% Function: preform the second stage of the Linear Prediction Speech Analysis.
%           ==> locate the GCIs

function [gci,ir]=lpana2(signal,basic,voicetype,rsd,ntotal,nframe);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%                                %%
%%%%%%%  Find Glottal closure index    %%
%%%%%%%                                %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%retrieve the basic specification
F_len=basic(5);
O_lap=basic(6); 
Order=basic(4);
M_len=F_len-O_lap;

 % This step will create the following variables :
 %
 % 	ploc = glottal closure instant located in each frame.
 %	gci = glottal closre index
 %	ppa = estimated pitch period for each voiced frame.

% allocate the vector space
ploc=zeros(nframe,10);
ppa=zeros(1,nframe);  

% Smooth integrated residue.
ir1=integ(rsd);  
ir1=filtfilt([1 -1],[1 -.99],ir1);
ir=filtfilt([1 -1],conv([1 -.9],[1 -.7]),ir1); 

% Perform frame-based pitch estimation.
% For details, please refer to Section 4.1.3 of Hu's disseration.
for kf=2:nframe-1
    if voicetype(kf) 
	rr=ir(M_len*(kf-1)+Order+1+M_len/2-256:M_len*kf+Order-M_len/2+256);
	ss=signal(M_len*(kf-1)+Order+1+M_len/2-256:M_len*kf+Order-M_len/2+256);
	r=rr((1:256)*2);
	r=hanning(256)'.*r;
	rc=real(ifft(abs(fft(r))));
	rc=[zeros(1,15) rc(16:128)];
	[mm nn]=max(rc(1:125));
	[m n]=max(rc(1:nn-12));
	if m>.7*mm
		nn=n;
	end;		
	b=(rc(nn+1)-rc(nn-1))/2;
	a=(rc(nn+1)+rc(nn-1))/2-rc(nn);
	devia=-b/(2*a);
	ppa(kf)=round((nn-1+devia)*2);
    end;
end;

% Smooth the pitch contour to prevent unreasonable excursion.
ppsmooth=1;
if ppsmooth
	for kf=3:nframe-1
		if ppa(kf-1)&ppa(kf)
			ppa1=ppa(kf-2:kf-1);
			pplast=round(mean(ppa1(ppa1~=0)));
			ppno=max([round(ppa(kf)/pplast) 1]);
			if ppno~=1
				ppa(kf)=round(ppa(kf)/ppno);
			end;
		end;
	end;
end;

%disp('Coarse pitch estimation is done!');

for kf=2:nframe-1
    if voicetype(kf)
	ss=signal(M_len*(kf-1)+Order+1-M_len/2: M_len*kf+M_len/2+Order);
	rr=rsd(M_len*(kf-1)+Order+1-M_len/2: M_len*kf+M_len/2+Order);
	irload=ir(M_len*(kf-1)+Order+1-M_len/2: M_len*kf+M_len/2+Order);
	[m nrr]=min(irload(M_len/2+1:3*M_len/2));
	% Search for the minimum. The location of this minimum is considered
	% the first GCI.
	nrr=nrr+M_len/2;
	if nrr==M_len/2+1
		while irload(nrr-1)<m
			m=irload(nrr-1);
			nrr=nrr-1;
		end;
	elseif nrr==3*(M_len/2)
		while irload(nrr+1)<m
			m=irload(nrr+1);
			nrr=nrr+1;
		end;
	end;
	
	% Approximate the smoothed curve by two straight lines. 
	% Choose the intersection of these two lines as the GCI. 
	irthd=irload(nrr)/3;	
	nadd=1;
	while irload(nrr+nadd)<irthd
		nadd=nadd+1;
	end;
	xn1=0:nadd;
	yn1=irload(nrr:nrr+nadd);
	aa=polyfit(xn1,yn1,1);
	nsub=1;
	while irload(nrr-nsub)<irthd
		nsub=nsub+1;
	end;
	x2=-nsub:0;
	y2=irload(nrr-nsub:nrr);
	bb=polyfit(x2,y2,1);
	nrr=nrr+round((bb(2)-aa(2))/(aa(1)-bb(1)));	
	
	startp=M_len*(kf-1)+Order;
	ptype=ir(startp+nrr-M_len/2-15:startp+nrr+30-M_len/2); % prototype
	
	% Search for the other GCI's by examining the corresponding
	% crosscorrelation.
	for k=1-10:260
		r1=ir(startp+(k-15:k+30));
%		rs(k+10)=ptype*r1'/sqrt(r1*r1');
		rs(k+10)=ptype*r1';
	end;	
	nrs=nrr+10-M_len/2;
	pp=ppa(kf);
	pph=floor(pp*.35);
	mid=nrs;
	pk=mid+pp;
	pptmp1=[];
	k1=1;
	% Search forward and backward within the range from 10 to 260.
	while pk < (F_len+10)
		 range=pk-pph:min([pk+pph F_len+20]);
		[m n]=max(rs(range));
		pptmp1(k1)=n-1+pk-pph;
		pk= pptmp1(k1)+pp;
		k1=k1+1;
	end;
	pk=mid-pp;
	pptmp2=[];
	k2=1;
	while pk > 10
		range=max([1 pk-pph]):pk+pph;
		[m n]=max(rs(range));
		pptmp2(k2)=n-1+max([1 pk-pph]);
		pk= pptmp2(k2)-pp;
		k2=k2+1;
	end;
	pptmp=[rev(pptmp2) mid pptmp1]-10;
	pptmp=pptmp(pptmp>=1 & pptmp<=230);
	pptmp=pptmp(rs(pptmp+10)>rs(nrs)*.2);
	ploc(kf,1:length(pptmp))=pptmp;
    end;
end;

% Eliminate the redundancy occuring at frame boundaries and smooth the
% pitch period if necessary.

ploc1=zeros(size(ploc));
gci=[];
for kf=2:nframe-1
	if ploc(kf,1)
	    	if ploc(kf-1,1)==0
			ppup=ploc(kf,:);
	    	else
			ppup=ppdown;
		end;
		ppdown=ploc(kf+1,:);
		ppmid=min(ppup(ppup>M_len)-M_len);
		ppup=ppup(ppup>0 & ppup<=M_len);
		if length(ppmid)==0
			ppadd=max(ppup)-M_len;
		else
			ppadd=ppmid;
		end;
		ppdown=ppdown(ppdown>0);
		ppmid1=max(ppdown(ppdown <= ppadd+25));
		if length(ppmid1) & length(ppmid)
			ppmid=round((ppmid+ppmid1)/2);
		end;
		ppdown=[ppmid ppdown(ppdown>ppadd+25)];
		if length(ppup)
			ploc1(kf,1:length(ppup))=ppup;
		end;
		Element=ploc1(kf,find(ploc1(kf,:)~=0));
		gci=[gci Element+(kf-1)*M_len+Order];
	end;
end;

%disp('GCI identification is done!');