pkpk_1c.m

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%	pkpk_1c.m
%
%	jmw and Dr. Hu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


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

%disp('loading ./signal.mat from hard disk ...');
%load signal;
%=sprintf('     signal name is %s',name);
%disp(s);


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

% start timer
t_start = clock;

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);

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

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


wait_counter = 2;
set(wait_window_slider,'value',wait_counter);
set(wait_window_display,'String','Now calculating the GCI using coarse analysis');
pause(0.01);

for kf=2:nframe-2
   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;
   wait_counter = wait_counter+1;
   set(wait_window_slider,'value',wait_counter);
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;

set(wait_window_display,'String','GCI coarse analysis done!..Now calculating the GCIs using fine analysis ');
pause(0.01);

% 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.

wait_counter = 2;
set(wait_window_slider,'value',wait_counter);
for kf=2:nframe-2
   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;
   wait_counter = wait_counter+1;
set(wait_window_slider,'value',wait_counter);
end;


% Adjust the position of each GCI under the criteria that no two GCI's
% are located within 25 samples.
set(wait_window_display,'String','Peak Picking');
wait_counter = 2;
set(wait_window_slider,'value',wait_counter);
pause(0.01);
ploc1=zeros(size(ploc));
gci=[];

for kf=2:nframe-2
	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;
   wait_counter = wait_counter+1; 
   set(wait_window_slider,'value',wait_counter);
end;

set(wait_window_display,'String','Saving variables to hard disk');
pause(0.01);


s=sprintf('saving ./track/out/%s_g.mat file to hard disk ...',name);
s=sprintf('save ./track/out/%s_g.mat gci',name);
eval(s);
set(wait_window_display,'String','GCI detection completed...Thank you for being patient!');
pause(0.01);

close (wait_window);
clear wait_window;
clear wait_counter;