//====================================================================

//      scilab version windows scilab-4.1.2 is used.

//      this is still draft, out of guarantee !

//-------------------------------------------------------------------------------------------

// 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 should be done by your own risk.

// Thanks.

//====================================================================

// waveform data read

//exec("wavedata0.sci",-1);      // old 1 line sample

exec("wavedata2.sci",-1);    // 2 lines sample

//

function [x]=readx()

 x=input("#input *.sci  ","string");

 disp(x);

endfunction

//--------------------------------------------------------------------

fs=fs0;  // sampling frequency

dt0=1./fs;

dph0= %pi / 360.;

maxlevel=32768.;

// initial parameter sample

wx0 = [  1.190476E-02 ;  2.022400E+04 ;  1.495997E-01 ;  0.000000E+00 ;  1.190476E-02 ;  1.011200E+04 ;  1.795196E-01 ;  5.235988E-01 ;  4.754286E+02];

//====================================================================

// global variables 

sp1=1;

ep1=m0;

m=1;

ym=zeros(1,1);

tm=zeros(1,1);

wm=ones(1,1);

f=0.;

xopt=zeros(1,1);

xopt2=zeros(1,1);

gropt=zeros(1,1);

x0=zeros(1,1);

binf0=zeros(1,1);

bsup0=zeros(1,1);

//------ type definition ----------

type0=1;   // when type0=1, expand

           // when type0=2, reduce

           // when type0=3, transit 

//---------------------------------

//           portion patterns, if only 2 notoki

ptype0=1;   // when ptype0=1, expand->reduce

            // when ptype0=2, reduce->expand

//           portion patterns, 3 notoki, 3C3

//

//---------------------------------

kcands=zeros(1,1);

scan_xopt=zeros(1,1);

scan_f=zeros(1,1);

scan_type=zeros(1,1);

//====================================================================

MAX_LOOP_COUNTER0=20;

MAX_LOOP_COUNTERB0=20;  

FITABLE_LEVEL0=75.;

ONE_STEP0=10;

MINIMUM_LEGTH0=2;  // MINIMUM LEGTH is MINIMUM_LEGTH0 x sizeof(wx0)

//====================================================================

// scilab global variables

global eh_var1    // qty of graphic windows

eh_var1=0;

//--------------------------------------------------------------------

// estimated function

function y = yth(t, x)

y  = exp(x(1)*t). * ( x(2) * sin(  x(3) * t +  x(4) ) )  +  exp(x(5)*t). * ( x(6) * sin(  (x(3)+x(7)) * t +  (x(4)+x(8)) ) )+ x(9)

endfunction

function e = myfun0(x, tm, ym, wm)

e = wm.*( yth(tm, x) - ym )

endfunction

function [f,xopt, gropt] = simplecall1()

// 1- the simplest call

 [f,xopt, gropt] = leastsq(list(myfun0,tm,ym,wm),x0)

//... call plot 

 xset('window',0); 

 clf();

 plot(tm,ym,'+',tm,yth(tm,xopt),'-',tm,(ym-yth(tm,xopt)),'--r')

endfunction

function [f,xopt, gropt] = simplecall1b(tm,ym,wm,x0,disp0)

// 1- the simplest call

 [f,xopt, gropt] = leastsq(list(myfun0,tm,ym,wm),x0)

//... call plot 

 if disp0 >= 0 then

  xset('window',disp0); 

  clf();

  plot(tm,ym,'+',tm,yth(tm,xopt),'-',tm,(ym-yth(tm,xopt)),'--r')

  //plot(tm,ym,'+',tm,yth(tm,xopt),'-')

 end

endfunction

function [ gosa1, gosa1abs ] = gosa(tm,ym,xopt)

  gosa1=ym-yth(tm,xopt);

  gosa1abs=abs(gosa1);

endfunction

//

//=========================================================================================

// This needs more more time than leastsq

// no bounds no kaidemo, bounds zi ha detekonaide, husei na kai ni naru

function f = myfun0b(x, varagin)

 e = wm.*( yth(tm, x) - ym )

 f =  norm(e)^2

endfunction

function [f,xopt] = simplecall2(binf,bsup)

 //[f,xopt]=optim(list(NDcost,myfun0b,0) ,x0)      //Simplest call

 [f,xopt]=optim(list(NDcost,myfun0b,0), 'b',binf,bsup ,x0)    //bounds on x0 call

//... call plot 

 xset('window',1); 

 clf();

 plot(tm,ym,'+',tm,yth(tm,xopt),'-')

endfunction

//===================================================================

