agif_go.m

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%				
% Automatic GIF program		
% Author : Minkyu Lee		
% Date : 12-Oct-1994		
%				
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Modified by Karthik on May 28 1997
% Modified by D. G. Childers

% Constants for general use

delta=0.25;
beta=0.75;
eta=1e-5;
Nf=256;			% Frame length for pitch Asynchronous analysis
No=56;			% Number of overlapping points for asyn. analysis

% Some extra points before and after the analysis frame seleted
ext_points=20;

LPF3=fir1(3,0.3);
HPF4=fir1(8,0.2,'high');

% Check the flags for various analysis method

asyn_flag=get(agif_cb_asyn_h,'Value');
iter_flag=get(agif_cb_iter_h,'Value');
arma_flag=get(agif_cb_arma_h,'Value');
clpc_flag=get(agif_cb_clpc_h,'Value');
cwrls_flag=get(agif_cb_cwrls_h,'Value');

% Get the analysis data 

xdata=speech(xsmin:xsmax);
ydata=xdata;
ydata=filter(LPF3,1,ydata);

% For Asynchronous analysis, find out the analysis frame
if asyn_flag == 1
   tlen=round(xsmax-xsmin);
   if tlen < Nf
      disp('The analysis frame is too short for pitch asynchronous analysis');
      asyn_flag=0;
      set(agif_cb_asyn_h,'Value',asyn_flag);
   end
   num_frame=round(tlen/(Nf-No));
end

% For pitch synchronous analysis, get the glottal closure instants 
% And find out the number of analysis frames

if iter_flag == 1 | arma_flag == 1 | clpc_flag == 1 | cwrls_flag == 1
   S=sprintf('./output/%s_GCI.mat',name);
   if exist(S) == 0
      disp('Glottal closure detection')
      B1;
   else
      disp('Glottal closure detection is done already');
   end
   disp(' ');
   disp('Loading glottal closure instance data')
   S=sprintf('load ./output/%s_GCI.mat',name);
   eval(S);
   
   num_gci=length(gci);
   for ii=1:num_gci,
      if gci(ii) > round(xsmin)
         break;
      end
   end
   
   for jj=ii:num_gci,
      if gci(jj) > round(xsmax);
         break;	
      end
   end
   
   clear gci_a;
   for kk=1:jj-ii+1,
      gci_a(kk)=gci(ii+kk)-round(xsmin);
   end
   num_gci=length(gci_a);
end

% initialize plotting windows


if exist('pzplot_flag') == 1 & pzplot_flag == 1 
   PV=[150 20 400 500];
   s2 = 'Pole Zero Plot';
   while exist('agif_pzplot_win_h')==1
      try1 = 'get(agif_pzplot_win_h,''position'');';
      eval(try1,catch2);
      if check ==0
         clear agif_pzplot_win_h;
         check = 1;
         break;
      end
      s1 = get(agif_pzplot_win_h,'Name');
      if ~strcmp(s1,s2)
         clear agif_pzplot_win_h;
         break;
      end
      figure(agif_pzplot_win_h);
      break;
   end;
   
   if exist('agif_pzplot_win_h')~=1;
      agif_pzplot_win_h=figure('Position',PV,...
         'Resize','on',...
         'Numbertitle','off',...
         'Name',s2);
   end;
end;

if exist('spec_flag') == 1 & spec_flag == 1 
   PV=[150 20 400 500];
   s2 = 'Spectrogram - Glottal Inverse Filtering';
   % Open analysis window
   while exist('agif_spec_win_h')==1
      try1 = 'get(agif_spec_win_h ,''position'');';
      eval(try1,catch2);
      if check ==0
         clear agif_spec_win_h;
         check = 1;
         break;
      end
      s1 = get(agif_spec_win_h,'Name');
      if ~strcmp(s1,s2)
         clear agif_spec_win_h;
         break;
      end
      figure(agif_spec_win_h);
      break;
   end;
   
   if exist('agif_spec_win_h')~=1;
      agif_spec_win_h=figure('Position',PV,...
         'Resize','on',...
         'Numbertitle','off',...
         'Name',s2);
   end
end;

if asyn_flag == 1 
   PV=[270 44 515 268];
   s2 = 'Output Asynchronous LPC';
   
   % Open analysis window
   while exist('agif_asyn_win_h')==1
      try1 = 'get(agif_asyn_win_h,''position'');';
      eval(try1,catch2);
      if check ==0
         clear agif_asyn_win_h;
         check = 1;
         break;
      end
      s1 = get(agif_asyn_win_h,'Name');
      if ~strcmp(s1,s2)
         clear agif_asyn_win_h;
         break;
      end
      figure(agif_asyn_win_h);
      break;
   end;

   if exist('agif_asyn_win_h')~=1;
      agif_asyn_win_h=figure('Position',PV,...
         '	Resize','on',...
         'Numbertitle','off',...
         'Name',s2);
   end
