This Tubes Model Waveform Generation is a sample program just for reference, of which name is TubeMode53for412c.sci. This is for windows scilab-4.1.2., but it maybe function on scilab-4.1.2 of some linux os.  Scilab is an open source platform for numerical computation.  If you calculate with certain versions of scilab,  this sample program will not work well , because some functions have changed  other than some older scilab  version. For example, if  scilab-3.1.1, wavwrite to make  .wav file  and xarrows to display tube length and area  does not work  right.  If it's scilab-3.1.1, wavwrite, real vector shall  be inversed row or column and  xarrows, arsize shall be adjusted.

//------------------------------------------------------------------------------------------
//
// Two Tubes Model for Vocal Tract and Three Tubes Model for Vocal Tract
// with Noise Source of which frequency band is limited
//
// .........................................................................................
// PLEASE PAY ATTENTION that this program may have some bugs and also you may modify program
// to fit your circumstances. If you use this program, everything done by your own risk.
// Thanks.
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
xmess='This program is for windows scilab-4.1.2 . If windows scilab-3.1.1, wavwrite and xarrows doesnot work well. Old scilab-2.7.2 plot doesnot work in this program. '
xmess1='1 is /a/.  2 is /ae/. 3 is /o/. 4 is nasal/u/. 5 is /i/. '
xmode2=input('Please Input number 1 or 2 or 3 or 4 or 5.  And  Enter !')
//
//
// Physical Parameters (length and area) Setup
//
D_QTY=7;            // set dimension of list. This defines length of generated wave
                    // Generated wave's length will be around D_QTY * (N1+N2+N3).
//-----------------------------
if xmode2 == 5 then
 T_QTY=5;
elseif xmode2 == 4 then
 T_QTY=4;
elseif xmode2 == 3 then
 T_QTY=3;
else
 T_QTY=2;
end
//----------------------------
//T_QTY=5;            //  T_QTY defines which Tubes Model  is used.
                    // set T_QTY 2 or other than below values when Two Tubes Model       
                    // set T_QTY 3 when Three Tubes Model of serial type
                    // set T_QTY 4 when Three Tubes Model of "T" type
                    // set T_QTY 5 when Two Tubes Model with "rl" noise

if xmode2==1 then                
// sample length & area for /a/ from problems 3.8 in Digital Processing of Speech Signals by L.R.Rabiner and R.W.Schafer
// set T_QTY=2 and use Two Tubes Model
L1=[9, 9, 9, 9, 9, 9, 9];     // set list of 1st tube's length by unit is [cm]
A1=[1, 1, 1, 1, 1, 1, 1];      // set list of 1st tube's area by unit is [cm^2]
L2=[8, 8, 8, 8, 8, 8, 8];     // set list of 2nd tube's length by unit is [cm]
A2=[7, 7, 7, 7, 7, 7, 7];      // set list of 2nd tube's area by unit is [cm^2]
end

if xmode2==2 then
// sample length & area for /ae/ from problems 3.8 in Digital Processing of Speech Signals by L.R.Rabiner and R.W.Schafer
// set T_QTY=2 Two Tubes Model
L1=[4, 4, 4, 4, 4, 4, 4];     // set list of 1st tube's length by unit is [cm]
A1=[1, 1, 1, 1, 1, 1, 1];      // set list of 1st tube's area by unit is [cm^2]
L2=[13, 13, 13, 13, 13, 13, 13];     // set list of 2nd tube's length by unit is [cm]
A2=[8, 8, 8, 8, 8, 8, 8];      // set list of 2nd tube's area by unit is [cm^2]
end

if xmode2==5 then
// sample length & area for /i/ from problems 3.8 in Digital Processing of Speech Signals by L.R.Rabiner and R.W.Schafer
// (((set T_QTY=2 and use Two Tubes Model)))
// set T_QTY=5 and use Two Tubes Model with "rl" noise
L1=[9, 9, 9, 9, 9, 9, 9];     // set list of 1st tube's length by unit is [cm]
A1=[8, 8, 8, 8, 8, 8, 8];      // set list of 1st tube's area by unit is [cm^2]
L2=[6, 6, 6, 6, 6, 6, 6];     // set list of 2nd tube's length by unit is [cm]
A2=[1, 1, 1, 1, 1, 1, 1];      // set list of 2nd tube's area by unit is [cm^2]
end

// sample for /u/
// set T_QTY=2 and use Two Tubes Model
//L1=[10, 10, 10, 10, 10, 10, 10];     // set list of 1st tube's length by unit is [cm]
//A1=[7, 7, 7, 7, 7, 7, 7];      // set list of 1st tube's area by unit is [cm^2]
//L2=[7, 7, 7, 7, 7, 7, 7];     // set list of 2nd tube's length by unit is [cm]
//A2=[3, 3, 3, 3, 3, 3, 3];      // set list of 2nd tube's area by unit is [cm^2]

if xmode2==3 then
// sample length & area for /o/ modoki
// set T_QTY=3 and use Three Tubes Model of serial type
L1=[3, 3, 3, 3, 3, 3, 3];     // set list of 1st tube's length by unit is [cm]
A1=[6, 6, 6, 6, 6, 6, 6];      // set list of 1st tube's area by unit is [cm^2]
L2=[7, 7, 7, 7, 7, 7, 7];     // set list of 2nd tube's length by unit is [cm]
A2=[1, 1, 1, 1, 1, 1, 1];      // set list of 2nd tube's area by unit is [cm^2]
L3=[8, 8, 8, 8, 8, 8, 8];     // set list of 3rd tube's length by unit is [cm]
A3=[7, 7, 7, 7, 7, 7, 7];      // set list of 3rd tube's area by unit is [cm^2]
end

