garea_m2.m

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

%garea_m2
%This use to be a function
%
% ***********************************************************************
%modified by D. G. Childes 5/20/98

global catch2;
global check;

clear D_POSL d_posl D_POSR d_posr D_Ag dpar par_r par
clear car dcar degg ddegg egg


modlpar=MODEL_PARAMETERS;

 L=modlpar(1)/10;
 W=modlpar(2)/10;
 T=modlpar(3)/10;
 vpd=modlpar(4);
 hpd=modlpar(5);
 xim=modlpar(6)/10;
 f0=modlpar(7);
 qa=modlpar(8);
 qs=modlpar(9);
 factor=100/W;
 fs=10000;
 xi2=(qa*xim)*factor;
 xi1=(qa+qs)*xim*factor;
 nz=15;
 ny=15;
 ccons=0.1;
 acons=0.1;
 tcons=2*pi*f0*(1/fs);
 npoint=round(fs/f0);

dy=200/ny;
dz=50/nz;
da=dy*dz;
sy=0:(200/ny):200;
sz=0:(50/nz):50;
dyw=L/ny;
dzw=T/nz;
daw=dyw*dzw;
[notisldr,noteh]=Wait_message(0);

for nt=1:npoint
  parea=0;
  gcar=0;
  set(notisldr,'value',[nt/npoint]);
  drawnow;
  for i=1:(ny+1)
     sdis=10^8;
     for j=1:(nz+1)
         sdisp(i,j)=( xi1-(xi1-xi2)*(sz(j)/50))*(1-(sy(i)/200));
         dispr(i,j)=sdisp(i,j)+100;
         displ(i,j)=100-sdisp(i,j);
         SINE1=sin(tcons*nt-vpd*(sz(j)/50)/(pi)-hpd*(sy(i)/200)/(pi)+3*pi/4);
         SINE=sin(pi*sy(i)/200)*SINE1;
         posr(i,j)=dispr(i,j)+xim*SINE*factor;
         posl(i,j)=displ(i,j)-xim*SINE*factor;
         cdis=(posr(i,j)-posl(i,j))/factor;
         if  (cdis<=0)
             gcar=gcar+1;
         end
         if  (cdis<sdis)
             sdis=cdis;
         end 

         if (posr(i,j)<=100)
            posr(i,j)=100;
         end
        
         if (posl(i,j)>=100)
            posl(i,j)=100;
         end
        
     end

     if  (sdis>0)
         parea=parea+sdis;
     end
    
  end
  
  if nt==1
     d_posl=[posl];
     d_posr=[posr];
  else
     d_posl=[d_posl;posl];
     d_posr=[d_posr;posr];
  end
  dpar(nt)=parea*dyw;
  dcar(nt)=gcar*daw;
  degg(nt)=ccons/(acons+(gcar*daw));
  
end

delete(noteh);

ddegg=diff(degg);
par=[dpar dpar dpar dpar]*100;
car=[dcar dcar dcar dcar]*100;
egg=[degg degg degg degg];

%
D_POSL=d_posl;
D_POSR=d_posr;
D_Ag=dpar;
par_r=par;
%


PV =[8 29 366 519];
s2 = 'Glottal Area Waveforms';
% Open analysis window
while exist('rbn1_output_window_3')==1
   try1 = 'get(rbn1_output_window_3,''position'');';
   eval(try1,catch2);
   if check ==0
      clear rbn1_output_window_3;
      check = 1;
      break;
   end
   s1 = get(rbn1_output_window_3,'Name');
   if ~strcmp(s1,s2)
      clear rbn1_output_window_3;
      break;
   end
   figure(rbn1_output_window_3);
   break;
end;

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

tscal=0:1000/fs:(length(par)-1)*1000/fs;   
subplot(3,1,1);
plot(tscal,par);
grid
ht=max(par)*1.1;
VV=axis;
axis([0  VV(2)  0  ht])
qa1=round(qa*1000)/1000;
qs1=round(qs*1000)/1000;
xim1=round(xim*10*1000)/1000;
na=num2str(qa1);
ns=num2str(qs1);
nim=num2str(xim1);
pv=num2str(vpd);
ph=num2str(hpd);
f0num=num2str(f0);

title( 'Projected glottal area' )
xlabel('Time (ms)')
ylabel('Unit=mm*mm')
subplot(3,1,2)
plot(tscal,car)
grid
ht=max(car)*1.1;
VV=axis;
axis([0  VV(2)  0  ht])
sL=num2str(L*10);
sW=num2str(W*10);
sT=num2str(T*10);
title('Vocal fold contact area');
xlabel('Time (ms)')
ylabel('Unit=mm*mm')
subplot(3,1,3)
plot(tscal,egg)
grid
ht=max(egg)*1.1;
VV=axis;
axis([0  VV(2)  0  ht])
title('Simulated EGG waveform ')
xlabel('Time (ms)')