end

   
if iter_flag == 1 
   PV=[270 44 515 268];
   s2 = 'Output Iterative GIF';
   
   % Open analysis window
   while exist('agif_iter_win_h')==1
      try1 = 'get(agif_iter_win_h,''position'');';
      eval(try1,catch2);
      if check ==0
         clear agif_iter_win_h;
         check = 1;
         break;
      end
      s1 = get(agif_iter_win_h,'Name');
      if ~strcmp(s1,s2)
         clear agif_iter_win_h;
         break;
      end
      figure(agif_iter_win_h);
      break;
   end;
   
   if exist('agif_iter_win_h')~=1;
      agif_iter_win_h=figure('Position',PV,...
         'Resize','on',...
         'Numbertitle','off',...
         'Name',s2);
   end
end


if arma_flag == 1 
   PV=[270 44 515 268];
   s2 = 'Output(ARMA)';

   % Open analysis window
   while exist('agif_arma_win_h')==1
      try1 = 'get(agif_arma_win_h,''position'');';
      eval(try1,catch2);
      if check ==0
         clear agif_arma_win_h;
         check = 1;
         break;
      end
      s1 = get(agif_arma_win_h,'Name');
      if ~strcmp(s1,s2)
         clear agif_arma_win_h;
         break;
      end
      figure(agif_arma_win_h);
      break;
   end;
   
   if exist('agif_arma_win_h')~=1;
      agif_arma_win_h=figure('Position',PV,...
         'Resize','on',...
         'Numbertitle','off',...
         'Name',s2);
   end
end

if clpc_flag == 1 
   PV=[270 44 515 268];
   s2 = 'Output (Closed LP)';
   
   % Open analysis window
   while exist('agif_clpc_win_h')==1
      try1 = 'get(agif_clpc_win_h,''position'');';
      eval(try1,catch2);
      if check ==0
         clear agif_clpc_win_h;
         check = 1;
         break;
      end
      s1 = get(agif_clpc_win_h,'Name');
      if ~strcmp(s1,s2)
         clear agif_clpc_win_h;
         break;
      end
      figure(agif_clpc_win_h);
      break;
   end;
   
   if exist('agif_clpc_win_h')~=1;
      agif_clpc_win_h=figure('Position',PV,...
         'Resize','on',...
         'Numbertitle','off',...
         'Name',s2);
   end
end

pma=ones(1,num_zeros+1);
par=ones(1,num_poles);

disp(' ');
disp('Now, doing glottal inverse filtering using specified method(s) ');

% 
% Pitch asynchronous analysis
%

if asyn_flag == 1
	clear dg_final g_final;
   for iter=1:num_frame,
      bp=(iter-1)*(Nf-No)+1;
      ep=bp+Nf-1;
      if ep > tlen
         ep=tlen;
      end
      
      s=ydata(bp+ext_points:ep+ext_points);
      
      % do a fixed pre-emphasis
      ds=filter([1 -0.9], 1, s);

      % apply a hamming window
      ds_w=hamming(length(ds)).*ds;
      
      % do LP analysis on the frame
      A=lpc(ds_w,num_poles);
      [H,W]=freqz( 1, A, 512);
      
      if exist('pzplot_flag') == 1 & pzplot_flag == 1
         figure(agif_pzplot_win_h);
         zplane(1, A);
         title('Pole-Zero plot for Asynchronous LPC analysis');
      end
      
      if exist('spec_flag') == 1 & spec_flag == 1
         figure(agif_spec_win_h);
         plot(20*log(abs(H)));
			title('Spectrum for Asynchronous LPC analysis');
		end
	
		
		%Get a differentiated glottal wave estimate
		dg=filter(A, 1, s);

		% get a glottal wave estimate
		isum=0;
		for i=1:length(dg),
			isum=isum+dg(i);
			g(i)=isum;
		end
      dg_out1=filter([1 -1],[1 -0.99],dg);
      dg_final(bp+No:ep)=dg_out1(No+1:length(dg_out1));
   end
   
	V=axis;
	axis([1 length(dg_final) V(3) V(4)]);

	% Remove components around d.c.
	dg_final=filter([1 -1],[1 -.99],dg_final);
   dg_final=filter(LPF3,1,dg_final);
   figure(agif_asyn_win_h);
	plot(dg_final) ;
	title('dgvv by pitch asynchronous analysis (AR model)') 
end

%
% iterative GIF analysis
%

