TubeModel53for412b.sci for windows scilab-4.1.2

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