if xmode2==4 then
// sample for  nasal /u/   This is not enough for nasal sound effect.
// set T_QTY=4 and use Three Tubes Model of "T" type
L1=[5, 5, 5, 5, 5, 5, 5];     // set list of 1st tube's length by unit is [cm]
A1=[4, 4, 4, 4, 4, 4, 4];      // set list of 1st tube's area by unit is [cm^2]
L2=[9, 9, 9, 9, 9, 9, 9];     // set list of 2nd tube's length by unit is [cm]
A2=[2, 2, 2, 2, 2, 2, 2];      // set list of 2nd tube's area by unit is [cm^2]
L3=[6, 6, 6, 6, 6, 6, 6];     // set list of 3rd tube's length by unit is [cm]
A3=[2, 2, 2, 2, 2, 2, 2];      // set list of 3rd tube's area by unit is [cm^2]
end


c0=35000;   // constant.  the speed of sound is round 35000 cm/second

if T_QTY == 3 then   //Three Tubes Model of serial type
 tu1=L1./c0;   // delay time in 1st tube
 tu2=L2./c0;   // delay time in 2nd tube
 tu3=L3./c0;   // delay time in 3rd tube
 r1=(A2-A1)./(A2+A1);  // reflection coefficient between 1st tube and 2nd tube
 r2=(A3-A2)./(A3+A2);  // reflection coefficient between 2nd tube and 3rd tube
elseif T_QTY == 4 then   //Three Tubes Model of T type
 tu1=L1./c0;   // delay time in 1st tube
 tu2=L2./c0;   // delay time in 2nd tube
 tu3=L3./c0;   // delay time in 3rd tube
 r12=((A2+A3)-A1)./(A3+A2+A1);  // reflection coefficient between 1st tube and others
 r21=(A2-(A1+A3))./(A3+A2+A1);  // reflection coefficient between 2nd tube and others
 r31=(A3-(A1+A2))./(A3+A2+A1);  // reflection coefficient between 3rd tube and others
else                //Two Tubes Model
 tu1=L1./c0;   // delay time in 1st tube
 tu2=L2./c0;   // delay time in 2nd tube
 r1=(A2-A1)./(A2+A1);  // reflection coefficient between 1st tube and 2nd tube
end


//---------------------------------------------------------------------------
// Input Glottal Volume Velocity for Two Tubes Model for Vocal Tract
// A.E.Rosenberg's formula as Glottal Volume Velocity
// Physical Parameters (duration) Setup
//
//D_QTY         // set dimension of list. This defines length of generated wave
tclosed=[5, 5, 5, 5, 5, 5, 5] ;      // set glottis closed duration (off time) by unit is millisecond [ms]
trise=  [6, 6, 6, 6, 6, 6, 6];     // set glottis opening duration (rise time) by unit is millisecond [ms]
tfall=  [2, 2, 2, 2, 2, 2, 2]; // set glottis closing duration (fall time) by unit is millisecond [ms]

// below (tfall is smaller) is  for /i/ that needs higher frequency signals as input
//tclosed=[4, 4, 4, 4, 4, 4, 4] ;      // set glottis closed duration (off time) by unit is millisecond [ms]
//trise=  [6, 6, 6, 6, 6, 6, 6];     // set glottis opening duration (rise time) by unit is millisecond [ms]
//tfall=  [0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7]; // set glottis closing duration (fall time) by unit is millisecond [ms]

amplitude=[1, 1, 1, 1, 1, 1, 1];   // set value of multiplier

// reflection coefficient between glottis and 1st tube
rg_rise=0.9;    //  set some value (1. or less) when glottis is opening
rg_fall=0.95;   //  set some value (1. or less) when glottis is closing
rg_closed=1;   // When glottis is closed,
               // reflection coefficient between glottis and 1st tube is 1
              //  because glottis's gate is closed and vocal tract is free from glottis's influence
//---------------------------------------------------------------------------
// Output Sound Pressure for Two Tubes Model for Vocal Tract
// Convert from Volume Velocity to Sound Pressure using High Pass Filter Model for Mouth Radiation
// Parameters ( cut off frequecny ) Setup
//
fc=1000;   // set cut off frequency of High Pass Filter by unit is [Hz]
fnew=fc;

// reflection coefficient between 2nd tube and mouth
rl=0.9;   // set adjust reduction response. rl is variable, however use one constant in this.
//
//---------------------------------------------------------------------------
//  "rl" noise is,        When Volume Velocity of the flow is above certain value,
//  Noise occured by turblent flow at narrow space near outside (mouth)
//  and  as new sound source, it returns through Tubes Model.
//  On calculation, this noise will be added to at "rl" reflection inside.
//  In this, one easy model is used,
//  Noise on/off is decided value of smoothing Volume Velocity by Low Pass Filter.
//  Addional Noise level is constant, beside actual depeds on value of Volume Velocity.
//
//  Property of noise  Setup
threshold_occur_minus = -1.1; // set threshold of Volume Velocity of minus side when noise occurs
threshold_occur_plus = 1.1; // set threshold of Volume Velocity of plus side when noise occurs
                            // two side threshold, minus side and plus side, is kunikuno-saku.
mix_amplitude=0.05;   // set mix amplitude at "rl" point
frq_center= 3000;  // set center frequency of noise by unit is [Hz]
freq_band= 2000;   // set frequency band of noise by unit is [Hz]
nstep_fir=129;    // set step number of band pass filter of linear-phase FIR filter
                  //   to transfer random noise to noise with limited frequency band

// Low Pass Filter to smooth Volume Velocity
fsm=1000;       // set cut off frequency of Low Pass Filter by unit is [Hz]
//
//
// reflection coefficient of 3rd tube's terminal side. 3rd tube is supposed that one side is closed.
//  This parameter is used only when T_QTY is set to 4 as Three Tubes Model of "T" type
rl3=-0.97;   // set certain a value, beside ideal value is -1


//---------------------------------------------------------------------------
// Other preparation: Sampling Frequency Setup
//
fs=44100;   // set sampling frequency by unit is [Hz]
ts=1/fs;
//
Max_Freq=7500;   // set plot maximum frequency by unit is [Hz]
//
//==========================================================================
//
// Step 1   Input Glottal Volume Velocity Generation
//
mess1='+++enter step 1  Input Glottal Volume Velocity Generation';
mess1
// length0 is wave length
length0=0; 
N1= zeros(1, D_QTY);
N2= zeros(1, D_QTY);
N3= zeros(1, D_QTY);
for loop=1:D_QTY
 N1(loop)= int( tclosed(loop) / 1000 / ts );
 N2(loop)= int( trise(loop) / 1000 / ts);
 N3(loop)= int( tfall(loop) / 1000 / ts);
 LL(loop)= N1(loop)+N2(loop)+N3(loop);
 length0 = length0 + LL(loop);
