その３のページにもどる
 目次
 その５の ページにすすむ



< FFT分析スペクトルと２ 管模型からのスペクトルの比較 その4 >

rl の推定



口の内での反響の強さを示す　rl 　を推定することを考えてみましょう。

rl は管模型の出口側の反射係数であり、定義域は



である。　rlが０であると管模型の出口での反射はなくなり　最後の管を素通して出てくるようなイメージになる。逆に　rl が１だと１００％　管の内部に反射して戻って、出口から出てこなくなる（はず理論的には）。


下記は２管模型の「あ」において　ｒｌ　を変化させたときの　スペクトルの様子を示す。
rl　が１に近いほど、２つの山頂をもつ山は険しくなっていく。　
逆に、ｒｌ　が０にちかいと山頂が一つの山になっていく。第１管のスペクトルに近くなることが予想される。






ちなみに、この「あ」の例の場合だと、　ｒｌ　は下記のように限度一杯（上限は１）になっている。








参考に、SCILABのデモ　プログラム　a_wavfile_edit_278b.sci　を載 せておきます。　入力の音声サンプルは、その１のページにあります。


このプログラムを実 際に動作させるためには、 ご使用されるプログラミング環境に合わせて修正・変更・追加などが必要かもしれません。 また、バグが含まれている可能性があります。 万一、このデモンストレーション用プログラムを動作させる場合は、あなたの責任でおこなってくださいね。

//-----------------------------------------------------------------------------------------
//   A .wav Edit  and   Compare with 2 tubes model's wave or 3 tubes model's wave.  
//
//      a trial of leastsq method to estimate tube model parameters including rl            
//      Using focus weight to compare (local)part of object, evaluation by siguma |wmf(i) *(x(i)-y(i))|^2 between spectrum
//      And  hoping into limit by addition of (1 + exp( fact0 * x0)) * (1 + exp(fact0 * x1))
//     
//
//   for window scilab-4.1.2
//
// .........................................................................................
// PLEASE PAY ATTENSION that this program may have some bugs and also you may adjust program 
// to fit your circumstances. If you use this program, everything should be done at your own risk.
// Thanks.
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
T_DEMO=1;   // set T_DEMO=1 for demonstration, beside set T_DEMO=0 when normal
SEL_CODE=1; // dummy set
//
//--------------------------------------------------------------------
// Get Scilab Version
// global
   scilab_version_number=4;  // dummy set
//
   ver0=getversion();
   p0=strindex(ver0,'-');
   p1=strindex(ver0,'.');
   ver1=part(ver0,(p0(1)+1):(p1(1)-1));
   scilab_version_number=str2code(ver1)
//---------------------------------------------------------------------
// some functions to replace from scilab-4.1.2 to scilab-5.1.1
function [xmode00]=x_choose0(a,b,c)
  if scilab_version_number >= 5 then
    xmode00=x_choose(a,['Please select one and double-click it'],c);
  else
    xmode00=x_choose(a,b,c);
  end
endfunction
//----------------------------------------------------------------------
// Sound format has a difference among scilab version.
// for sound wav file
if scilab_version_number == 4 then
 f412=1;  // flag to detect scilab 4.1.2
 f511=0;
elseif scilab_version_number >= 5 then
 f412=1;
 f511=1;  // flag to detect scilab 5.1.1
else
 f412=0;
 f511=0;
end
//---------------------------------------------------------------------
// get Scilab current directory
getcwd
//===================================================================
// global variables   No. 1
// for .wav 1
chs1=0;    // channels
qty1=0;    // data quantity,  samples
size1=0;   // actual loaded data quantity
fs1=0;     // sampling frequency
bit1=0;    // bits
wavfile1=''; //  file name,  this will be overwritten.
wdat1=zeros(1,10);// data of the .wav,  only one first channel, this will be overwritten.
// for fft
fft_point1=512;
db_fft1=zeros(1,512);   // this will be overwritten.
phi_fft1=zeros(1,512);  // this will be overwritten.
Min_Freq=150;   // set plot maximum frequency by unit is [Hz]
Max_Freq=7500;   // set plot maximum frequency by unit is [Hz]

// for display
width0=400;   // default display width
ydat0=zeros(1,width0+1);  // for all data
ydat1=zeros(1,width0+1);  // for portion
xziku0=[0:1:width0];  // for all data
xziku1=[0:1:width0];  // for portion
sp1=1;       // start point for actual display
ep1=1;       // end point for actual display

// for others
defch1=1;  // default channel 1
f_win=1;  // flag if windows
//---------------------------------------
//
function save_wavprop()
    save('wavprop.dat',fs1,size1,bit1,wdat1,sp1,ep1,ydat0,ydat1,xziku0,xziku1);
endfunction
function [fs1,size1,bit1,wdat1,sp1,ep1,ydat0,ydat1,xziku0,xziku1]=load_wavprop()
   load('wavprop.dat','fs1','size1','bit1','wdat1','sp1','ep1','ydat0','ydat1','xziku0','xziku1');
endfunction
//
function save_wfft()
    save('wavfft.dat',fft_point1,db_fft1,phi_fft1);
endfunction
function [fft_point1,db_fft1,phi_fft1]=load_wfft()
    load('wavfft.dat','fft_point1','db_fft1','phi_fft1');
endfunction
//
//--------------------------------------
function [wavfile1,chs1,qty1,fs1,bit1,f412,SEL_CODE]=get_wavfile_pro()
// choose wave file
if T_DEMO == 1 then
 //+ select
 SEL_CODE=x_choose0(['a_sample';'o_sample'; 'i_sample'; 'manual select'],['Please select one'],'default')
 if SEL_CODE == 4 then  // when manual select
    SEL_CODE = 999;
 elseif SEL_CODE <= 0 then
   SEL_CODE=1;
 end   
 //- select
 if SEL_CODE == 1 then  // a_sample
   disp(' This is demonstration. Copy this file and a_sample.wav to your scilab''s   bin or home   directory.');
   [x,ierr]=fileinfo('a_sample.wav');
   xs=size(x);
   if xs(1)==0 then
     disp('Choose a_sample.wav file');
     if scilab_version_number >= 5 then
       wavfile1=uigetfile(["*.wav"],"","Choose  a_sample.wav file");
     else
       wavfile1=tk_getfile(Title="Choose a_sample.wav file");
     end
   else
     wavfile1='a_sample.wav';
   end
 elseif SEL_CODE == 2 then   // o_sample
   disp(' This is demonstration. Copy this file and o_sample.wav to your scilab''s   bin or home   directory.');
   [x,ierr]=fileinfo('o_sample.wav');
   xs=size(x);
   if xs(1)==0 then
     disp('Choose o_sample.wav file');
     if scilab_version_number >= 5 then
       wavfile1=uigetfile(["*.wav"],"","Choose  o_sample.wav file");
     else
       wavfile1=tk_getfile(Title="Choose o_sample.wav file");
     end
   else
     wavfile1='o_sample.wav';
   end
 elseif SEL_CODE == 3 then   // i_sample
   disp(' This is demonstration. Copy this file and i_sample.wav to your scilab''s   bin or home   directory.');
   [x,ierr]=fileinfo('i_sample.wav');
   xs=size(x);
   if xs(1)==0 then
     disp('Choose i_sample.wav file');
     if scilab_version_number >= 5 then
       wavfile1=uigetfile(["*.wav"],"","Choose  i_sample.wav file");
     else
       wavfile1=tk_getfile(Title="Choose i_sample.wav file");
     end
   else
     wavfile1='i_sample.wav';
   end
 else  // include SEL_CODE == 999
   if scilab_version_number >= 5 then
    wavfile1=uigetfile(["*.wav"],"","Choose  a .wav file");
   else
    wavfile1=tk_getfile(Title="Choose a .wav file name");
   end
 end

else
  if scilab_version_number >= 5 then
    wavfile1=uigetfile(["*.wav"],"","Choose  a .wav file");
  else
    wavfile1=tk_getfile(Title="Choose a .wav file name");
 end
 SEL_CODE = 999;
end
// 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 display the .wav property
function disp_pro_wav()
 disp(qty1,'data quantity',bit1, 'bits',chs1,'channels',fs1,'sampling frequency');
endfunction
//--------------------------------------
function [wdat1, size1]=read_one_ch_of_wav(wavfile1)
 y=wavread(wavfile1);
 ysize=size(y);
 if ysize(2) == chs1 then
   wdat1=zeros(1,ysize(1));
   if chs1 == 1 then
//     for v=1:ysize(1)
//       wdat1(v)=y(v);
//    end
       wdat1=y';
   else
    disp('Channels are more than two, but, only 1st channel is stored.');
    for v=1:ysize(1)
      wdat1(v)=y(v,defch1);
    end
   end
   size1=ysize(1);
   disp(size1,'stored data size is');
 elseif ysize(1) == chs1 then  // for scilab-4.1.2
 wdat1=zeros(1,ysize(2));
 if chs1 == 1 then
//     for v=1:ysize(2)
//       wdat1(v)=y(v);
//     end
    wdat1=y;
   else
    disp('Channels are more than two, but, only 1st channel is stored.');
    for v=1:ysize(2)
      wdat1(v)=y(defch1,v);
    end
   end
   size1=ysize(2);
   disp(size1,'stored data size is');
 end
endfunction
//-----------------------------------------------
function plot_wave1(disp0)
 wb0=xget('window');  // stack old window
 xset('window',disp0);   // create new windows
 clf();
 subplot(211);
 //plot(xziku0,ydat0,'b');
 plot2d(xziku0,ydat0,style=[color("blue")]);
 xtitle('waveform all');
 // When linux scilab-3.1 2nd subplot does not work well.
 subplot(212);
 //plot(xziku1,ydat1,'b');
 plot2d(xziku1,ydat1,style=[color("blue")]);
 xtitle('waveform selected');
 xset('window',wb0);   // push old windows
endfunction
//--------------------------------------
function [ydat1,xziku1]=make_width_data( work1, wsize1)
wsize2=wsize1 - sp1 + 1;
wsize3=ep1-sp1+1;
if wsize3  > (width0+1) then
 wstep1= wsize3 / (width0+1);
 for v=1:(width0+1)
  ydat1(v)= work1( sp1 + int(wstep1 * (v - 1)));
  xziku1(v)= sp1 + int(wstep1 * (v - 1));
 end
else
 for v=1:(width0+1)
  if v <= wsize3 then
   ydat1(v)=work1((sp1-1)+v);
  else
   ydat1(v)=0.;
  end
   xziku1(v)=(sp1-1)+v;
 end
end
//
// plot_wave1(0);
endfunction
//----------------------------------------------------------------
function [wsp1, wep1]= reset_sp1_ep1(wsize1)
 wsp1=1;
 wep1=wsize1;