function [wsp1, wep1, m ,tm, ym, wm, wtype]= set_sp1_ep1(wsize1,wtm,wym,disp0)

 wb0=xget('window');  // stack old window

 xset('window',disp0);   // create new windows

 clf();

 subplot(211);

 plot(wtm,wym);

 txt1=['start point(sec.1)';'end point(sec.1)';'type'];

 wstr1=sprintf('%d',sp1);

 wstr2=sprintf('%d',ep1);

 wstr3=sprintf('%d',type0);

 sig1=x_mdialog('set start point, end point, and type',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

 wtype=arg3;

 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

 [m,tm,ym,wm]= set_mtmymwm(wsp1,wep1,wym);

 subplot(212);

 plot(tm,ym);

 xset('window',wb0);   // push old windows

endfunction

//--------------------------------------------------------------------

function [m,tm,ym,wm]= set_mtmymwm(wsp1,wep1,wym)

 m=wep1-wsp1+1;

 ci0=1;

 cr0=0.0;

 tm=ones(m,1);

 ym=ones(m,1);

 for v=wsp1:wep1

  tm(ci0)= cr0;

  ym(ci0)= wym(v);

  ci0=ci0+1;

  cr0=cr0+1.;

 end

 wm=ones(m,1);

endfunction

//-------------------------------------------------------------------

function [wmean, kmax, wmax, durmax, kmin, wmin, durmin, maxabs, maxpoint, kyokusei ]=detectp(wym)

 // cal mean

 s0=size(wym);

 s1=s0(1);

 wmean=0.;

 maxpoint=1;

 maxabs=abs( wym(1));

 for v=1:s1

  wmean = wmean + wym(v);

  if abs(wym(v)) > maxabs then

   maxpoint=v;

   maxabs=abs(wym(v));

  end

 end

 wmean= wmean / s1;

 durmax=s1;

 durmin=s1;

 // call diff

 wmax=zeros(1,1);

 wmin=zeros(1,1);

 kmax=0;

 kmin=0;

 dym=diff(wym);

 s0=size(dym);

 s1=s0(1);

 // detect zero-cross points

 for v=2:s1

    if ( ( dym(v-1) > 0.) & ( dym(v) <= 0.) ) then

       kmax=kmax+1;

       wmax(kmax)=v;

    elseif ( ( dym(v-1) < 0.) & ( dym(v) >= 0.) ) then

       kmin=kmin+1;

       wmin(kmin)=v;

    end

 end

 // cal peak-peak duration mean

 if kmax >= 2 then

   durmax=0.;

   for v=2:kmax

     durmax=durmax + ( wmax(v) - wmax(v-1));

   end

   durmax=durmax / (kmax-1);

 end

 if kmin >= 2 then

   durmin=0.;

   for v=2:kmin

     durmin=durmin + ( wmin(v) - wmin(v-1));

   end

   durmin=durmin / (kmin-1);

 end

 //

 if dym(1) >= 0. then

  kyokusei =0.;

 else

  kyokusei =1.;

 end

endfunction

//--------------------------------------------------------------------

function [wmean, kmax, wmax, durmax, kmin, wmin, durmin, maxabs, maxpoint, kyokusei ]=detectzc(wym)

// detect zero cross of wym

// cal mean

 s0=size(wym);

 s1=s0(1);

 wmean=0.;

 maxpoint=1;

 maxabs=abs( wym(1));

 for v=1:s1

  wmean = wmean + wym(v);

  if abs(wym(v)) > maxabs then

   maxpoint=v;

   maxabs=abs(wym(v));

  end

 end

 wmean= wmean / s1;

 durmax=s1;

 durmin=s1;

 // call diff

 wmax=zeros(1,1);

 wmin=zeros(1,1);

 kmax=0;

 kmin=0;

 dym=wym;

 s0=size(dym);

 s1=s0(1);

 // detect zero-cross points

 for v=2:s1

    if ( ( dym(v-1) > 0.) & ( dym(v) <= 0.) ) then

       kmax=kmax+1;

       wmax(kmax)=v;

    elseif ( ( dym(v-1) < 0.) & ( dym(v) >= 0.) ) then

       kmin=kmin+1;

       wmin(kmin)=v;

    end

 end

 // cal peak-peak duration mean

 if kmax >= 2 then

   durmax=0.;

   for v=2:kmax

     durmax=durmax + ( wmax(v) - wmax(v-1));

   end

   durmax=durmax / (kmax-1);

 end

 if kmin >= 2 then

   durmin=0.;

   for v=2:kmin

     durmin=durmin + ( wmin(v) - wmin(v-1));

   end

   durmin=durmin / (kmin-1);

 end

 //

 if (dym(2) - dym(1)) >= 0. then

  kyokusei =0.;

 else

  kyokusei =1.;

 end

endfunction

//--------------------------------------------------------------------

function [x0,binf,bsup]=set_x0( wx0, wym, wtype0)

 s0=size(wx0);

 s1=s0(1);

 x0=zeros(s1,1);

 binf0=zeros(s1,1);

 bsup0=zeros(s1,1);

 for v=1:s1

  x0(v)=wx0(v);

 end

 [wmean, kmax, wmax, durmax, kmin, wmin, durmin, maxabs, maxpoint, kyokusei ]=detectp(wym);

 // set x0 according to TYPE

 s0=size(wym);

 s1=s0(1);

 

 if wtype0 == 2 then       // reduce

   mix0=0.5;

   expand0=exp(0.);

   x0(1)= -1. / s1;

   x0(2)= maxabs * mix0 / expand0;

   x0(3)= (2. * %pi) / durmax; 

   x0(4)= %pi * kyokusei;

   x0(5)= -1. / s1;

   x0(6)= maxabs * ( 1. - mix0) / expand0;

   x0(7)= x0(3); // frequency double

   x0(8)= 0.;    // phase same

   x0(9)=wmean;

  

   // DC

   binf(9)= -1. * maxlevel;

   bsup(9)= maxlevel;

   // phase

   binf(8)= -1. * %pi;

   bsup(8)= %pi;

   binf(4)= -2. * %pi;

   bsup(4)= 2. * %pi;

   // frequency

   binf(3)= 0.;

   bsup(3)= %pi / 2.;

   // ??? what is good setting ?

   binf(7)=  0.;   

   bsup(7)=   %pi / 4.;  

   // level

   binf(2)= 0.;

   bsup(2)= maxlevel;

   binf(6)= 0.;

   bsup(6)= maxlevel;

   // exp( a< 0

   binf(1)= -1.0 * log(maxlevel);

   bsup(1)= 0.;

   binf(5)= -1.0 * log(maxlevel);

   bsup(5)= 0.;

 elseif wtype0 ==3 then    // change: this is mikansei !

   expand0=exp(0.);

   expand1=exp(1.0);

   x0(1)= -1. / s1;   // reduce

   x0(2)= maxabs *  expand0;

   x0(3)= (2. * %pi) / durmax; 

   x0(4)= %pi * kyokusei;

   x0(5)= 1. / s1;    // expand

   x0(6)= maxabs *  expand1;

   x0(7)= 0.;    // frequency same

   x0(8)= 0.;    // phase same

   x0(9)=wmean;

  

   // DC

   binf(9)= -1. * maxlevel;

   bsup(9)= maxlevel;

   // phase

   binf(8)= -1. * %pi;

   bsup(8)= %pi;

   binf(4)= -2. * %pi;

   bsup(4)= 2. * %pi;

   // frequency

   binf(3)= 0.;

   bsup(3)= %pi / 2.;

   // ??? what is good setting ?

   binf(7)=   -1.0 * %pi / 4.;   

   bsup(7)=   %pi / 4.;  

   // level

   binf(2)= -1.0 * maxlevel;

   bsup(2)= 0.;

   binf(6)= 0.;

   bsup(6)= maxlevel;

   // exp( a< 0

   binf(1)= 0.;

   bsup(1)= log(maxlevel);

   // exp( a> 0

   binf(5)= 0.;

   bsup(5)= log(maxlevel);





 else   // default wtype = 1  expand

   mix0=0.5;

   expand0=exp(1.0);

   x0(1)= 1. / s1;

   x0(2)= maxabs * mix0 / expand0;

   x0(3)= (2. * %pi) / durmax; 

   x0(4)= %pi * kyokusei;

   x0(5)= 1. / s1;

   x0(6)= maxabs * ( 1. - mix0) / expand0;

   x0(7)= x0(3); // frequency double

   x0(8)= 0.;    // phase same

   x0(9)=wmean;

  

   // DC

   binf(9)= -1. * maxlevel;

   bsup(9)= maxlevel;

   // phase

   binf(8)= -1. * %pi;

   bsup(8)= %pi;

   binf(4)= -2. * %pi;

   bsup(4)= 2. * %pi;

   // frequency

   binf(3)= 0.;

   bsup(3)= %pi / 2.;

   // ??? what is good setting ?

   binf(7)=  0.;   

   bsup(7)=   %pi / 4.;  

   // level

   binf(2)= 0.;

   bsup(2)= maxlevel;

   binf(6)= 0.;

   bsup(6)= maxlevel;

   // exp( a> 0

   binf(1)= 0.;

   bsup(1)= log(maxlevel);

   binf(5)= 0.;

   bsup(5)= log(maxlevel);

 end

endfunction

//--------------------------------------------------------------------

function [wstr1,wstr2,wstr3,wstr0]=disp_xopt(f,xopt, xflag)



if xflag == 0 then

 wstr1=sprintf('%f',f);

 wstr1='f= ' + wstr1;

 disp(wstr1);

 //

 disp('unit converted to frequecny Hz and phase deg');

 wstr0='unit converted to frequecny Hz and phase deg';

 wstr1='exp( ' + sprintf('%f', xopt(1)) + ' * t ) * ' + sprintf('%f', xopt(2)) + ' * sin( ' + sprintf('%f', xopt(3) * fs / (%pi * 2)) + ' * t + ' + sprintf('%f', xopt(4) / (%pi * 2) * 360.) + ' )';

 wstr1=wstr1 + ' + exp( ' + sprintf('%f', xopt(5)) + ' * t ) * ' + sprintf('%f', xopt(6)) + ' * sin( (' + sprintf('%f', xopt(3) * fs / (%pi * 2)) + ' + ' + sprintf('%f', xopt(7) * fs / (%pi * 2)) + ') * t + ('  + sprintf('%f', xopt(4) / (%pi * 2) * 360.)  + ' + ' + sprintf('%f', xopt(8) / (%pi * 2) * 360.) + ') )';

 wstr1=wstr1 + ' + ' + sprintf('%f', xopt(9));

 disp(wstr1);

elseif xflag==2 then

 wstr0='unit converted to frequecny Hz and phase deg';

 wstr1='exp( ' + sprintf('%f', xopt(1)) + ' * t ) * ' + sprintf('%f', xopt(2)) + ' * sin( ' + sprintf('%f', xopt(3) * fs / (%pi * 2)) + ' * t + ' + sprintf('%f', xopt(4) / (%pi * 2) * 360.) + ' )';

 wstr2='+ exp( ' + sprintf('%f', xopt(5)) + ' * t ) * ' + sprintf('%f', xopt(6)) + ' * sin( (' + sprintf('%f', xopt(3) * fs / (%pi * 2)) + ' + ' + sprintf('%f', xopt(7) * fs / (%pi * 2)) + ') * t + ('  + sprintf('%f', xopt(4) / (%pi * 2) * 360.)  + ' + ' + sprintf('%f', xopt(8) / (%pi * 2) * 360.) + ') )';

 wstr3='+ ' + sprintf('%f', xopt(9));

else    // xflag == 1

 disp( f, 'f');

 disp( xopt(1), 'exp(');

 disp( xopt(2), ' * sin(');

 disp( xopt(3) * fs / (%pi * 2), 'base freq');

 disp( xopt(4) / (%pi * 2) * 360. , 'base phase');

 disp( xopt(5), 'exp(');

 disp( xopt(6), ' * sin(');

 disp( xopt(7) * fs / (%pi * 2), 'sabun freq');

 disp( xopt(8) / (%pi * 2) * 360. , 'sabun phase');

 disp( xopt(9) , 'DC');

end

endfunction

//---------------------------------------------------------------------------------

function [hatyou]=get_hatyou( wxopt)

 // return wave length of basic one

 if wxopt(3) > (wxopt(3)+wxopt(7)) then

   basic0= wxopt(3) + wxopt(7);

 else

   basic0= wxopt(3);

 end

 //

 hatyou= ceil(%pi * 2 / basic0);

 //

endfunction

//---------------------------------------------------------------------------------

function [rtcode]=is_check_type1( wx0, wym )

// when rtcode is 0,  the wx0 is type 1,

// else it isnot.  in rtcode,  2~N shows  negative point. 

 

 //--- checking sign ---

 rtcode=0; 

 if ( (wx0(1) > 0.) & (wx0(5) > 0.) )  then // (1)check exp's sign 

 else

   rtcode=rtcode + 1;

 end  // (1)



 if ( (wx0(2) > 0.) & (wx0(6) > 0.) )  then // (2)check amplitude's sign

 else

   rtcode=rtcode + 2;

 end  // (2)



 if (wx0(3) > 0.)    then   // (3)check base frequency  sign

 else

   rtcode=rtcode + 4;

 end  // (3)



 if ( (wx0(3) + wx0(7)) > 0.)    then // (4)check secondary frequency  sign

 else

   rtcode=rtcode + 8;

 end  // (4)



 //--- checking similarity

 s0=size(wym);

 s1=s0(1);

 paramslop=0.5;          // (5) limit slop angle

 e1=exp(wx0(1) * s1);

 e5=exp(wx0(5) * s1);

 if (((e1 >= 1.0) & (e5 >= 1.0)) | ((e1 < 1.0) & (e5 < 1.0))) then // slop direction same desuka.

   if e1 >= e5 then

      if e5 > (e1 * paramslop) then

      else

        rtcode=rtcode + 16;

      end

   else

      if e1 > (e5 * paramslop) then

      else

        rtcode=rtcode + 16;

      end

   end

 else

  rtcode=rtcode + 16;

 end

 





 //--- checking mix ratio of two signals

 parammix=0.1;          // (6) minimum mix ratio

 if( abs(wx0(2)) > abs(wx0(6))) then

    if ( abs(wx0(6)) > ( abs(wx0(2)) *  parammix )) then

    else

      rtcode=rtcode + 32;

    end

 else

    if ( abs(wx0(2)) > ( abs(wx0(6)) *  parammix )) then

    else

      rtcode=rtcode + 32;

    end

 end   // (6)

endfunction

//---------------------------------------------------------------------------------

function [rtcode]=is_check_type2( wx0, wym )

// when rtcode is 0,  the wx0 is type 1,

// else it isnot.  rtcode's 2~N shows  negative point. 

 

 //--- checking sign ---

 rtcode=0; 

 if ( (wx0(1) < 0.) & (wx0(5) < 0.) )  then // (1)check exp's sign 

 else

   rtcode=rtcode + 1;

 end  // (1)



 if ( (wx0(2) > 0.) & (wx0(6) > 0.) )  then // (2)check amplitude's sign

 else

   rtcode=rtcode + 2;

 end  // (2)



 if (wx0(3) > 0.)    then   // (3)check base frequency  sign

 else

   rtcode=rtcode + 4;

 end  // (3)



 if ( (wx0(3) + wx0(7)) > 0.)    then // (4)check secondary frequency  sign

 else

   rtcode=rtcode + 8;

 end  // (4)



 //--- checking similarity

s0=size(wym);

 s1=s0(1);

 paramslop=0.5;          // (5) limit slop angle

 e1=exp(wx0(1) * s1);

 e5=exp(wx0(5) * s1);

 if (((e1 >= 1.0) & (e5 >= 1.0)) | ((e1 < 1.0) & (e5 < 1.0))) then // slop direction same desuka.

   if e1 >= e5 then

      if e5 > (e1 * paramslop) then

      else

        rtcode=rtcode + 16;

      end

   else

      if e1 > (e5 * paramslop) then

      else

        rtcode=rtcode + 16;

      end

   end

 else

  rtcode=rtcode + 16;

 end

 

 //--- checking mix ratio of two signals

 parammix=0.1;          // (6) minimum mix ratio

 if( abs(wx0(2)) > abs(wx0(6)) ) then

    if ( abs(wx0(6)) > ( abs(wx0(2)) *  parammix )) then

    else

      rtcode=rtcode + 32;

    end

 else

    if ( abs(wx0(2)) > ( abs(wx0(6)) *  parammix )) then

    else

      rtcode=rtcode + 32;

    end

 end   // (6)

endfunction

//--------------------------------------------------------------------

function [yopt]=irekae(wopt)

  s0=size(wopt);

  s1=s0(1);

  yopt=zeros(s1,1);

  for v=1:s1

    yopt(v)=wopt(v);

  end

  //  if frequency diffrence is negative, exchange 1st and 2nd kou

  // y  = exp(x(1)*t). * ( x(2) * sin(  x(3) * t +  x(4) ) )  +  exp(x(5)*t). * ( x(6) * sin(  (x(3)+x(7)) * t +  (x(4)+x(8)) ) )+ x(9)

  if wopt(7) < 0. then

    yopt(1)=wopt(5);

    yopt(2)=wopt(6);

    yopt(3)=wopt(3)+wopt(7);

    yopt(4)=wopt(4)+wopt(8);

    //

    yopt(5)=wopt(1);

    yopt(6)=wopt(2);

    yopt(7)= -1. * wopt(7);

    yopt(8)= -1. * wopt(8);

  end

  // 

endfunction

//--------------------------------------------------------------------

function [kcands]=get_candidates(wym)

 //  kcands(2*n)  term start point

 //  kcands(2*n+1)  term end point



 [wmean, kmax, wmax, durmax, kmin, wmin, durmin, maxabs, maxpoint, kyokusei ]=detectp(wym);

 /////[wmean, kmax, wmax, durmax, kmin, wmin, durmin, maxabs, maxpoint, kyokusei ]=detectzc(wym);

 s0=size(wym);

 s1=s0(1);

 kcands=zeros(1,1);

 // default full term

 kcands(1)=1;

 kcands(2)=s1;

 if (( kyokusei <= 0.) & ( kmin > 0)) then   // rise -> peak up peak minus

   for v=1:kmin

     if v==1 then

       kcands(2*v-1)=1;

     else

       kcands(2*v-1)=wmin(v-1);

     end

     kcands(2*v)=wmin(v);

   end

   // what to do last term ?

   if s1 > wmin(kmin) then

     ////kcands(2*(kmin+1)-1)= wmin(kmin);

     ////kcands(2*(kmin+1))=s1;

   end

 elseif (( kyokusei >= 1.) & ( kmax > 0))                    // fall -> peak up peak plus

   for v=1:kmax

     if v==1 then

       kcands(2*v-1)=1;

     else

       kcands(2*v-1)=wmax(v-1);

     end

     kcands(2*v)=wmax(v);

   end

   // what to do last term ?

   if s1 > wmax(kmax) then

     ////kcands(2*(kmax+1)-1)= wmax(kmax);

     ////kcands(2*(kmax+1))=s1;

   end

 end

endfunction

//--------------------------------------------------------------------

function [scan_xopt,scan_f,scan_type]=scan_candidates(kcands,wym,display0)

 s0=size(kcands);

 s1=s0(1)/2;

 ws0=size(wx0);

 ws1=ws0(1);

 scan_xopt=zeros(s1,ws1);

 scan_f=zeros(s1);

 scan_type=zeros(s1);

 //s1=2;  // this line comment out !

 //

 //+display control flag

 disp0=-1;

 pdisp0=0;

 fdisp0=-1;   // fdisp0 >= 10 disp_xopt 

 //-display control flag

 ////+one more search flag, using [hatyou]=get_hatyou( wxopt) and set length >= hatyou

 ////use_hatyou0=0; 

 ////-one more search flag

 for loop=1:s1

   if display0 >= 1 then

    wstr1='loop= ' + sprintf('%d',loop) ;

    disp(wstr1);

   end

   // set failure pattern

   scan_type(loop)=-1.;

   scan_f(loop)=-1.;

   //

   wsp1=kcands(2*loop-1);

   wep1=kcands(2*loop);

   [m,wtm,wym2,wwm]= set_mtmymwm(wsp1,wep1,wym);

   //(1) try type1

   wtype0=1

   [wx0,binf0,bsup0]=set_x0( wx0, wym2, wtype0);

   disp0=disp0+pdisp0;

   [f,xopt2, gropt] = simplecall1b(wtm,wym2,wwm,wx0,disp0);

   [xopt]=irekae(xopt2);





   [rtcode]=is_check_type1(xopt,wym2); 

   if rtcode <= 0. then

     scan_type(loop)=wtype0;

     scan_f(loop)=f;

     for v=1:ws1

      scan_xopt(loop,v)=xopt(v);

     end

   end

    if fdisp0 >= 1 then

     wstr1='loop= ' + sprintf('%d',loop) + '  type= ' + sprintf('%d',wtype0) + ' rtcode= ' + sprintf('%d',rtcode);

     disp(wstr1);

     if fdisp0 >= 10 then

       disp_xopt(f,xopt, 0);

     end

    end

   //

   //(2) try type2

   wtype0=2

   [wx0,binf0,bsup0]=set_x0( wx0, wym2, wtype0);

   disp0=disp0+pdisp0;

   [f,xopt2, gropt] = simplecall1b(wtm,wym2,wwm,wx0,disp0);

   [xopt]=irekae(xopt2);

   [rtcode]=is_check_type2(xopt,wym2);   // disp(rtcode);

   if rtcode <= 0. then

    if  ( (scan_f(loop) <= -1.) | ((scan_f(loop) > 0.) & ( f < scan_f(loop))) ) then

     scan_type(loop)=wtype0;

     scan_f(loop)=f;

     for v=1:ws1

      scan_xopt(loop,v)=xopt(v);

     end  // for v=1:ws1

    end

   end  // if rtcode <= 0. then

 

    if fdisp0 >= 1 then

     wstr1='loop= ' + sprintf('%d',loop) + '  type= ' + sprintf('%d',wtype0) + ' rtcode= ' + sprintf('%d',rtcode);

     disp(wstr1);

     if fdisp0 >= 10 then

      disp_xopt(f,xopt, 0);

     end    

    end

 end  // for loop=1:s1

 //

 //

 // kekka display

 wxopt=zeros(ws1,1);

 if display0 >= 2 then

   

   for loop=1:s1

      wstr1='loop= ' + sprintf('%d',loop) + '  type= ' + sprintf('%d',scan_type(loop));

     disp(wstr1);

     if scan_type(loop) >= 1 then

        for v=1:ws1

         wxopt(v)=scan_xopt(loop,v);

        end  

       disp_xopt(scan_f(loop),wxopt, 0);

     end

   end  // for loop=1:s1

 end  // if display0 >= 1 then

 //



endfunction

//--------------------------------------------------------------------

function plot_result(wym,kcands,scan_xopt,scan_f,scan_type)

 xset('window',0); 

 clf();

 // get entire waveform

 s0=size(wym);

 s1=s0(1);

 wsp1=1;

 wep1=s1;

 [m,wtm,wym2,wwm]= set_mtmymwm(wsp1,wep1,wym);

 // cal waveform from xopt

 wyth=zeros(s1,1);

 wfreq=zeros(s1,1);

 wdfreq=zeros(s1,1);

 wtype=zeros(s1,1);

 s2=size(kcands);

 s3=s2(1)/2;

 ws0=size(wx0);

 ws1=ws0(1);

 wxopt=zeros(ws1,1);

 for loop=1:s3

  if scan_type(loop) >= 1 then

    wsp1=kcands(2*loop-1);

    wep1=kcands(2*loop);

    t0=0.;

    for v=wsp1:wep1

      wtype(v)=scan_type(loop);

      wfreq(v)=scan_xopt(loop,3) * fs / (%pi * 2);

      wdfreq(v)=scan_xopt(loop,7) * fs / (%pi * 2);

      for v2=1:ws1

        wxopt(v2)=scan_xopt(loop,v2);

      end

      wyth(v)=yth(t0, wxopt);

      t0=t0+1.;

    end

  end   

 end

 subplot(411);

 plot(wtm,wym2,'+',wtm,wyth,'-');

 wstr1 = '+:raw waveform   and  estimated waveform by xopt'

 xtitle( wstr1,'','' );

 subplot(412);

 plot(wtm,wfreq,'-');

 wstr1 = 'basic frequency    Hz'

 xtitle( wstr1,'','' );

 subplot(413);

 plot(wtm,wdfreq,'-');

  wstr1 = 'frequency difference    Hz'

 xtitle( wstr1,'','' );

 subplot(414);

 plot(wtm,wtype,'-');

 wstr1 = 'type 1:expand 2:reduce 0:undefined'

 xtitle( wstr1,'','' );

endfunction

//--------------------------------------------------------------------

function test0_show(wym0,PLAN_SP0EP0,PLAN_SP1EP1,PLAN_XOPT,PLAN_TYPE,rtcode_xloop)

 global eh_var1

 if eh_var1 > 0 then

    xset('window',(eh_var1 - 1)); 

    clf();

    //

   if rtcode_xloop > 0 then

     swx0=size(wx0);

     swx1=swx0(1);

     xopt=zeros(swx1,1);

     for v=1:swx1

       xopt(v)=PLAN_XOPT(rtcode_xloop,v);

     end

     type0=PLAN_TYPE(rtcode_xloop);

     swym0=size(wym0);

     swym1=swym0(1);

     sp1=PLAN_SP0EP0(rtcode_xloop,1);

     ep1=PLAN_SP0EP0(rtcode_xloop,2);

     [m,tm,ym,wm]= set_mtmymwm(sp1,ep1,ym0);

     sp2=PLAN_SP1EP1(rtcode_xloop,1);

     ep2=PLAN_SP1EP1(rtcode_xloop,2);

     tm2=zeros(ep2-sp2+1);

     ym2=zeros(ep2-sp2+1);

     for v=1:(ep2-sp2+1)

      tm2(v)=sp2+v-1;

      ym2(v)=ym(sp2-sp1+v);

     end

     plot((tm+sp1),ym,'+');

     plot( tm2,yth((tm2-sp1),xopt),'-');

     plot(tm2,(  ym2- yth(  (tm2-sp1),xopt  )   ),'--r');

     f=0; // dummy value set

     [wstr1,wstr2,wstr3,wstr0]=disp_xopt(f,xopt, 2);

     wstr4='type= ' + sprintf('%d',type0);

     xstring(sp1,1,[ wstr0 ; wstr1 ; wstr2 ; wstr3 ; wstr4 ]);





  //

  else // if rtcode_xloop <= 0 then

 

  end  // if rtcode_xloop > 0 then

 end  // if eh_var1 > 0 then

endfunction

//--------------------------------------------------------------------

function test0(wym0,wtm0,MAX_LOOP_COUNTER,MAX_LOOP_COUNTERB,ONE_STEP)

 //(0)clear all graphic windows

 global global eh_var1

 if eh_var1 > 1 then

  for v=1:eh_var1

   xdel( (v-1) );

  end

 end

 eh_var1=1;

 //(1)

 [kcands]=get_candidates(wym0); 

 disp( kcands);

 s0=size(kcands);

 s1=s0(1);

 index0=1;

 swx0=size(wx0);

 swx1=swx0(1);

 swym0=size(wym0);

 swym1=swym0(1);

 //(2) 1st step

 sp1=kcands(2*index0-1);

 ep1=kcands(2*index0);

 //(3)search loop

 for xloop=1:MAX_LOOP_COUNTERB

   [PLAN_SP0EP0,PLAN_SP1EP1,PLAN_XOPT,PLAN_TYPE,rtcode_xloop]=test1(wym0,sp1,ep1,MAX_LOOP_COUNTER,ONE_STEP);

   // result shows

   eh_var1=eh_var1+1;

   test0_show(wym0,PLAN_SP0EP0,PLAN_SP1EP1,PLAN_XOPT,PLAN_TYPE,rtcode_xloop);

   //

   if rtcode_xloop >0 then

      sp1=PLAN_SP1EP1(rtcode_xloop,2);

      // set index0 as -1

      index0=-1;

      for v=1:(s1/2)

       // In which range is sp1

       if ( (sp1 >= kcands(2*v-1)) & (sp1 < kcands(2*v) )) then

         index0=v;

         break;

       end

      end

      if index0 < 0 then  // out of kcands range

       break;

      else

       ep1=kcands(2*index0);

        // make it more than minimum leghth 

       if( (ep1-sp1) < (MINIMUM_LEGTH0 * swx1) ) then

         ep1= sp1 + MINIMUM_LEGTH0 * swx1;

        end

        // check out of data length

        if ep1 > swym1 then

         break;

        end

      end

   else

     break;

   end

   if xloop == MAX_LOOP_COUNTERB then

     disp('waring: xloop == MAX_LOOP_COUNTERB in function test0');

   end

 end  // for xloop=1:2

      //(3)search loop

 // delete working window

 xdel(0);

  //

 xset('window',0);   

 clf();

 plot(wtm0,wym0);

endfunction

//---------------------------------------------------------------------

function [PLAN_SP0EP0,PLAN_SP1EP1,PLAN_XOPT,PLAN_TYPE,rtcode_xloop]=test1(wym0,sp1,ep1,MAX_LOOP_COUNTER,ONE_STEP)

//

//MAX_LOOP_COUNTER=2;

PLAN_SP0EP0=zeros(MAX_LOOP_COUNTER,2);

PLAN_SP1EP1=zeros(MAX_LOOP_COUNTER,2);  // fitable area

xs0=size(wx0);

xs1=xs0(1);

PLAN_XOPT=zeros(MAX_LOOP_COUNTER,xs1);

PLAN_TYPE=zeros(MAX_LOOP_COUNTER,1);



  ws0=size(wym0);

  ws1=ws0(1);



  //(3) search loop

  for xloop=1:MAX_LOOP_COUNTER

   type0=1;

   [m,tm,ym,wm]= set_mtmymwm(sp1,ep1,ym0);

   [x0,binf0,bsup0]=set_x0( wx0, ym, type0);

   //(4)simple leastsq call

   [f1,xopt1, gropt1] = simplecall1(); 



   type0=2;

   [m,tm,ym,wm]= set_mtmymwm(sp1,ep1,ym0);

   [x0,binf0,bsup0]=set_x0( wx0, ym, type0);

   //(4)simple leastsq call

   [f2,xopt2, gropt2] = simplecall1(); 



   if f1 < f2 then

    xoptc=xopt1;

    f=f1;

    gropt=gropt1;

   else

    xoptc=xopt2;

    f=f2;

    gropt=gropt2;

   end



   [xopt]=irekae(xoptc);

   disp_xopt(f,xopt,0);

    //

    wtype=-1;

    [rtcode]=is_check_type1( xopt, wm );

    if rtcode==0 then

     wtype=1;

    else

      [rtcode]=is_check_type2( xopt, wm );

      if rtcode==0 then

       wtype=2;

      end

    end

    for v35=1:xs1

      PLAN_XOPT(xloop,v35)=xopt(v35);

    end

    PLAN_TYPE(xloop)=wtype;

   //(5)

    [fsp1,fep1,dsp1,dep1]=find_fitable_area(tm,ym,xopt,FITABLE_LEVEL0);

    PLAN_SP1EP1(xloop,1)=sp1+dsp1;

    PLAN_SP1EP1(xloop,2)=ep1-dep1;

    PLAN_SP0EP0(xloop,1)=sp1;

    PLAN_SP0EP0(xloop,2)=ep1;

   wstr1='sp1 fsp= ' + sprintf('%d %d',sp1,(sp1+dsp1)) + '  ep1 fep= ' + sprintf('%d %d',ep1,(ep1-dep1)) + '  type= ' + sprintf('%d',wtype) ;

   disp(wstr1);

   //(6)

   if xloop == 1 then

     if ((fsp1== -1) | (fep1== -1)) then

       rtcode_xloop=-1;

       break;

     end

   else 

      if (PLAN_SP1EP1(xloop,2) - PLAN_SP1EP1(xloop,1)) < (PLAN_SP1EP1(xloop-1,2) - PLAN_SP1EP1(xloop-1,1)) then

        rtcode_xloop=xloop-1;

        break;

      end

   end

   //(7)

    if ep1+ONE_STEP <=ws1 then

      ep1=ep1+ONE_STEP;

    else

      rtcode_xloop=xloop;

      break;

    end

    if xloop == MAX_LOOP_COUNTER then

     disp('waring: xloop == MAX_LOOP_COUNTER in function test1');

   end

  end  //for xloop=1:MAX_LOOP_COUNTER0

      // (3) search loop

  wstr1='rtcode_xloop= ' + sprintf('%d',rtcode_xloop);

  disp(wstr1);

endfunction

//------------------------------------------------------------

function [wsp1,wep1,dsp1,dep1]=find_fitable_area(wtm,wym,wxopt,FITABLE_LEVEL)

// this finds only bottom and its round area

//FITABLE_LEVEL=75.;

//

 s0=size(wym);

 s1=s0(1);

 [ gosa1, gosa1abs ] = gosa(wtm,wym,wxopt);

 // find bottom points

 [ min0, k ]= min( gosa1abs);

 //

 if min0 > FITABLE_LEVEL then

   // there is no fitable area.

   wsp1=-1;

   wep1=-1;

   dsp1=-1;

   dep1=-1;

 else

   wep1=s1;

   for v=k:s1

    if gosa1abs(v) > FITABLE_LEVEL then

      wep1=v-1;

      break;

    end

   end

   wsp1=1;

   for v=k:-1:1

    if gosa1abs(v) > FITABLE_LEVEL then

      wsp1=v+1;

      break;

    end

   end

   dsp1=wsp1-1;

   dep1=s1-wep1;

 end

endfunction

//----------------------------------------------------------------

function disp_fitable_area()

   [fsp1,fep1,dsp1,dep1]=find_fitable_area(tm,ym,xopt,FITABLE_LEVEL0);

   wstr1='sp1 fsp= ' + sprintf('%d %d',sp1,(sp1+dsp1)) + '  ep1 fep= ' + sprintf('%d %d',ep1,(ep1-dep1))  ;

   disp(wstr1);

endfunction

//

//--------------------------------------------------------------------

//

// ++ MENU MENU MENU =================================================

//

//  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_set');

//////addmenu('step1_set',[ '(1) set_point()','(2) set_init','(3) call leastsq','(4) call optim']);

//////step1_set=[ '[sp1, ep1, m, tm, ym, wm, type0 ]= set_sp1_ep1(m0,tm0,ym0,10);'; '[x0,binf0,bsup0]=set_x0( wx0, ym, type0);' ; '[f,xopt, gropt] = simplecall1(); disp_xopt(f,xopt,0);';'[f,xopt] = simplecall2(binf0,bsup0); disp_xopt(f,xopt,0);' ] ;

addmenu('step1_set',[ '(1) set_point()','(2) set_init','(2b) Or x0=xopt','(3) call leastsq','(4-1) check xopt as type 1','(4-2) check xopt as type 2','(5) hatyou of basic wave ','(6) disp_fitable_area','(9) read *.sci file and exce']);

step1_set=[ '[sp1, ep1, m, tm, ym, wm, type0 ]= set_sp1_ep1(m0,tm0,ym0,0);'; '[x0,binf0,bsup0]=set_x0( wx0, ym, type0);'; ' x0=xopt;'  ; '[f,xopt2, gropt] = simplecall1(); [xopt]=irekae(xopt2); disp_xopt(f,xopt,0);';'[rtcode]=is_check_type1(xopt,ym); disp(rtcode);' ; '[rtcode]=is_check_type2(xopt,ym); disp(rtcode);' ; '[hatyou]=get_hatyou( xopt); disp(hatyou);' ; 'disp_fitable_area();' ; '[rf0]=readx(); exec( rf0,-1);'] ;

//

end //if ADD_MENU_PLOT == 1 then

//

//

ADD_MENU_PLOT2=-1;       // set ADD_MENU 1 to add menu button of some functions

//

//

if ADD_MENU_PLOT2 == 1 then

delmenu('step2_set');

addmenu('step2_set', ['(1) get_candidate term','(2) scan_candidate','(3) plot result' ]);

step2_set=['[kcands]=get_candidates(ym0); disp(kcands);','[scan_xopt,scan_f,scan_type]=scan_candidates(kcands,ym0,2)','plot_result(ym0,kcands,scan_xopt,scan_f,scan_type)'];

//

end //if ADD_MENU_PLOT2 == 1 then

//

ADD_MENU_PLOT3=1;       // set ADD_MENU 1 to add menu button of some functions

//

//

if ADD_MENU_PLOT3 == 1 then

delmenu('step3_set');

addmenu('step3_set', ['(1) test']);

step3_set=['test0(ym0,tm0,MAX_LOOP_COUNTER0,MAX_LOOP_COUNTERB0,ONE_STEP0);'];

//

end //if ADD_MENU_PLOT2 == 1 then



//

//----------------------------------------------------------------------

T_DEMO=1;   // set T_DEMO=1 for demonstration, beside set T_DEMO=0 when normal 

//

//

//-------------------------------

T_DEMO_1=0;

if T_DEMO_1 == 1 then

 // (1) set_et_plot()

 sp1=1;

 ep1=35;

 type0=1;

 wstr1='sp1= ' + sprintf('%d',sp1) + '  ep1= ' + sprintf('%d',ep1) + '  type= ' + sprintf('%d',type0) ;

 disp(wstr1);

 [m,tm,ym,wm]= set_mtmymwm(sp1,ep1,ym0);

 //(2) set_init

 [x0,binf0,bsup0]=set_x0( wx0, ym, type0);

 //(3)simple leastsq call

 [f,xopt2, gropt] = simplecall1(); 

 [xopt]=irekae(xopt2);

 disp_xopt(f,xopt,0);

end

//-------------------------------

T_DEMO_2=0;

if T_DEMO_2 == 1 then

 //(1)get candidates

 [kcands]=get_candidates(ym0); 

 kcands

 //(2)scan candidates

 [scan_xopt,scan_f,scan_type]=scan_candidates(kcands,ym0,2);

 //(3)plot

 plot_result(ym0,kcands,scan_xopt,scan_f,scan_type);

end

//-------------------------------

if T_DEMO == 1 then

 test0(ym0,tm0,MAX_LOOP_COUNTER0,MAX_LOOP_COUNTERB0,ONE_STEP0);

end