end
// yg is Glottal Volume Velocity.  rg is reflection coefficient between glottis and 1st tube
 yg= zeros(1,length0+1);
 rg= zeros(1,length0+1);
 tc0=1;
 yg(tc0)=0.;
 yg(tc0)=amplitude(1) * yg(tc0);
 for loop=1:D_QTY
 for t0=1:int(LL(loop))
 tc0=tc0+1;
 if t0 < N1(loop) then
   yg(tc0) =0;
   rg(tc0)=rg_closed;
 elseif t0 <= (N2(loop) + N1(loop)) then
   yg(tc0)= 0.5 * ( 1 - cos( ( %pi / N2(loop) ) * (t0 - N1(loop)) ) );
   rg(tc0)=rg_rise;
 elseif t0 <= (N3(loop)+N2(loop)+N1(loop)) then
   yg(tc0)= cos( ( %pi / ( 2 * N3(loop) )) * ( t0 - (N2(loop)+N1(loop)) ) );
   rg(tc0)=rg_fall;
 else
  yg(tc0) =0;
  rg(tc0)=rg_closed;
 end
 yg(tc0)=amplitude(loop) * yg(tc0);
end  // t0
end  // loop
// function save yg (Glottal Volume Velocity) as one sound wave file
function save_as_sound_file_yg( filename)
max01=0.;
for loop=1:length0
 if abs( yg(loop) ) > max01 then
    max01 = abs( yg(loop));
 end
end
 y3out2=zeros(1,length0+1);
if max01 > 0. then
 for loop=1:length0
  y3out2(loop)= yg(loop) / max01;
 end
end
 wavwrite(y3out2, fs, 16, filename );
endfunction
// Function for yg (Glottal Volume Velocity) frequency-phase response
function yg_bode(ni)
frq=[100:25:Max_Freq];
frqs=size(frq);
resp=zeros(1,frqs(2));
stp=0;
if ni > 1 then
for v=1:(ni-1)
 stp=stp + LL(v);
end
end
stp=stp+N1(ni);
y00=0.;
for v=1:(N2(ni)+N3(ni)+1)
 y00= y00 + yg(stp + v);
end
for w=1:frqs(2)
xw=2 * %pi * frq(w) * ts;
yi=0;
for v=1:(N2(ni)+N3(ni)+1)
  yi= yi + yg(stp + v) * %e^(-1. * xw * (v-1) * %i);
end  // for v
resp(w)=yi / y00;
end  // for w

[db,phi] =dbphi(resp);
 clf();
 bode(frq,db,phi);
endfunction
//
// End of Step 1   Input Glottal Volume Velocity Generation
//
//==========================================================================
//
// Step 2   Three Tubes/Two Tubes Model for Vocal Tract Calculation
//
if T_QTY == 3 then   //Three Tubes Model of serial type
 mess2='+++enter step 2  Three Tubes Model of serial type for Vocal Tract Calculation';
 mess2
// yt is output of Two Tubes Model for Vocal Tract
M1= zeros(1, D_QTY);
M2= zeros(1, D_QTY);
M3= zeros(1, D_QTY);
y2tm= zeros(1,length0+1);
for loop=1:D_QTY
M1(loop)= int( tu1(loop) / ts ) + 1;
M2(loop)= int( tu2(loop) / ts ) + 1;
M3(loop)= int( tu3(loop) / ts ) + 1;
end
MAX_M1=M1(1);
MAX_M2=M2(1);
MAX_M3=M3(1);
for loop=1:D_QTY
 if M1(loop) > MAX_M1 then
   MAX_M1= M1(loop)
 end
 if M2(loop) > MAX_M2 then
   MAX_M2= M2(loop)
 end
 if M3(loop) > MAX_M3 then
   MAX_M3= M3(loop)
 end
end
ya1= zeros(1,MAX_M1);
ya2= zeros(1,MAX_M1);
yb1= zeros(1,MAX_M2);
yb2= zeros(1,MAX_M2);
yc1= zeros(1,MAX_M3);
yc2= zeros(1,MAX_M3);
 tc0=1;
 ya1(tc0)=0.;
 ya2(tc0)=0.;
 yb1(tc0)=0.;
 yb2(tc0)=0.;
 yc1(tc0)=0.;
 yc2(tc0)=0.;
 y2tm(tc0)=0.;

 for loop=1:D_QTY
 for t0=1:int(LL(loop))

 tc0=tc0+1;

 for jl=MAX_M1:-1:2
  ya1(jl)=ya1(jl-1);
  ya2(jl)=ya2(jl-1);
 end
 for jl=MAX_M2:-1:2
  yb1(jl)=yb1(jl-1);
  yb2(jl)=yb2(jl-1);
 end
 for jl=MAX_M3:-1:2
  yc1(jl)=yc1(jl-1);
  yc2(jl)=yc2(jl-1);
 end

 ya1(1)= ((1. + rg(tc0) ) / 2.) * yg(tc0) + rg(tc0) * ya2(M1(loop));
 ya2(1)= -1. * r1(loop) * ya1(M1(loop))  +  ( 1. - r1(loop) ) * yb2(M2(loop));
 yb1(1)= ( 1 + r1(loop) ) * ya1(M1(loop)) + r1(loop) * yb2(M2(loop));
 yb2(1)= -1. * r2(loop) * yb1(M2(loop))  +  ( 1. - r2(loop) ) * yc2(M3(loop));
 yc1(1)= ( 1 + r2(loop) ) * yb1(M2(loop)) + r2(loop) * yc2(M3(loop));
 yc2(1)=  -1. * rl  * yc1(M3(loop));
 y2tm(tc0)= (1 + rl) * yc1(M3(loop));
