sugain.c

seismic software,very useful

		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	*/
	static float *epowfac;		/* exponent stretchs	*/

	if (first) {
		epowfac = ealloc1float(nt);
		for (i = 0; i < nt; i++) 
			epowfac[i] = exp(epow * (tmin + i * dt));

		first = false;
	}

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


/* Zero out outliers */
void do_trap(
	register float trap,	/* zero if magnitude > trap	*/
	register int nt		/* number of samples		*/
)
{
	int i;

	for(i=0;i<nt;i++) {
		if (ABS(tr.data[i]) > trap) tr.data[i] = 0.0;
	}

}


/* Hard clip outliers */
void do_clip(
	register float clip,	/* hard clip if magnitude > clip	*/
	register int nt		/* number of samples			*/
)
{
	int i;
	for(i=0;i<nt;i++) {
		if(tr.data[i]>clip) { 
			tr.data[i] = clip;
		} else if(tr.data[i]<-clip) {
			tr.data[i] = -clip;
		}
	}
}


/* Quantile clip on magnitudes of trace values */
void do_qclip(
	float qclip,	/* quantile at which to clip	*/
	int nt		/* number of sample points	*/
)
{
	register int i;
	register float *dataptr = tr.data;	/* ptr to trace data	*/
	static bool first = true;	/* first entry flag		*/
	static float *absdata;		/* absolute value trace		*/
	static int iq;			/* index of qclipth quantile	*/
	float clip;			/* ... value of rank[iq]	*/

	if (first) {
		absdata = ealloc1float(nt);
 		iq = (int) (qclip * nt - 0.5); /* round, don't truncate */
		first = false;
	}

	/* Clip on value corresponding to qth quantile */
	for (i = 0; i < nt; ++i)  absdata[i] = ABS(tr.data[i]);
 	clip = quant(absdata, iq, nt);
	do_clip(clip, nt);
}


/* Quantile balance */
void do_qbal(
	float qclip,	/* quantile at which to clip	*/
	int nt		/* number of sample points	*/
)
{
	register int i;
	register float *dataptr = tr.data;	/* ptr to trace data	*/
	static bool first = true;	/* first entry flag		*/
	static float *absdata;		/* absolute value trace		*/
	static int iq;			/* index of qclipth quantile	*/
	float bal;			/* value used to balance trace	*/

	if (CLOSETO(qclip, 1.0)) { /* balance by max magnitude on trace */
		bal = ABS(tr.data[0]);
		for (i = 1; i < nt; ++i)  bal = MAX(bal, tr.data[i]);

		if (CLOSETO(bal, 0.0)) {
			return;
		} else {
			for (i = 1; i < nt; ++i)  tr.data[i] /= bal;
			return;
		}
	} else if (first) {
		absdata = ealloc1float(nt);
 		iq = (int) (qclip * nt - 0.5); /* round, don't truncate */
		first = false;
	}

	/* Balance by quantile value (qclip < 1.0) */
	for (i = 0; i < nt; ++i)  absdata[i] = ABS(tr.data[i]);
	bal = quant(absdata, iq, nt);

	if (CLOSETO(bal, 0.0)) {
		return;
	} else {
		for (i = 1; i < nt; ++i)  tr.data[i] /= bal;
		do_clip(1.0, nt);
		return;
	}
}


