次のページ     目次



<音声の波形の 生成を理解す る>

その4：　2管声道モデルを使って音をつなぐ
　「あえ」と「うあ」の生成と２つの母音の間のつながり方をさぐるこころみ


　音と音のつながりの部分はどの様につなげるのであろうか？ ２つの管（チューブ）をつなぎ合わせた模型を使って、もっとも簡単と思われる変形を管（チューブ）にすることで２つの音をつないでみた。２つの音からなる 生成した、「あえ」の音と「うあ」の .wavサンプルを参考にリンクしておこう。


右図は、この模型を使った「あ」から「え」の変化の例である。 管のつなぎ目の位置（赤い枠で示される管の断面積が変化している場所）を、真ん中から右の方向に移動することで、音を「あ」から「え」に変化させている。 周波数特性を見ると、周波数の低い側にあった揃った２つの波の一方がより高い周波数側に移動して、低い側は１つになり、より高い周波数側が２つの波になっ ているように移行していく。

もう一つ、右図に、「う」から「え」の例を示す。 「う」の音は汎用的でいろいろな形状がとりえるのであるが、ここでは、一本の管（チューブ）からはじめて、その真ん中を中心に外側（右の）の管をじょじょ に開いて行くことで、「う」から「あ」の音の方向に変化させている。周波数特性を見ると、（等しい間隔の）ばらばらの波が、近接していく２つの揃った波に 変化していくのがわかる。



また、フランス発の数値演算ライブラリ SCILAB　を使って波形を生成 するサンプルプログラムを作ってみたので紹介しよう。

ピッチ単位毎に管（チューブ）の変形を設定すことができる、声門波形と2管声道モデルによる音声波形の生成のサンプルプログラム ContTubeModel52.sciを以下に示す。
　
このサンプルでは、T_DEMOを１に設定して、半自動動作のデモ　モードで起動するようになっています。はじめに、「あえ」なら /ae/ を、 「うあ」なら /ua/ を選択します。

もしも、手動で操作するときの操作の流れとしては、

  1. メニュー　step1_edit　の中の
    　(1)set duration 　ピッチ単位で生成する数（長さ）を設定する。 管（チューブ）の個数は2に設定する。

  2.   (2)edit lrpa　　　管（チューブ）の断面積と長さを　r1 とl1を使って設定する。　

  3.   (3)make pareters　rとlから断面積と長さを計算する。

  4. メニュー　step2_wave_generate　の中の
    　(1)wave generate　音の生成の計算をする。

  5. メニュー　step3_check_area　の中の
    　(1)preparation bodebode　　周波数特性を計算する。

  6.   (2) set event handler window(1)  上図の様な、管（チューブ）の断面積と周波数特性の図を表示させる。

のようになる。

また、OSのサウンドの設定が上手くいっていて、scilabのバージョン(ここではwindowsのscilab-4.1.2を想定している)が上手く 合えば、 メニュー　menu_sound の中の　(2)play sound　 を使って生成した音を直接再生できる可能性もあります。または、 (3)save sound　で、　.wavの形式のファイルで一旦保存して、他のソフトウエア等を使って再生することを想定しています。

管（チューブ）の断面積と長さの設定は、ｒ1とｌ1の値を使って設定します。


　左から順番に、番号（この例では一番上の行が1番目を示している。次の行が2番目）、 ｌ１、ｒ１、そして、デフォルトである値に設定してあるピッチ間隔時間からの変 化をパーセントで設定（あまり大きな変化量は設定できない。符号はマイナスだと周波数が低くなる。）、振幅の大きさ（最大で１まで設定可能。波形の大きさ は正数で０から１の間で設定する）。
各番号が１ピッチ区間に相当していて、個々の要素を手動で設定した後、最後に、OK　ボタンを押す。
 
この設定画面の方法は、ちょっと使い勝手がよくないですが。 他に、編集した情報を、scilab の保存形式　または、テキスト形式のファイルで保存することもできます。



//------------------------------------------------------------------------------------------
//
// Two Tubes Model for Vocal Tract and Three Tubes Model for Vocal Tract
// with Noise Source of which frequency band is limited
// Continuous transformation of Tube, every pitch unit
// .........................................................................................
// 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.
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
T_DEMO=1;   // set T_DEMO=1 for demonstration, beside set T_DEMO=0 when normal
xmess='This program is only 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. '
//
//
//===================================================================
// scilab global variables
//===================================================================
// scilab global variables
global eh_var1
global eh_var2
global eh_var3
global eh_var4
global eh_var5
//
//==========================================================================
//
// Step -1
//
//------------------------------------------------------
// global variables
T_QTY=2;   // default set
D_QTY=7;  //  default set
memo1=''; // dummy set
memo2=''; // dummy set
// dummy set
A1=zeros(1,1);  // dummy set
A2=zeros(1,1);  // dummy set
A3=zeros(1,1);  // dummy set
L1=zeros(1,1);  // dummy set
L2=zeros(1,1);  // dummy set
L3=zeros(1,1);  // dummy set
tclosed=zeros(1,1);  // dummy set
trise=zeros(1,1);  // dummy set
tfall=zeros(1,1);  // dummy set
amplitude=zeros(1,1);  // dummy set
//--------------------------------------------------------------------------
if T_DEMO==1 then
  xmode2=x_choose(['/ae/';'/ua/'],['Please select one'],'default')
  if xmode2 <= 0 then
   sel_demo1=1;
  else
   sel_demo1=xmode2;
  end
  xmode2=0;