end  // t0
end  //loop
elseif T_QTY == 4 then   //Three Tubes Model of T type
 mess2='+++enter step 2  Three Tubes Model of T type for Vocal Tract Calculation';
 mess2
// yt is output of Two Tubes Model for Vocal Tract
M1= zeros(1, D_QTY);
M2= zeros(1, D_QTY);
M3= zeros(1, D_QTY);
y2tm= zeros(1,length0+1);
for loop=1:D_QTY
M1(loop)= int( tu1(loop) / ts ) + 1;
M2(loop)= int( tu2(loop) / ts ) + 1;
M3(loop)= int( tu3(loop) / ts ) + 1;
end
MAX_M1=M1(1);
MAX_M2=M2(1);
MAX_M3=M3(1);
for loop=1:D_QTY
 if M1(loop) > MAX_M1 then
   MAX_M1= M1(loop)
 end
 if M2(loop) > MAX_M2 then
   MAX_M2= M2(loop)
 end
 if M3(loop) > MAX_M3 then
   MAX_M3= M3(loop)
 end
end
ya1= zeros(1,MAX_M1);
ya2= zeros(1,MAX_M1);
yb1= zeros(1,MAX_M2);
yb2= zeros(1,MAX_M2);
yc1= zeros(1,MAX_M3);
yc2= zeros(1,MAX_M3);
 tc0=1;
 ya1(tc0)=0.;
 ya2(tc0)=0.;
 yb1(tc0)=0.;
 yb2(tc0)=0.;
 yc1(tc0)=0.;
 yc2(tc0)=0.;
 y2tm(tc0)=0.;

 for loop=1:D_QTY
 for t0=1:int(LL(loop))

 tc0=tc0+1;

 for jl=MAX_M1:-1:2
  ya1(jl)=ya1(jl-1);
  ya2(jl)=ya2(jl-1);
 end
 for jl=MAX_M2:-1:2
  yb1(jl)=yb1(jl-1);
  yb2(jl)=yb2(jl-1);
 end
 for jl=MAX_M3:-1:2
  yc1(jl)=yc1(jl-1);
  yc2(jl)=yc2(jl-1);
 end
 ya1(1)= ((1. + rg(tc0) ) / 2.) * yg(tc0) + rg(tc0) * ya2(M1(loop));
 ya2(1)= -1. * r12(loop) * ya1(M1(loop))  +  ( 1. - r12(loop) ) * yb2(M2(loop)) + ( 1.  - r12(loop)) * yc2(M3(loop));
 yb1(1)= ( 1 + r21(loop) ) * ya1(M1(loop)) + r21(loop) * yb2(M2(loop)) + ( 1. + r21(loop)) * yc2(M3(loop));
 yb2(1)=  -1. * rl  * yb1(M2(loop));
 y2tm(tc0)= (1 + rl) * yb1(M2(loop));
 yc1(1)= ( 1. + r31(loop) ) * ya1(M1(loop)) + r31(loop) * yc2(M3(loop)) + ( 1. + r31(loop) ) * yb2(M2(loop)) ;
 yc2(1)=  -1. * rl3  * yc1(M3(loop));
end  // t0
end  //loop
else   //Two Tubes Model
 mess2='+++enter step 2  Two Tubes Model for Vocal Tract Calculation';
 mess2
// yt is output of Two Tubes Model for Vocal Tract
M1= zeros(1, D_QTY);
M2= zeros(1, D_QTY);
y2tm= zeros(1,length0+1);
for loop=1:D_QTY
M1(loop)= int( tu1(loop) / ts ) + 1;
M2(loop)= int( tu2(loop) / ts ) + 1;
end
MAX_M1=M1(1);
MAX_M2=M2(1);
for loop=1:D_QTY
 if M1(loop) > MAX_M1 then
   MAX_M1= M1(loop)
 end
 if M2(loop) > MAX_M2 then
   MAX_M2= M2(loop)
 end
end
ya1= zeros(1,MAX_M1);
ya2= zeros(1,MAX_M1);
yb1= zeros(1,MAX_M2);
yb2= zeros(1,MAX_M2);

 tc0=1;
 ya1(tc0)=0.;
 ya2(tc0)=0.;
 yb1(tc0)=0.;
 yb2(tc0)=0.;
 y2tm(tc0)=0.;

 for loop=1:D_QTY
 for t0=1:int(LL(loop))

 tc0=tc0+1;

 for jl=MAX_M1:-1:2
  ya1(jl)=ya1(jl-1);
  ya2(jl)=ya2(jl-1);
 end
 for jl=MAX_M2:-1:2
  yb1(jl)=yb1(jl-1);
  yb2(jl)=yb2(jl-1);
 end

 ya1(1)= ((1. + rg(tc0) ) / 2.) * yg(tc0) + rg(tc0) * ya2(M1(loop));
 ya2(1)= -1. * r1(loop) * ya1(M1(loop))  +  ( 1. - r1(loop) ) * yb2(M2(loop));
 yb1(1)= ( 1 + r1(loop) ) * ya1(M1(loop)) + r1(loop) * yb2(M2(loop));
 yb2(1)=  -1. * rl  * yb1(M2(loop));
 y2tm(tc0)= (1 + rl) * yb1(M2(loop));
end  // t0
end  //loop
end  //End of Two Tubes model
// Function for disply two tubes model's frequency-phase response when glottis is closed
function y = two_tubes( xw, ni, rg0, rl0 )
 yi  = 0.5 * ( 1 + rg0 ) * ( 1 + r1(ni))  * ( 1 + rl0 ) * %e^( -1 * ( tu1(ni) + tu2(ni) ) .* xw * %i);
 yb =  1 + r1(ni) * rg0 * %e^( -2 * tu1(ni) .* xw * %i) + r1(ni) * rl0 * %e^( -2 * tu2(ni) .* xw * %i) + rl0 * rg0 * %e^( -2 * (tu1(ni) + tu2(ni)) .* xw * %i);
 y = yi ./ yb;
