sugain.c

seismic software,very useful

/* SUGAIN: $Revision: 1.14 $ ; $Date: 90/12/22 08:14:32 $		*/

/*----------------------------------------------------------------------
 * Copyright (c) Colorado School of Mines, 1990.
 * All rights reserved.
 *
 * This code is part of SU.  SU stands for Seismic Unix, a processing line
 * developed at the Colorado School of Mines, partially based on Stanford
 * Exploration Project (SEP) software.  Inquiries should be addressed to:
 *
 *  Jack K. Cohen, Center for Wave Phenomena, Colorado School of Mines,
 *  Golden, CO 80401  ([e-mail])
 *----------------------------------------------------------------------
 */

#include "su.h"
#include "segy.h"
#include "subc.h"

/*********************** self documentation *****************************/
string sdoc =
" 									\n"
" SUGAIN - apply various types of gain to display traces		\n"
" 									\n"
" sugain <stdin >stdout [optional parameters]				\n"
" 							        	\n"
" Required parameters:							\n"
" 	none (no-op)							\n"
" 							        	\n"
" Optional parameters: 							\n"
" 	tpow=0.0	multiply data by t^tpow		    		\n"
" 	epow=0.0	multiply data by exp(epow*t)	       		\n"
" 	gpow=1.0	take signed gpowth root of scaled data		\n"
" 	agc=0		flag; 1 = do automatic gain control		\n"
" 	gagc=0		flag; 1 = ... with gaussian taper		\n"
" 	wagc=0.5	agc window in seconds (use if agc=1 or gagc=1)	\n"
" 	trap=0.0	zero any value whose magnitude exceeds trapval 	\n"
" 	clip=0.0	clip any value whose magnitude exceeds clipval	\n"
" 	qclip=1.0	clip by quantile on absolute values on trace	\n"
" 	qbal=0		flag; 1 = balance traces by qclip and scale	\n"
" 	pbal=0		flag; 1 = bal traces by dividing by rms	value	\n"
" 	jon=0		flag; 1 means tpow=2, gpow=.5, qclip=.95	\n"
"       t=0.0           time (in ms) where scale is to be applied	\n" 
" 	scale=1.0	scale applied at t      			\n"
"                       (e.g., t=0.,1000,4000,6000 scale=1.,1.5,1.7,1.0)\n"
" 							        	\n"
" Operation order:							\n"
" 							        	\n"
" out(t) = scale(t)*BAL{CLIP[AGC{[t^tpow * exp(epow * t) * in(t)]^gpow}]}\n"
" 							        	\n"
" Notes:								\n"
" 	The jon flag selects the parameter choices discussed in		\n"
" 	Claerbout's Imaging the Earth, pp 233-236.   			\n"
" 	Selected traces can be marked as not to be gained by		\n"
" 	using sumark to set the tr.mark header field to 1.		\n"
" 							        	\n"
;
/**************** end self doc *******************************************/

/* Credits:
 *	SEP: Jon
 *	CWP: Jack, Brian, Dave
 *            : Zhiming Li (add time varying scale)
 *
 * Technical Reference:
 *	Jon's second book, pages 233-236.
 */

/* subroutine prototypes */
void do_tpow(float tpow, register float tmin, register float dt, int nt);
void do_epow(float epow, register float tmin, register float dt, int nt);
void do_trap(register float trap, register int nt);
void do_clip(register float clip, register int nt);
void do_qclip(float qclip, int nt);
void do_qbal(float qclip, int nt);
void do_agc(int iwagc, int nt);
void do_gagc(int iwagc, int nt);
float quant(float *a, int k, int n);


#define	TPOW     0.0
#define	EPOW     0.0
#define	GPOW     1.0
#define	TRAP     0.0
#define	CLIP     0.0
#define	QCLIP    1.0
#define	SCALE    1.0
#define	WAGC     0.5
#define	AGC      0
#define	GAGC     0
#define	PBAL 	 0
#define	QBAL 	 0
#define	JON 	 0


