pitch_determine.m

这是一个用于语音信号处理的工具箱

%Program from the analysis toolboxes
%Originally by jmw and Dr. Hu
%
%Modified by Karthik to fit the formant synthesis program

%Defining some constants that are used often

function [ppa,voicetype,gci] = Pitch_determine(excitation);

num_poles = 14;

wait_window = figure('Numbertitle','off',...
   'Color','white',...
   'Name','Please wait',...
   'Position',[117 248 560 150]);
pause(0.01)

signal=excitation';
wait_counter = 0;

% Normalize the amplitude of the analyzing signal

signal=14000/max(abs(signal))*signal;

% Remove components around d.c.
signal=filtfilt([1 -1],[1 -.99],signal);

% Segmentize the underlying signal

ntotal=length(signal);
old_ntotal = ntotal;
nframe=floor(ntotal/200);
if ntotal-nframe*200-13 > 0
	nframe=nframe+1;
end;
ntotal=nframe*200+13;
signal(ntotal)=0;

for i=old_ntotal+1:ntotal,
	signal(i)=rand(1) - 0.5;
end;

cofa=zeros(nframe,14);  
curmaxp=zeros(1,nframe);
compa=1e7;

wait_window_slider = uicontrol('Style','Slider',...
   'Position',[50 100 410 30],...
   'Max',nframe,...
   'Min',0,...
   'Value',wait_counter);

wait_window_display = uicontrol('Style','Text',...
   'Position',[50 50 410 30],...
   'Backgroundcolor','white',...
   'Foregroundcolor','blue',...
   'String','Please wait..This program is now calculating pitch periods from your signal');

pause(0.01);

% Calculate the lpc coefficients and get residue;
% Please see "lpc_1a.m" for the I/O arguments.

% ********************* do for FIRST segment *************************
sso=signal(1:250+13);
[cofa1,emp(1),energy(1),residue1]=lpc_1a(sso,num_poles-1);
residue1=sqrt(energy(1)/(residue1*residue1'))*residue1;
cofa(1,:)=cofa1;
residue(14:250+13)=residue1;
ss1=residue1(201:250);
wait_counter =1;
set(wait_window_slider,'value',wait_counter);
drawnow;
% ******************** do for middle segments ************************
for k=2:nframe-1

	sso=signal((k-1)*200+1:k*200+13+50);
	[cofa1,emp(k),energy(k),residue1]=lpc_1a(sso,num_poles-1);
	residue1=sqrt(energy(k)/(residue1*residue1'))*residue1;
	energy1=sso(1+13:100+13)*sso(1+13:100+13)'; 
	energy2=sso(101+13:200+13)*sso(101+13:200+13)';
   energyidx(k)=sqrt(energy1*energy2);
	ss2=residue1(1:50);
	ss=(ss1.*(50:-1:1)+ss2.*(1:50))/51;
	residue((k-1)*200+14:k*200+50+13)=[ss residue1(51:250)];
	ss1=residue1(201:250);
   cofa(k,1:length(cofa1))=cofa1;
   wait_counter = wait_counter+1;
   set(wait_window_slider,'value',wait_counter);
end;

% ***************** do for LAST segment *******************************
k=nframe;
sso=signal((k-1)*200+1:k*200+13);
[cofa1,emp(k),energy(k),residue1]=lpc_1a(sso,num_poles-1);
residue1=sqrt(energy(k)/(residue1*residue1'))*residue1;
energy1=sso(1+13:100+13)*sso(1+13:100+13)'; 
energy2=sso(101+14:200+13)*sso(101+14:200+13)';
energyidx(k)=sqrt(energy1*energy2);
ss2=residue1(1:50);
ss=(ss1.*(50:-1:1)+ss2.*(1:50))/51;
residue((k-1)*200+14:k*200+13)=[ss residue1(51:200)];
cofa(k,1:length(cofa1))=cofa1;
wait_counter = nframe;
set(wait_window_slider,'value',wait_counter);

set(wait_window_display,'String','LP Analysis over...Doing V/UV classification');
pause(0.01);
% *************************** V/U classification **************
voicetype(1)=0;
for kf=2:nframe-1
    if emp(kf)>.2 & energyidx(kf)>3*compa
    %if emp(kf)>.4 & energyidx(kf)>3*compa  % hu's fix to elim false V detection
	voicetype(kf)=1;
    else
	voicetype(kf)=0;
    end;
end;
voicetype(nframe)=0;

% ************************* eliminate UVU and VUV possibilities *******
for k=2:nframe-1
	if sum(voicetype(k-1:k+1)) > 1
		voicetype(k)=1;
	else
		voicetype(k)=0;
	end;
end;

set(wait_window_display,'String','Saving variables to temp.mat file');
pause(0.01);
save temp cofa energy residue emp voicetype nframe ntotal;

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


%s=sprintf('clear %s',name);
%eval(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 cepstrum.
% 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','GCI detection completed...Thank you for being patient!');
drawnow;
pause(1);

close (wait_window);