endfunction
function two_tubes_model_bode(ni)
 frq=[100:25:Max_Freq]';
 [db,phi] =dbphi(two_tubes(2 * %pi * frq, ni,rg_closed ,rl));
 clf();
 bode(frq,db,phi);
endfunction
// Function for disply three tubes model's frequency-phase response when glottis is closed
function y = three_tubes( xw, ni, rg0, rl0 )
if T_QTY == 4 then  // "T" type
 yi  = 0.5 * ( 1. + rg0 ) * ( 1. + r21(ni))  * ( 1. + rl0 ) * %e^( -1. * ( tu1(ni) + tu2(ni)) * xw * %i);
 yb =  1.0;
 yb = yb  + r12(ni) * rg0 * %e^( -2 * tu1(ni) * xw * %i) + r21(ni) * rl0 * %e^( -2 * tu2(ni) * xw * %i);
 yb = yb  + rg0 * rl0 * ( 1. - r12(ni) + r21(ni) ) * %e^( -2 * (tu1(ni) + tu2(ni)) * xw * %i);
 yc = rl3 * ( 1. + r31(ni));
 yc = yc * (1. + rg0 * %e^( -2 * tu1(ni) * xw * %i));
 yc = yc * (1. - rl0 * %e^( -2 * tu2(ni) * xw * %i));
 yc = yc * %e^( -1 * tu3(ni) * xw * %i);
 yd = 1. - rl3 * %e^( -2 * tu3(ni) * xw * %i);
 if yd == 0. then    // There is room for reconsideration in this limit value.
  y=0.;   // because yc/yd is infinity when yd = 0
 else
  yc = yc / yd;
  y = yi / ( yb + yc);
 end
else   // serial type
 yi  = 0.5 * ( 1. + rg0 ) * ( 1. + r1(ni)) * ( 1. + r2(ni)) * ( 1. + rl0 ) * %e^( -1. * ( tu1(ni) + tu2(ni)  + tu3(ni)) * xw * %i);
 yb =  1 + rg0 * r1(ni) * %e^( -2 * tu1(ni) * xw * %i) + r1(ni) * r2(ni) * %e^( -2 * tu2(ni) * xw * %i) + r2(ni) * rl0 * %e^( -2 * tu3(ni) * xw * %i);
 yb = yb + rg0 * r2(ni) * %e^( -2 * (tu1(ni) + tu2(ni)) * xw * %i) + r1(ni) * rl0 * %e^( -2 * (tu2(ni) + tu3(ni)) * xw * %i);
 yb = yb + rg0 * r1(ni) * r2(ni) * rl0 *  %e^( -2 * (tu1(ni) + tu3(ni)) * xw * %i);
 yb = yb + rg0 * rl0 *  %e^( -2 * (tu1(ni) + tu2(ni) + tu3(ni)) * xw * %i);
 y = yi / yb;
end
endfunction
function three_tubes_model_bode(ni)
 frq3=[100:25:Max_Freq];
 frq3s=size(frq3);
 resp3=zeros(1,frq3s(2));
 for v=1:frq3s(2)
   resp3(v)= three_tubes(2 * %pi * frq3(v), ni,rg_closed ,rl);
 end
 [db,phi] =dbphi(resp3);
 clf();
 bode(frq3,db,phi);
endfunction
// Function plot sqrt(area) vs length
function plot_area(ni)
clf();
plot2d(0,0,1,"010","",[-5,-10,20,10],[5,5,4,5]);  // wake wo egaku
if T_QTY == 4 then  //Three Tubes Model of T type
ya1=sqrt(A1(ni) );
ya2=sqrt(A2(ni) );
ya3=sqrt(A3(ni) );
if ya2 > ya1 then
 yy0=ya2;
else
 yy0=ya1;
end
yy0=10- (yy0 * 2);
// Area scale is double than other Tube Models.
yy1=ya1  + yy0;
yy2=ya2  + yy0;
yy3=yy0 - L3(ni);
xx1=L1(ni) - (ya3 / 2);
xx2=L1(ni) + (ya3 / 2);
xp=[0    0    xx1  xx1   xx2   xx2  L1(ni)+L2(ni) L1(ni)+L2(ni)  L1(ni)  L1(ni)   0 ]';
yp=[yy1 yy0  yy0   yy3   yy3   yy0   yy0          yy2            yy2     yy1     yy1]';
xpolys( xp, yp, [5 5 5 5 5 5 5 5 5 5 5] );
xar=[-2 -0.5];
yar=[yy0+ ( ya1 / 2 ) yy0+ ( ya1 / 2 ) ];
xarrows(xar, yar, 20, 5);
yar=[yy0 + ( ya2 / 2 )   yy0+ ( ya2 / 2 ) ];
xar=[ L1(ni)+L2(ni)+0.5 L1(ni)+L2(ni)+2];
xarrows(xar, yar, 20, 5);
xstring(xx1, yy3 - 1, "CLOSED");
xstring(0,-9,"Attension: Area scale is double than others Tube Models");
//
elseif T_QTY == 3 then  //Three Tubes Model of serial type
ya1=sqrt(A1(ni) );
ya2=sqrt(A2(ni) );
ya3=sqrt(A3(ni) );
xp=[0    0    L1(ni)  L1(ni) L1(ni)+L2(ni) L1(ni)+L2(ni) L1(ni)+L2(ni)+L3(ni) L1(ni)+L2(ni)+L3(ni) L1(ni)+L2(ni) L1(ni)+L2(ni)  L1(ni)  L1(ni) 0   ]';
yp=[ya1 -ya1  -ya1    -ya2       -ya2      &nnbsp;   -ya3       -ya3      &nnbsp;         ya3                  ya3             ya2          ya2     ya1    ya1 ]';
xpolys( xp, yp, [5 5 5 5 5 5 5 5 5 5 5 5] );
xar=[-2 -0.5];
yar=[0 0];
xarrows(xar, yar, 20, 5);
xar=[ L1(ni)+L2(ni)+L3(ni)+0.5 L1(ni)+L2(ni)+L3(ni)+2];
xarrows(xar, yar, 20, 5);
else  //Two Tubes Model
ya1=sqrt(A1(ni) );
ya2=sqrt(A2(ni) );
xp=[0    0    L1(ni)  L1(ni)   L1(ni)+L2(ni)   L1(ni)+L2(ni)   L1(ni)  L1(ni) 0   ]';
yp=[ya1 -ya1  -ya1    -ya2       -ya2      &nnbsp;   ya2             ya2     ya1    ya1 ]';
xpolys( xp, yp, [5 5 5 5 5 5 5 5] );
xar=[-2 -0.5];
yar=[0 0];
xarrows(xar, yar, 20, 5);
xar=[ L1(ni)+L2(ni)+0.5 L1(ni)+L2(ni)+2];
xarrows(xar, yar, 20, 5);
 if T_QTY==5 then
   xar=[ L1(ni)+L2(ni)-0.5 L1(ni)+L2(ni)-2];
   xarrows(xar, yar, 20, 3);
 end
