garea_calculate.m

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

%Function to calculate Glottal Area, Glottal Flow and displacements 
%of the two masses for Ishizaka and Flanagans, 1972 two mass model

%Author : Karthik
% modified by D. G. Childers 5/19-20/98

%The difference equations in Flanagans paper are implemented using this model.  
%All notation is as used in the original paper by Ishizaka and Flanagan
%For the constants, the values suggested by Ishizaka and Flanagan are used. 

function [UG,garea1,garea2,x1,x2] = garea_calculate(Parameters); 


%Get values of parameters
Length = Parameters(1);
Pressure = Parameters(4);
m1 = Parameters(5);
m2 = Parameters(6);
d1 = Parameters(7);
d2 = Parameters(8);

excursion1 = Parameters(9);
excursion2 = Parameters(10);
Q = Parameters(11);
No_of_Segments = Parameters(13);

%Initialize some variables
UG(1) = 0;
U = zeros(1,10);
Synth_Speech(1)= 0;
x1(1) = 0;
x2(1) = 0;
%Initialize difference variables
diff = zeros(1,No_of_Segments);

A = Parameters(12)*ones(No_of_Segments,1);  %Assume constant vocal tract area function


T = 1/100000;
Tinv = 1/T;

%Find static glottal areas
AG0_1=2*excursion1*Length;
AG0_2=2*excursion2*Length;

%Initialize window to display progress
wait_window = figure('Numbertitle','off',...
   'Color',[1 1 1],...
   'Name','Please wait',...
   'Position',[117 248 560 150]);
pause(0.01);

wait_counter = 0;

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

wait_window_display = uicontrol('Style','Text',...
   'Position',[50 50 410 30],...
   'Backgroundcolor','white',...
   'Foregroundcolor','blue',...
   'String','Please wait while glottal area waveforms are calculated');

drawnow;

%Set the values of the various constants in the model to the values suggested by
%Ishizaka and Flanagan. all stiffness values are in dyne/cm. 

%To account for fundamental frequency, scale the masses by Q

norm_m1 = m1/Q;
norm_m2 = m2/Q;

k1 = 80000;                    %Linear Stiffness of spring1(k1)
k2 = 2000;   						  %Linear Stiffness of spring2(k2)

norm_k1 = k1*Q;
norm_k2 = k2*Q;

kc = 25000;						  %Linear stiffness of coupled spring(kc)
norm_kc = kc*Q^2;

%A nonlinear characteristic is given to spring 1 and 2 for realistic simulation of
%vocal fold characteristics. Define these non linear stiffness constants. 

eta_k1 = 100;						%Same for spring1 and 2(etak1)
norm_eta_k1 = eta_k1*Q^2;

%Contact spring constants
h1 = 3*k1;         %Linear stiffness of contact spring(h1)
h2 = 24000;           
etah = 500;		    %Non linear stiffness - same for both contact springs(etah1)

%Define the various damping ratios

%Damping ratios for open phase (si1)
si1_open = 0.1;     
si2_open = 0.6;

%Closed phase damping ratios = saturation values of damping ratios
si1_close = 1.1;                              
si2_close = 1.6;


%Some more general constants
ro =  1.14e-03;          %Density of air
mu = 1.864e-04;          %Viscosity of air 
C = 3.53e04;             %Velocity of sound in air

%Call function to evaluate equivalent circuit parameter values
%This function uses the vocal tract area function stored in a file
%For more details about this file, refer to Wu's dissertation.
[R L C] = eq_circuit(No_of_Segments,A);


%NOTE :
%HERE THE ITERATION IS REPEATED A FIXED NUMBER OF TIMES ... IF NEED BE
%YOU CAN CHANGE IT TO BE REPEATED BY USING THE LENGTH OF UTTERANCE USING
%A DOUBLE LOOP, AND CHANGE THE PARAMTERS IN BETWEEN LOOPS.
%Number of times to repeat the iteration
%Divide the length of the utterance into intervals of 0.001 seconds each 
%The vocal tract parameters are assumed to be constant in this interval
%Note that for this program the pressure and initial glottal areas have been assumed 
%to be constant throughout the length of the utterance. If necessary, these parameters can
%be changed after each 0.001 sec interval....this is the reason for using a double loop 
%in the calculations

