sutxtaper.c

su 的源代码库

      float fac=1. ;		/* trace weighting factor */
        
	efread (&tr, 1, HDRBYTES, headerfp);
	efread (trbuf, FSIZE, nt, tracefp);
				
	/* factor for first traces */
	
	
	fac*= ( itr < tr1 ? *taperv++ : 1. );

	/* add (tr2+1)-1 to the array pointer if assymetric tapering is required 
	 * tr1=tr2  : ntap-tr1-2 = -1
	 * tr1!=tr2 : ntap-tr1-2 = tr2 */
	if ( itr == tr1 ) taperv+=ntap-tr1-2; 
	
	/* factor for last traces */
	fac*=( itr > (ntr-tr2-1) ? *taperv-- : 1. );
	
      if ( verbose && 
          (tr1+tr2) && 
	  (fac < 1.) ) fprintf(stderr,"trace %i factor %g \n",itr,fac);
       
       { register int i;       
           for (i=0;i<nt;i++)
	        tr.data[i]=fac*trbuf[i];
       }

        
	tr.unscale= 1./fac ;
	
       puttr(&tr); 	 	
  } 
  
  
  efclose(tracefp);
  efclose(headerfp);
  efclose(x2fp);
  return(CWP_Exit());
}


#define EPS     3.8090232       /* exp(-EPS*EPS) = 5e-7, "noise" level  */
				/* see sugain.c				*/
void taper ( float t1, float t2, int tap_type, float T, float dt, 
	     float *trace )
/*********************************************************************
sweep taper - tapers the sweep
**********************************************************************
Input: 
t1	  start taper in ms
t2	  end taper in ms
tap_type  type of taper to apply: 1 linear, 2 sine, 3 cosine
T	  trace duration in ms
dt	  sample rate in ms
 
Output:
trace	  array of tapered samples
*********************************************************************
This subroutine tapers a sweep mainly to reduce Gibbs phenomena.
Taper coulld be one of the specified above.
*********************************************************************
References:
Any book on Vibroseis.
*********************************************************************
Author: Tagir Galikeev				  Date:7 Oct 1994
Rewrite: Tagir Galikeev				  Date:  Oct 2002
*********************************************************************/
{
	int nt, i, nt1, nt2;
	float env=0.0, f, x;

	nt = (int)(T / dt + 1);
	nt1 = (int)(t1 / dt + 1);
	nt2 = (int)(t2 / dt + 1);
	/* apply start taper */
	if( nt1 > 1 ) {
		for (i=0; i<nt1; i++) {
	  		f = (float)i / (float)nt1;
	  		switch ((char) tap_type)	{
	  			case 1: env=f;
	  				break;
	  			case 2: env=sin(PI*f/2.);
	  				break;
	  			case 3: env=0.5*(1.0-cos(PI*f));
	  				break;
	  			case 4: x=EPS*(1-f);
	  				env=exp(-(x*x));
	  				break;
	  			case 5: x=2.0*(1-f);
	  				env=exp(-(x*x));
	  				break;
	  			default:err (" taper ?!");
	  		}
	  		trace[i] *= env;
		}
	}
	/* apply end taper */
	if( nt2 > 1 ) {
		for (i=0; i<nt2; i++) {
	  		f = (float)i / (float)nt2;
			switch ((char) tap_type)	{
	  			case 1: env=f;
	  				break;
	  			case 2: env=sin(PI*f/2.);
	  				break;
	  			case 3: env=0.5*(1.0-cos(PI*f));
	  				break;
	  			case 4: x=EPS*(1-f);
	  				env=exp(-(x*x));
	  				break;
	  			case 5: x=2.0*(1-f);
	  				env=exp(-(x*x));
	  				break;
	  			default:err (" taper ?!");
	  		}
			trace[nt-i]  *= env;
		}
	}
}