end
endfunction
//
// End of Step 2   Two Tubes Model for Vocal Tract Calculation
//
//==========================================================================
//
// Sub-Step 2-2    "rl" noise Calculation
//
if T_QTY == 5 then
mess22='+++enter sub-step 2-2 rl noise Calculation' ;
mess22
//  band pass filter of linear-phase FIR filter
//  0<cfreq(1),cfreq(2)<.5   sampling frequency = 1
cfreq(1) = (frq_center - ( freq_band / 2)) / fs;
cfreq(2) = (frq_center + ( freq_band / 2)) / fs;
fpar(1)=1;  // dummy data
fpar(2)=0;  // dummy data
[wft,wfm,fr]=wfir( 'bp' , nstep_fir ,cfreq, 'hm' ,fpar);
//
rand('normal');   // random generator is set to a Gaussian
length1= length0 + 1 + nstep_fir;
y_ran1=rand(1, length1);
y_ran2=zeros(1, (length0 + 1));
sc0=0;
mess221='---please wait for some time.  noise calculating' ;
mess221
sizex=1024;
for loop=1:(length0+1)
 if int(loop / sizex) > sc0 then
  sc0=int(loop / sizex);
  sc0 * sizex   // counter display
 end
 ydz = 0.;
 for loop2=1:nstep_fir
  ydz = ydz + wft(loop2) * y_ran1(1 + nstep_fir - loop2 + loop);
 end
 y_ran2(loop)=ydz;
end
//
mess23='+++enter sub-step 2-3 Re-Calculation of Two Tubes Model for Vocal Tract ' ;
mess23
//
// low pass filter's coefficients
wcsm=2. * %pi * fsm;
al1sm= -1. * %e^( -1. * wcsm * ts);
bl0sm= wcsm * ts;
bl1sm= 0.;
// IIR type low pass filter
 LI1sm=2;
 LI2sm=2;
 ya1sm= zeros(1,LI1sm);
 yb1sm= zeros(1,LI2sm);
//
// yt is output of Two Tubes Model for Vocal Tract
M1= zeros(1, D_QTY);
M2= zeros(1, D_QTY);
y2tm= zeros(1,length0+1);
y2tm_lpf= zeros(1,length0+1);
for loop=1:D_QTY
M1(loop)= int( tu1(loop) / ts ) + 1;
M2(loop)= int( tu2(loop) / ts ) + 1;
end
MAX_M1=M1(1);
MAX_M2=M2(1);
for loop=1:D_QTY
 if M1(loop) > MAX_M1 then
   MAX_M1= M1(loop)
 end
 if M2(loop) > MAX_M2 then
   MAX_M2= M2(loop)
 end
end
ya1= zeros(1,MAX_M1);
ya2= zeros(1,MAX_M1);
yb1= zeros(1,MAX_M2);
yb2= zeros(1,MAX_M2);

y_ran3=zeros(1, (length0 + 1));

 tc0=1;
 ya1(tc0)=0.;
 ya2(tc0)=0.;
 yb1(tc0)=0.;
 yb2(tc0)=0.;
 y2tm(tc0)=0.;
 y2tm_lpf(tc0)=0.;

 for loop=1:D_QTY
 for t0=1:int(LL(loop))

 tc0=tc0+1;

 for jl=MAX_M1:-1:2
  ya1(jl)=ya1(jl-1);
  ya2(jl)=ya2(jl-1);
 end
 for jl=MAX_M2:-1:2
  yb1(jl)=yb1(jl-1);
  yb2(jl)=yb2(jl-1);
 end

 ya1(1)= ((1. + rg(tc0) ) / 2.) * yg(tc0) + rg(tc0) * ya2(M1(loop));
 ya2(1)= -1. * r1(loop) * ya1(M1(loop))  +  ( 1. - r1(loop) ) * yb2(M2(loop));
 yb1(1)= ( 1 + r1(loop) ) * ya1(M1(loop)) + r1(loop) * yb2(M2(loop));
 //
 y2tm(tc0)= (1 + rl) * yb1(M2(loop));

 // smoothing y2tm
 // low pass filter
 ya1sm(LI1sm)=ya1sm(1);
 yb1sm(LI2sm)=yb1sm(1);
 ya1sm(1)=y2tm(tc0);
 yb1sm(1)= bl0sm * ya1sm(1) + bl1sm * ya1sm(LI1sm) - al1sm * yb1sm(LI2sm);
 y2tm_lpf(tc0)=yb1sm(1);
 if y2tm_lpf(tc0) > threshold_occur_plus then
  ynoise = y_ran2(tc0);
 elseif y2tm_lpf(tc0) < threshold_occur_minus then
  ynoise = y_ran2(tc0);
 else
  ynoise = 0.;
 end

 y_ran3(tc0)=ynoise;
 yb2(1)=  -1. * rl  * yb1(M2(loop)) + mix_amplitude * ynoise;  // mix with noise

end  // t0
end  //loop

