function[x1r,x2r,x3r,x4r]=wgen(num2,t,residual,duration,sr,x1r,x2r,x3r,x4r,fl,fu,nt,ie,ffmin,ffmax);
%
%  wgen.m   version 2.1   March 30, 2011
%
%
    fu=ffmax;
    fl=ffmin;
%
tp=2*pi;
ave=0.;
sq=0.;
%
jki=0;
jflag=1;
%
ave=mean(residual);   
sd=std(residual);
%
for(i=1:num2)
    asd(i)=(residual(i)^2.);
end
%
out1=sprintf(' ave=%12.4g  sd=%12.4g  nt=%d',ave,sd,nt);
disp(out1);
%
am=2.*sd;
n=num2;
%	       
errormax=1.0e+53;
%
disp(' ');
disp('  Trial     Error      Amplitude   Freq(Hz)   NHS    delay(sec) ');
%
noct=log(fu/fl)/log(2)
%
hnt=0.5*nt;
%
    tp=2*pi;               % needed
    tp_ffmax=tp*ffmax;     % needed
%
for(j=1:nt)
%
    ran1=rand;
    ran2=rand;	
    ran3=rand;
%				
    x1=rand;
    x2=rand;
    x33=rand;
    x4=rand;
%
    x1=2.*am*(x1-0.5);		% amplitude
%
    if(rand>0.7)
        x2=((fu-fl)*x2+fl);	% freq  
%%%                    out1=sprintf('lin  x2=%8.4g ',x2);
%%%            disp(out1)
    else
        x2=x2*noct;
        x2=fl*exp(x2*log(2));
%%%                   out1=sprintf('log  x2=%8.4g ',x2);
%%%            disp(out1)
    end
%
    x2=x2*tp;
%
    x3= 3+(2*round(x33*30));	% nhs
%
    x4=x4*0.8*duration;		% delay
%
    itype=1;
%
	if(ran1>0.4 && ran1<=0.5 && j>hnt)
%						
	    x2=x2r(ie)*(0.99+0.02*rand);
	    x4=x4r(ie)*(0.99+0.02*rand);
%
	    if(ran2<=0.25)
		    x3=x3r(ie)-4;
        end    
		if(ran2>0.25 && ran2<=0.5)
		    x3=x3r(ie)-2;
        end
	    if(ran2>0.50 && ran2<=0.75)
			x3=x3r(ie)+2;
        end
		if(ran2>0.75 && ran2<=1.0)
		    x3=x3r(ie)+4;
        end
%
		x1=x1r(ie)*(0.95+0.10*rand);
%
	    if(ran3>0.5)
			x1=-x1;
        end
%
		if(x3<3)
            x3=3;
        end
%
		itype=2;  % mainly NHS
%
    end
%   
	if(ran1>0.5 && ran1<=0.6 && j>hnt)  
%			
				x1=x1r(ie)*(0.98+0.04*rand);
				x2=x2r(ie);
				x3=x3r(ie);
				x4=x4r(ie);
%
				itype=3; % amp  
    end
	if(ran1>0.6 && ran1<=0.7 && j>hnt)
%			
				x1=x1r(ie);
				x2=x2r(ie)*(0.99+0.02*rand);	
				x3=x3r(ie);
				x4=x4r(ie);
%				
				itype=4;  % freq
    end
	if(ran1>0.8 && ran1<=0.9 && j>hnt)
%			
				x1=x1r(ie);
				x2=x2r(ie);
				x3=x3r(ie);
				x4=x4r(ie)*(0.99+0.02*rand);
%
				itype=5;  % delay
    end
	if(ran1>0.9 && ran1<=1. && j>hnt)
%			
				x1=x1r(ie)*(0.999+0.002*rand);
				x2=x2r(ie)*(0.999+0.002*rand);
				x3=x3r(ie);
				x4=x4r(ie)*(0.999+0.002*rand);
%
				itype=6;  % all but NHS
    end    
%
	while(1)
%
%    x2=fn
%    x3=nhs
%    x4=delay
%
	    if( tp*x3/(2.*x2) + x4 < duration )
			break;
		else
			x3=x3-2;
        end
    end      
%            
	if(x3==1)
				x1=0;
				x2=fu*tp;
				x3=3;
				x4=0;
                disp(' ref 1');
%
    end
%
    if(j==1 || x3 < 3)
				x1=0.;
				x2=fu*tp;
				x3=3;
				x4=0.;
               disp(' ref 2');            
    end
%
    error=0.;
%
    t1=x4 + t(1);
    t2=t1 + tp*x3/(2.*x2);
%
    if(j==1)
            x1=0.;
    end
%
    fmin=(1/duration)*1.75;
    if(fmin<ffmin)
        fmin=ffmin;
    end
%
    if( (x2/tp)<fmin || x3 < 3)
        out1=sprintf('(x2/tp)= %8.4g  fmin=%8.4g  x3=%d ',(x2/tp),fmin,x3);
        disp(out1);
 %
        x2=fmin*tp;
        x3=3;
        x4=0.;
                   disp(' ref 3');    
    end
%
    if( x2>tp_ffmax)
            out1=sprintf('R4 x2/tp=%8.4g  ffmax=%8.4g',x2/tp,ffmax);
            disp(out1)
        x2=((ffmax-ffmin)*rand + ffmin)*tp;
        x3=3;
        x4=0.;
       disp(' ref 4');    
    end
%
	for(i=1:num2)
%		
		    tt=t(i);
%
			y=0.;
%
		    if( tt>= t1 && tt <= t2)
%				
			    arg=x2*(tt-t1);  
%
				y=x1*sin(arg/double(x3))*sin(arg);
%
				error=error+((residual(i)-y)^2.);
%
			else
			    error=error+asd(i);
            end
    end          
%            
    error=sqrt(error);
%
    if(error<errormax && x2>=(fmin*tp))
%
				x1r(ie)=x1;
				x2r(ie)=x2;
				x3r(ie)=x3;
				x4r(ie)=x4;
%
				out1=sprintf(' %d  %11.4e  %9.4f  %9.4f  %d  %9.4f ',j,error,x1,x2/tp,x3,x4);
                disp(out1);
%
				errormax=error;
%
%%                out1=sprintf(' x3r[%d]=%d \n',ie,x3r(ie));
%%                disp(out1);
    end
%
    if(x2r(ie)<ffmin*tp)
        x2r(ie)=ffmin*tp;
        x1r(ie)=1.0e-20;
        x3r(ie)=3;
        x4r(ie)=0;
    end
%
end