suop.c

su 的源代码库

/* Copyright (c) Colorado School of Mines, 2006.*/
/* All rights reserved.                       */

/* SUOP: $Revision: 1.24 $ ; $Date: 2006/07/10 19:16:24 $		*/

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

/*********************** self documentation **********************/
char *sdoc[] = {
" 								",
" SUOP - do unary arithmetic operation on segys 		",
" 								",
" suop <stdin >stdout op=abs					",
" 								",
" Required parameters:						",
"	none							",
"								",
" Optional parameter:						",
"	op=abs		operation flag				",
"			abs   : absolute value			",
"			avg   : remove average value		",
"			ssqrt : signed square root		",
"			sqr   : square				",
"			ssqr  : signed square			",
"			sgn   : signum function			",
"			exp   : exponentiate			",
"			slog  : signed natural log		",
"			slog10: signed common log		",
"			cos   : cosine				",
"			sin   : sine				",
"			tan   : tangent				",
"			cosh  : hyperbolic cosine		",
"			sinh  : hyperbolic sine			",
"			tanh  : hyperbolic tangent		",
"			norm  : divide trace by Max. Value	",
"			db    : 20 * slog10 (data)		",
"			neg   : negate value			",
"			posonly : pass only positive values	",
"			negonly : pass only negative values	",
"                       sum   : running sum trace integration   ",
"                       diff  : running diff trace differentiation",
"                       refl  : (v[i+1] - v[i])/(v[i+1] + v[i]) ",
"			mod2pi : modulo 2 pi			",
"			inv   : inverse				",
"			nop   : no operation			",
"								",
" Note:	Binary ops are provided by suop2.			",
" Operations inv, slog and slog10 are \"punctuated\",		", 
" meaning that if, the input contains 0 values,			",
" 0 values are returned.					",	
"								",
" For file operations on non-SU format binary files use:  farith",
NULL};

/* Credits:
 *
 *	CWP: Shuki, Jack
 *	Toralf Foerster: norm and db operations, 10/95.
 *
 * Notes:
 *	If efficiency becomes important consider inverting main loop
 *      and repeating operation code within the branches of the switch.
 *
 *	Note on db option.  The following are equivalent:
 *	... | sufft | suamp | suop op=norm | suop op=slog10 |\
 *		sugain scale=20| suxgraph style=normal
 *
 *	... | sufft | suamp | suop op=db | suxgraph style=normal
 */
/**************** end self doc ***********************************/

#ifndef	TWOPI
#define	TWOPI (2.0 * PI)
#endif


#define	FABS	1
#define	SSQRT	2
#define	SQR	3
#define	SSQR	4
#define EXP	5
#define SLOG	6
#define SLOG10	7
#define COS	8
#define SIN	9
#define TAN	10
#define COSH	11
#define SINH	12
#define TANH	13
#define NORM    14	/* normalize at maximal value   */
#define DB      15	/* for frequency plot in db     */
#define SIGN	16	/* signum function		*/

#define NEG	17
#define NOP	18
#define ONLY_POS 19
#define ONLY_NEG 20
#define SUM    21
#define REFL   22
#define DIFF   23
#define MOD_2PI 24
#define INV	25
#define AVG	26

segy tr;

