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.2   October 28, 2011
%
nth=length(residual);
%
tp=2*pi;
ave=0.;
sq=0.;
%
jki=0;
jflag=1;
%
ave=mean(residual);   
sd=std(residual);
%
for(i=1:nth)
    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);
%
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.4)
        nq=noct*rand();
        x2=fl*(2^nq);
    else
        x2=((fu-fl)*x2+fl);	% freq 
    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>100)
%						
	    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>100)  
%			
				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>100)
%			
				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>100)
%			
				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>100)
%			
				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)
%
	    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;
%
    end
%
    if(j==1 || x3 < 3)
				x1=0.;
				x2=fu*tp;
				x3=3;
				x4=0.;
    end
%
    error=0.;
%
    t1=x4 + t(1);
    t2=t1 + tp*x3/(2.*x2);
%
    if(j==1)
            x1=0.;
    end
%
    fmin=duration/1.75;
    if(fmin<ffmin)
        fmin=ffmin;
    end
%
    if( (x2/tp)<fmin || x3 < 3)
        x2=fmin*tp;
        x3=3;
        x4=0.;
    end
%
    if( (x2/tp)>ffmax)
        x2=((ffmax-ffmin)*rand + ffmin)*tp;
        x3=3;
        x4=0.;
    end
%
%   if( tt>= t1 && tt <= t2)
%
    T=t(nth)-t(1);
%    
    x=t1-t(1);    
    c2=x/T;
    c1=1-c2;
%   
    Rn1=1*c1+nth*c2;
    n1=ceil(Rn1); 
%
    x=t2-t(1);    
    c2=x/T;
    c1=1-c2;
    Rn2=1*c1+nth*c2;
    n2=floor(Rn2);
    if( (n2+1)>nth)
        n2=nth-1;
    end
%
    error=sum(asd(1:(n1-1)))+sum(asd((n2+1):nth));
%
    y=zeros(nth,1);
    arg=zeros(nth,1);
    arg(n1:n2)=x2*(t(n1:n2)-t1);
    beta=arg/double(x3);
%    
    y(n1:n2)=x1*sin(beta(n1:n2)).*sin(arg(n1:n2));
%
	error=error+sum((residual(n1:n2)-y(n1:n2)).^2.);
%            
    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