/* Automatic Gain Control--standard box */
void do_agc(int iwagc, int nt)
{
	static bool first = true;
	static float *agcdata;
	register int i;
	double val;
	double sum;
	register int nwin;
	register float rms;


	/* allocate room for agc'd data */
	if (first) {
		first = false;
		agcdata = ealloc1float(nt);
	}


	/* compute initial window for first datum */
	sum = 0.0;
	for (i = 0; i <= iwagc; ++i) {
		val = tr.data[i];
		sum += val*val;
	}
	nwin = iwagc + 1;
	rms = sum/nwin;
	if(rms>0.) {
		agcdata[0] = tr.data[0]/sqrt(rms);
	} else {
		agcdata[0] = 0.;
	}

	/* ramping on */
	for (i = 1; i <= iwagc; ++i) {
		val = tr.data[i+iwagc];
		sum += val*val;
		++nwin;
		rms = sum/nwin;
		if(rms>0.) {
			agcdata[i] = tr.data[i]/sqrt(rms);
		} else {
			agcdata[i] = 0.;
		}
	}

	/* middle range -- full rms window */
	for (i = iwagc + 1; i <= nt - iwagc; ++i) {
		val = tr.data[i+iwagc];
		sum += val*val;
		val = tr.data[i-iwagc-1];
		sum -= val*val;
		rms = sum/nwin;
		if(rms>0.) {
			agcdata[i] = tr.data[i]/sqrt(rms);
		} else {
			agcdata[i] = 0.;
		}
	}

	/* ramping off */
	for (i = nt - iwagc + 1; i < nt; ++i) {
		val = tr.data[i-iwagc-1];
		sum -= val*val;
		--nwin;
		rms = sum/nwin;
		if(rms>0.) {
			agcdata[i] = tr.data[i]/sqrt(rms);
		} else {
			agcdata[i] = 0.;
		}
	}

	/* copy data back into trace */
	memcpy((char*)tr.data, (char*)agcdata, nt*FSIZE);

	return;
}


#define EPS	3.8090232	/* exp(-EPS*EPS) = 5e-7, "noise" level	*/

/* Automatic Gain Control--gaussian taper */
void do_gagc(int iwagc, int nt)
{
	static bool first=true;	/* first entry flag			*/
	static float *agcdata;	/* agc'd data				*/
	static float *w;	/* Gaussian window weights		*/
	static float *d2;	/* square of input data			*/
	static float *s;	/* weighted sum of squares of the data	*/
	float u;		/* related to reciprocal of std dev	*/
	float usq;		/* u*u					*/


	if (first) {
		first = false;

		/* Allocate room for agc'd data */
		agcdata = ealloc1float(nt);

		/* Allocate and compute Gaussian window weights */
		w = ealloc1float(iwagc);  /* recall iwagc is HALF window */
		u = EPS / ((float) iwagc);
		usq = u*u;
		{
			register int i;
			float floati;

			for (i = 1; i < iwagc; ++i) {
				floati = (float) i;
				w[i] = exp(-(usq*floati*floati));
			}
		}

		/* Allocate sum of squares and weighted sum of squares */
		d2 = ealloc1float(nt);
		s  = ealloc1float(nt);
	}


	/* Agc the trace */
	{
		register int i, j, k;
		register float val;
		register float wtmp;
		register float stmp;

		/* Put sum of squares of data in d2 and
		/* initialize s to d2 to get center point set */
		for (i = 0; i < nt; ++i) {
			val = tr.data[i];
			s[i] = d2[i] = val * val;
		}

		/* Compute weighted sum s; use symmetry of Gaussian */
		for (j = 1; j < iwagc; ++j) {
			wtmp = w[j];
			for (i = j; i < nt; ++i)  s[i] += wtmp*d2[i-j]; 
			k = nt - j;
			for (i = 0; i < k; ++i)   s[i] += wtmp*d2[i+j]; 
		}

		for (i = 0; i < nt; ++i) {
			stmp = s[i];
			agcdata[i] = (CLOSETO(stmp, 0.0) ?
					0.0 : tr.data[i]/sqrt(stmp));
		}

		/* Copy data back into trace */
		memcpy((char*)tr.data, (char*)agcdata, nt*FSIZE);
	}


	return;
}


/*
 * QUANT - find k/n th quantile of a[]
 *
 * Works by reordering a so a[j] < a[k] if j < k.
 *
 * Parameters:
 *    a		- data
 *    k		- indicates quantile
 *    n		- number of points in data
 *
 * This is Hoare's algorithm worked over by SEP (#10, p100) and Brian.
 */

float quant(float *a, int k, int n)
{
	register int i, j;
	int low, hi;
	register float ak, aa;

	low = 0; hi = n-1;

	while (low < hi) {
		ak = a[k];
		i = low;
		j = hi;
		do {
			while (a[i] < ak) i++;
			while (a[j] > ak) j--;
			if (i <= j) {
				aa = a[i]; a[i] = a[j]; a[j] = aa;
				i++;
				j--;
			}
		} while (i <= j);

		if (j < k) low = i;

		if (k < i) hi = j;
	}

	return(a[k]);
}