int
main(int argc, char **argv)
{
	cwp_String op="abs";	/* operation: abs, exp, ..., 		*/
	int iop=FABS;		/* integer abbrev. for op in switch	*/
	int nt;			/* number of samples on input trace	*/
	float dt;		/* time sampling interval */
	float twobydt;		/* 2*dt */


	/* Initialize */
	initargs(argc, argv);
	requestdoc(1);


	/* Get information from first trace */
	if (!gettr(&tr)) err("can't get first trace");
	nt = tr.ns;
	dt = tr.dt/1000000.0;
	twobydt=2*dt;


	/* Get operation, recall iop initialized to the default FABS */
	getparstring("op", &op);
	if      (STREQ(op, "ssqrt"))	iop = SSQRT;
	else if (STREQ(op, "sqr"))	iop = SQR;
	else if (STREQ(op, "ssqr"))	iop = SSQR;
	else if (STREQ(op, "sgn"))	iop = SIGN;
	else if (STREQ(op, "exp"))	iop = EXP;
	else if (STREQ(op, "slog"))	iop = SLOG;
	else if (STREQ(op, "slog10"))	iop = SLOG10;
	else if (STREQ(op, "cos"))	iop = COS;
	else if (STREQ(op, "sin"))	iop = SIN;
 	else if (STREQ(op, "tan"))	iop = TAN;
	else if (STREQ(op, "cosh"))	iop = COSH;
	else if (STREQ(op, "sinh"))	iop = SINH;
	else if (STREQ(op, "tanh"))	iop = TANH;
	else if (STREQ(op, "norm"))     iop = NORM;
	else if (STREQ(op, "db"))       iop = DB;
	else if (STREQ(op, "neg"))      iop = NEG;
	else if (STREQ(op, "nop"))      iop = NOP;
	else if (STREQ(op, "posonly"))  iop = ONLY_POS;
	else if (STREQ(op, "negonly"))  iop = ONLY_NEG;
	else if (STREQ(op, "sum"))      iop = SUM;
	else if (STREQ(op, "refl"))     iop = REFL;
	else if (STREQ(op, "diff"))     iop = DIFF;
	else if (STREQ(op, "mod2pi"))  	iop = MOD_2PI;
	else if (STREQ(op, "inv"))  	iop = INV;
 	else if (STREQ(op, "avg"))  	iop = AVG;
	else if (!STREQ(op, "abs"))
		err("unknown operation=\"%s\", see self-doc", op);


	/* Main loop over traces */
	do {

		switch(iop) { register int i;
		case FABS:
			for (i = 0; i < nt; ++i)
				tr.data[i] = ABS(tr.data[i]);
		break;
		case SSQRT:
			for (i = 0; i < nt; ++i) {
				float x = tr.data[i];
				tr.data[i] = SGN(x) * sqrt(ABS(x));
			}
		break;
		case SQR:
			for (i = 0; i < nt; ++i) {
				float x = tr.data[i];
				tr.data[i] = x * x;
			}
		break;
		case SSQR:
			for (i = 0; i < nt; ++i) {
				float x = tr.data[i];
				tr.data[i] = SGN(x) * x * x;
			}
		break;
		case SIGN:
			for (i = 0; i < nt; ++i) {
				float x = tr.data[i];
				tr.data[i] = SGN(x);
			}
		break;
		case EXP:
			for (i = 0; i < nt; ++i)
				tr.data[i] = exp(tr.data[i]);
		break;
		case SLOG:
			for (i = 0; i < nt; ++i) {
				float x = tr.data[i];

				if (ABS(x) > 0) {
					tr.data[i] = SGN(x) * log(ABS(x));
				} else {
					tr.data[i] = 0;
				}
			}
		break;
		case SLOG10:
			for (i = 0; i < nt; ++i) {
				float x = tr.data[i];
				if (ABS(x) > 0) {
					tr.data[i] = SGN(x) * log10(ABS(x));
				} else {
					tr.data[i] = 0;
				}
					
			}
		break;
		case COS:
			for (i = 0; i < nt; ++i)
				tr.data[i] = cos(tr.data[i]);
		break;
		case SIN:
			for (i = 0; i < nt; ++i)
				tr.data[i] = sin(tr.data[i]);
		break;
		case TAN:
			for (i = 0; i < nt; ++i)
				tr.data[i] = tan(tr.data[i]);
		break;
		case COSH:
			for (i = 0; i < nt; ++i)
				tr.data[i] = cosh(tr.data[i]);
		break;
		case SINH:
			for (i = 0; i < nt; ++i)
				tr.data[i] = sinh(tr.data[i]);
		break;
		case TANH:
			for (i = 0; i < nt; ++i)
				tr.data[i] = tanh(tr.data[i]);
		break;
		case NORM:
		      { float x, max;
			max = 0.0;
			for (i = 0; i < nt; ++i) {
				x = ABS (tr.data [i]);
				if (max < x) max = x;
			}
			if (max != 0.0)
				for (i = 0; i < nt; ++i) tr.data [i] /= max;
		      }	/* end scope of x, max */
		break;
 
		case DB:
			if (iop == DB)  {
				float x;
				for (i = 0; i < nt; ++i) {
					x = tr.data[i];
					tr.data[i] = (ABS(x) > 0) ?
					   20.0*SGN(x)*log10(ABS(x)) : 0;
				}
			}
		break;
		
		case NEG:
			for (i = 0; i < nt; ++i)
				tr.data[i] = -tr.data[i];
		break;
		case NOP:
		break;
		case ONLY_POS:
		{
			float x;		     
	       		for (i = 0; i < nt; ++i) {   
				x = tr.data[i];      
				tr.data[i] = (x > 0) ? x : 0;  
			}
		}			      
		break; 
		case ONLY_NEG:
		{
			float x;		     
	       		for (i = 0; i < nt; ++i) {   
				x = tr.data[i];      
				tr.data[i] = (x < 0) ? x : 0;  
			}
		}			      
		break;
		case SUM:
		{
	       		for (i = 1; i < nt; ++i) {   
				tr.data[i] += tr.data[i-1];
			}
		}			      
		break; 
		case REFL:
		{
	       		for (i = nt-1; i > 0; --i) {   
				float numer=(tr.data[i] - tr.data[i-1]);
				float denom=(tr.data[i] + tr.data[i-1]);
				if (denom == 0.0) denom = 1.0;
				tr.data[i] = ( numer/denom ) ;
			}
			tr.data[0] = 0.0;
		}			      
		break;
		case DIFF:
		{
			float *datatmp=NULL; /* temporary data storage */

			/* allocate space for temporary data */
			datatmp = ealloc1float(nt);

			/* copy data from tr.data */
			memcpy((void *)datatmp, (const void *) tr.data, nt*FSIZE);
			/* do centered differences for the rest */
	       		for (i = 2; i < nt-2; ++i) {   
				float numer;
				numer=(datatmp[i+1] - datatmp[i-1]);
				tr.data[i] = ( numer/twobydt ) ;
			}
	
			/* simple diffrence for tr.data[0] */
			tr.data[0] = (datatmp[1] - datatmp[0])/dt;
			tr.data[nt-1] = (datatmp[nt-1] - datatmp[nt-2])/dt;

			/* centered difference for tr.data[1] */
			tr.data[1] = (datatmp[2] - datatmp[0])/twobydt;
			tr.data[nt-2] = (datatmp[nt-1] - datatmp[nt-3])/twobydt;
		}			      
		break;

 		case AVG:
 		        {
 			float sum = 0;float avg = 0 ;
 		        for (i = 0; i < nt; ++i)
 			        sum = sum + tr.data[i];
 			avg = sum / nt ;
 			for (i = 0; i < nt; ++i)
 				tr.data[i] = tr.data[i] - avg;
 			}
 		break;

		case MOD_2PI:
			for (i = 0; i < nt; ++i)
			{	while(tr.data[i]<0) tr.data[i] += TWOPI;
				while(tr.data[i]>=TWOPI) tr.data[i] -= TWOPI;
		        }
		break;
		case INV:
			for (i = 0; i < nt; ++i)
				if (tr.data[i]) tr.data[i] = 1.0 / tr.data[i];
		break;

		default:  /* defensive programming */
			err("mysterious operation=\"%s\"", op);
		} /* end scope of i */
		
		
		
		puttr(&tr);

	} while (gettr(&tr));


	return(CWP_Exit());
}