SCILAB のバージョンがwindows用のScilab-4.1.2向けのTubeModel53for412ver2.sciを以下に示す。
（Windows Scilab-5.4.0 (32bit), Linux Scilab-5.5.2(64bit) でも動作させることができた。）
（使い方）
起動して exec ファイルとして読み込ませると
1 is /a/. 2 is /ae/. 3 is /o/. 4 is nasal/u/. 5 is /i/. 6 or 7 is manual-input.
Please Input number 1 or 2 or 3 or 4 or 5 or 6 or 7. And Enter !-->
が表示されるので、 6 を入れて Ｅｎｔｅｒキーを押す。
次に、２管声道モデルの各パラメーターを手動で入力していく。 各パラメーターの説明については ここを見てください。
+input l1 (example -0.0588235) -->
+input r1 (exampe 0.75) -->
+input ttl_Length[cm] (example 17) -->
+input ttl_Area[cm^2] (example 8) -->
何も入力しないで Enterキーを押すと exampleの値が設定される。
今回は、 l1 と r1 の値はexampleに固定にして、長さttl_Lengthと面積ttl_Areaを変化させて 正規化された周波数応答がどのように変化するかを計算する。
最後に、上部右端にある functions の中の
save_as_sound_file('y3out.wav')
を呼び出すと、合成した音を wavファイルの y3out_xmode*_Lng*_Ara*.wav として保存できて、 メディアプレーヤーなどで音を聞くことができる。
（警告）
このサンプルを実際に動作させるためには、
ご使用されるプログラミング環境に合わせて修正・変更・追加などが必要かもしれません。
また、バグが含まれている可能性があります。
万一、このサンプルを動作させる場合は、あなたの責任でおこなってください。
//------------------------------------------------------------------------------------------
//
// 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/. 6 or 7 is manual-input.'
xmode2=input('Please Input number 1 or 2 or 3 or 4 or 5 or 6 or 7. 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;
elseif (xmode2 == 6) | (xmode2 == 7) then // manual input is use Two Tubes Model
T_QTY=2;
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==6 then // manual input mode : example is
/a/
// set T_QTY=2 and use Two Tubes Model
l1in=input('+input l1 (example -0.0588235) ');
if l1in == [] then
l1in= -0.0588235;
disp(l1in);
end
r1in=input('+input r1 (exampe 0.75) ');
if r1in == [] then
r1in= 0.75;
disp(r1in);
end
ttl_Lengthx=input('+input ttl_Length[cm] (example 17) ');
if ttl_Lengthx == [] then
ttl_Lengthx= 17;
disp(ttl_Lengthx);
end
ttl_Areax=input('+input ttl_Area[cm^2] (example 8) ');
if ttl_Areax == [] then
ttl_Areax= 8;
disp(ttl_Areax);
end
L2x= ( 1. + l1in ) * ttl_Lengthx / 2.;
L1x= ttl_Lengthx - L2x;
A2x= ( 1. + r1in ) * ttl_Areax/ 2.;
A1x= ttl_Areax - A2x;
L1=[L1x, L1x, L1x, L1x, L1x, L1x, L1x ]; // set list of 1st tube's length by unit is [cm]
A1=[A1x, A1x, A1x, A1x, A1x, A1x, A1x ]; // set list of 1st tube's area by unit is [cm^2]
L2=[L2x, L2x, L2x, L2x, L2x, L2x, L2x ]; // set list of 2nd tube's length by unit is [cm]
A2=[A2x, A2x, A2x, A2x, A2x, A2x, A2x ]; // set list of 2nd tube's area by unit is [cm^2]
end
if xmode2==7 then // manual input mode example
is
/ae/
// set T_QTY=2 and use Two Tubes Model
l1in=input('+input l1 (example 0.53) ');
if l1in == [] then
l1in= 0.53;
disp(l1in);
end
r1in=input('+input r1 (exampe 0.78) ');
if r1in == [] then
r1in= 0.78;
disp(r1in);
end
ttl_Lengthx=input('+input ttl_Length[cm] (example 17) ');
if ttl_Lengthx == [] then
ttl_Lengthx= 17;
disp(ttl_Lengthx);
end
ttl_Areax=input('+input ttl_Area[cm^2] (example 9) ');
if ttl_Areax == [] then
ttl_Areax= 9;
disp(ttl_Areax);
end
L2x= ( 1. + l1in ) * ttl_Lengthx / 2.;
L1x= ttl_Lengthx - L2x;
A2x= ( 1. + r1in ) * ttl_Areax/ 2.;
A1x= ttl_Areax - A2x;
L1=[L1x, L1x, L1x, L1x, L1x, L1x, L1x ]; // set list of 1st tube's length by unit is [cm]
A1=[A1x, A1x, A1x, A1x, A1x, A1x, A1x ]; // set list of 1st tube's area by unit is [cm^2]
L2=[L2x, L2x, L2x, L2x, L2x, L2x, L2x ]; // set list of 2nd tube's length by unit is [cm]
A2=[A2x, A2x, A2x, A2x, A2x, A2x, A2x ]; // set list of 2nd tube's area by unit is [cm^2]
end
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]
//
frq_com=[100:25:Max_Freq];
frq_com_s=size(frq_com);
resp_com=zeros(1,frq_com_s(2));
resp_yg=ones(1,frq_com_s(2));
resp_hpf=ones(1,frq_com_s(2));
//
//==========================================================================
//
// 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,30,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
-ya3
-ya3
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
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
if xmode2 == 6 then
wstr1=sprintf('%d', int(L1(1)+L2(1))); // ttl_Length
wstr2=sprintf('%d', int(A1(1)+A2(1))); // ttl_Area
wstr0='y3out_xmode6_Lng' + wstr1 + '_Ara' + wstr2 + '.wav';
disp(wstr0);
wavwrite(y3out2, fs, 16, wstr0 );
elseif xmode2 == 7 then
wstr1=sprintf('%d', int(L1(1)+L2(1))); // ttl_Length
wstr2=sprintf('%d', int(A1(1)+A2(1))); // ttl_Area
wstr0='y3out_xmode7_Lng' + wstr1 + '_Ara' + wstr2 + '.wav';
disp(wstr0);
wavwrite(y3out2, fs, 16, wstr0 );
else
wavwrite(y3out2, fs, 16, filename );
end
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
//
//
// End of Step 4
//
//================================================================================
//
// Step 5 Plot amplitude[dB] vs normalized frequecny by 1st peak (f1)
//
//
//-------------------------------------------------
mess5='+++enter step 5 Plot';
mess5
// +yg response
function [res0]=yg_response()
stp=0;
stp=stp+N1(1);
y00=0.;
for v=1:(N2(1)+N3(1)+1)
y00= y00 + yg(stp + v);
end
for w=1:frq_com_s(2)
xw=2 * %pi * frq_com(w) * ts;
yi=0;
for v=1:(N2(1)+N3(1)+1)
yi= yi + yg(stp + v) * %e^(-1. * xw * (v-1) * %i);
end // for v
res0(w)=yi / y00; //resp_yg(w)=yi / y00;
end // for w
//
//[db,phi] =dbphi(res0);
//clf();
//bode(frq_com,db,phi);
endfunction
// -yg response
//-------------------------------------------------
// +hpf
function [res1]=hpf_response()
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);
// 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
res1=hpf1(2 * %pi * frq_com ); //resp_hpf=hpf1(2 * %pi * frq_com );
//[db,phi] =dbphi(resp_hpf);
//clf();
//bode(frq_com,db,phi);
endfunction
// -hpf
//----------------------------------------------------------
//----------------------------------------------------------
function [peak_db,peak_frq, peak_count0, avr_pwr, two_db_pwr, def_frq]=two_tubes_model_res2(ni, xr_flag, show_bode, nth)
for v=1:frq_com_s(2)
dummy0= two_tubes(2 * %pi * frq_com(v), ni, rg_closed ,rl);
resp_com(v)= resp_hpf(v) * resp_yg(v) * dummy0;
end
[db,phi] =dbphi(resp_com);
if show_bode==1 then
clf();
bode(frq_com,db,phi);
end
//
peak_count0=0;
peak_db=ones(1,frq_com_s(2));
peak_frq=ones(1,frq_com_s(2));
peak_index=ones(1,frq_com_s(2)); // added
for v=3:(frq_com_s(2)-2)
//+peak detect using 4points around the point
if db(v) > db(v-2) then
if db(v) > db(v-1) then
if db(v) > db(v+1) then
if db(v) > db(v+2) then
//+detect peak of which Q is more than q0, check only higher frequency
side
q0=2.;
dfrq0=frq_com(v) / q0;
frq_w1=frq_com(v)- ( dfrq0 / 2.0 ); // lower frequency side
frq_w2=frq_com(v)+ ( dfrq0 / 2.0 ); // higher frequency side
w2_v=-1.;
for v2=(v+1):frq_com_s(2)
if (db(v) - db(v2)) > 3.0 then // diffrence is
more than 3dB
w2_v=v2;
break;
end
if frq_com(v2) > frq_w2 then
w2_v=-3.;
break;
end
if db(v2) > db(v) then // not peak in Q area
w2_v=-2.;
break;
end
end //-v2
if w2_v > 0 then
//
peak_count0=peak_count0+1;
peak_db(peak_count0)=db(v);
peak_frq(peak_count0)=frq_com(v);
peak_index(peak_count0)=v; // added
//
end //if w2_v > 0 then
end
end
end
end
//-peak detect
end
//
//+++ added
if peak_count0 > 0 then
db_list=zeros(1,frq_com_s(2) - (peak_index(1)-1) );
phi_list=zeros(1,frq_com_s(2) - (peak_index(1)-1) );
frq_com_list=zeros(1,frq_com_s(2) - (peak_index(1)-1) );
for i=1:(frq_com_s(2) - (peak_index(1)-1))
db_list(i)=db( peak_index(1) +(i-1));
phi_list(i)=phi( peak_index(1) +(i-1));
frq_com_list(i)=frq_com(peak_index(1) +(i-1)) / frq_com(peak_index(1));
end
end
if show_bode==3 then
clf();
//subplot(311);
//plot_area(1);
subplot(211);
gainplot(frq_com,db',phi');
wstr1=sprintf('%f', L1(1)+L2(1)); // ttl_Length
wstr2=sprintf('%f', A1(1)+A2(1)); // ttl_Area
wstr3=sprintf('%f', r1(1)); // r1
wstr4=sprintf('%f', (L2(1)-L1(1)) / (L2(1)+L1(1)) ); //l1
wstr0='l1= ' + wstr4 + ' r1= ' + wstr3 + ' ttl_Length[cm]= ' + wstr1 + ' ttl_Area[cm^2]= ' + wstr2;
xtitle(wstr0,'Frequency[Hz]', 'Amplitude[dB]');
subplot(212);
if nth < 2 then // what until n-th plot ?
nth_new = ceil(frq_com_list(frq_com_s(2) - (peak_index(1)-1)) );
else
nth_new=nth;
end
plot2d(frq_com_list, db_list,style=[color("red")],rect=[1,min(db),nth_new,max(db)]);
xtitle('amplitude[dB] vs normalized frequecny by 1st peak (f1)','Normalized frequency', 'Amplitude[dB]');
end
//--- added
//
if xr_flag == 1 then
mess7z=' peak_count0 ';
messsp0=' ';
disp(peak_count0,mess7z);
for v=1:peak_count0
disp(peak_frq(v),messsp0, peak_db(v),messsp0, v);
end
end
// max two peak points under condition that they are neighborhood (next to each other)
pointn=1;
if peak_count0 > 1 then
v=pointn;
if peak_db(v) < 0. then
pk1= 1.0 / 10^(-1.0 * peak_db(v) / 20.0);
else
pk1= 10^( peak_db(v) / 20.0);
end
if peak_db(v+1) < 0. then
pk2= 1.0 / 10^(-1.0 * peak_db(v+1) / 20.0);
else
pk2= 10^( peak_db(v+1) / 20.0);
end
pk_ttl=pk1+pk2;
peak_db_ttl=20. * log10( pk_ttl);
//
if peak_count0 > 2 then
for v=2:(peak_count0-1)
if peak_db(v) < 0. then
pk1= 1.0 / 10^(-1.0 * peak_db(v) / 20.0);
else
pk1= 10^( peak_db(v) / 20.0);
end
if peak_db(v+1) < 0. then
pk2= 1.0 / 10^(-1.0 * peak_db(v+1) / 20.0);
else
pk2= 10^( peak_db(v+1) / 20.0);
end
pk_ttl2=pk1+pk2;
peak_db_ttl2=20. * log10( pk_ttl2);
if peak_db_ttl2 > peak_db_ttl then
peak_db_ttl=peak_db_ttl2;
pointn=v;
end
//
end
end //-if peak_count0 > 2 then
end //-if peak_count0 > 1 then
if xr_flag == 1 then
mess7z2=' peak_ttl of two points ';
disp(pointn,mess7z2);
disp(peak_db_ttl);
end
//+ average power
pk2=0.;
for v=1:frq_com_s(2)
if db(v) < 0. then
pk1= 1.0 / 10^(-1.0 * db(v) / 20.0);
else
pk1= 10^( db(v) / 20.0);
end
pk2=pk2+pk1;
end
pk2=pk2 / frq_com_s(2);
avr_pwr=10. * log10( pk2);
if xr_flag == 1 then
mess7z3=' average power ';
disp(avr_pwr,mess7z3);
end
//-------------------------------------------------------------
//+ two db power: strength scale is dB not linear
pointn2=1;
if peak_count0 > 1 then
two_db_pwr=( peak_db(1)+peak_db(2) ) /2.;
//
if peak_count0 > 2 then
for v=2:(peak_count0-1)
two_db_pwr2= ( peak_db(v)+peak_db(v+1)) / 2.;
if two_db_pwr2 > two_db_pwr then
two_db_pwr = two_db_pwr2;
pointn2=v;
end
//
end
end //-if peak_count0 > 2 then
end //-if peak_count0 > 1 then
//
//----------------------------------------------------------------------
//
// define def_frq= (w2-w0)/w0
//
//
def_frq= (peak_frq(pointn2+1)-peak_frq(pointn2))/peak_frq(pointn2);
if xr_flag == 1 then
mess7z4=' peak_db_pwr of two points, def_frq[rate] ';
disp(pointn2,mess7z4);
disp(def_frq , two_db_pwr);
end
//+ average power
pk2=0.;
for v=1:frq_com_s(2)
pk2=pk2+db(v);
end
pk2=pk2 / frq_com_s(2);
if xr_flag == 1 then
mess7z5=' average db power ';
disp(pk2,mess7z3);
end
//- two db power
//-----------------------------------------------------------------
endfunction
//
[resp_yg]=yg_response();
[resp_hpf]=hpf_response();
// +SWITCH
if (xmode2 == 6) | (xmode2==7) then // when manual input
xset('window',1); // create new windows
clf();
plot_area(1);
xset('window',2); // create new windows
clf();
[peak_db0,peak_frq0,peak_count0,avr_pwr0,two_db_pwr0,def_frq0]=two_tubes_model_res2(1,0,3,5);
xset('window',0); // back to old window
end
// -SWITCH
//
//================================================================================
//
// 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()'
; 'plot normailzed response'] );
functions=['two_tubes_model_bode(1)'; 'fft_plot(1024)';
'save_as_sound_file(''y3out.wav'')'; 'plot_area(1) '; 'plot_yg_yout()'
; ' two_tubes_model_res2(1,0,3,5)'] ;
end
end //if ADD_MENU == 1 then/
//
//===========================================================================================
//