endfunction
//----------------------------------------------------------------
function [wsp1, wep1]= set_sp1_ep1(wsize1)
 txt1=['start point';'end point'];
 wstr1=sprintf('%d',sp1);
 wstr2=sprintf('%d',ep1);
 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
 //
//disp(arg2,'arg2',arg1,'arg1');
//
 wsp1=sp1;
 wep1=ep1;
 if arg1 >= 1 then
  if arg2 <= wsize1 then
   if arg1 < arg2 then
     wsp1=arg1;
     wep1=arg2;
     disp(wep1,'ep1',wsp1,'sp1');
   end
  end
 end
endfunction
//----------------------------------------------------------------
function snd_play1()
// 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!
wdatw=zeros(ep1-sp1+1);
for v=sp1:ep1
 wdatw(v-sp1+1)=wdat1(v);
end

if f_win == 1 then
if f412 == 1 then    // for windows scilab-4.1.2
 sound(wdatw' ,fs1,bit1);
else                 // for windows scilab-3.1.1
 if fs1 == 22050 then
  sound(wdatw,fs1,bit1);
  disp('This function may work, if luckily.');
 elseif fs1 == 44100 then    // down-sampling from 44100 to 22050
    wdatw2=zeros((ep1-sp1)/2+1);
    for v=1:((ep1-sp1)/2+1)
    wdatw2(v)=wdat1(sp1 + 2 * (v -1) );
    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('Sorry, but, this function does not work.');
end
endfunction
//
//----------------------------------------------------------------
function snd_save1()
//
wavefilename= input(' + enter file name for saved .wav file =>',["string"]);
//
wdatw=zeros(ep1-sp1+1);
for v=sp1:ep1
 wdatw(v-sp1+1)=wdat1(v);
end
if f412 == 1 then    // for windows scilab-4.1.2
 wavwrite(wdatw' ,fs1,bit1,wavefilename );
else                 // for windows scilab-3.1.1
 wavwrite(wdatw,fs1,bit1, wavefilename);
end
endfunction
//--- check platform is windows -----------------------------------
if isdir('c:\\') then
 f_win=1;
else
 f_win=0;
end
//
//=================================================================
//
//   +MAIN (1)program starts
//
[wavfile1,chs1,qty1,fs1,bit1,f412,SEL_CODE]=get_wavfile_pro();
disp_pro_wav();
[wdat1, size1]=read_one_ch_of_wav(wavfile1);
sp1=1;
ep1=size1;
[ydat0,xziku0]=make_width_data( wdat1, size1);
if T_DEMO==1 then
 if SEL_CODE == 1 then // a_sample
  sp1=1900;
  ep1=2413;
 elseif SEL_CODE == 2 then // o_sample
  sp1=2700;
  ep1=3300;
 elseif SEL_CODE == 3 then // i_sample
  sp1=4100;
  ep1=4700;
 else   // include SEL_CODE == 999
  [sp1, ep1]= set_sp1_ep1(size1);
 end
end
[ydat1,xziku1]=make_width_data( wdat1, size1);
plot_wave1(0);
//
// -MAIN (1)program starts
// select .wav and set portion to make it fft
//==============================================================++=
//
//
//
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 agwru tame
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 plot_fft1(disp0)
 wb0=xget('window');  // stack old window
 xset('window',disp0);   // create new windows
 [frq1,sfrq1,is1,ie1]=set_frq(0);
 clf();
 subplot(211);
 gainplot(frq1,db_fft1,phi_fft1);
 xtitle('frequency response of waveform selected by fft analysis');
 subplot(212);
 //plot(xziku1,ydat1,'b');
 plot2d(xziku1,ydat1,style=[color("blue")]);
 xtitle('waveform selected');
 xset('window',wb0);   // push old windows
endfunction

//----------------------------------------------------
function  [db_fft1,phi_fft1]=do_fft_wav()
 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);
 end   // -if size1 <= ( sp1 + fft_point1 - 1 ) then
endfunction
//=================================================================
//  peak detect on FFT spectrum
// 
//    (1)2-3 smoothing on FFT data / nizi sanzi takousiki tekigou niyoru smoothing
//    (2)differentiation (def)     / heikatsuka bibun
//    (3)detect peak as the point from plus to minus
//
sn=0;  // data quantity.  this will be overwritten.
sm1=2;  // set order of 2-3 smoothing
sm2=2;  // set order of 5 points  differentiation (def)
        // if sm2=1, more peak candidate may appear
        // if sm2=2, sometime gentle slop (nadaraka) peak will be lost

swnd= zeros(1,128);  // swnd(1),swnd(2),....  this will be overwritten.

sm1out=zeros(1,512);   // this will be overwritten. calculate weighted average
sm2out=zeros(1,512);   // this will be overwritten. def of sm1out

npk=0;
pklist=zeros(9,512);  // nth, freq, its value of peak candinates
                      // nth, freq, its value of estimated left edge
                      // nth, freq, its value of estimated right edge
//------------------------------------------------------
function plot_fft_sm1(disp0)
 [frq1,sfrq1,is1,ie1]=set_frq(0);

 w0=xget('window');
 xset('window',disp0);
 clf();
 subplot(211);
 gainplot(frq1,db_fft1,phi_fft1);
 //plot2d(frq1,db_fft1,color("black"),logflag="ln");
 xtitle('frequency response of waveform selected by fft analysis');
 plot2d(frq1,sm1out,color("green"),logflag="ln");
 legend(['org'; 'smoothed'],3);
 //xtitle('smoothed frequency response ');
 subplot(212);
 plot2d(frq1,sm2out,color("red"),logflag="ln");
 xtitle('def of smoothed frequency response ');
 xset('window',w0);
endfunction
//----------------------------------------------------------
function [sn,swnd,sm1out,sm2out,npk,pklist]=smooth1(disp_code)
 // get db_fft1 data size
 ax=size(db_fft1);
 sn=ax(2);
 swnd= zeros(1,128);
 // preparetion smoothing 2-3
 sm=sm1;
 snorm=(4.0 * sm * sm - 1.0) * (2.0 * sm + 3.0 ) / 3.0;
 sl=3.0 * sm * (sm + 1.0) -1.0;
 for loop=1:(sm+1)
  swnd(loop)= sl - 5.0 * (loop -1.0) * (loop - 1.0);
 end
 // calculate weighted average
 sm1out=zeros(sn);
 dmax0= 0.;
 sigma0=0.;
 ep0=0.;
 for v=1:sn
   if v < (sm+1) then
     sm1out(v)=db_fft1(v);
   elseif v > (sn-sm) then
     sm1out(v)=db_fft1(v);
   else
     sum0 = db_fft1(v) * swnd(1);
     for v2=1:sm
       sum0 = sum0 + db_fft1(v+v2) * swnd(1+v2);
       sum0 = sum0 + db_fft1(v-v2) * swnd(1+v2);
     end
     sm1out(v)= sum0 / snorm;
     //
     sa0=db_fft1(v) - sm1out(v);
     sigma0=sigma0+sa0*sa0;
   end
     ep0=ep0+db_fft1(v);
     if sm1out(v) > dmax0 then
       dmax0= sm1out(v);
     end
  end
 //
 // calculate  def of smoothed
 sm=sm2;
 sm2out=zeros(sn);
 for v=1:sn
  if v<(sm+1) then
   sm2out(v)=0.;
  elseif v>(sn-sm) then
   sm2out(v)=0.;
  else
   sm2out(v)=0.;
   for v2=1:sm
     sm2out(v)=sm2out(v)+ ( sm1out(v+v2) - sm1out(v-v2)) * v2;
   end
  end
 end
// find from + to - in def
pklist=zeros(9,sn);
[frq1,sfrq1,is1,ie1]=set_frq(0);
npk=0;
for v=2:sn
  if sm2out(v-1) > 0 then
    if sm2out(v) < 0 then
      dmax1=-9999.0;
      stack_v=-1;
      //.. ..
      for v2=(v-sm):(v+sm)
        if v2 < 1 then
        elseif v2 > sn then
        else
          if db_fft1(v2) > dmax1 then
           dmax1=db_fft1(v2);
           stack_v=v2;
          end
        end
      end
      //.. ..
      if stack_v > 0 then
        npk=npk+1;
        pklist(1,npk)=stack_v;
        pklist(2,npk)=frq1(stack_v);
        pklist(3,npk)=dmax1;
      end
    end
  end
end
// find left edge estimated based on smoothed signal
 for loop=1:npk
   if pklist(1,loop) == 1 then
      pklist(4,loop) = 1;
      pklist(5,loop) = frq1(1);
      pklist(6,loop) = sm1out(1);
   else
     pklist(4,loop) = 1;
     pklist(5,loop) = frq1(1);
     pklist(6,loop) = sm1out(1);
    for v=pklist(1,loop):-1:2
     if sm1out(v-1) > sm1out(v) then
      pklist(4,loop) = v;
      pklist(5,loop) = frq1(v);
      pklist(6,loop) = sm1out(v);
      break;
     end
    end
 
   end
 end  // for loop=1:npk
//..
// find right edge estimated based on smoothed signal
 for loop=1:npk
   if pklist(1,loop) == sn then
      pklist(7,loop) = sn;
      pklist(8,loop) = frq1(sn);
      pklist(9,loop) = sm1out(sn);
   else
     pklist(7,loop) = sn;
     pklist(8,loop) = frq1(sn);
     pklist(9,loop) = sm1out(sn);
    for v=pklist(1,loop):1:(sn-1)
     if sm1out(v+1) > sm1out(v) then
      pklist(7,loop) = v;
      pklist(8,loop) = frq1(v);
      pklist(9,loop) = sm1out(v);
      break;
     end
    end
 
   end
 end  // for loop=1:npk
//...
 if disp_code == 1 then
  disp('');
  //-- peak candinate list out
  disp('-- peak candinate list out --');
  disp('-- peak [Hz] [dB]  left edge [Hz] [dB]  right edge [Hz] [dB] --');
  for v=1:npk
  // disp([pklist(1,v) pklist(2,v) pklist(3,v)]);
    disp([ pklist(2,v) pklist(3,v) pklist(5,v) pklist(6,v) pklist(8,v) pklist(9,v)]);
  end
 //-- peak candinate list out
 disp('');
 end
endfunction
//
//=================================================================
//
//   +MAIN (2)program starts
//    about fft
if T_DEMO==1 then
 [db_fft1, phi_fft1]=do_fft_wav();
 plot_fft1(1);
 //
 [sn,swnd,sm1out,sm2out,npk,pklist]=smooth1(1);
 plot_fft_sm1(2);
end
//
//  -MAIN (2)program starts
//
//==============================================================++=
//
//
// global variables  No. 2
//
c0=35000;   // constant.  the speed of sound is round 35000 cm/second
ts=1.0 / fs1;
//
L1=1;  // this will be overwritten.
L2=1;  // this will be overwritten.
L3=1;  // this will be overwritten.
A1=1;  // this will be overwritten.
A2=1;  // this will be overwritten.
A3=1;  // this will be overwritten.
//
ttl_Length=18.5;  // set  total length of combined tubes, Two Tubes, and etc               
ttl_Area=10;    // a constant in this program.  set total area of combined tubes, Two Tubes, and etc  for dummy display
l1=0; // this will be overwritten.
l2=0; // this will be overwritten.
rg=1; //  a constant in this program.  When glottis
r1=0.8; // this will be overwritten.
r2=0.8; // this will be overwritten.
rl=0.9;   // set adjust reduction response. rl is variable, however use one constant in this.
fc=1000;   // this will be overwritten. set cut off frequency of High Pass Filter by unit is [Hz]
trise=6.0; // this will be overwritten.
//tfall=0.7; // this will be overwritten.
tfall=2.0; // this will be overwritten.
tclosed=5.0;  // use ONLY in waveform making
//
yg_res1= ones(1,1024);  // this will be overwritten.
tube2_res1= ones(1,1024);  // this will be overwritten.
hpf_res1= ones(1,1024);  // this will be overwritten.
overall_res1=zeros(1,1024); // this will be overwritten.
db_2tube=zeros(1,1024); // this will be overwritten.
phi_2tube=zeros(1,1024); // this will be overwritten.
//
// others
WT_QTY=2;  // WT_QTY=2: standard 2 tubes model
N_REPEAT=7; // how many section to generate waveform ?
Wtubepropdat= 'tubeprop.dat';
//
//
//--------------------------------------------------------------------------
function save_2tube()
    save(Wtubepropdat,L1,L2,L3,A1,A2,A2,ttl_Length,ttl_Area,l1,l2,rg,r1,r2,fc,trise,tfall,tclosed,yg_res1,tube2_res1,hpf_res1,overall_res1,db_2tube,phi_2tube,WT_QTY);
endfunction
function [L1,L2,L3,A1,A2,A3,ttl_Length,ttl_Area,l1,l2,rg,r1,r2,fc,trise,tfall,tclosed,yg_res1,tube2_res1,hpf_res1,overall_res1,db_2tube,phi_2tube,WT_QTY]=load_2tube()
   load(Wtubepropdat,'L1','L2','L3','A1','A2','A3','ttl_Length','ttl_Area','l1','l2','rg','r1','r2','fc','trise','tfall','tclosed','yg_res1','tube2_res1','hpf_res1','overall_res1','db_2tube','phi_2tube','WT_QTY');
endfunction
//
function [WT_QTY]= set_tube_model()
 l1=list('model',(WT_QTY-1),['2 tubes','3 tubes  serial type']);
 //l1=list('model',1,['2 tubes']);
 wrep=x_choices('Select tube model',list(l1));
//
 if wrep== 1 then
  WT_QTY=2;
 elseif  wrep== 2 then
  WT_QTY=3;
 end
endfunction
//
//-------------------------------------------------
function [r1,r2,l1,l2,ttl_Length,rl,fc,trise,tfall,tclosed]= set_2tube_para()
 wstr1=sprintf('%f',r1);
 wstr2=sprintf('%f',l1);
 wstr3=sprintf('%f',ttl_Length);
 wstr4=sprintf('%f',fc);
 wstr5=sprintf('%f',trise);
 wstr6=sprintf('%f',tfall);
 wstr7=sprintf('%f',tclosed);
 wstr8=sprintf('%f',r2);
 wstr9=sprintf('%f',l2);
 wstr10=sprintf('%f',rl);
 //
 if WT_QTY == 3 then
  txt1=['r1';'l1'; 'r2'; 'l2'; 'total_tube_Length'; 'rl' ; 'fc'; 'trise'; 'tfall' ; 'tclosed (use only for generation)'];
  sig1=x_mdialog('Input parameters',txt1, [wstr1 ; wstr2 ; wstr8 ; wstr9; wstr3 ; wstr10 ; wstr4; wstr5; wstr6; wstr7 ]);
  if sig1==[] then
   r1=evstr(wstr1);
   l1=evstr(wstr2);
   r2=evstr(wstr8);
   l2=evstr(wstr9);
   ttl_Length=evstr(wstr3);
   rl=evstr(wstr10);
   fc=evstr(wstr4);
   trise=evstr(wstr5);
   tfall=evstr(wstr6);
   tclosed=evstr(wstr7);
  else
   r1=evstr(sig1(1));
   l1=evstr(sig1(2));
   r2=evstr(sig1(3));
   l2=evstr(sig1(4));
   ttl_Length=evstr(sig1(5));
   rl=evstr(sig1(6));
   fc=evstr(sig1(7));
   trise=evstr(sig1(8));
   tfall=evstr(sig1(9));
   tclosed=evstr(sig1(10));
  end
 else
  txt1=['r1';'l1'; 'total_tube_Length'; 'rl' ; 'fc'; 'trise'; 'tfall' ; 'tclosed (use only for generation)'];
  sig1=x_mdialog('Input parameters',txt1, [wstr1 ; wstr2 ; wstr3 ; wstr10; wstr4; wstr5; wstr6; wstr7 ]);
  if sig1==[] then
    r1=evstr(wstr1);
    l1=evstr(wstr2);
    ttl_Length=evstr(wstr3);
    rl=evstr(wstr10);
    fc=evstr(wstr4);
    trise=evstr(wstr5);
    tfall=evstr(wstr6);
    tclosed=evstr(wstr7);
    r2=evstr(wstr8);
    l2=evstr(wstr9);
  else
    r1=evstr(sig1(1));
    l1=evstr(sig1(2));
    ttl_Length=evstr(sig1(3));
    rl=evstr(sig1(4));
    fc=evstr(sig1(5));
    trise=evstr(sig1(6));
    tfall=evstr(sig1(7));
    tclosed=evstr(sig1(8));
    r2=evstr(wstr8);
    l2=evstr(wstr9);
   end
 end
endfunction
//
// Function plot sqrt(area) vs length
function plot_area(disp0)
 wb0=xget('window');  // stack old window
 xset('window',disp0);   // create new windows
//
if f412 == 1 then
 arrowsize1=20;
else
 arrowsize1=1;
end
clf();
plot2d(0,0,1,"010","",[-5,-10,30,10],[5,5,4,5]);  // wake wo egaku
if WT_QTY == 4 then  //Three Tubes Model of T type
ya1=sqrt(A1);
ya2=sqrt(A2);
ya3=sqrt(A3);
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;
xx1=L1 - (ya3 / 2);
xx2=L1 + (ya3 / 2);
xp=[0    0    xx1  xx1   xx2   xx2  L1+L2 L1+L2  L1  L1   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, arrowsize1, 5);
yar=[yy0 + ( ya2 / 2 )   yy0+ ( ya2 / 2 ) ];
xar=[ L1+L2+0.5 L1+L2+2];
xarrows(xar, yar, arrowsize1, 5);
xstring(xx1, yy3 - 1, "CLOSED");
xstring(0,-9,"Attension: Area scale is double than others Tube Models");
//
elseif WT_QTY == 3 then  //Three Tubes Model of serial type
ya1=sqrt(A1);
ya2=sqrt(A2);
ya3=sqrt(A3);
xp=[0    0    L1  L1 L1+L2 L1+L2 L1+L2+L3 L1+L2+L3 L1+L2 L1+L2  L1  L1 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, arrowsize1, 5);
xar=[ L1+L2+L3+0.5 L1+L2+L3+2];
xarrows(xar, yar, arrowsize1, 5);
else  //Two Tubes Model
ya1=sqrt(A1 );
ya2=sqrt(A2 );
xp=[0    0    L1  L1   L1+L2   L1+L2   L1  L1 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, arrowsize1, 5);
xar=[ L1+L2+0.5 L1+L2+2];
xarrows(xar, yar, arrowsize1, 5);
 if WT_QTY==5 then
   xar=[ L1+L2-0.5 L1+L2-2];
   xarrows(xar, yar, arrowsize1, 3);
 end
end

 xset('window',wb0);   // push old windows
endfunction
//
//-------------------------------------------------
// +hpf
function [hpf_res1]=hpf_response( point0)
wc=2. * %pi * fc;
al1= -1. * %e^( -1. * wc * ts);
bl0= wc * ts;
// Then, convert them to high pass filter's coefficients
fnew=fc;
wcnew=2. * %pi * fnew;
alfa= -1. *  cos( (wc + wcnew) * ts /2. ) / cos( (wc - wcnew) * 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);
// Function for disply high pass filter's frequency-phase response
[frq1,sfrq1,is1,ie1]=set_frq(point0);
hpf_res1=zeros(1,sfrq1);
for v=1:sfrq1
 xw=2 * %pi * frq1(v);
 yi= bh0  + bh1 * cos( xw * ts) - bh1 * sin( xw * ts ) * %i;
 yb=1.0  + ah1 * cos( xw * ts ) -  ah1 * sin( xw  * ts ) * %i;
  hpf_res1(v)= yi /yb;
end
endfunction
// -hpf
//----------------------------------------------------------
//
//
function [ tube2_res1, L1, L2, L3,  A1, A2, A3]= two_tubes2_resp()
if WT_QTY==3 then
//
// 3 tubes model definition
//
//    L1+L2+L3 = ttl_Length (total tubes length)
//
//   l1= (L2-L1)/(L2+L1)
//   l2= (L3-L2)/(L3+L2)
//
//    A1+A2+A3= ttl_Area 
//
//    r1= (A2-A1)/(A2+A1)
//    r2= (A3-A2)/(A3+A2)
//
  bunbo= l1 * l2 + ( l1 - l2 ) + 3.0;
  L1= ttl_Length * ( l2 - 1.0 ) * (l1 - 1.0 ) / bunbo ;
  L2= ttl_Length * ( l2 - 1.0 ) * ( -1.0 * l1 - 1.0 ) / bunbo ;
  L3= ttl_Length * ( l2 + 1.0 ) * (l1 + 1.0 ) / bunbo ;
  tu1=L1/c0;
  tu2=L2/c0;
  tu3=L3/c0;
  bunbo= r1 * r2 + ( r1 - r2 ) + 3.0;
  A1= ttl_Area * ( r2 - 1.0 ) * (r1 - 1.0 ) / bunbo ;
  A2= ttl_Area * ( r2 - 1.0 ) * ( -1.0 * r1 - 1.0 ) / bunbo ;
  A3= ttl_Area * ( r2 + 1.0 ) * (r1 + 1.0 ) / bunbo ;
  [frq1,sfrq1,is1,ie1]=set_frq(1024);
  tube2_res1=zeros(1,sfrq1);
  for v=1:sfrq1
   xw= 2 * %pi * frq1(v);
   //+
   yi  = 0.5 * ( 1. + rg ) * ( 1. + r1) * ( 1. + r2) * ( 1. + rl ) * %e^( -1. * ( tu1 + tu2  + tu3) * xw * %i);
   yb =  1 + rg * r1 * %e^( -2 * tu1 * xw * %i) + r1 * r2 * %e^( -2 * tu2 * xw * %i) + r2 * rl * %e^( -2 * tu3 * xw * %i);
   yb = yb + rg * r2 * %e^( -2 * (tu1 + tu2) * xw * %i) + r1 * rl * %e^( -2 * (tu2 + tu3) * xw * %i);
   yb = yb + rg * r1 * r2 * rl *  %e^( -2 * (tu1 + tu3) * xw * %i);
   yb = yb + rg * rl *  %e^( -2 * (tu1 + tu2 + tu3) * xw * %i);
   tube2_res1(v) = yi / yb;
   //-
  end
else
 L2= ( 1. + l1 ) * ttl_Length / 2.;
 L1= ttl_Length - L2;
 tu1=L1/c0;
 tu2=L2/c0;
 A2= ( 1. + r1 ) * ttl_Area/ 2.;
 A1= ttl_Area - A2;
 //
 [frq1,sfrq1,is1,ie1]=set_frq(1024);
 tube2_res1=zeros(1,sfrq1);
 for v=1:sfrq1
  xw= 2 * %pi * frq1(v);
  yi  = 0.5 * ( 1 + rg ) * ( 1 + r1)  * ( 1 + rl ) * %e^( -1 * ( tu1 + tu2 ) * xw * %i);
  yb =  1 + r1 * rg * %e^( -2 * tu1 * xw * %i) + r1 * rl * %e^( -2 * tu2 * xw * %i) + rl * rg * %e^( -2 * (tu1 + tu2) * xw * %i);
  tube2_res1(v) = yi / yb;
 end
 // +dummy data set
 L3=1;
 A3=1;
 // -dummy data set
end
endfunction
//
// +yg response
function [yg_res1]=yg_resp(points)
 N2= int( trise / 1000 / ts);
 N3= int( tfall / 1000 / ts);
 sizew=max((N2+N3+1),points);
 yg=zeros(1,sizew);
 //
 for tc0=1:sizew
  if  tc0 <= (N2+1)  then
    yg(tc0)= 0.5 * ( 1 - cos( ( %pi / N2 ) * (tc0 - 1) ) );
  elseif tc0<= (N2+N3+1) then
   yg(tc0)= cos( ( %pi / ( 2 * N3 )) * ( tc0 - (N2+1) ) );
  else
   yg(tc0)=0.0;
  end
 end

 [frq1,sfrq1,is1,ie1]=set_frq(points);
 yg_res1=zeros(1,sfrq1);
 // +do points(1024) fft and get portion of data
 if sizew == points then
    fftwout=fft( yg, -1 );
    v=1;
    is2= int(is1 * points / fft_point1);
    ie2= int(ie1 * points / fft_point1);
    for w=is2:ie2
      yg_res1(v)=fftwout(is2+v) / fftwout(1);
      v=v+1;
    end
 else   // -do 1024 fft and get portion of data
//
 y00=0.;
 for v=1:(N2+N3+1)
  y00= y00 + yg(v);
 end
//
 for w=1:sfrq1
 xw=2 * %pi * frq1(w) * ts;
 yi=0;
 for v=1:(N2+N3+1)
  yi= yi + yg(v) * %e^(-1. * xw * (v-1) * %i);
 end  // for v
 yg_res1(w)=yi / y00; 
 end  // for w
end
endfunction
// -yg response
//-------------------------------------------------
//
function [overall_res1,db_2tube, phi_2tube ]=overall_response()
[frq1,sfrq1,is1,ie1]=set_frq(1024);
 for v=1:sfrq1
   overall_res1(v)= hpf_res1(v) * tube2_res1(v) * yg_res1(v) ;
 end
 [db_2tube,phi_2tube] =dbphi(overall_res1);
endfunction
//-----------------------------------------------
function plot_2tube(disp0, kcode)
//                         When kcode =0, normal, beside kcode=1, leastsq result
 wb0=xget('window');  // stack old window
 xset('window',disp0);   // create new windows
 [frq1,sfrq1,is1,ie1]=set_frq(0);
 clf();
 subplot(211);
 gainplot(frq1,db_fft1,phi_fft1);
 xtitle('frequency response of waveform selected by fft analysis');
  [frq1,sfrq1,is1,ie1]=set_frq(1024);
 subplot(212);
 gainplot(frq1,db_2tube',phi_2tube');
 if kcode == 0 then
  xtitle('frequency response of tubes model by initial value for comparison to one of waveform selected by fft analysis ');
 elseif kcode == 1 then
  xtitle('frequency response of tubes model by leastsq method result from the initial value');
 elseif kcode == 3 then
 wstr0='   frequency response of tubes model';
 wstr1=sprintf('%f',l1);
 wstr2=sprintf('%f',r1);
 wstr3=sprintf('%f',l2);
 wstr4=sprintf('%f',r2);
 if WT_QTY == 3 then
   wstr5='l1=' + wstr1 + '  r1=' + wstr2 + '  l2=' + wstr3 + '  r2=' + wstr4 + wstr0;
 else  // when WT_QTY == 2 then
   wstr5='l1=' + wstr1 + '  r1=' + wstr2  + wstr0;
 end
   xtitle(wstr5);
end
 xset('window',wb0);   // push old windows
endfunction

//
//=================================================================
//
//   +MAIN (3)program starts
//    about matching to 2 tubes model
if T_DEMO==1 then
 trise=6.0; // just fit for 1st demo version
 tfall=0.7; // just fit for 1st demo version
 if SEL_CODE == 1 then // a_sample
   WT_QTY=2;
   r1=0.8;
   l1=0.;
   ttl_Length=18.5;
   rl=0.9;
 elseif SEL_CODE == 2 then // o_sample
   WT_QTY=3;
   r1=-0.5;
   l1=0.5;
   r2=0.8;
   l2=0.;
   ttl_Length=23.;
   rl=0.9;
 elseif SEL_CODE == 3 then // i_sample
   WT_QTY=2;
   r1=-0.2;
   l1=0.;
   ttl_Length=18.5;
   rl=0.9;
 else   // include SEL_CODE == 999
  [WT_QTY]= set_tube_model();
  [r1,r2,l1,l2,ttl_Length,rl,fc,trise,tfall,tclosed]= set_2tube_para();
 end
 [yg_res1]=yg_resp(1024);
 [ tube2_res1, L1, L2, L3, A1, A2, A3 ]= two_tubes2_resp();
 [hpf_res1]=hpf_response(1024);
 [overall_res1,db_2tube, phi_2tube ]=overall_response();
 plot_2tube(3,0);
end
//
//  -MAIN (3)program starts
//
//==============================================================++=
//
//
//-----------------------------------------------------------------
// global variables   No. 3
//
yg= zeros(1,1024);  // this will be overwritten.
y2tm= zeros(1,1024);  // this will be overwritten.
y3out= zeros(1,1024);  // this will be overwritten.
// other
amplitude= ones(1,2); // this will maybe  be overwritten.
Wy3out_filename='dummy.wav';

//---------------------------------------------------------
function [N_REPEAT]= set_section_qty()
 txt1=['sections quantity'];
 wstr1=sprintf('%f',N_REPEAT);
 sig1=x_mdialog('How many sections to generate',txt1, [wstr1]);
 N_REPEAT=int(evstr(sig1(1)));
endfunction
//---------------------------------------------------------
function [yg,y2tm,y3out]= make_waveform()
 L1z=zeros(1,N_REPEAT);
 L2z=zeros(1,N_REPEAT);
 L3z=zeros(1,N_REPEAT);
 A1z=zeros(1,N_REPEAT);
 A2z=zeros(1,N_REPEAT);
 A3z=zeros(1,N_REPEAT);
 trisez=zeros(1,N_REPEAT);
 tfallz=zeros(1,N_REPEAT);
 tclosedz=zeros(1,N_REPEAT);
 amplitudez=ones(1,N_REPEAT);
 for v=1:N_REPEAT
  L1z(v)=L1;
  L2z(v)=L2;
  L3z(v)=L3;
  A1z(v)=A1;
  A2z(v)=A2;
  A3z(v)=A3;
  trisez(v)=trise;
  tfallz(v)=tfall;
  tclosedz(v)=tclosed;
 end
 //
 if WT_QTY == 3 then   //Three Tubes Model of serial type
 tu1=L1z./c0;   // delay time in 1st tube
 tu2=L2z./c0;   // delay time in 2nd tube
 tu3=L3z./c0;   // delay time in 3rd tube
 r1z=(A2z-A1z)./(A2z+A1z);  // reflection coefficient between 1st tube and 2nd tube
 r2z=(A3z-A2z)./(A3z+A2z);  // reflection coefficient between 2nd tube and 3rd tube
elseif WT_QTY == 4 then   //Three Tubes Model of T type
 tu1=L1z./c0;   // delay time in 1st tube
 tu2=L2z./c0;   // delay time in 2nd tube
 tu3=L3z./c0;   // delay time in 3rd tube
 r12z=((A2z+A3z)-A1z)./(A3z+A2z+A1z);  // reflection coefficient between 1st tube and others
 r21z=(A2z-(A1z+A3z))./(A3z+A2z+A1z);  // reflection coefficient between 2nd tube and others
 r31z=(A3z-(A1z+A2z))./(A3z+A2z+A1z);  // reflection coefficient between 3rd tube and others
else                //Two Tubes Model
 tu1=L1z./c0;   // delay time in 1st tube
 tu2=L2z./c0;   // delay time in 2nd tube
 r1z=(A2z-A1z)./(A2z+A1z);  // reflection coefficient between 1st tube and 2nd tube
end
 //
 disp('+++enter step 1  Input Glottal Volume Velocity Generation');
 D_QTY=N_REPEAT;
 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( tclosedz(loop) / 1000 / ts );
  N2(loop)= int( trisez(loop) / 1000 / ts);
  N3(loop)= int( tfallz(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
 // 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=rg;   // 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
 yg= zeros(1,length0+1);
 rgz= zeros(1,length0+1);
 tc0=1;
 yg(tc0)=0.;
 yg(tc0)=amplitudez(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;
   rgz(tc0)=rg_closed;
 elseif t0 <= (N2(loop) + N1(loop)) then
   yg(tc0)= 0.5 * ( 1 - cos( ( %pi / N2(loop) ) * (t0 - N1(loop)) ) );
   rgz(tc0)=rg_rise;
 elseif t0 <= (N3(loop)+N2(loop)+N1(loop)) then
   yg(tc0)= cos( ( %pi / ( 2 * N3(loop) )) * ( t0 - (N2(loop)+N1(loop)) ) );
   rgz(tc0)=rg_fall; 
 else
  yg(tc0) =0;
  rgz(tc0)=rg_closed;
 end
  yg(tc0)=amplitudez(loop) * yg(tc0);
 end  // t0
 end  // loop
 //
 //
if WT_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. + rgz(tc0) ) / 2.) * yg(tc0) + rgz(tc0) * ya2(M1(loop));
 ya2(1)= -1. * r1z(loop) * ya1(M1(loop))  +  ( 1. - r1z(loop) ) * yb2(M2(loop));
 yb1(1)= ( 1 + r1z(loop) ) * ya1(M1(loop)) + r1z(loop) * yb2(M2(loop));
 yb2(1)= -1. * r2z(loop) * yb1(M2(loop))  +  ( 1. - r2z(loop) ) * yc2(M3(loop));
 yc1(1)= ( 1 + r2z(loop) ) * yb1(M2(loop)) + r2z(loop) * yc2(M3(loop));
 yc2(1)=  -1. * rl  * yc1(M3(loop));
 y2tm(tc0)= (1 + rl) * 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. + rgz(tc0) ) / 2.) * yg(tc0) + rgz(tc0) * ya2(M1(loop));
 ya2(1)= -1. * r1z(loop) * ya1(M1(loop))  +  ( 1. - r1z(loop) ) * yb2(M2(loop));
 yb1(1)= ( 1 + r1z(loop) ) * ya1(M1(loop)) + r1z(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
//
disp('+++enter step 3  Output Sound Pressure for Two Tubes Model for Vocal Tract Calculation');
// At first, low pass filter's coefficients
 wcz=2. * %pi * fc;
 al1= -1. * %e^( -1. * wcz * ts);
 bl0= wcz * ts;
 wc_org=wcz;
 wc_new=2. * %pi * fc;
 alfa=-1. *  cos( (wc_org + wc_new) * ts /2. ) / cos( (wc_org - wc_new) * 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

 disp('--- Waveform Generation finished.');
endfunction     // end ofmake_waveform()
//----------------------------------------------------------
function y3out_snd_play1()
// 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!
sy3out=size(y3out);
y3outz=zeros(sy3out(2));
vm0=abs(y3out(1));
for v=1:sy3out(2)
  if abs(y3out(v)) > vm0 then
    vm0=abs(y3out(v));
  end
end
if vm0 > 0 then
 for v=1:sy3out(2)
  y3outz(v)= y3out(v) / vm0;
 end
end

if f_win == 1 then
if f412 == 1 then    // for windows scilab-4.1.2
 sound(y3outz' ,fs1,bit1);
else                 // for windows scilab-3.1.1
 if fs1 == 22050 then
  sound(y3outz,fs1,bit1);
  disp('This function may work, if luckily.');
 elseif fs1 == 44100 then    // down-sampling from 44100 to 22050
    wdatw2=zeros((sy3out(2)/2));
    for v=1:(sy3out(2)/2)
    wdatw2(v)=y3outz(2 * (v -1) +1);
    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(y3outz' ,fs1,bit1);
 else                 // for scilab-3.1.1
  if f511== 1 then  // for windows scilab-5.1.1
         sound(wdat2',fs1,bit1);
  else
    if fs1 == 22050 then
     sound(y3outz,fs1,bit1);
     disp('This function may work, if luckily.');
    elseif fs1 == 44100 then    // down-sampling from 44100 to 22050
      wdatw2=zeros((sy3out(2)/2));
      for v=1:(sy3out(2)/2)
       wdatw2(v)=y3outz(2 * (v -1) +1);
      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
end
endfunction
//---------------------------------------------------------


function y3out_snd_save1()
//
 wavefilename= input(' + enter file name for saved .wav file =>',["string"]);
//
sy3out=size(y3out);
y3outz=zeros(sy3out(2));
vm0=abs(y3out(1));
for v=1:sy3out(2)
  if abs(y3out(v)) > vm0 then
    vm0=abs(y3out(v));
  end
end
if vm0 > 0 then
 for v=1:sy3out(2)
  y3outz(v)= y3out(v) / vm0;
 end
end
//
if f412 == 1 then    // for windows scilab-4.1.2
 wavwrite(y3outz' ,fs1,bit1,wavefilename );
else                 // for windows scilab-3.1.1
 wavwrite(y3outz, fs1,bit1, wavefilename);
end

endfunction
//---------------------------------------------------------
function plot_waveform(disp0)
 wb0=xget('window');  // stack old window
 xset('window',disp0);   // create new windows
 clf();
 if f412 == 1 then
   plot(yg,'b'); 
   plot(y3out,'r');
   xtitle('generated wave by tubes model','Samples (Time)', 'Amplitude');
  legend('Glottal Volume Velocity','Sound Pressure',2);
 elseif f_win == 1 then
   plot(yg,'b'); 
   plot(y3out,'r');
   xtitle('generated wave by tubes model','Samples (Time)', 'Amplitude');
  legend('Glottal Volume Velocity','Sound Pressure',2);
 else
  ygs=size(yg);
  xzikur=[1:1:ygs(2)];
  //plot2d(xzikur,yg,style=[color("blue")]);
  plot2d(xzikur,y3out,style=[color("red")]);
  xtitle('generated wave by tubes model','Samples (Time)', 'Amplitude');
  //legend('Glottal Volume Velocity','Sound Pressure',2);
 end
 xset('window',wb0);   // push old windows
endfunction
//
//=================================================================
//
//   +MAIN (5)program starts
//    generation waveform
if T_DEMO==1 then
 [yg,y2tm,y3out]= make_waveform();
 plot_waveform(5);
end
//
//  -MAIN (5)program starts
//
//==============================================================++=
//
//
//F_DUMP=0;
//if F_DUMP == 1 then
//// r1m=[-0.9:0.1:0.9];
//// l1m=[-0.8:0.1:0.8];
// r1m=[-0.5:0.5:0.5];
// l1m=[-0.5:0.5:0.5];
// sr1m=size(r1m);
// sl1m=size(l1m);
// //
// for v1=1:sr1m(2)
//  for v2=1:sl1m(2)
//    noz= (r1m(v1)+1.0) * 100.;
//    noz= noz + (l1m(v2)+1.0) * 10000.;
//    wstr1=sprintf('./wavdat/%i.wav',noz);
//    wstr2=sprintf('./wavdat/%i.dat',noz);
//    disp( wstr1);
//    //
//    r1=r1m(v1);
//    l1=l1m(v2);
//    [yg_res1]=yg_resp(1024);
//    [ tube2_res1, L1, L2, L3, A1, A2, A3 ]= two_tubes2_resp();
//    [hpf_res1]=hpf_response(1024);
//    [overall_res1,db_2tube, phi_2tube ]=overall_response();
//    Wtubepropdat= wstr2;
//    save_2tube();
//    //
//    [yg,y2tm,y3out]= make_waveform();
//    Wy3out_filename=wstr1;
//    y3out_snd_save1();
//  end
// end
//end  // if F_DUMP == 1 then
//
//-----------------------------------------------------------------
//=================================================================
//
//  refrence:  db_fft1   ,portion of fft point fft_point1(=default value is 512)
//  tube output: db_2tube    ,portion of fft point 1024
//
//-----------------------------
// global
yg_res9=zeros(2,1);    // dummy set
hpf_res9=zeros(2,1);   // dummy set
wm8=ones(2,1);          // dummy set
wm9=ones(2,1);          // dummy set
tm8=ones(2,1);          // dummy set
x_init8=zeros(2,1);      // dummy set
xopt=zeros(2,1);        // dummy set
ym8=zeros(2,1);         // dummy set

ma0=['ttl_Length ' ; 'l1         ' ; 'r1         ' ; 'l2         ' ; 'r2         ' ; 'rl         '];
ma1=['ttl_Length ' ; 'l1         ' ; 'r1         ' ; 'rl         '];
//----------------------------------------------------------------
function [tm8, ym8, x_init8, wm9]=prepare8()
  [frq1,sfrq1,is1,ie1]=set_frq(fft_point1);
 wm9=ones(sfrq1,1);
 ym8=zeros(sfrq1,1);
 tm8=zeros(sfrq1,1);
 for v=1:sfrq1
  tm8(v)= 2 * %pi * frq1(v);
  ym8(v)=db_fft1(v);
 end
 //.
 if  WT_QTY == 3 then
  // when 3 tube model serial
  x_init8=zeros(6,1);
  x_init8(1,1)=ttl_Length;
  x_init8(2,1)=l1;
  x_init8(3,1)=r1;
  x_init8(4,1)=l2;
  x_init8(5,1)=r2;
  x_init8(6,1)=rl;
 else
  // when 2 tube model
  x_init8=zeros(4,1);
  x_init8(1,1)=ttl_Length;
  x_init8(2,1)=l1;
  x_init8(3,1)=r1;
  x_init8(4,1)=rl;
 end
 disp(' + prepare8 done');
endfunction
//----------------------------------------------------------------
function [yg_res9, hpf_res9]=prepare9()
 yg_res9=yg_resp(fft_point1)';
 hpf_res9=hpf_response(fft_point1)';

 disp(' + prepare9 done'); 
endfunction
//-----------------------------------------------------------------
function [r1,r2,l1,l2,ttl_Length,rl]=exchange8()  // set leastsq result as tube paras
  ttl_Length=xopt(1);
  l1=xopt(2);
  r1=xopt(3);
  rl=xopt(4);
  if  WT_QTY == 3 then
   l2=xopt(4);
   r2=xopt(5);
   rl=xopt(6);
  else
   l2=0;   // dummy set
   r2=0.8;  // dummy set
   rl=0.9;  // dummy set
  end
endfunction
function [r1,r2,l1,l2,ttl_Length,rl]=exchange9()  // set initial value as tube paras
  ttl_Length=x_init8(1);
  l1=x_init8(2);
  r1=x_init8(3);
  rl=x_init8(4);
  if WT_QTY == 3 then
   l2=x_init8(4);
   r2=x_init8(5);
   rl=x_init8(6);
  else
   l2=0;   // dummy set
   r2=0.8;  // dummy set
   rl=0.9;  // dummy set
  end
endfunction
//-----------------------------------------------------------------
//(1)test1
//   (1 + exp( fact0 * x0)) * (1 + exp(fact0 * x1))
//    x0=x - fact1, x1= fact2 -x
//
//   ex: (1 + exp( 50 * (x - 0.9)) * (1 + exp( 50 * (-0.9 -x))
//        When fact0 = 50, fact1 = 0.9, fact2 = -0.9
//        hoping  that -1 < x < 1    (  fact2 ~< x ~< fact1 )
//
limit_switch0=0.;  // set 1 when r1r2,l1,l2 limit is enable, set 0 when limit is disable
fact0=20.;
// for limit of 3 tubes model
fact1_3r1=-0.1;
fact2_3r1=-0.9;
fact1_3r2=0.9;
fact2_3r2=0.1;
fact1_3l1=0.9;
fact2_3l1=-0.9;
fact1_3l2=0.9;
fact2_3l2=-0.9;
// for limit of 2 tubes model
fact1_2r1=0.9;
fact2_2r1=-0.9;
fact1_2l1=0.9;
fact2_2l1=-0.9;
//
//
limit_switch2=1.0;  // set 1 when rl limit is enable, set 0 when limit is disable
// for limit of rl
fact1_rl=1.;
fact2_rl=0.5;
//
//=====================================================
//
//   +MAIN (L)program starts
//
if T_DEMO==1 then
 if SEL_CODE == 1 then // a_sample
   fact1_2r1=0.9;
   fact2_2r1=-0.9;
   fact1_2l1=0.9;
   fact2_2l1=-0.9;
   fact1_rl=1.;
   fact2_rl=0.5;
   limit_switch0=0.0;   // r1,r2,l1,l2 off
   limit_switch2=1.0;   // rl on
 elseif SEL_CODE == 2 then // o_sample
   fact1_3r1=-0.1;
   fact2_3r1=-0.9;
   fact1_3r2=0.9;
   fact2_3r2=0.1;
   fact1_3l1=0.9;
   fact2_3l1=-0.9;
   fact1_3l2=0.9;
   fact2_3l2=-0.9;
   fact1_rl=1.;
   fact2_rl=0.5;
   limit_switch0=1.0;   // r1,r2,l1,l2 on
   limit_switch2=1.0;   // rl on
 elseif SEL_CODE == 3 then // i_sample
   fact1_2r1=-0.1;
   fact2_2r1=-0.9;
   fact1_2l1=0.9;
   fact2_2l1=-0.9;
   fact1_rl=1.;
   fact2_rl=0.5;
   limit_switch0=0.0;   // r1,r2,l1,l2 off
   limit_switch2=1.0;   // rl on
 else   // include SEL_CODE == 999
   // as same as initial values
 end
 //
 // +for only debug
 ////limit_switch0=1.0;   // r1,r2,l1,l2 on
 ////limit_switch2=1.0;   // rl on
 // -for only debug

 

end
//
// -MAIN (L)program starts
//
//
//=====================================================
//
//---------------------------------------------
function [switch0,act0,act1_3l1,act2_3l1,act1_3r1,act2_3r1,act1_3l2,act2_3l2,act1_3r2,act2_3r2,act1_2l1,act2_2l1,act1_2r1,act2_2r1]= edit_limit0()
 txt1=['switch on(1)/off(0)';'exp(a0*x),a0=';'l1 plus';'l1 minus';'r1 plus';'r1 minus';'l2 plus';'l2 minus';'r2 plus';'r2 minus'];
 txt2=['switch on(1)/off(0)';'exp(a0*x),a0=';'l1 plus';'l1 minus';'r1 plus';'r1 minus'];
 wstr0=sprintf('%f',limit_switch0);
 wstr1=sprintf('%f',fact0);
 wstr2=sprintf('%f',fact1_3l1);
 wstr3=sprintf('%f',fact2_3l1);
 wstr4=sprintf('%f',fact1_3r1);
 wstr5=sprintf('%f',fact2_3r1);
 wstr6=sprintf('%f',fact1_3l2);
 wstr7=sprintf('%f',fact2_3l2);
 wstr8=sprintf('%f',fact1_3r2);
 wstr9=sprintf('%f',fact2_3r2);
 wstr10=sprintf('%f',fact1_2l1);
 wstr11=sprintf('%f',fact2_2l1);
 wstr12=sprintf('%f',fact1_2r1);
 wstr13=sprintf('%f',fact2_2r1);

 if WT_QTY == 3 then
   sig1=x_mdialog('Set limit',txt1, [wstr0 ; wstr1 ; wstr2  ; wstr3 ; wstr4 ; wstr5 ; wstr6 ; wstr7 ; wstr8 ; wstr9 ]);
  if sig1==[] then
   switch0=evstr(wstr0);
   act0=evstr(wstr1);
   act1_3l1=evstr(wstr2);
   act2_3l1=evstr(wstr3);
   act1_3r1=evstr(wstr4);
   act2_3r1=evstr(wstr5);
   act1_3l2=evstr(wstr6);
   act2_3l2=evstr(wstr7);
   act1_3r2=evstr(wstr8);
   act2_3r2=evstr(wstr9);
   act1_2l1=fact1_2l1;
   act2_2l1=fact2_2l1;
   act1_2r1=fact1_2r1;
   act2_2r1=fact2_2r1;
  else
   switch0=evstr(sig1(1));
   act0=evstr(sig1(2));
   act1_3l1=evstr(sig1(3));
   act2_3l1=evstr(sig1(4));
   act1_3r1=evstr(sig1(5));
   act2_3r1=evstr(sig1(6));
   act1_3l2=evstr(sig1(7));
   act2_3l2=evstr(sig1(8));
   act1_3r2=evstr(sig1(9));
   act2_3r2=evstr(sig1(10));
   act1_2l1=fact1_2l1;
   act2_2l1=fact2_2l1;
   act1_2r1=fact1_2r1;
   act2_2r1=fact2_2r1;
  end
 else  // if WT_QTY == 2 then
   sig1=x_mdialog('Set limit',txt2, [wstr0 ; wstr1 ; wstr10  ; wstr11 ; wstr12 ; wstr13 ]);
  if sig1==[] then
   switch0=evstr(wstr0);
   act0=evstr(wstr1);
   act1_2l1=evstr(wstr10);
   act2_2l1=evstr(wstr11);
   act1_2r1=evstr(wstr12);
   act2_2r1=evstr(wstr13);
   act1_3l1=fact1_3l1;
   act2_3l1=fact2_3l1;
   act1_3r1=fact1_3r1;
   act2_3r1=fact2_3r1;
   act1_3l2=fact1_3l2;
   act2_3l2=fact2_3l2;
   act1_3r2=fact1_3r2;
   act2_3r2=fact2_3r2;
  else
   switch0=evstr(sig1(1));
   act0=evstr(sig1(2));
   act1_2l1=evstr(sig1(3));
   act2_2l1=evstr(sig1(4));
   act1_2r1=evstr(sig1(5));
   act2_2r1=evstr(sig1(6));
   act1_3l1=fact1_3l1;
   act2_3l1=fact2_3l1;
   act1_3r1=fact1_3r1;
   act2_3r1=fact2_3r1;
   act1_3l2=fact1_3l2;
   act2_3l2=fact2_3l2;
   act1_3r2=fact1_3r2;
   act2_3r2=fact2_3r2;
  end
 end
endfunction
//-----------------------------------------------------------------
function [switch2,act1_rl,act2_rl]= edit_limit2()
 txt1=['switch on(1)/off(0)';'rl plus';'rl minus'];
 wstr0=sprintf('%f',limit_switch2);
 wstr2=sprintf('%f',fact1_rl);
 wstr3=sprintf('%f',fact2_rl);
 sig1=x_mdialog('Set limit',txt1, [wstr0 ; wstr2 ; wstr3]);
 if sig1==[] then
   switch2=evstr(wstr0);
   act1_rl=evstr(wstr2);
   act2_rl=evstr(wstr3);
 else
   switch2=evstr(sig1(1));
   act1_rl=evstr(sig1(2));
   act2_rl=evstr(sig1(3));
 end
endfunction
//-----------------------------------------------------------------
function test_plot1exp(disp0)

 // limit_switch0=0.0;   // r1,r2,l1,l2 off
 if limit_switch0 >= 1. then
   disp(' r1,r2,l1,l2 limit control is ON');
 else
   disp(' r1,r2,l1,l2 limit control is OFF');
 end
 // limit_switch2=1.0;   // rl on
 if limit_switch2 >= 1. then
   disp(' rl limit control is ON');
 else
   disp(' rl limit control is OFF');
 end

 wb0=xget('window');  // stack old window
 xset('window',disp0);   // create new windows
 // dx0=[-1:0.1:1];
 clf();
 //subplot(311);
 //dx0=[-1:0.1:(fact1_3r1+0.2)];
 //plot2d( dx0,(1 + exp( fact0 * ( dx0 - fact1_3r1))),style=[color("red")]  );
 //subplot(312);
 //dx0=[(fact2_3r1-0.2):0.1:1];
 //plot2d( dx0,(1 + exp( fact0  * ( fact2_3r1 - dx0))),style=[color("blue")]  );
 //subplot(313);

if WT_QTY == 3 then
 if limit_switch0 >= 1 then
 subplot(511);
 dx0=[(fact2_3l1-0.2):0.1:(fact1_3l1+0.2)];
 plot2d( dx0,(1 + exp( fact0 * ( dx0 - fact1_3l1))) .* (1 + exp( fact0  * ( fact2_3l1 - dx0))),style=[color("red")]  );
 wstr0='decrease value depend on l1, hoping l1 inside the limit ';
 xtitle(wstr0);

 subplot(512);
 dx0=[(fact2_3r1-0.2):0.1:(fact1_3r1+0.2)];
 plot2d( dx0,(1 + exp( fact0 * ( dx0 - fact1_3r1))) .* (1 + exp( fact0  * ( fact2_3r1 - dx0))),style=[color("red")]  );
 wstr0='decrease value depend on r1, hoping r1 inside the limit ';
 xtitle(wstr0);

 subplot(513);
 dx0=[(fact2_3l2-0.2):0.1:(fact1_3l2+0.2)];
 plot2d( dx0,(1 + exp( fact0 * ( dx0 - fact1_3l2))) .* (1 + exp( fact0  * ( fact2_3l2 - dx0))) ,style=[color("red")] );
 wstr0='decrease value depend on l2, hoping l2 inside the limit ';
 xtitle(wstr0);

 subplot(514);
 dx0=[(fact2_3r2-0.2):0.1:(fact1_3r2+0.2)];
 plot2d( dx0,(1 + exp( fact0 * ( dx0 - fact1_3r2))) .* (1 + exp( fact0  * ( fact2_3r2 - dx0))),style=[color("red")]  );
 wstr0='decrease value depend on r2, hoping r2 inside the limit ';
 xtitle(wstr0);

 end // if limit_switch0 >= 1 then

 if limit_switch2 >= 1 then
 subplot(515);
 dx0=[(fact2_rl-0.2):0.1:(fact1_rl+0.2)];
 plot2d( dx0,(1 + exp( fact0 * ( dx0 - fact1_rl))) .* (1 + exp( fact0  * ( fact2_rl - dx0))),style=[color("red")]  );
 wstr0='decrease value depend on rl, hoping rl inside the limit ';
 xtitle(wstr0);
 end // if limit_switch2 >= 1 then

else   //if WT_QTY == 2 then
 if limit_switch0 >= 1 then
 subplot(311);
 dx0=[(fact2_2l1-0.2):0.1:(fact1_2l1+0.2)];
 plot2d( dx0,(1 + exp( fact0 * ( dx0 - fact1_2l1))) .* (1 + exp( fact0  * ( fact2_2l1 - dx0))),style=[color("red")]  );
 wstr0='decrease value depend on l1, hoping l1 inside the limit ';
 xtitle(wstr0);

 subplot(312);
 dx0=[(fact2_2r1-0.2):0.1:(fact1_2r1+0.2)];
 plot2d( dx0,(1 + exp( fact0 * ( dx0 - fact1_2r1))) .* (1 + exp( fact0  * ( fact2_2r1 - dx0))),style=[color("red")]  );
 wstr0='decrease value depend on r1, hoping r1 inside the limit ';
 xtitle(wstr0);
 end //if limit_switch0 >= 1 then

 if limit_switch2 >= 1 then
 subplot(313);
 dx0=[(fact2_rl-0.2):0.1:(fact1_rl+0.2)];
 plot2d( dx0,(1 + exp( fact0 * ( dx0 - fact1_rl))) .* (1 + exp( fact0  * ( fact2_rl - dx0))),style=[color("red")]  );
 wstr0='decrease value depend on rl, hoping rl inside the limit ';
 xtitle(wstr0);
 end //if limit_switch2 >= 1 then

end

 xset('window',wb0);   // push old windows
endfunction
//
//
//-----------------------------------------------------------------
// y = tube2_resp9(tm8, x_init8)
//
// [frq1,sfrq1,is1,ie1]=set_frq(fft_point1);
// clf();
// subplot(211);
// gainplot(frq1,y,phi_fft1);
//
function y = tube2_resp9(xw, pa)  // when 2 tube model
   L2= ( 1. + pa(2)) * pa(1) / 2.;   // L2= ( 1. + l1 ) * ttl_Length / 2.;
   L1= pa(1) - L2;                   // L1= ttl_Length - L2;
   tu1=L1/c0;
   tu2=L2/c0;
   A2= ( 1. + pa(3)) * ttl_Area/ 2.; // A2= ( 1. + r1 ) * ttl_Area/ 2.;
   A1= ttl_Area - A2;
   yi  = 0.5 * ( 1 + rg ) * ( 1 + pa(3))  * ( 1 + pa(4) ) * %e^( -1 * ( tu1 + tu2 ) .* xw * %i);
   yb =  1 + pa(3) * rg * %e^( -2 * tu1 .* xw * %i) + pa(3) * pa(4) * %e^( -2 * tu2 .* xw * %i) + pa(4) * rg * %e^( -2 * (tu1 + tu2) .* xw * %i);
   yc =  yg_res9 .* hpf_res9 .* ( yi ./ yb) ;
   y = 20. * log10(abs(yc));

   y= y - limit_switch0 * ((1 + exp( fact0 * ( pa(3) - fact1_2r1))) *  (1 + exp( fact0  * ( fact2_2r1 - pa(3) )   )  )) * wm9;
   y= y - limit_switch0 * ((1 + exp( fact0 * ( pa(2) - fact1_2l1))) *  (1 + exp( fact0  * ( fact2_2l1 - pa(2) )   )  )) * wm9;
   y= y - limit_switch2 * ((1 + exp( fact0 * ( pa(4) - fact1_rl))) *  (1 + exp( fact0  * ( fact2_rl - pa(4) )   )  )) * wm9;



endfunction
//--------------------------------
// y = tube3_resp9(tm8, x_init8)
//
// [frq1,sfrq1,is1,ie1]=set_frq(fft_point1);
// clf();
// subplot(211);
// gainplot(frq1,y,phi_fft1);
function y = tube3_resp9(xw, pa)  // when 3 tube model serial
  bunbo= pa(2) * pa(4) + ( pa(2) - pa(4) ) + 3.0;  // bunbo= l1 * l2 + ( l1 - l2 ) + 3.0;
  L1= pa(1) * ( pa(4) - 1.0 ) * (pa(2) - 1.0 ) / bunbo ; // L1= ttl_Length * ( l2 - 1.0 ) * (l1 - 1.0 ) / bunbo ;
  L2= pa(1) * ( pa(4) - 1.0 ) * ( -1.0 * pa(2) - 1.0 ) / bunbo ; // L2= ttl_Length * ( l2 - 1.0 ) * ( -1.0 * l1 - 1.0 ) / bunbo ;
  L3= pa(1) * ( pa(4) + 1.0 ) * (pa(2) + 1.0 ) / bunbo ; // L3= ttl_Length * ( l2 + 1.0 ) * (l1 + 1.0 ) / bunbo ;
  tu1=L1/c0;
  tu2=L2/c0;
  tu3=L3/c0;
  bunbo= pa(3) * pa(5) + ( pa(3) - pa(5) ) + 3.0; // bunbo= r1 * r2 + ( r1 - r2 ) + 3.0;
  A1= ttl_Area * ( pa(5) - 1.0 ) * (pa(3) - 1.0 ) / bunbo ; // A1= ttl_Area * ( r2 - 1.0 ) * (r1 - 1.0 ) / bunbo ;
  A2= ttl_Area * ( pa(5) - 1.0 ) * ( -1.0 * pa(3) - 1.0 ) / bunbo ; // A2= ttl_Area * ( r2 - 1.0 ) * ( -1.0 * r1 - 1.0 ) / bunbo ;
  A3= ttl_Area * ( pa(5) + 1.0 ) * (pa(3) + 1.0 ) / bunbo ; // A3= ttl_Area * ( r2 + 1.0 ) * (r1 + 1.0 ) / bunbo ;

   yi  = 0.5 * ( 1. + rg ) * ( 1. + pa(3)) * ( 1. + pa(5)) * ( 1. + pa(6) ) * %e^( -1. * ( tu1 + tu2  + tu3) .* xw * %i);
   yb =  1 + rg * pa(3) * %e^( -2 * tu1 .* xw * %i) + pa(3) * pa(5) * %e^( -2 * tu2 .* xw * %i) + pa(5) * pa(6) * %e^( -2 * tu3 .* xw * %i);
   yb = yb + rg * pa(5) * %e^( -2 * (tu1 + tu2) .* xw * %i) + pa(3) * pa(6) * %e^( -2 * (tu2 + tu3) .* xw * %i);
   yb = yb + rg * pa(3) * pa(5) * pa(6) *  %e^( -2 * (tu1 + tu3) .* xw * %i);
   yb = yb + rg * pa(6) *  %e^( -2 * (tu1 + tu2 + tu3) .* xw * %i);
   yc =  yg_res9 .* hpf_res9 .* ( yi ./ yb) ;

    //(1)test1
   // for r1
   // yc = ( 1. / ((1 + exp( fact0 * ( pa(3) - fact1_3r1))) * (1 + exp( fact0  * ( fact2_3r1 - pa(3)))))   ) * yc;

   y = 20. * log10(abs(yc));
   y= y - limit_switch0 * ((1 + exp( fact0 * ( pa(3) - fact1_3r1))) *  (1 + exp( fact0  * ( fact2_3r1 - pa(3) )   )  )) * wm9;
   y= y - limit_switch0 * ((1 + exp( fact0 * ( pa(5) - fact1_3r2))) *  (1 + exp( fact0  * ( fact2_3r2 - pa(5) )   )  )) * wm9;
   y= y - limit_switch0 * ((1 + exp( fact0 * ( pa(2) - fact1_3l1))) *  (1 + exp( fact0  * ( fact2_3l1 - pa(2) )   )  )) * wm9;
   y= y - limit_switch0 * ((1 + exp( fact0 * ( pa(4) - fact1_3l2))) *  (1 + exp( fact0  * ( fact2_3l2 - pa(4) )   )  )) * wm9;
   y= y - limit_switch2 * ((1 + exp( fact0 * ( pa(6) - fact1_rl))) *  (1 + exp( fact0  * ( fact2_rl - pa(6) )   )  )) * wm9;


endfunction
//-----------------------------------------------------------------
//  e=myfun(x_init8, tm8, ym8, wm8) 
function e = myfun(x, tm, ym, wm)   // when 2 tube model
e = wm .*( tube2_resp9(tm, x) - ym )
endfunction
function e = myfun2(x, tm, ym, wm)   // when 3 tube model serial
e = wm .*( tube3_resp9(tm, x) - ym )
endfunction
//-----------------------------------------------------------------
function [fopt,xopt, grdopt]=leastsq_main1()
  //
  // check weight
  a=0.;
  s0=size(wm8);
  wm8c=ones(s0(1),1);
  for v=1:s0(1)
    a=a+abs( wm8(v));
  end
  if a == 0. then
   disp('+WARNING: set weight as all 1, because wm8 are all 0');
  else
   for v=1:s0(1)
     wm8c(v)=wm8(v);
   end
  end


  // 1- the simplest call
  if  WT_QTY == 3 then
   // when 3 tube model serial
   [fopt,xopt, grdopt] = leastsq(list(myfun2,tm8,ym8,wm8c),x_init8);
  else
   // when 2 tube model
   [fopt,xopt, grdopt] = leastsq(list(myfun,tm8,ym8,wm8c),x_init8);
  end
  //... call
  s0=size(xopt);
  s0s=s0(1);
  disp( 'initial value -> result value');
  for v=1:s0s
    wstr1=sprintf('%f',xopt(v));
    wstr2=sprintf('%f',x_init8(v));
    if  WT_QTY == 3 then
       wstr3='  ' + ma0(v) + wstr2 + ' -> ' + wstr1;
    else
       wstr3='  ' + ma1(v) + wstr2 + ' -> ' + wstr1;
    end
    disp(wstr3);
  end
  disp(' ');

endfunction
//-----------------------+++ ---- ++ ------------------------------
function [wm8]=reset_wm8()
  [frq1,sfrq1,is1,ie1]=set_frq(fft_point1);
  wm8=zeros(sfrq1,1);
  // draw all 1
  disp('+set weight all 0');
endfunction
//----------------
function plot_wm8(disp0)
 wb0=xget('window');  // stack old window
 xset('window',disp0);   // create new windows
 [frq1,sfrq1,is1,ie1]=set_frq(0);
 clf();
 subplot(211);
 gainplot(frq1,db_fft1,phi_fft1);
 xtitle('frequency response of waveform selected by fft analysis');
  [frq1,sfrq1,is1,ie1]=set_frq(fft_point1);
 subplot(212);
 plot2d(frq1,wm8',style=[color("red")],logflag="ln",rect=[frq1(1),-1,frq1(sfrq1),2]);
 wstr0=' wmf: weight for leastsq evalution ';
 xtitle(wstr0);
 xset('window',wb0);   // push old windows
endfunction
//----------------------------------------------------------------
function [wm8c]= set_wm8( arg1, arg2, arg3)
 [frq1,sfrq1,is1,ie1]=set_frq(fft_point1);
 wm8c=zeros(sfrq1,1);
 for v=1:sfrq1
  wm8c(v)=wm8(v);
 end
 for v=1:sfrq1
   if (frq1(v) >= arg1) & (frq1(v) <= arg2) then
     wm8c(v)=arg3;
   end
 end
endfunction
//---------------------------------------------
function [wm8c]= edit_wm8()
 txt1=['from (frequency)';'to (frequency)';'weight value'];
 sp1=0.;
 ep1=0.;
 vl0=1.;
 wstr1=sprintf('%f',sp1);
 wstr2=sprintf('%f',ep1);
 wstr3=sprintf('%f',vl0);
 sig1=x_mdialog('Set weight value',txt1, [wstr1 ; wstr2  ; wstr3]);
 if sig1==[] then
  arg1=evstr(wstr1);
  arg2=evstr(wstr2);
  arg3=evstr(wstr3);
 else
  arg1=evstr(sig1(1));
  arg2=evstr(sig1(2));
  arg3=evstr(sig1(3));
 end

 [wm8c]= set_wm8( arg1, arg2, arg3);

endfunction
//
//=================================================================
//
//   +MAIN (4)program starts
//    leastsq method to estimate tube model parameters
if T_DEMO==1 then
 
  disp(' ');
  disp('a trial of leastsq method to estimate tube model parameters including rl  using focuus weight.');
  disp(' ');

  // once set wm8 as all ones
  [wm8]=reset_wm8();

  if SEL_CODE == 1 then // a_sample
    [wm8]= set_wm8( 516.79688 ,1808.7891, 1.0 );
    [wm8]= set_wm8( 2153.3203 ,3703.7109, 1.0);
    plot_wm8(6);
  elseif SEL_CODE == 2 then // o_sample
    [wm8]= set_wm8( 172.26562, 1200. , 1.0 );
    [wm8]= set_wm8( 1894.9219, 3703.7109, 1.0);
    plot_wm8(6);
  elseif SEL_CODE == 3 then // i_sample
    [wm8]= set_wm8( 200., 5000. , 1.0 );
    plot_wm8(6);
  else   // include SEL_CODE == 999
    plot_wm8(6);
    [wm8]= edit_wm8();
    plot_wm8(6);
    [wm8]= edit_wm8();
    plot_wm8(6);
  end

  test_plot1exp(11);

  [tm8, ym8, x_init8, wm9]=prepare8();
  [yg_res9, hpf_res9]=prepare9();
  [fopt,xopt, grdopt]=leastsq_main1();

  // set result value
  [r1,r2,l1,l2,ttl_Length,rl]=exchange8();
  [yg_res1]=yg_resp(1024);
  [ tube2_res1, L1, L2, L3, A1, A2, A3 ]= two_tubes2_resp();
  [hpf_res1]=hpf_response(1024);
  [overall_res1,db_2tube, phi_2tube ]=overall_response();
  plot_2tube(4,1);
  // back to initial value
  [r1,r2,l1,l2,ttl_Length,rl]=exchange9();
  [yg_res1]=yg_resp(1024);
  [ tube2_res1, L1, L2, L3, A1, A2, A3 ]= two_tubes2_resp();
  [hpf_res1]=hpf_response(1024);
  [overall_res1,db_2tube, phi_2tube ]=overall_response();

end
//
//  -MAIN (4)program starts
//
//==============================================================++=
//
//
//
//
//=======$======&==========================#=====================
//
//  Add menu buttons of which name are 'a_plotwav' and ''
//
ADD_MENU_PLOT=1;       // set ADD_MENU 1 to add menu button of some functions
//
//
if ADD_MENU_PLOT == 1 then
delmenu('step1_plotwav');
addmenu('step1_plotwav',[ '(0)read_wav' ; '(1)set_et_plot()'; 'reset_plot()'; 'load_wavprop()'; 'save_wavprop()']);
step1_plotwav=[ '[wavfile1,chs1,qty1,fs1,bit1,f412,SEL_CODE]=get_wavfile_pro(); disp_pro_wav(); [wdat1, size1]=read_one_ch_of_wav(wavfile1); sp1=1; ep1=size1; [ydat0,xziku0]=make_width_data( wdat1, size1);' ; '[sp1, ep1]= set_sp1_ep1(size1); [ydat1,xziku1]=make_width_data( wdat1, size1);  plot_wave1(0);' ; '[sp1, ep1]= reset_sp1_ep1(size1); [ydat1,xziku1]=make_width_data( wdat1, size1);' ; '[fs1,size1,bit1,wdat1,sp1,ep1,ydat0,ydat1,xziku0,xziku1]=load_wavprop(); plot_wave1(0)' ; 'save_wavprop()'] ;
end //if ADD_MENU_PLOT == 1 then
//
//-----------------------------------------------------------------
//
ADD_MENU_FFT=1;       // set ADD_MENU 1 to add menu button of some functions
//
//
if ADD_MENU_FFT == 1 then
delmenu('step2_fft');
addmenu('step2_fft',[ '(1)set points'; '(2)fft cal'; '(3)plot smooth, peak candinate' ; 'load_wfft()'; 'save_wfft()']);
step2_fft=[ '[fft_point1]=set_fft_points()' ; '[db_fft1, phi_fft1]=do_fft_wav(); plot_fft1(1)'; '[sn,swnd,sm1out,sm2out,npk,pklist]=smooth1(1); plot_fft_sm1(2)'  ;'[fft_point1,db_fft1,phi_fft1]=load_wfft(); plot_fft1(1)' ; 'save_wfft()'] ;
end //if ADD_MENU_FFT == 1 then
//
//-----------------------------------------------------------------
//
ADD_MENU_TUBE2=1;       // set ADD_MENU 1 to add menu button of some functions
//
//
if ADD_MENU_TUBE2 == 1 then
delmenu('step3_tubes');
addmenu('step3_tubes',[ '(1)set_tube_model()'; '(2)set_tube_initial_para()'; '(3)plot_tube_freq(3)'; 'plot_area(10)'; 'load_tube()'; 'save_tube()'; 'plot_tube_freq_l1_r1']);
step3_tubes=[ '[WT_QTY]= set_tube_model()' ; '[r1,r2,l1,l2,ttl_Length,rl,fc,trise,tfall,tclosed]= set_2tube_para()'; '[yg_res1]=yg_resp(1024); [ tube2_res1, L1, L2, L3, A1, A2, A3 ]= two_tubes2_resp(); [hpf_res1]=hpf_response(1024); [overall_res1,db_2tube, phi_2tube ]=overall_response(); plot_2tube(3,0);' ; 'plot_area(10)'; '[L1,L2,L3,A1,A2,A3,ttl_Length,ttl_Area,l1,l2,rg,r1,r2,fc,trise,tfall,tclosed,yg_res1,tube2_res1,hpf_res1,overall_res1,db_2tube,phi_2tube,WT_QTY]=load_2tube(); plot_2tube(3,0);' ; 'save_2tube();' ; ' plot_2tube(3,3)'] ;
end //if ADD_MENU_TUBE2 == 1 then
//-----------------------------------------------------------------
//
ADD_MENU_LEASTSQURE=1;       // set ADD_MENU 1 to add menu button of some functions
//
//
if ADD_MENU_LEASTSQURE == 1 then
delmenu('step4_leastsq');
addmenu('step4_leastsq',[ '(0)set_weight_all_0' ; '(1)edit_plot_weight' ; '(2)edit_limit' ; '(3)do_leastsq' ; '(4)plot_freq(4)' ; '(4)set_result_as_initial_para' ; 'plot_1+exp_ 1+exp(11)']);
step4_leastsq=['[wm8]=reset_wm8();' ; 'plot_wm8(6); [wm8]= edit_wm8(); plot_wm8(6);'; '[limit_switch0,fact0,fact1_3l1,fact2_3l1,fact1_3r1,fact2_3r1,fact1_3l2,fact2_3l2,fact1_3r2,fact2_3r2,fact1_2l1,fact2_2l1,fact1_2r1,fact2_2r1]= edit_limit0();  [limit_switch2,fact1_rl,fact2_rl]= edit_limit2(); test_plot1exp(11);'    ;' [tm8, ym8, x_init8, wm9 ]=prepare8(); [yg_res9, hpf_res9]=prepare9(); [fopt,xopt, grdopt]=leastsq_main1();' ; '[r1,r2,l1,l2,ttl_Length,rl]=exchange8(); [yg_res1]=yg_resp(1024); [ tube2_res1, L1, L2, L3, A1, A2, A3 ]= two_tubes2_resp(); [hpf_res1]=hpf_response(1024); [overall_res1,db_2tube, phi_2tube ]=overall_response(); plot_2tube(4,1); [r1,r2,l1,l2,ttl_Length,rl]=exchange9(); [yg_res1]=yg_resp(1024); [ tube2_res1, L1, L2, L3, A1, A2, A3 ]= two_tubes2_resp(); [hpf_res1]=hpf_response(1024); [overall_res1,db_2tube, phi_2tube ]=overall_response();';'[r1,r2,l1,l2,ttl_Length,rl]=exchange8();' ; 'test_plot1exp(11);'];
end
//
//-----------------------------------------------------------------
//
ADD_MENU_GEN=1;       // set ADD_MENU 1 to add menu button of some functions
//
//
if ADD_MENU_FFT == 1 then
delmenu('step5_generation');
if (f_win == 1) | (scilab_version_number >= 4) then
 addmenu('step5_generation',[ '(1)set_section_qty()' ; '(2)make_waveform()' ; 'plot_waveform()'; 'y3out_snd_play1()'; 'y3out_snd_save1()' ]);
 step5_generation=[ '[N_REPEAT]= set_section_qty()'; '[yg,y2tm,y3out]= make_waveform();' ; 'plot_waveform(5)'; 'y3out_snd_play1()'; 'y3out_snd_save1()'] ;
else
 addmenu('step5_generation',[ '(1)set_section_qty()' ; '(2)make_waveform()' ; 'plot_waveform()';  'y3out_snd_save1()' ]);
 step5_generation=[ '[N_REPEAT]= set_section_qty()'; '[yg,y2tm,y3out]= make_waveform();' ; 'plot_waveform(5)';  'y3out_snd_save1()'] ;
end
end //if ADD_MENU_GEN == 1 then
//
//
//-----------------------------------------------------------------
//
if f_win == 1 then
 ADD_MENU_SND=0; 
else
 ADD_MENU_SND=0;       // set ADD_MENU 1 to add menu button of some functions
end
//
//
if ADD_MENU_SND == 1 then
delmenu('sound_play');
addmenu('sound_play',[ 'snd_play1()' ; 'snd_save1()' ]);
sound_play=[ 'snd_play1()' ; 'snd_save1()' ] ;
end //if ADD_MENU_SND == 1 then
//
//==============================================================++=
//
//////addmenu('&go');   // only windows shortcut can  alt+()



No.2   2009年７月18日