else
  xmode2=x_choose(['edit mode';'/a/';'/ae/';'/o/';'nasal /u/';'/i/'],['Please select one'],'default')
 if xmode2 <= 0 then
   xmode2=1;
 else
   xmode2=xmode2-1;   // for compatibility to 412b.sci
                      // xmode2==0, edit mode
 end       
end
//
//
// Physical Parameters (length and area) Setup
//
if xmode2 == 0 then
 //
else
 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).
end
//-----------------------------
if xmode2 == 5 then
 T_QTY=5;
elseif xmode2 == 4 then
 T_QTY=4;
elseif xmode2 == 3 then
 T_QTY=3;
elseif xmode2 == 1 then
 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

//----------------------------------------------------------------
function [tu1,tu2,tu3,r1,r2,r12,r21,r31]=make_tu_r()
 tu3=zeros(1,1);  // dummy set
 r1=zeros(1,1);   // dummy set
 r2=zeros(1,1);   // dummy set
 r12=zeros(1,1);   // dummy set
 r21=zeros(1,1);   // dummy set
 r31=zeros(1,1);   // dummy set
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
elseif T_QTY == 2 then                //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
elseif T_QTY == 5 then                //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
endfunction

//---------------------------------------------------------------------------
// Input Glottal Volume Velocity for Two Tubes Model for Vocal Tract
// A.E.Rosenberg's formula as Glottal Volume Velocity
// Physical Parameters (duration) Setup
//
 // default tclosed, trise, tfall
 tclosed0=5.;
 trise0=6.;
 tfall0=2.;
if xmode2 == 0 then
//
else
 //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
end

// 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;
//
Min_Freq=100;   // set plot maximum frequency by unit is [Hz]
Max_Freq=7500;   // set plot maximum frequency by unit is [Hz]
//
//==========================================================================
//
// Step 0   Edit
//
//------------------------------------------------------
// global variables
ttl_Length=17;  // set  total length of combined tubes, Two Tubes, Three Tubes, and etc               
ttl_Area2=10;    // set  total area of combined tubes, Two Tubes,  for dummy display
ttl_Area3=14;    // set  total area of combined tubes, Three Tubes,  for dummy display
lrpa=zeros(1,1); // rlap is matrix of n-th,l1,r1,pitch(+-%),amplitude(0-1),(noise on/off)
//
//
//-------------------------------------------------
function [D_QTY,T_QTY,lrpa,memo1x]= set_QTY()
 txt1=['duration, pitch cycle number'; 'number of Tubes (2 or 3)';'memo'];
  wstr1=sprintf('%d',D_QTY);
  wstr2=sprintf('%d',T_QTY);
  wstr3=memo1;
  sig1=x_mdialog('set duration and type of Tube Model',txt1, [wstr1; wstr2; wstr3]);
 if sig1==[] then
  arg1=evstr(wstr1);
  arg2=evstr(wstr2);
  memo1x=wstr3;
 else
  arg1=evstr(sig1(1));
  arg2=evstr(sig1(2));
  memo1x=sig1(3);
 end
 D_QTY=int(arg1);
 T_QTY=int(arg2);
 //
 if T_QTY == 2 then
  lrpa=zeros(D_QTY,5);
  for v=1:D_QTY
    lrpa(v,1)=v;
    lrpa(v,2)=0.;
    lrpa(v,3)=0.8;
    lrpa(v,4)=0.;
    lrpa(v,5)=1.;
  end
 elseif T_QTY == 3 then
  lrpa=zeros(D_QTY,7);
  for v=1:D_QTY
    lrpa(v,1)=v;
    lrpa(v,2)=0.4;
    lrpa(v,3)=-0.7;
    lrpa(v,4)=0.07;
    lrpa(v,5)=0.75;
    lrpa(v,6)=0.;
    lrpa(v,7)=1.;
  end
 end
endfunction
//--------------------------------------------------------------------------
function [A1,A2,L1,L2]=cal_AL2(r1,l1)  
// for Two Tubes Model
 r1s=size(r1);   
 size_r1=r1s(2);   // size of array r1
l1s=size(l1);
 size_l1=l1s(2);   // size of array l1
 //
 L1=zeros(1,size_l1);
 L2=zeros(1,size_l1);
 for v=1:size_l1
  L2(v)= ( 1. + l1(v) ) * ttl_Length / 2.;
  L1(v)= ttl_Length - L2(v);
 end
 A1=zeros(1,size_r1);
 A2=zeros(1,size_r1);
 for v=1:size_r1
  A2(v)= ( 1. + r1(v) ) * ttl_Area2/ 2.;
  A1(v)= ttl_Area2 - A2(v);
 end
endfunction
//--------------------------------------------------------------------------
function [A1,A2,A3,L1,L2,L3]=cal_AL3(r1,l1,r2,l2)  
// for Three Tubes Model
//         tube  #1  #2  #3
//              ___ ___ ___
//  glottis=>  |___|___|___| => out
//               L1  L2  L3
 r1s=size(r1);   
 size_r1=r1s(2);   // size of array r1