//
end  // if T_QTY == 5 then
// Function for disply two tubes model with "rl" noise frequency-phase response when glottis is closed
//
function y = two_tubes_rl_noise( xw, ni, rg0, rl0 )
 //yi  = 0.5 * ( 1. + rg0 ) * ( 1. + r1(ni)) * ( 1. + rl0 ) * %e^( -1. * ( tu1(ni) + tu2(ni)) * xw * %i);  // this is for ug
 yi  = r1(ni) * r1(ni) * rg0 *  %e^( -2. * tu1(ni) * xw * %i) + r1(ni) *  %e^( -2. * tu2(ni) * xw * %i) + (1. - r1(ni) * r1(ni)) * rg0 * %e^( -2 * (tu1(ni) + tu2(ni)) * xw * %i);
 yi = (1. + rl0 ) * yi;
 yb =  1 + r1(ni) * rg0 * %e^( -2 * tu1(ni) * xw * %i) + r1(ni) * rl0 * %e^( -2 * tu2(ni) * xw * %i) + rl0 * rg0 * %e^( -2 * (tu1(ni) + tu2(ni)) * xw * %i);
 y = yi / yb;
endfunction
function two_tubes_rl_noise_bode(ni)
 frq3=[100:25:Max_Freq];
 frq3s=size(frq3);
 resp3=zeros(1,frq3s(2));
 for v=1:frq3s(2)
   resp3(v)= two_tubes_rl_noise(2 * %pi * frq3(v), ni,rg_closed ,rl);
 end
 [db,phi] =dbphi(resp3);
 clf();
 bode(frq3,db,phi);
endfunction
// Function FFT y_ran2(noise with limited frequency band) and plot frequency-phase
// input argument lng is fft points 1024 or others
// points 1024, when nzyou=10
// points 2048, when nzyou=11
// points 4096, when nzyou=12
function ntz=fft_plot_ran2( lng )
nzyou= int(log2( lng ));
ntz=2^nzyou;
// check if data length is enough for fft
ct0=0;
for loop=2:D_QTY
 ct0=ct0+LL(loop);
end
if ct0 > ntz then    // when actual data > ntz
xin1=zeros(1,ntz);
for loop=1:ntz
 xin1(loop) = y_ran2( loop);
end
win_hn=window('hn',ntz);   // make hanning windows data
xin2= xin1 .* win_hn;
fftout=fft(xin2,-1);
dfs= fs / ntz;
p100= int(100. / dfs);
p10000= int(Max_Freq / dfs);
frq=[(dfs * p100):dfs:(dfs * p10000)];
frqs=size(frq);
resp=zeros(1,frqs(2));
ct0=1;
for loop=p100:p10000
 resp(ct0)=fftout( loop + 1) / fftout(1);
 ct0=ct0+1;
end
[db,phi] =dbphi(resp);
clf();
bode(frq,db,phi);
end  // when actual data > ntz
endfunction
// Function for disply low pass filter's frequency-phase response
function y = lpf1( xw )
 yi= bl0sm  + bl1sm * cos( xw * ts) - bl1sm * sin( xw * ts ) * %i;
 yb=1.0  + al1sm * cos( xw * ts ) -  al1sm * sin( xw  * ts ) * %i;
 y=yi ./yb;
endfunction
function lpf_bode()
 frq=[100:25:Max_Freq]';
 [db,phi] =dbphi(lpf1(2 * %pi * frq ));
 clf();
 bode(frq,db,phi);
endfunction
//
// End of Sub-Step 2-2  "rl" noise Calculation
//
//==========================================================================
//
// Step 3   Output Sound Pressure for Two Tubes Model for Vocal Tract Calculation
//
mess3='+++enter step 3  Output Sound Pressure for Two Tubes Model for Vocal Tract Calculation';
mess3
// At first, low pass filter's coefficients
wc=2. * %pi * fc;
al1= -1. * %e^( -1. * wc * ts);
bl0= wc * ts;
// Then, convert them to high pass filter's coefficients
function y=f_alfa( wc_org,  wc_new)
 y= -1. *  cos( (wc_org + wc_new) * ts /2. ) / cos( (wc_org - wc_new) * ts / 2. );
endfunction
wnew=2. * %pi * fnew;
alfa=f_alfa( wc, wnew );
// high pass filter's coefficients
 bh0= bl0 / (1.0 -  al1 * alfa);
 bh1= alfa * bl0 / (1.0 -  al1 * alfa);
 ah1= (alfa - al1) / (1.0 -  al1 * alfa);
// IIR type high pass filter
y3out= zeros(1,length0+1);
 LI1=2;
 LI2=2;
 ya1= zeros(1,LI1);
 yb1= zeros(1,LI2);
//loop
 for loop=1:length0
  ya1(LI1)=ya1(1);
  yb1(LI2)=yb1(1);

 ya1(1)=y2tm(loop);
 yb1(1)= bh0 * ya1(1) + bh1 * ya1(LI1) - ah1 * yb1(LI2);
 y3out(loop)=yb1(1);

end   //loop
// Function for disply high pass filter's frequency-phase response
function y = hpf1( xw )
 yi= bh0  + bh1 * cos( xw * ts) - bh1 * sin( xw * ts ) * %i;
 yb=1.0  + ah1 * cos( xw * ts ) -  ah1 * sin( xw  * ts ) * %i;
 y=yi ./yb;
endfunction
function hpf_bode()
 frq=[100:25:Max_Freq]';
 [db,phi] =dbphi(hpf1(2 * %pi * frq ));
 clf();
 bode(frq,db,phi);
endfunction
// function save y3out (Sound Pressure) as one sound wave file
function save_as_sound_file( filename)
max01=0.;
for loop=1:length0
 if abs( y3out(loop) ) > max01 then
    max01 = abs( y3out(loop));
 end
end
 y3out2=zeros(1,length0+1);
if max01 > 0. then
 for loop=1:length0
  y3out2(loop)= y3out(loop) / max01;
 end
end
 wavwrite(y3out2, fs, 16, filename );