void weights ( int tr1, int tr2, float amp, float low, int type, float *w)
/*********************************************************************
sweep taper - tapers the sweep
**********************************************************************
Input: 
tr1	  number of traces to apply begin taper
tr2	  number of traces to apply end taper 
amp       maximum amplitude factor (=1.)
low	  minimum amplitude factor (=0.)
type	  type of taper to apply: 1 linear, 2 sine, 3 cosine
 
Output:
w	  array of taper weights 
*********************************************************************
This subroutine computes the taper weights 
*********************************************************************
Author: Tagir Galikeev				  Date:7 Oct 1994
Rewriten: Gerald Klein				  Date:31 Mar 2004
*********************************************************************/
{ 
   if ( tr2 && (tr1-tr2) ) { /* end taper differs from begin taper */
                register int i;
		float env=0.0, f, x;
		/* set taper weights for last traces; fill array from end */
                for (i = 0; i <= tr2; ++i) {
			f = (float) (i)/tr2;
			switch ((char) type)	{
	  			case 1: env = low + (amp - low) * f;
	  				break;
	  			case 2: env=sin(PI*f/2.);
	  				break;
	  			case 3: env=0.5*(1.0-cos(PI*f));
	  				break;
	  			case 4: x=EPS*(1-f);
	  				env=exp(-(x*x));
	  				break;
	  			case 5: x=2.0*(1-f);
	  				env=exp(-(x*x));
	  				break;
	  			default:err (" taper ?!");
	  		}
                        w[1+tr1+i] = env ;
	        } 		
   } 	
   if (tr1) { 	/* set taper weights for first traces */
                register int i;
		float env=0.0, f, x;
	           for (i = 0; i <= tr1; i++) {
			f = (float) (i)/tr1;
			switch ((char) type)	{
	  			case 1: env = low + (amp - low) * f;
	  				break;
	  			case 2: env=sin(PI*f/2.);
	  				break;
	  			case 3: env=0.5*(1.0-cos(PI*f));
	  				break;
	  			case 4: x=EPS*(1-f);
	  				env=exp(-(x*x));
	  				break;
	  			case 5: x=2.0*(1-f);
	  				env=exp(-(x*x));
	  				break;
	  			default:err (" taper ?!");
	  		}
                          w[i] = env ;
	        }
   }
   return;
}

void xweights( const float *y, float ymax, float ymin, float dy, 
               float low, float amp, int ntap, int taper_type, float *w)
/*********************************************************************
trace tapers from x2 header valus
**********************************************************************
Input: 
type	  type of taper to apply: 1 linear, 2 sine, 3 cosine
 
Output:
w	  array of taper weights 

*********************************************************************
This subroutine computes the taper weights 
*********************************************************************
Author: Gerald Klein				  Date:31 Mar 2004
*********************************************************************/
{
  register int i;
  float env=0.0, f=0., x;
  
  for (i=0;i<ntap;i++) {

      /*  Reduce weights in the range of  min-dy  to  min 
       *                        and       max     to  max+dy
       *  set to  0. for y < min-dy and to 1. for y > max+dy
       */

     if  (y[i] < ymin) {
        if (y[i] > (ymin-dy)) {
		             f = 1. - ( ymin - y[i] ) / dy ;
        } else f = 0. ;
     } else 
     if (y[i] > ymax)  {
        if (y[i] < (ymax+dy)) {
	                     f = ( ymax + dy - y[i] ) / dy  ;
        } else  f=0. ;
     } else  f=1. ;

    /* limit weight range from  low  to  amp */
    /*f = low + (amp - low) * f;*/
    
    /* apply shape of the taper corresponding to taper_type */
      switch ((char) taper_type)	{
	  			case 1: env = low + (amp - low) * f;
	  				break;
	  			case 2: env=sin(PI*f/2.);
	  				break;
	  			case 3: env=0.5*(1.0-cos(PI*f));
	  				break;
	  			case 4: x=EPS*(1-f);
	  				env=exp(-(x*x));
	  				break;
	  			case 5: x=2.0*(1-f);
	  				env=exp(-(x*x));
	  				break;
	  			default:err (" taper ?!");
	  		}
    
    w[i] = env ;

   }
  return ;
}