if iter_flag == 1

	dg_final=[];
	g_final=[];
	ext_left=200;
	ext_right=200;
	LPF3=fir1(3,0.5);
	shift=0;

   for cur_gci=1:num_gci-1
      ydata_a=speech(gci_a(cur_gci)+xsmin-ext_left+shift:...
         gci_a(cur_gci+1)+xsmin+ext_right+shift);
      ydata_p=filter([1 -0.9],1,ydata_a);
      ydata_lp=filter(LPF3,1,ydata_a);
      
      % Estimate glottal effect to speech
      ydata_w=hamming(length(ydata_p)).*ydata_p;
      %Hg1=lpc_cov(ydata_a,1);
      Hg1=lpc(ydata_w,1);
      
      % Eliminate the estimated glottal contribution	
      tmp=filter(Hg1, 1, ydata_lp);
      tmp_w=hamming(length(tmp)).*tmp;
      
      % The first estimate for the vocal tract
      Hvt1=lpc(tmp_w,num_poles);
      % Eliminate the effect of vocal tract
      tmp=filter(Hvt1,1,ydata_lp);
      % The first estimate for the glottal excitation
      g1=integ_1a(tmp);
      g1_w=hamming(length(g1)).*g1;
      %Second iteration begins here
      Hg2=lpc(g1_w,4);
      % Eliminate the estimated glottal contribution
      tmp=filter(Hg2,1,ydata_lp);
      tmp_w=hamming(length(tmp)).*tmp;
      % The final model for the vocal tract
      Hvt2=lpc(tmp_w,num_poles);

      %Fine tune the first formant and bandwidth
      dg2=filter(Hvt2,1,ydata_lp);
      dg2=filter([1 -1],[1 -.99],dg2);
      % The final estimate for the glottal excitation
      g2=integ_1a(dg2);
      g2=filter([1 -1],[1 -.99],g2);
      
      if exist('pzplot_flag') == 1 & pzplot_flag == 1
         figure(agif_pzplot_win_h);
         zplane(1, Hvt2);
         title('Pole-Zero plot for ARMA analysis');
      end
      
      if exist('spec_flag') == 1 & spec_flag == 1
         Fs_a=fft(ydata_p,1024);
         [h_vt w_vt]=freqz(1,Hvt2,512);
         figure(agif_spec_win_h);
         plot(20*log(abs(h_vt))+20*log(abs(Fs_a(1))),'r');
         hold on;
         plot(20*log(abs(Fs_a(1:512))));
         hold off
         title('Spectrum for iterative analysis');
      end
      
      % Connect each analysis frame
      dg_final=[dg_final ; dg2(ext_left+1:length(dg2)-ext_right)];
      g_final=[g_final ; g2(ext_left+1:length(g2)-ext_right)];
      % Plot
      figure(agif_iter_win_h);
      subplot(311)
      plot(dg2(ext_left+1:length(dg2)-ext_right));
		title('dgvv by iteratice analysis (AR model)');
      subplot(312)
      plot(g2(ext_left+1:length(g2)-ext_right));
		title('gvv by iterative analysis');
		pause(0.01);
	end
	dg_final = filter([1 -1],[1 -.99],dg_final);
   dg_final = filter(LPF3,1,dg_final);
   figure(agif_iter_win_h);
   subplot(313);
   plot(dg_final);
   title('Low pass filtered dgvv');
	clear dg_final g_final;
end

%
% ARMA model analysis
%

if arma_flag == 1

	dg_final=zeros(size(xdata));
	dg_mfinal=zeros(size(xdata));
	g_final=zeros(size(xdata));
	g_mfinal=zeros(size(xdata));
	nshift=20;
	
        Zi=ydata(gci_a(1)-1:-1:gci_a(1)-num_poles);
	for iter=1:num_gci-2
		bp_ext=gci_a(iter)+1-ext_points;
		bp_n=gci_a(iter)+1;
		ep_n=gci_a(iter+1);
		Nf=ep_n-bp_n+1;
		ext_sample=bp_n-bp_ext;
		
		s_n=ydata(bp_ext:ep_n);
		s_n1=s_n(ext_points+1:length(s_n));

		% do a fixed pre-emphasis
		ds_n=filter([1 -0.9], 1, s_n);
		% do LP analysis on the frame
		A=lpc_cov(ds_n(ext_points-num_poles+1:length(ds_n)),num_poles);
      % get a differentiated glottal wave estimate
      dg_n=filter(A, 1, s_n);
      dg_n1=dg_n(ext_points+1:length(dg_n));
      [minv mindex]=min(dg_n1);
      mindex=mindex+ext_points;
      g_n=integ_1a(dg_n);
      g_n1=g_n(ext_points+1:length(g_n));
   end
end