l1s=size(l1);
 size_l1=l1s(2);   // size of array l1
 //
 L1=zeros(1,size_l1);
 L2=zeros(1,size_l1);
 L3=zeros(1,size_l1);
 for v=1:size_l1
  y0=l1(v) * l2(v) + ( l1(v) - l2(v) ) + 3.;
  L1(v)= ( l2(v) - 1. ) * ( l1(v) - 1. ) / y0 * ttl_Length;
  L2(v)=  -1.0 * ( 1. + l1(v) ) * ( l2(v) - 1. ) / y0 * ttl_Length;
  L3(v)= ( l2(v) + 1. ) * ( l1(v) + 1. ) / y0 * ttl_Length;
 end
 A1=zeros(1,size_r1);
 A2=zeros(1,size_r1);
 A3=zeros(1,size_r1);
 for v=1:size_r1
  y0=r1(v) * r2(v) + ( r1(v) - r2(v) ) + 3.;
  A1(v)= ( r2(v) - 1. ) * ( r1(v) - 1. ) / y0 * ttl_Area3;
  A2(v)=  -1.0 * ( 1. + r1(v) ) * ( r2(v) - 1. ) / y0 * ttl_Area3;
  A3(v)= ( r2(v) + 1. ) * ( r1(v) + 1. ) / y0 * ttl_Area3;
 end
endfunction
//------------------------------------------------------------------
function [lrpar,memo2x]=edit_lrpa()
 //
 if T_QTY == 2 then
   [lrpar]=x_matrix('Edit no, l1, r1, pitch(+-%), amplitude(0-1)',lrpa);
   memo2x=' no, l1, r1, pitch(+-%), amplitude(0-1)';
 elseif T_QTY == 3 then
   [lrpar]=x_matrix('Edit no, l1, r1, l2, r2 ,pitch(+-%), amplitude(0-1)',lrpa);
   memo2x=' no, l1, r1, l2, r2 ,pitch(+-%), amplitude(0-1)';
 end
 //
 if lrpar==[] then   // set old setting
   s0=size(lrpa);
   lrpar=zeros(s0(1),s0(2));
   for v1=1:s0(1)
     for v2=1:s0(2)
       lrpar(v1,v2)=lrpa(v1,v2);
     end
   end
 else
 // 
 end
 //...
endfunction 
//--------------------------------------------------------------------------
function [A1,A2,A3,L1,L2,L3,tclosed,trise,tfall,amplitude]=make_paras()
 A1=zeros(1,1);  // dummy set
 A2=zeros(1,1);  // dummy set
 A3=zeros(1,1);  // dummy set
 L1=zeros(1,1);  // dummy set
 L2=zeros(1,1);  // dummy set
 L3=zeros(1,1);  // dummy set
 //
 r1=zeros(1,D_QTY);
 r2=zeros(1,D_QTY);
 l1=zeros(1,D_QTY);
 l2=zeros(1,D_QTY);
 //
 s0=size(lrpa);
 if T_QTY == 2 then
   vp=4;
   va=5;
   for v=1:D_QTY
     l1(v)=lrpa(v,2);
     r1(v)=lrpa(v,3);
   end
   [A1,A2,L1,L2]=cal_AL2(r1,l1);  

 elseif T_QTY == 3 then
   vp=6;
   va=7;
   for v=1:D_QTY
     l1(v)=lrpa(v,2);
     r1(v)=lrpa(v,3);
     l2(v)=lrpa(v,4);
     r2(v)=lrpa(v,5);
   end
   [A1,A2,A3,L1,L2,L3]=cal_AL3(r1,l1,r2,l2);  

 end
 
 // common
 tclosed=zeros(1,D_QTY);
 trise=zeros(1,D_QTY);
 tfall=zeros(1,D_QTY);
 amplitude=zeros(1,D_QTY);
 pl0=tclosed0+trise0+tfall0;
 for v=1:D_QTY
    amplitude(v)=lrpa(v,va);
    //
    adjp=lrpa(v,vp);
    delta= pl0 * adjp / 100.;
    if delta > 0 then   // plus, more high speed
       if delta < tclosed0 then 
         tclosed(v) = tclosed0 - delta;
        else
         tclosed(v) = 0.0;
         disp('waring: cannot pitch control because tclosed0 is not enough');
        end
    else   // when minus, more low speed
       tclosed(v)=tclosed0+delta;  
    end
    trise(v)=trise0;
    tfall(v)=tfall0;
 end
 //
 disp(' - makeing parameters was done. ');
 //
endfunction
//--------------------------------------------------------------------------
function save_lrpa()
  filo=input(' + enter file name for save =>',["string"]);
  save( filo,D_QTY,T_QTY,lrpa,memo1);
  disp(' - save was done. ');
endfunction
function [D_QTY,T_QTY,lrpa,memo1]=load_lrpa()
  filo=tk_getfile(Title="Choose load file");
  load( filo,'D_QTY','T_QTY','lrpa');
  memo1='';
  load( filo,'memo1');
  disp(' - load was done. ');
