%
%    sf_engine.m  ver 1.6  August 23, 2013
%
function[a,x1r,x2r,x3r,x4r,x5r,cnew] = ...
    sf_engine(a,t,num2,flow,fup,ie,nt,dur,x1r,x2r,x3r,x4r,...
                                             x5r,running_sum,original,istr)
%
rand('state',0)
tp=2*pi;
%
ave=mean(a);
sd=std(a);
%
am=2.*sd;
n=num2;
%
dt=(t(n)-t(1))/(n-1);
sr=1/dt;
%
pu=1.;
pl=0.;
%
errormax=1.0e+53;
%
jk=0;
%
delta=0.001;
nnn=fix(0.1*nt);
%
zcf=0;
fft_freq=0;
freq_est=0;
%
if(istr==2)
    a=a-mean(a);
%
    iband=1;
    fl=sr/8;
    fh=0.;
    iphase=1;
    [a,mu,sd,rms]=...
                  Butterworth_filter_function_alt(a,dt,iband,fl,fh,iphase);
%
    nzc=0;
    pzc=0;
    for i=1:length(a)-1
%
       if(a(i) <=0 && a(i+1) > 0)
           pzc=pzc+1;
       end
%       
       if(a(i) >=0 && a(i+1) < 0)
           nzc=nzc+1;      
       end 
    end
%
    zcf=((pzc+nzc)/2)/dur;  
%
    m_choice=1;  % mean removal
    h_choice=1;  % rectangual window
%    
    N=2^floor(log(n)/log(2));
    [freq,full,phase,complex_FFT]=full_FFT_core(m_choice,h_choice,a,N,dt);    
%
    magnitude_FFT=[freq full];
%    
    [~,fft_freq]=find_max(magnitude_FFT);
%
    out1=sprintf('\n zero crossing frequency = %8.4g Hz ',zcf);
    disp(out1);
%    
    out1=sprintf(' FFT frequency = %8.4g Hz ',fft_freq);
    disp(out1);
%
end
%
disp(' ');
disp('  Trial     Error      Amplitude   Freq(Hz)   Phase(rad)  damp   delay(sec) ');
%
for j=1:nt
%
    if(istr==2)
        if(rand<0.5)
            freq_est=zcf;
        else    
            freq_est=fft_freq;
        end    
    end
%
    if(j>nnn)
%
      fr=(x2r(ie)/tp);
%
      if( (fup-fr)/fup < delta)
          fup=fup*(1+delta);
      end
 %
      if(flow>1.0e-06)
        if( (fr-flow)/flow < delta)
            flow=flow*(1-delta);
        end  
      end
%      
    end
%
	jk=jk+1;
 %           
    if(jk==10000)
		out4=sprintf('\n %ld ',j);
        disp(out4);
		jk=0;
    end
%
	x1=rand;
	x2=rand;
	x3=rand;
	x4=rand*rand;
	x5=rand;
%
    if(rand<0.2)
        x5=0.;
    end
%
    n1=fix(2.*nt/3.);
    n2=fix(4.*nt/5.);
	if( j < n1 )
				x1=2.*am*x1;
				x2=((fup-flow)*x2+flow)*tp;
				x3=((pu-pl)*x3+pl)*tp;
				x4=(x4^2.5);
				x5=dur*(x5^2.);
    end
    if( j >= n1 && j < n2)
				x1=x1r(ie)*(0.98+0.04*x1);
				x2=x2r(ie)*(0.98+0.04*x2);
				x3=x3r(ie)*(0.98+0.04*x3);
%				x4=x4r(ie)*(0.80+0.40*x4);
				x5=x5r(ie)*(0.98+0.04*x5);
    end
    if( j >= n2 )
				x1=x1r(ie)*(0.99+0.02*x1);
				x2=x2r(ie)*(0.99+0.02*x2);
				x3=x3r(ie)*(0.99+0.02*x3);
				x4=x4r(ie)*(0.90+0.20*x4);
				x5=x5r(ie)*(0.99+0.02*x5);
    end
%
    if(istr==2 && freq_est < sr/6)
        if(rand<0.3)
            x2=freq_est*(0.98+0.04*rand());
        else
            if(rand<0.5)
                x2=freq_est*(0.95+0.10*rand());
            else
                x2=freq_est*(0.90+0.20*rand()); 
            end    
        end
        x2=x2*tp;
%
        if(j>nnn && rand<0.3 )
            if(rand<0.5)
                x2=x2r(ie)*(0.99+0.02*x2);
            else
                x2=x2r(ie)*(0.995+0.01*x2);                
            end
        end    
%
    end
%
%%%
	if(x2 > fup*tp || freq_est > sr/6 )
                x2 = fup*tp;
    end    
%
	if(x2 < flow*tp)
                x2 = flow*tp;
    end 
%
	if(x4 >= 0.85 )
                x4=0.05;
    end    
%		
	if(j==1)
                x1=0.;
    else
        if(x1<1.0e-12)
            x1=am*rand;
        end
    end
%
	error=0.;
%        
    domegan=x4*x2;
    omegad=x2*sqrt(1.-x4^2);
%
    yr=zeros(n,1);
    for i=1:n
%
		tt=t(i)-t(1);
		y=0.;      
%
		if(tt>x5)
            ta=tt-x5;
			y=x1*exp(-domegan*ta)*sin(omegad*ta-x3);
            yr(i)=y;
        end
%				
		error=error+(abs(a(i)-y));
%   
    end
	error=sqrt(error);
%
	if(error<errormax)
%
				x1r(ie)=x1;
				x2r(ie)=x2;
				x3r(ie)=x3;
				x4r(ie)=x4;
				x5r(ie)=x5;
%
       out4 = sprintf(' %6ld  %13.4e  %10.4g %9.4f %9.4f %10.5f %9.4f ',j,error,x1,x2/tp,x3,x4,x5);
       disp(out4)    
%
	   errormax=error;
%
       figure(2)
       cnew=running_sum+yr';
       plot(t,original,t,cnew);
       legend ('Input Data','Synthesis',1);
       xlabel('Time(sec)');
%			
     end
%				
end
%
ave=0.;
sq=0.;
%   
for i=1:n
%
    tt=t(i)-t(1);
%		 
    domegan=x4r(ie)*x2r(ie);
    omegad=x2r(ie)*sqrt(1.-x4r(ie)^2);
%	
    if(tt>x5r(ie))
        ttt=tt-x5r(ie);
        a(i)=a(i)-x1r(ie)*exp( -domegan*ttt )*sin( omegad*ttt - x3r(ie) );
        ave=ave+a(i);
	    sq= sq+(a(i)^2.);
    end
end
%
ave=ave/n;   
sd=sqrt( (sq/n) -(ave^2.));
%	
out4=sprintf('\n  ave=%12.4g  sd=%12.4g \n',ave,sd);
disp(out4)    