endfunction
//
//
// End of Step 3   Output Sound Pressure for Two Tubes Model for Vocal Tract
//
//================================================================================
//
// Step 4  Plot Glottal Volume Velocity and Sound Pressure
//
mess4='+++enter step 4  Plot generated sound wave';
mess4
function plot_yg_yout()
clf();
plot(yg,'b');
plot(y3out,'r');
xtitle('generated wave','Samples (Time)', 'Amplitude');
legend('Glottal Volume Velocity','Sound Pressure',2);
endfunction
function plot_yg_y2tm_y2tm_lpf_yo()
clf();
plot(yg,'b');
plot(y2tm,'y');
plot(y2tm_lpf,'g');
plot(y3out,'r');
//plot(y2tm,'r');
//plot(y_ran3,'g');  // noise plot
xtitle('generated wave','Samples (Time)', 'Amplitude');
legend('Glottal Volume Velocity','2nd Tube Edge Volume Velocity','Smoothed Volume Velocity','Sound Pressure',2);
endfunction

plot_yg_yout();

// Function FFT y3out(Sound Pressure) and plot frequency-phase
// FFT starts from 2nd section
// input argument lng is fft points 1024 or others
// points 1024, when nzyou=10
// points 2048, when nzyou=11
// points 4096, when nzyou=12
function ntz=fft_plot( lng )
nzyou= int(log2( lng ));
ntz=2^nzyou;
// check if data length is enough for fft
ct0=0;
for loop=2:D_QTY
 ct0=ct0+LL(loop);
end
if ct0 > ntz then    // when actual data > ntz
xin1=zeros(1,ntz);
ofs=LL(1);  // FFT starts from 2nd section
for loop=1:ntz
 xin1(loop) = y3out( loop + ofs);
end

win_hn=window('hn',ntz);   // make hanning windows data
xin2= xin1 .* win_hn;
fftout=fft(xin2,-1);
dfs= fs / ntz;
p100= int(100. / dfs);
p10000= int(Max_Freq / dfs);
frq=[(dfs * p100):dfs:(dfs * p10000)];
frqs=size(frq);
resp=zeros(1,frqs(2));
ct0=1;
for loop=p100:p10000
 resp(ct0)=fftout( loop + 1) / fftout(1);
 ct0=ct0+1;
end
[db,phi] =dbphi(resp);
clf();
bode(frq,db,phi);
end  // when actual data > ntz
endfunction
//
//================================================================================
//
// other functions
//
//yg_bode(1)    // plot yg (Glottal Volume Velocity) frequency-phase response
//two_tubes_model_bode(1)  // plot Two Tubes Model's frequency-phase response
//three_tubes_model_bode(1)  // plot Three Tubes Model's frequency-phase response
//hpf_bode()     // plot high pass filter's frequency-phase response
//save_as_sound_file_yg( 'yg.wav' )  // save yg (Glottal Volume Velocity) as one .wav file
//save_as_sound_file( 'y3out.wav' )  // save y3out (Sound Pressure) as one .wav file
//plot_area(1)  // plot sqrt(area) vs length
//plot_yg_yout() // plot yg (Glottal Volume Velocity) and y3out (Sound Pressure)
//fft_plot( 1024 )  // FFT y3out(Sound Pressure)  (1024 points) and plot frequency-phase
//
//two_tubes_rl_noise_bode(1) // plot Two Tubes Model with "rl" noise frequency-phase response
//fft_plot_ran2( 1024 )  //  FFT y_ran2(noise with limited frequency band) and plot frequency-phase
//plot_yg_y2tm_y2tm_lpf_yo() // plot yg (Glottal Volume Velocity) y2tm(Volume Velocity) y2tm_lpf(smoothed Volume Velocity) and y3out (Sound Pressure)
//lpf_bode()     // plot low pass (smoothing) filter's frequency-phase response
//
//
//===========================================================================================
//
//  Add menu button of which name is 'functions'
//
ADD_MENU=1;       // set ADD_MENU 1 to add menu button of some functions
//
//
if ADD_MENU == 1 then
delmenu('functions');
if T_QTY == 4 then
addmenu('functions',['three_tubes_model_bode(1)'; 'fft_plot(1024)'; 'save_as_sound_file(''y3out.wav'')'; 'plot_area(1) '; 'plot_yg_yout()'] );
functions=['three_tubes_model_bode(1)'; 'fft_plot(1024)'; 'save_as_sound_file(''y3out.wav'')'; 'plot_area(1) '; 'plot_yg_yout()'] ;
elseif T_QTY == 3 then
addmenu('functions',['three_tubes_model_bode(1)'; 'fft_plot(1024)'; 'save_as_sound_file(''y3out.wav'')'; 'plot_area(1) '; 'plot_yg_yout()'] );
functions=['three_tubes_model_bode(1)'; 'fft_plot(1024)'; 'save_as_sound_file(''y3out.wav'')'; 'plot_area(1) '; 'plot_yg_yout()'] ;
elseif T_QTY == 5 then
addmenu('functions',['two_tubes_model_bode(1)'; 'fft_plot(1024)'; 'save_as_sound_file(''y3out.wav'')'; 'plot_area(1) '; 'plot_yg_yout()'; 'two_tubes_rl_noise_bode(1)'; 'fft_plot_ran2(1024)'; 'plot_yg_y2tm_y2tm_lpf_yo()' ]);
functions=['two_tubes_model_bode(1)';  'fft_plot(1024)'; 'save_as_sound_file(''y3out.wav'')'; 'plot_area(1) '; 'plot_yg_yout()' ; 'two_tubes_rl_noise_bode(1)'; 'fft_plot_ran2(1024)'; 'plot_yg_y2tm_y2tm_lpf_yo()'] ;
else
addmenu('functions',['two_tubes_model_bode(1)'; 'fft_plot(1024)'; 'save_as_sound_file(''y3out.wav'')'; 'plot_area(1) '; 'plot_yg_yout()'] );
functions=['two_tubes_model_bode(1)';  'fft_plot(1024)'; 'save_as_sound_file(''y3out.wav'')'; 'plot_area(1) '; 'plot_yg_yout()'] ;
end
end //if ADD_MENU == 1 then/
//
//===========================================================================================
//
// date: 7th December 2007


No.2, 12  December 2007