endfunction
//--------------------------------------------------------------------------
//
//
//==========================================================================
//
// Step 1   Input Glottal Volume Velocity Generation
//
//------------------------------------------------------
// global variables
rg=zeros(1,1); // dummy set
yg=zeros(1,1); // dummy set
N1=zeros(1,1); // dummy set
N2=zeros(1,1); // dummy set
N3=zeros(1,1); // dummy set
LL=zeros(1,1); // dummy set
length0=0;
//
function [rg,yg,N1,N2,N3,LL,length0]=step1()
disp('+++enter step 1  Input Glottal Volume Velocity Generation');
// length0 is wave length
length0=0; 
N1= zeros(1, D_QTY);
N2= zeros(1, D_QTY);
N3= zeros(1, D_QTY);
LL= 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
endfunction // step1
// 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
//
//------------------------------------------------------
// global variables
M1= zeros(1, 1);  // dummy set
M2= zeros(1, 1); // dummy set
M3= zeros(1, 1); // dummy set
y2tm= zeros(1,1); // dummy set
//
function [y2tm,M1,M2,M3]=step2()
if T_QTY == 3 then   //Three Tubes Model of serial type
 disp('+++enter step 2  Three Tubes Model of serial type for Vocal Tract Calculation');
// 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
 disp('+++enter step 2  Three Tubes Model of T type for Vocal Tract Calculation');
// 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
 disp('+++enter step 2  Two Tubes Model for Vocal Tract Calculation');