segy tr;


main(int argc, char **argv)
{
	int jon;	/* flag to get Claerbout values			*/
	int agc;	/* agc flag					*/
	int gagc;	/* gaussian agc flag				*/
	int pbal;	/* power balance flag				*/
	int qbal;	/* quantile balance flag			*/
	float tpow;	/* exponent of t 				*/
	float epow;	/* deattenutation coefficient			*/
	float gpow;	/* dynamic compression power			*/
	float rmsq;	/* root mean square of a trace			*/
	float trap;	/* zero any larger value magnitude than trapval	*/
	float clip;	/* clip any larger value magnitude than clipval	*/
	float qclip;	/* clip at qth quantile (100qth percentile)	*/
	float wagc;	/* size of agc window in seconds		*/
	int iwagc;	/* ... half window in samples			*/
	int nt;		/* number of samples on trace			*/
	float tmin;	/* delay recording time in secs			*/
	float dt;	/* sample rate in secs				*/
	float *scale;	/* scale factor					*/
	float *t;	/* time to apply scale 				*/
	int nscale; 	/* number of t,scale pairs			*/
	float *scalet;   /* scale at every time sample 			*/
	int itmp, one=1, it;
	float tmp;
	FILE *infp=stdin, *outfp=stdout;


	/* Initialize */
	initargs(argc, argv);
	askdoc(1);
	
	/* large file extension */
	file64(infp);
	file64(outfp);	

	/* Get nt from first trace */
	if (!gettr(&tr)) err("can't get first trace");
	nt   = (int) tr.ns;
	dt   = (float) tr.dt/1000000.0;   /* microsecs to secs */
	tmin = (float) tr.delrt/1000.0;   /* millisecs to secs */
	if (!dt) getparfloat("dt", &dt);
	if (!dt) MUSTFGETPAR("dt", &dt);


	/* Get parameters */
	if (!getparfloat ("tpow" , &tpow))	tpow     = TPOW;
	if (!getparfloat ("epow" , &epow))	epow     = EPOW;
	if (!getparfloat ("gpow" , &gpow))	gpow     = GPOW;
	if (!getparfloat ("trap" , &trap))	trap     = TRAP;
	if (!getparfloat ("clip" , &clip))	clip     = CLIP;
	if (!getparfloat ("qclip", &qclip))	qclip    = QCLIP;
	if (!getparfloat ("wagc" , &wagc))	wagc     = WAGC;
	if (!getparint   ("agc"  , &agc))	agc      = AGC;
	if (!getparint   ("gagc" , &gagc))	gagc     = GAGC;
	if (!getparint   ("pbal" , &pbal))	pbal 	 = PBAL;
	if (!getparint   ("qbal" , &qbal))	qbal 	 = QBAL;
	if (!getparint   ("jon"  , &jon))	jon 	 = JON;

	nscale = countparval("scale");
	itmp = countparval("t");
	if(itmp!=nscale) err("check t and scale ");
	if(nscale>0) {
		scale = (float*) malloc(nscale*sizeof(float));
		t = (float*) malloc(nscale*sizeof(float));
		getparfloat("scale", scale);
		getparfloat("t", t);
		scalet = (float*) malloc(nt*sizeof(float));	
		for(it=0;it<nscale;it++) t[it] = t[it] * 0.001;
		for(it=0;it<nt;it++) {
			if(nscale==1) {
				scalet[it] = scale[0];
			} else {
				tmp = tmin + it*dt;
				lin1d_(t,scale,&nscale,&tmp,
					scalet+it,&one,&itmp);
			}
		}
	}


	/* Data validation */
	if (trap < 0.0) err("trap = %f, must be positive", trap);
	if (clip < 0.0) err("clip = %f, must be positive", clip);
	if (qclip < 0.0 || qclip > 1.0) 
		err("qclip = %f, must be between 0 and 1", qclip);
	if (agc || gagc) {
		iwagc = NINT(wagc/dt);
		if (iwagc < 1) err("wagc=%g must be positive", wagc);
		if (iwagc > nt) err("wagc=%g too long for trace", wagc);
		iwagc >>= 1;  /* windows are symmetric, so work with half */
	}
	if (jon) { 
		tpow  = 2.0;
		gpow  = 0.5;
		qclip = 0.95;
	}


	/* Main loop over traces */
	do {
		if (!tr.mark) {
			if (!CLOSETO(tpow, 0.0)) {
				do_tpow(tpow, tmin, dt, nt);
			}
			if (!CLOSETO(epow, 0.0)) {
				do_epow(epow, tmin, dt, nt);
			}
			if (!CLOSETO(gpow, 1.0)) {
				register int i;
				register float val;

				if (CLOSETO(gpow, 0.5)) {
					for (i = 0; i < nt; ++i) {
						val = tr.data[i];
						tr.data[i] = (val >= 0.0) ?
							sqrt(val) : -sqrt(-val);
					}
				} else if (CLOSETO(gpow, 2.0)) {
					for (i = 0; i < nt; ++i) {
						val = tr.data[i];
						tr.data[i] = val * ABS(val);
					}
				} else {
					for (i = 0; i < nt; ++i) {
						val = tr.data[i];
						tr.data[i] = (val >= 0.0) ?
							 pow(val, gpow) :
							-pow(-val, gpow);
					}
				}
			}
			if (agc)		   do_agc(iwagc, nt);
			if (gagc)		   do_gagc(iwagc, nt);
			if (trap > 0.0)		   do_trap(trap, nt);
			if (clip > 0.0)		   do_clip(clip, nt);
			if (qclip < 1.0 && !qbal)  do_qclip(qclip, nt);
			if (qbal)		   do_qbal(qclip, nt);
			if (pbal) {
				register int i;
				register float val;

				/* rmsq = sqrt (SUM( a()*a() ) / nt) */
				rmsq = 0.0;
				for (i = 0; i < nt; ++i) {
					val = tr.data[i];
					rmsq += val * val;
				}
				rmsq = sqrt(rmsq / nt);

				if (!CLOSETO(rmsq, 0.0)) {
					for (i = 0; i < nt; ++i)
						tr.data[i] /= rmsq;
				}
			}
			if (nscale>0) {
				register int i;
				for (i=0;i< nt;++i)  tr.data[i] *= scalet[i];
			}
		}
		puttr(&tr);

	} while(gettr(&tr));


	return EXIT_SUCCESS;
}