%interval = l_utter/0.001;
%interval = 250;
%Each 0.001 second interval corresponds to the sampling time. Therefore the number of
%sampling periods in each interval is given by
%samp_periods = 0.001/sampling_time;

%Iteration begins here
%for i = 1:interval
%If necessary change Ps, glottal areas etc. here
%for j = 1:samp_periods
%From the values of the constants, calculate the values for each component
%Based on Flanagan and Ishizaka
%Initial Calculation for the glottal volume velocity
%The notation used is the same as in Ishizaka and Flanagan 
%For more details refer to Ishizaka and Flanagan
for n = 1:3000
   %n = (i-1)*samp_periods + j;
   %Glottal areas as seen by the two masses
   garea1(n) = AG0_1 + 2*Length*x1(n);
   garea2(n) = AG0_2 + 2*Length*x2(n);
   
   R_k1 = (0.19*ro)/(garea1(n)^2);
    
   temp = garea2(n)/A(1);
   R_k2 = ro*(0.5 - temp*(1-temp))/garea1(n)^2;
   
   L_g1 = (ro*d1)/(garea1(n));
   L_g2 = (ro*d2)/(garea2(n));
   
   R_v1 = (12*mu*Length^2*d1)/(garea1(n)^3);
   R_v2 = (12*mu*Length^2*d2)/(garea2(n)^3);
   
   %Find the sum of differences in glottal areas. This quantity is used
   %in the equation to calculate glottal flow    
   
   for s = 1:No_of_Segments
      if s==1
         diff(s) = diff(s) + T/C(s)*(UG(n) - U(n,1));
      else
         diff(s) = diff(s) + T/C(s)*(U(n,s) - U(n,s-1));
      end
   end
   %The equation for the new value of UG can be written out as a quadratic equation 
   %of the form ax^2 + bx + c = 0.  Defining the values for a, b and c
   L_temp = (L_g1 + L_g2 + L(1))*Tinv;
   A_UG = R_k1 + R_k2;
   B_UG = L_temp + R(1)+R_v1+R_v2;
   C_UG =  -L_temp*UG(n)+(T/C(1)*diff(1))-Pressure;
   
   if C_UG >= 0
      UG(n+1) = 0;
   else
      UG(n+1) = (-B_UG + sqrt(B_UG^2 - 4*A_UG*C_UG))/(2*A_UG);
   end
   
   %For subsequent vocal tract segments we have a first order equation, 
   % Calculate a and b in the equation ax+b = 0
   %Calculate U for the next segments. This can be done in a loop
   %For ease the radiation resistance, inductance etc. are a part of the R and L vectors
   %A summation that requrires addition of the differences between glottal flows of 
   % adjacent segments occurs in these equations. The sum is accumulated in the variable
   %diff2 (The reason for the name diff2 is that this is the sum of successive differences)
   
   for s = 1:No_of_Segments
      if s ~= No_of_Segments
         a = Tinv*(L(s) + L(s+1))+ R(s) +R(s+1);
         %Split b into constants b1,b2 and b3 and add them for simplicity
         b = -Tinv*(L(s) + L(s+1))*U(n,s)+ diff(s)- diff(s+1);
      else
         %Calculate U at lips
         if n <= 1
            temp = 0;
         else
            temp = Synth_Speech(n-1);
         end
         a = Tinv*(L(s)+ L(s+1))+ R(s);
         b1 = -Tinv*(L(s) + L(s+1))*U(n,s);
         b2 = -Tinv*L(s+1)*(Synth_Speech(n) - temp);
         b = b1+b2+diff(s);
      end
      U(n+1,s) = -b/a;
   end
   
   %Now calculate radiation value of  U radiated(Synth_Speech)
   %Use the variable temp = Synth_Speech(n-1) from previous loop
   %This is done for computational efficiency
   a = (L(s+1)*Tinv)+ R(s+1);
   b = L(s+1)*Tinv*(-U(n+1,s)-Synth_Speech(n)+U(n,s));
   Synth_Speech(n+1) = -b/a;
      
   %Next step is to calculate the displacements
   %Define constants XMIN1 and XMIN2 that will be used often
   %A first order equation of the form ax+b=c is obtained
   %Determine the values of a and b
   XMIN1 = -(AG0_1/(2*Length));
   XMIN2 = -(AG0_2/(2*Length));
   %Calculate equivalent viscous resistances
   if x1(n) <= XMIN1
      r1 = 2*si1_close *sqrt(m1*k1);
   else
      r1 = 2*si1_open*sqrt(m1*k1);
   end
   if x2(n) <= XMIN2
      r2 = 2*si2_close * sqrt(m2*k2);
   else
      r2 = 2*si2_open*sqrt(m2*k2);
   end
   %Equations to calculate the displacements
   if x1(n)>XMIN1 
      Ax1 = norm_m1/T^2 + (r1/T) + norm_k1;
      if n <= 1
         temp = 0;
      else
         temp = x1(n-1);
      end
      Bx1 = (norm_m1/T^2)*(-2*x1(n)+temp)+ (r1/T)*(-x1(n))+norm_kc*(x1(n)- x2(n))...
         + norm_k1*eta_k1*(x1(n)^3);
      
   else
      Ax1 = norm_m1/T^2 + (r1/T) + norm_k1 + h1;
      if n <= 1
         temp = 0;
      else
         temp = x1(n-1);
      end
      b_temp  = (norm_m1/T^2)*(-2*x1(n)+temp)+(r1/T)*(-x1(n))+ norm_kc*(x1(n)- x2(n))...
         +norm_k1*eta_k1*x1(n)^3;
      Bx1  = h1*(AG0_1/(2*Length) + (etah*(x1(n)+ (AG0_1/(2*Length))^3)))+ b_temp;
   end
   if x2(n)>XMIN2
      Ax2 = norm_m2/T^2 + (r2/T) + norm_k2;
      if n <= 1
         temp = 0;
      else
         temp = x2(n-1);
      end
      Bx2 = (norm_m2/T^2)*(-2*x2(n)+temp)+ (r2/T)*(-x2(n))+ norm_kc*(x2(n)- x2(n))...
         + norm_k2*eta_k1*x2(n)^3;
   else
      Ax2  = norm_m2/T^2 + (r2/T) + norm_k2 + h2;
      if n <= 0
         temp = 0;
      else
         temp = x2(n-1);
      end
      b_temp2 = (norm_m2/T^2)*(-2*x2(n)+temp)+(r2/T)*(-x2(n))+norm_kc*(x2(n)- x2(n))...
         + norm_k2*eta_k1*x2(n)^3;
      Bx2 = h2*(AG0_2/(2*Length) + (etah*(x2(n)+ (AG0_2/(2*Length))^3)))+ b_temp2;
   end
   %Now use the forcing functions to calculate the  right hand side of the equation
   if (x1(n)>XMIN1 & x2(n)>XMIN2)
      PM1 = Pressure - 1.37*(ro/2)*((UG(n+1)/garea1(n))^2)...
         - (1/2)*(R_v1*UG(n+1)+ (L_g1/T)*(UG(n+1) - UG(n)));
      PM2 = PM1 -  (1/2)*(R_v1 + R_v2)*UG(n+1) ...
         + (L_g1+L_g2)*(UG(n+1) - UG(n))/(T*2) ...
         - (ro/2)*UG(n+1)^2*((1/garea2(n)^2) - (1/garea1(n)^2));
      Cx1 = PM1*d1*Length;
      Cx2 = PM1*d2*Length;
   elseif (x1(n) <= XMIN1 & x2(n)>XMIN2)
      Cx1 = Pressure*Length*d1;
      Cx2 = 0;
   elseif (x1(n) > XMIN1 & x2(n) <= XMIN2)
      Cx1 = Pressure*Length*d1;
      Cx2 = Pressure*Length*d2;
   else
      Cx1 = Pressure*(Length*d1);
      Cx2 = 0;
   end
   x1(n+1) = (Cx1 - Bx1)/Ax1;
   x2(n+1) = (Cx2 - Bx2)/Ax2;
   wait_counter = wait_counter + 1;
   set(wait_window_slider,'value',wait_counter);
   drawnow;
end
set(wait_window_display,'string','Thanks for your patience. Now plotting waveforms');
pause(1);

close(wait_window);
clear wait_counter wait_window wait_window_slider;