pkpk_1b.m

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%	pkpk_1b.m
%
%	jmw and Dr. Hu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp(' ');
disp('SCRIPT: pkpk_1b.m *************************************');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% load variables

file_string=sprintf('./temp/%s.mat',name);
s=sprintf('loading %s from hard disk ...',file_string);
disp(s);
s=sprintf('load %s',file_string);
eval(s);

signal=eval(name);

s=sprintf('clear %s',name);
eval(s);

disp('loading ./temp.mat from hard disk ...');
load temp/temp;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ploc=zeros(nframe,10);
ppa=zeros(1,nframe);

ir=integ_1a(residue);   % ir = integrated residue
ir=filtfilt([1 -1],[1 -.99],ir);

% Zero-phase filtering
ir=filtfilt([1 -1],conv([1 -.9],[1 -.7]),ir);

disp('Doing coarse analysis ...');

%  create message window and write message in that window also
   message_win3_f=figure('Unit','normalized',...
           'Position',[0.05 0.2 0.5 0.2],...
           'Resize','off',...
           'Color',BACK_COLOR,...
           'Numbertitle','off',...
           'Name','Message');
axis('off');        
ss=sprintf('Doing coarse analysis ...');
text(-0.05,1,ss,'color',[0 0 1],'FontSize',10);
pause(0.01);        


% First step : estimate the pitch period using a frame-based approach
% The method used here is similar to the cepstrum.
% For details, please see hu's dissertation.

% The result is recorded in the variable, "ppa".


for kf=2:nframe-1
    if voicetype(kf) 
	rr=ir(200*(kf-1)+14-156:200*kf+13+156);
	ss=signal(200*(kf-1)+14-156:200*kf+13+156);
	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;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Correct the erroneous pitch periods

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 analysis ok.')
disp('Finding glottal closure instants ...');

%print new message in message window
ss=sprintf('Now finding glottal closure instants ...');
text(-0.05,1,ss,'color',[0 0 1],'FontSize',10);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Identify the glottal closure intants (GCI's)
% Results are shown as a (nframe x 10) matrix, "ploc", in which each element 
% indicating the position of the GCI in the current frame.

for kf=2:nframe-1
   if voicetype(kf)
	ss=signal(200*(kf-1)+14-100:200*kf+113);
	rr=residue(200*(kf-1)+14-100:200*kf+113);
	irload=ir(200*(kf-1)+14-100:200*kf+113);
	[m nrr]=min(irload(101:300));
	nrr=nrr+100;
	if nrr==101
		while irload(nrr-1)<m
			m=irload(nrr-1);
			nrr=nrr-1;
		end;
	elseif nrr==300
		while irload(nrr+1)<m
			m=irload(nrr+1);
			nrr=nrr+1;
		end;
	end;
	irthd=irload(nrr)/3;	
	nadd=1;
	while irload(nrr+nadd)<irthd
		nadd=nadd+1;
	end;
	x1=0:nadd;
	y1=irload(nrr:nrr+nadd);
	aa=polyfit(x1,y1,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=200*(kf-1)+13;
	ptype=ir(startp+nrr-100-15:startp+nrr+30-100);
	for k=1-10:260
		r1=ir(startp+(k-15:k+30));
		rs(k+10)=ptype*r1';
	end;	
	nrs=nrr+10-100;
	pp=ppa(kf);
	pph=floor(pp*.35);
	mid=nrs;
	pk=mid+pp;
	pptmp1=[];
	k1=1;
	while pk < 260
		 p_range=pk-pph:min([pk+pph 270]);
		[m n]=max(rs(p_range));
		pptmp1(k1)=n-1+pk-pph;
		pk= pptmp1(k1)+pp;
		k1=k1+1;
	end;
	pk=mid-pp;
	pptmp2=[];
	k2=1;
	while pk > 10
		p_range=max([1 pk-pph]):pk+pph;
		[m n]=max(rs(p_range));
		pptmp2(k2)=n-1+max([1 pk-pph]);
		pk= pptmp2(k2)-pp;
		k2=k2+1;
	end;
	pptmp=[rev_1a(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;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Adjust the position of each GCI under the criteria that no two GCI's
% are located within 25 samples.

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>200)-200);
		ppup=ppup(ppup>0 & ppup<=200);
		if length(ppmid)==0
			ppadd=max(ppup)-200;
		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)*200+13-1];
	end;
end;
disp('     GCI peak picking ok.');

%print new message in message window
ss=sprintf('GCI peak picking ok');
text(-0.05,1,ss,'color',[0 0 1],'FontSize',10);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

s=sprintf('saving ./%s_GCI.mat file to hard disk ...',name);
disp(s);
s=sprintf('save ./temp/%s_GCI.mat gci',name);
eval(s);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% save ir
ir_save=ir;

% clean up time

clear a aa b bb cofa devia Element emp energy ir irload irthd k k1 k2 kf
clear m mid mm n nadd nframe nn nrr nrs nsub ntotal p_range pk ploc
clear ploc1 pp ppa ppa1 ppadd ppdown pph pplast ppmid ppmid1
clear ppno ppsmooth pptmp pptmp1 pptmp2 ppup ptype r r1 rc residue rr 
clear rs s ss startp t_start voicetype x1 x2 y1 y2 file_string 
clear signal gci 

close (message_win3_f);
clear message_win3_f;