// 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
endfunction // step2()
//-------------------------------------------------------------------
// 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 y = two_tubes2( 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          -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
//
//------------------------------------------------------
// global variables
y_ran2=zeros(1,1);    // dummy set
y2tm_lpf= zeros(1,1);   // dummy set
al1sm= -1.; // dummy set
bl0sm= 1.; // dummy set
bl1sm= 0.; // dummy set
//
function [y2tm,y2tm_lpf,y_ran2,al1sm,bl0sm,bl1sm]=step2_2()
if T_QTY == 5 then
disp('+++enter sub-step 2-2 rl noise Calculation') ;
//  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;
disp('---please wait for some time.  noise calculating') ;
sizex=1024;
for loop=1:(length0+1)
 if int(loop / sizex) > sc0 then
  sc0=int(loop / sizex);
  disp( 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
//
disp('+++enter sub-step 2-3 Re-Calculation of Two Tubes Model for Vocal Tract ') ;
//
// 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
endfunction  // step2_2()

// 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
//
//------------------------------------------------------
// global variables
y3out= zeros(1,1);  // dummy set
ah1=1.; // dummy set
bh0=1.; // dummy set
bh1=0.; // dummy set
//
function [y3out,ah1,bh0,bh1]=step3()
disp('+++enter step 3  Output Sound Pressure for Two Tubes Model for Vocal Tract Calculation');
// 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 );
alfa= -1. *  cos( (wc + wnew) * ts /2. ) / cos( (wc - wnew) * ts / 2. );
// 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
endfunction  // step3()
//
// 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
//
function step4()
disp('+++enter step 4  Plot generated sound wave');
plot_yg_yout();
endfunction
//
function plot_yg_yout()
clf();
plot(yg,'b');
plot(y3out,'r');
wstr1='generated wave' + ' : ' + memo1;
xtitle(wstr1,'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
// 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
//-------------------------------------------------------------------------
// + others function library
//
//------------------------------------------------------
// global variables
sp1=1;       // start point for actual display of section 1
ep1=1;       // end point for actual display of section 1
fft_point1=512;   // default set
db_fft1=zeros(1,1);   // dummy set
phi_fft1=zeros(1,1);  // dummy set
db_all=zeros(1,1);   // dummy set
phi_all=zeros(1,1);   // dummy set
//----------------------------------------------------------------
function [fft_point1]=set_fft_points()
 l1=list('points',3,['128','256','512','1024']);
 wrep=x_choices('Select FFT points',list(l1));
 fft_point1=512; // default
 if wrep== 1 then
  fft_point1=128;
 elseif  wrep== 2 then
  fft_point1=256;
 elseif  wrep== 3 then
  fft_point1=512;
 elseif  wrep== 4 then
  fft_point1=1024;
 end
endfunction
//-------------------------------------------------
function [frq1,sfrq1,is1,ie1]=set_frq( spec1024 )
// spec1024,  tube model response no syuhasuu bunkainou wo ageru tame
fs1=fs;  // to keep old source compatibility
dfsw= fs1 / fft_point1;
is1= ceil(Min_Freq / dfsw);
ie1= int(Max_Freq / dfsw);
if spec1024 == 1024 then
 dfsw2= fs1 / 1024.0;
else
 dfsw2=dfsw;
end
frq1=[(dfsw * is1):dfsw2:(dfsw * ie1)];
frqs=size(frq1);
sfrq1=frqs(2);
endfunction
//-----------------------------------------------
function  [db_fft1,phi_fft1,rt_code]=do_fft_wav(wdat1)
 // size1=length0;  // to keep old source compatibility
 [size1]=rt_size(wdat1);
 rt_code=-1;
 if size1 >= ( sp1 + fft_point1 - 1 ) then
  win_hn=window('hn',fft_point1);  // make hanning windows data
  wdatw=zeros(1, fft_point1);
  for v=1:fft_point1
   wdatw(v)=wdat1(sp1 + v - 1);
  end
  wdatw2= wdatw .* win_hn;
  fftwout=fft( wdatw2, -1 );
 //
 [frq1,sfrq1,is1,ie1]=set_frq(0);
 //
 respw=zeros(1,sfrq1);
 ct0=1;
 for loop=is1:ie1
  respw(ct0)=fftwout(loop+1);
  ct0=ct0+1;
 end
 [db_fft1,phi_fft1] =dbphi(respw);
 rt_code=0;
 end   // -if size1 <= ( sp1 + fft_point1 - 1 ) then
endfunction
//-----------------------------------------------
function [wsp1, wep1]= set_sp1_ep1(wdat1)
 // wsize1=length0;  // to keep old source compatibility
 [wsize1]=rt_size(wdat1);
 txt1=['start point(sec.1)';'end point(sec.1)'];
 wstr1=sprintf('%d',sp1);
 if ep1 > wsize1 then
   ep2 = wsize1;
 else
   ep2 = ep1;
 end
 if ep2 <= sp1 then
   ep2 = wsize1;
 end
 wstr2=sprintf('%d',ep2);
 sig1=x_mdialog('Input start point and end point for portion display.',txt1, [wstr1 ; wstr2 ]);
 if sig1==[] then
  arg1=evstr(wstr1);
  arg2=evstr(wstr2);
 else
  arg1=evstr(sig1(1));
  arg2=evstr(sig1(2));
 end
 wsp1=sp1;
 wep1=ep2;
 if arg1 >= 1 then
  if arg2 <= wsize1 then
   if arg1 < arg2 then
     wsp1=arg1;
     wep1=arg2;
   end
  end
 end
endfunction
//--------------------------------------------------
function [length1]=rt_size(wdat1)
 s0=size(wdat1);
 if s0(1) > s0(2) then
   length1 = s0(1);
 else
   length1 = s0(2);
 end
endfunction
//--------------------------------------------------
function display_portion(disp0,wdat0)
 wb0=xget('window');  // stack old window
 [length1]=rt_size(wdat0);
 if disp0 >= 1 then
  xset('window',disp0);   // create new windows
  clf();
  //(1) plot wave all
  subplot(311);
  xziku0=zeros(1,length1);
  for v=1:length1
   xziku0(v)=v;
  end
  plot(xziku0,wdat0,'b');
  wstr1='waveform all : ' + memo1;
  xtitle( wstr1 );
  //(2) plot wave portion
  subplot(312);
  length12=ep1-sp1+1;
  xziku1=zeros(1,length12);
  wdat1=zeros(1,length12);
  for v=1:length12
   xziku1(v)=v+sp1-1;
   wdat1(v)=wdat0(v+sp1-1);
  end
  plot(xziku1,wdat1,'r');
  wstr1='waveform portion : ' ;
  xtitle( wstr1 );
  subplot(311);
  plot(xziku1,wdat1,'r');
  //(3) plot fft of the portion
  subplot(313);
  [db_fft1,phi_fft1,rt_code]=do_fft_wav(wdat0);
  if rt_code >= 0 then
   [frq1,sfrq1,is1,ie1]=set_frq(0);
   gainplot(frq1,db_fft1,phi_fft1);
   wstr1 = 'frequency response of waveform part of which  ';
   wstr1 = wstr1  + ' start=' + string(sp1) + ' points=' + string( fft_point1);
   xtitle( wstr1 );
  end
 end  // disp0 >= 1 then
 xset('window',wb0);   // push old windows
endfunction
//--------------------------------------------------
// + event handler
function my_eventhandler(win,x,y,ibut)
 global eh_var1
 global eh_var2
 global eh_var3
 global eh_var4
 global eh_var5
 global eh_fm_disp0
      if ibut==-1000 then return,end
      [x,y]=xchange(x,y,'i2f')
       xinfo(msprintf('Event code %d at mouse position is (%f,%f)',ibut,x,y))
      if ibut== 100   // d
        if (eh_var2 + eh_var3 - 1) < D_QTY then
         eh_var2=eh_var2+1
        end
        plot_areas(eh_var1, eh_var2, eh_var3)
      elseif ibut== 117  // u
        if eh_var2 > 1 then
          eh_var2=eh_var2-1
        end
        plot_areas(eh_var1, eh_var2, eh_var3)
      elseif ibut==3 // Left mouse button has been clicked
        eh_var4=x
      elseif ibut==5 // Right mouse button has been clicked
        eh_var5=x
      end
// seteventhandler('my_eventhandler')
// seteventhandler('') //suppress the event handler
endfunction
//--------------------------------------------------
function plot_areas(disp_code, nstartline, nlines)
 w0=xget('window');  // stack old window
 if disp_code >= 1 then
  xset('window',disp_code);   // create new windows
  clf();
  //
     frq=[Min_Freq:25:Max_Freq];
     frqs=size(frq);
     db_fft=zeros(1,frqs(2));
     phi_fft=zeros(1,frqs(2));
     for vloop=1:nlines
       subplot(nlines ,2, 2*(vloop-1)+1);
       plot_area2(vloop+nstartline-1,vloop);

      subplot(nlines ,2, 2*(vloop-1)+2);
      for w=1:frqs(2)
        db_fft(w)= db_all(w,vloop+nstartline-1);
        phi_fft(w)=phi_all(w,vloop+nstartline-1);
      end
      gainplot(frq,db_fft,phi_fft);
      wstr1 = 'frequency response theorycal: ';
      xtitle( wstr1 );

    end   // for vloop=1:nlines
 end // if disp_code >= 1 then
 xset('window',w0);   // push old windows
endfunction
//----------------------------------------------------
function plot_area2(ni,nth)
plot2d(0,0,1,"010","",[-5,-10,20,10],[5,5,4,5]);  // wake wo egaku
if nth == 1 then
 wstr1 = ' No. ' + string(ni) + '   :memo(' + memo1 + ')';
else
 wstr1 = ' No. ' + string(ni)
end
xstring(-4,6, wstr1);
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
//------------------------------------------
function [db_all,phi_all]=make_bode()
//
 ni=1; // yg is only 1st calculate
 frq=[Min_Freq:25:Max_Freq];
 frqs=size(frq);
 resp=zeros(1,frqs(2));
 respt=zeros(1,frqs(2));
 resph=zeros(1,frqs(2));
 respa=zeros(1,frqs(2));
 db_all=zeros(frqs(2),D_QTY);
 phi_all=zeros(frqs(2),D_QTY);
 //
 disp(' + please wait for some time to cacluate bode ');
 //
 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;
 resph(w)=hpf1(2 * %pi * frq(w) );
end  // for w
//
for vloop=1:D_QTY
 for w=1:frqs(2)
  if T_QTY == 2 then
   respt(w)= two_tubes2(2 * %pi * frq(w), vloop,rg_closed ,rl);
  elseif (T_QTY == 3) | (T_QTY == 4) then
   respt(w)= three_tubes(2 * %pi * frq(w), vloop,rg_closed ,rl);
  end
  respa(w)=resp(w) * respt(w) * resph(w);
 end
 //
 [db_fft1,phi_fft1] =dbphi(respa);
 for w=1:frqs(2)
   db_all(w,vloop)=db_fft1(w);
   db_phi(w,vloop)=phi_fft1(w);
 end
//
end
disp(' - cacluation was done.  ');
endfunction
//----
// - event handler
//--------------------------------------------------
// - others function library
//-------------------------------------------------------------------------
function [D_QTY,T_QTY,lrpa,memo1x,memo2x]=sub_demo1( inum )
//  setting for demonstration
 if inum == 2 then   // when /ua/
 disp('+ a sample of /ua/ is setting. ');
memo1x='a sample of /ua/';
D_QTY=20;
T_QTY=2;
memo2x=' no, l1, r1, pitch(+-%), amplitude(0-1)';
lrpa=[ 1 , 0 , 0 , -10 , 0.1 ;
       2 , 0 , 0 , -7 , 0.2 ;
       3 , 0 , 0 , -4 , 0.3 ;
       4 , 0 , 0 , -3 , 0.4 ;
       5 , 0 , 0 , -2 , 0.5 ;
       6 , 0 , 0 , -1 , 0.6 ;
       7 , 0 , 0 , 0 , 0.6 ;
       8 , 0 , 0 , 0 , 0.6 ;
       9 , 0 , 0.1 , 0 , 0.6 ;
       10 , 0 , 0.2 , 0 , 0.6 ;
       11 , 0 , 0.3 , 0 , 0.6 ;
       12 , 0 , 0.4 , 0 , 0.7 ;
       13 , 0 , 0.5 , 0 , 0.7 ;
       14 , 0 , 0.6 , 0 , 0.8 ;
       15 , 0 , 0.7 , 0 , 0.8 ;
       16 , 0 , 0.8 , 0 , 0.9 ;
       17 , 0 , 0.8 , 0 , 0.9 ;
       18 , 0 , 0.8 , 0 , 0.7 ;
       19 , 0 , 0.8 , -3 , 0.4 ;
       20 , 0 , 0.8 , -5 , 0.1 ] ;
 elseif inum == 1 then // when /ae/
disp('+ a sample of /ae/ is setting. ');
memo1x='a sample of /ae/';
D_QTY=20;
T_QTY=2;
memo2x=' no, l1, r1, pitch(+-%), amplitude(0-1)';
lrpa=[ 1 , 0 , 0.75 , -10 , 0.1 ;
       2 , 0 , 0.75 , -7 , 0.2 ;
       3 , 0 , 0.75 , -4 , 0.3 ;
       4 , 0 , 0.75 , -3 , 0.4 ;
       5 , 0 , 0.75 , -2 , 0.5 ;
       6 , 0 , 0.75 , -1 , 0.6 ;
       7 , 0 , 0.75 , 0 , 0.6 ;
       8 , 0 , 0.75 , 0 , 0.6 ;
       9 , -0.1 , 0.75 , 0 , 0.6 ;
       10 , -0.1 , 0.75 , 0 , 0.6 ;
       11 , -0.25 , 0.75 , 1 , 0.6 ;
       12 , -0.35 , 0.75 , 2 , 0.6 ;
       13 , -0.45 , 0.75 , 3 , 0.6 ;
       14 , -0.45 , 0.75 , 3 , 0.6 ;
       15 , -0.55 , 0.75 , 3 , 0.6 ;
       16 , -0.55 , 0.75 , 3 , 0.6 ;
       17 , -0.55 , 0.75 , 3 , 0.6 ;
       18 , -0.55 , 0.75 , 2 , 0.6 ;
       19 , -0.55 , 0.75 , 0 , 0.5 ;
       20 , -0.55 , 0.75 , 0 , 0.1 ] ;
 end
endfunction
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// proccess step0
if xmode2==0 then
 if T_DEMO == 1 then
   [D_QTY,T_QTY,lrpa,memo1,memo2]=sub_demo1( sel_demo1 );
   [A1,A2,A3,L1,L2,L3,tclosed,trise,tfall,amplitude]=make_paras();
 else
   [D_QTY,T_QTY,lrpa,memo1]= set_QTY();
   [lrpar,memo2]=edit_lrpa();
   [A1,A2,A3,L1,L2,L3,tclosed,trise,tfall,amplitude]=make_paras();
 end
end // if xmode2==0 then
//
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//   proccess step1 to step 4
[rg,yg,N1,N2,N3,LL,length0]=step1();
[tu1,tu2,tu3,r1,r2,r12,r21,r31]=make_tu_r();
[y2tm,M1,M2,M3]=step2();
[y2tm,y2tm_lpf,y_ran2,al1sm,bl0sm,bl1sm]=step2_2();
[y3out,ah1,bh0,bh1]=step3();
step4();
//
if T_DEMO == 1 then
  [db_all,phi_all]=make_bode();
  eh_var1=1; eh_var2=1; eh_var3=3;
  w0x=xget('window');
  xset('window', eh_var1);
  seteventhandler('my_eventhandler');
  plot_areas(eh_var1, eh_var2, eh_var3);
  xset('window', w0x);
  //
  disp('+usage: key in d (=down scroll) or u (=up scroll) on Scilab Graphic (1) windows. ');
end
//
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//-------------------------------------------------------------------------
//================================================================================
//
// 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
//
//
//-----------------------------------------------------------------
//
//
//----------------------------------------------------------------
// for sound wav file
f412=1;  // flag to detect scilab 4.1.2
defch1=1;  // default channel 1
f_win=1;  // flag if windows
bit1=16;
fs1=fs;
//--- check platform is windows -----------------------------------
if isdir('c:\\') then
 f_win=1;
else
 f_win=0;
end
//-----------------------------------------------------------------
function [wavfile1,chs1,qty1,fs1,bit1,f412]=get_wavfile_pro()
// choose wave file
  wavfile1=tk_getfile(Title="Choose  any xxx.wav file name");
// read channels and samples
cs1=wavread(wavfile1,'size');
chs1=cs1(2);   // channel
qty1=cs1(1);   // samples
//...
if chs1 > 2 then   // invert data for scilab-4.1.2
 chs1=cs1(1);
 qty1=cs1(2);
 f412=1;
else
 f412=0;
end
//...
[y,fs1,bit1]=wavread(wavfile1,[1 1]);
endfunction
//-----------------------------------------------------------------
function snd_play2(wdat2)
// this sound doesnot work well on windows scilab-3.1.1 and linux scilab-3.1
// but, this sound works on windows scilab-4.1.2. But, linux scilab-4.1.2 doesnot!
// for wdat2
sz=size(wdat2);
size2=sz(2);

if f_win == 1 then
if f412 == 1 then    // for windows scilab-4.1.2
 sound(wdat2 ,fs1,bit1);
else                 // for windows scilab-3.1.1
 if fs1 == 22050 then
  sound(wdat2',fs1,bit1);
  disp('This function may work, if luckily.');
 elseif fs1 == 44100 then    // down-sampling from 44100 to 22050
    wdatw2=zeros(1,(int(size2/2)));
    for v=1:(int(size2/2))
    wdatw2(v)=wdat2(2 * v );
    end
    disp('down-sampling from 44100 to 22050');
    sound(wdatw2',22050,bit1);
    disp('This function may work, if luckily.');
 else
  //sound(wdatw,fs1,bit1);
  disp('This function does not work.');
 end
 end
else
 disp('If sound setting is good for scilab, this function maybe work.');
 if f412 == 1 then    // for scilab-4.1.2
  sound(wdat2 ,fs1,bit1);
 else                 // for scilab-3.1.1
  if fs1 == 22050 then
   sound(wdat2',fs1,bit1);
   disp('This function may work, if luckily.');
  elseif fs1 == 44100 then    // down-sampling from 44100 to 22050
    wdatw2=zeros(1,(int(size2/2)));
    for v=1:(int(size2/2))
    wdatw2(v)=wdat2(2 * v );
    end
    disp('down-sampling from 44100 to 22050');
    sound(wdatw2',22050,bit1);
    disp('This function may work, if luckily.');
  else
   //sound(wdatw,fs1,bit1);
   disp('This function does not work.');
  end
 end
end
endfunction
//
//----------------------------------------------------------------
function snd_save2(wdat2)
// for wdat2
//
wavefilename= input(' + enter file name for saved .wav file =>',["string"]);
//
if f412 == 1 then    // for windows scilab-4.1.2
 wavwrite(wdat2 ,fs1,bit1,wavefilename );
else                 // for windows scilab-3.1.1
 wavwrite(wdat2',fs1,bit1, wavefilename);
end
endfunction
//----------------------------------------------------------------
function save_text1()
// save lrpa as a text file
 u=file('open','lrpa.txt','unknown');
 wstr1=date();
 dt=getdate();
 wstr1=wstr1 + ' ' + string(dt(7)) + ':' + string(dt(8)) + ':' +string(dt(9));
 fprintf(u,'//   %s',wstr1);
 fprintf(u,'D_QTY=%d;',D_QTY);
 fprintf(u,'T_QTY=%d;',T_QTY);
 fprintf(u,'memo1=''%s'';',memo1);
 fprintf(u,'//  %s',memo2);
 s0=size(lrpa);
 for v=1:s0(1)
  if v == 1 then
    wstr1='lrpa=[ ';
    for w=1:s0(2)
      if w== s0(2) then
        wstr1 = wstr1 + string( lrpa(v,w) ) + ' ; ';
      else
        wstr1 = wstr1 + string( lrpa(v,w) ) + ' , ';
      end
    end
    fprintf(u,'%s',wstr1);
  elseif v == s0(1) then
    wstr1='       ';
    for w=1:s0(2)
      if w== s0(2) then
        wstr1 = wstr1 + string( lrpa(v,w) ) + ' ] ;';
      else
        wstr1 = wstr1 + string( lrpa(v,w) ) + ' , ';
      end
    end
    fprintf(u,'%s',wstr1);
  else
   wstr1='       ';
    for w=1:s0(2)
      if w== s0(2) then
        wstr1 = wstr1 + string( lrpa(v,w) ) + ' ; ';
      else
        wstr1 = wstr1 + string( lrpa(v,w) ) + ' , ';
      end
    end
    fprintf(u,'%s',wstr1);
  end
 end
 file('close',u);
endfunction
//
//----------------------------------------------------------------
//
//
ADD_MENU_EDIT=1;       // set ADD_MENU 1 to add menu button of some functions
//
//
if ADD_MENU_EDIT == 1 then
delmenu('step1_edit');
addmenu('step1_edit',['(1)set duration','(2)edit lrpa','(3)make parameters','(-)save lrpa as a file','(-)load a lrpa file','(--)save lrpa as a text file as lrpa.txt']);
step1_edit=['[D_QTY,T_QTY,lrpa,memo1]= set_QTY()','[lrpa,memo2]=edit_lrpa()','[A1,A2,A3,L1,L2,L3,tclosed,trise,tfall,amplitude]=make_paras()','save_lrpa()','[D_QTY,T_QTY,lrpa,memo1]=load_lrpa()','save_text1()'];
end //if ADD_MENU == 1 then/
//
//----------------------------------------------------------------
//
//
ADD_MENU_GEN=1;       // set ADD_MENU 1 to add menu button of some functions
//
//
if ADD_MENU_GEN == 1 then
delmenu('step2_wave_generate');
addmenu('step2_wave_generate',['(1)wave generate','(-)plot portion','(--)set fft points ']);
step2_wave_generate=['[rg,yg,N1,N2,N3,LL,length0]=step1(); [tu1,tu2,tu3,r1,r2,r12,r21,r31]=make_tu_r(); [y2tm,M1,M2,M3]=step2();[y2tm,y2tm_lpf,y_ran2,al1sm,bl0sm,bl1sm]=step2_2();[y3out,ah1,bh0,bh1]=step3();step4();disp('' - done. '');'  ,'[sp1, ep1]= set_sp1_ep1(y3out); display_portion(2,y3out)'   ,'[fft_point1]=set_fft_points()'];
end //if ADD_MENU == 1 then/
//-----------------------------------------------------------------
//
ADD_MENU_EH=1;   // set ADD_MENU 1 to add menu button of some functions
//
//
if ADD_MENU_EH == 1 then
delmenu('step3_check_area');
addmenu('step3_check_area',['(1)preparation bode','(2) set event handler window(1), key in d(=down) or u(=up)','(e) suppress handler']);
step3_check_area=[ '[db_all,phi_all]=make_bode();', 'eh_var1=1; eh_var2=1; eh_var3=3; w0x=xget(''window''); xset(''window'', eh_var1); seteventhandler(''my_eventhandler''); plot_areas(eh_var1, eh_var2, eh_var3); xset(''window'', w0x);' ,  'seteventhandler('''') ' ];
end // if ADD_MENU_EH == 1 then
//
//-----------------------------------------------------------------
//
ADD_MENU_PLAY=1;       // set ADD_MENU 1 to add menu button of some functions
//
//
if ADD_MENU_PLAY == 1 then
delmenu('menu_sound');
addmenu('menu_sound',['(1)check wav version','(2)play sound','(3)save sound']);
menu_sound=['[xwavfile1,xchs1,xqty1,xfs1,xbit1,f412]=get_wavfile_pro()','snd_play2(y3out)','snd_save2(y3out)'];
end //if ADD_MENU == 1 then/
//
//===========================================================================================
//
//  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/
//
//===========================================================================================
//
ADD_MENU_STATUS=1;
if ADD_MENU_STATUS == 1 then
function status1()
 ver0= getversion();
 wstr1=sprintf('scilab version is %s',ver0);
 disp(wstr1);
endfunction
function status2()
 sz=stacksize();
 gsz=gstacksize();
 wstr1=sprintf('current stack size is %d  of %d  ,(%f)', sz(2),sz(1), sz(2)/sz(1) );
 wstr2=sprintf('current global stack size is %d  of %d ,(%f)', gsz(2),gsz(1), gsz(2)/gsz(1));
 disp( wstr1);
 disp( wstr2);
endfunction
delmenu('status_');
addmenu('status_',[ 'scilab version' ; 'check stack memory size' ]);
status_=[ 'status1()' ; 'status2()' ] ;
end //if ADD_MENU_STATUS == 1 then
//
//----------------------------------------------------------------------------------------------------//


No.2c   2008年6月5日