/* Multiply by t^tpow */
void do_tpow(
	float tpow,		/* multiply data by t^tpow	*/
	register float tmin,	/* first time on record		*/
	register float dt,	/* sampling rate in seconds	*/
	int nt			/* number of samples		*/
)
{
	static bool first = true;	/* first entry flag	*/
	static float *tpowfac;		/* tpow values		*/
	register int i;			/* counter		*/

	if (first) { /* first entry, set up array of tpow factors */
		register float t;

		tpowfac = ealloc1float(nt);
		if(tpow==0.) {
			for (i = 0; i < nt; ++i)  tpowfac[i] = 1.;
		} else {
			for (i = 0; i < nt; ++i)  {
				t = tmin + i*dt;
				tpowfac[i] = (t!=0.)? pow(fabs(t),tpow): 0.;
			/*
			tpowfac[i] = (t >= 0.0) ? pow(t, tpow) : -pow(-t, tpow);
			*/
			}
		}
		first = false;
	} /* end first entry */

	for (i = 0; i < nt; ++i)  tr.data[i] *= tpowfac[i];
}


/* Exponential deattenuation  with deattenuation factor epow */
void do_epow(
	float epow,		/* coefficient of t in exponent	*/
	register float tmin,	/* first time on record		*/
	register float dt,	/* sampling rate in seconds	*/
	int nt			/* number of samples		*/
)
{
	register int i;			/* counter		*/
	static bool first = true;	/* first entry flag	*/