suattributes.c

su 的源代码库

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

/* SUATTRIBUTES:  $Revision: 1.27 $ ; $Date: 2005/03/09 17:29:26 $	*/

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

/*********************** self documentation **********************/
char *sdoc[] = {
" 									",
" SUATTRIBUTES - instantaneous trace ATTRIBUTES 			",
" 		 							",
" 									",
" suattributes <stdin >stdout mode=amp					",
" 									",
" Required parameters:							",
" 	none								",
" 									",
" Optional parameter:							",
" 	mode=amp	output flag 					",
" 	       		=amp envelope traces				",
" 	       		=phase phase traces				",
" 	       		=freq frequency traces				",
"			=freqw Frequency Weighted Envelope		",
"			=thin  Thin-Bed (inst. freq - average freq)	",
"			=bandwith Instantaneous bandwidth		",
"			=normamp Normalized Phase (Cosine Phase)	",
" 	       		=fdenv 1st envelope traces derivative		",
" 	       		=sdenv 2nd envelope traces derivative		",
" 	       		=q Ins. Q Factor				",
"	unwrap=		default unwrap=0 for mode=phase			",
" 			default unwrap=1 for mode=freq			",
" 			dphase_min=PI/unwrap				",
"	wint=		windowing for freqw				",
"			windowing for thin				",
"			default=1 					",
" 			o--------o--------o				",
" 			data-1	data	data+1				",
" 									",
" Notes:								",
" This program performs complex trace attribute analysis. The first three",
" attributes, amp,phase,freq are the classical Taner, Kohler, and	",
" Sheriff, 1979.							",
" 									",
" 									",
" The unwrap parameter is active only for mode=freq and mode=phase. The	",
" quantity dphase_min is the minimum change in the phase angle taken to be",
" the result of phase wrapping, rather than natural phase variation in the",
" data. Setting unwrap=0 turns off phase-unwrapping altogether. Choosing",
" unwrap > 1 makes the unwrapping function more sensitive to phase changes.",
" Setting unwrap > 1 may be necessary to resolve higher frequencies in	",
" data (or sample data more finely). The phase unwrapping is crude. The ",
" differentiation needed to compute the instantaneous frequency		",
" freq(t)= d(phase)/dt is a simple centered difference.			",
"	 					       			",
" Examples:								",
" suvibro f1=10 f2=50 t1=0 t2=0 tv=1 | suattributes2 mode=amp | ...	",
" suvibro f1=10 f2=50 t1=0 t2=0 tv=1 | suattributes2 mode=phase | ...	",
" suvibro f1=10 f2=50 t1=0 t2=0 tv=1 | suattributes2 mode=freq | ...	",
" suplane | suattributes mode=... | supswigb |...       		",

NULL};

/* Credits:
 *	CWP: Jack Cohen
 *      CWP: John Stockwell (added freq and unwrap features)
 *	UGM (Geophysics Students): Agung Wiyono
 *           email:[e-mail] (others)
 *					
 *
 * Algorithm:
 *	c(t) = hilbert_tranform_kernel(t) convolved with data(t)  
 *
 *  amp(t) = sqrt( c.re^2(t) + c.im^2(t))
 *  phase(t) = arctan( c.im(t)/c.re(t))
 *  freq(t) = d(phase)/dt
 *
 * Reference: Taner, M. T., Koehler, A. F., and  Sheriff R. E.
 * "Complex seismic trace analysis", Geophysics,  vol.44, p. 1041-1063, 1979
 *
 * Trace header fields accessed: ns, trid
 * Trace header fields modified: d1, trid

 */
/**************** end self doc ********************************/

#define	AMP		 1
#define	ARG		 2
#define	FREQ		 3
#define BANDWIDTH	 4
#define NORMAMP		 5
#define FREQW		 6
#define THIN		 7
#define FENV		 8
#define SENV		 9
#define Q		 10

/* function prototype of functions used internally */
void unwrap_phase(int n, float unwrap, float *phase);
void differentiate(int n, float h, float *f);
void twindow(int nt, int wtime, float *data);

segy tr;

int
main(int argc, char **argv)
{
	cwp_String mode;	/* display: real, imag, amp, arg	*/
	int imode=AMP;		/* integer abbrev. for mode in switch	*/
	register complex *ct;	/* complex trace			*/
	int nt;			/* number of points on input trace	*/
	float dt;		/* sample spacing			*/
	float *data;		/* array of data from each trace	*/
	float *hdata;		/* array of Hilbert transformed data	*/
	float unwrap;		/* PI/unwrap=min dphase assumed to by wrap*/
	int wint;		/* n time sampling to window */
	cwp_Bool seismic;	/* is this seismic data?		*/
	int ntout;
	/* Initialize */
	initargs(argc, argv);
	requestdoc(1);
	

	/* Get info from first trace */
	if (!gettr(&tr)) err("can't get first trace");
	nt = tr.ns;
	dt = ((double) tr.dt)/1000000.0;
	ntout = nt + nt -1;

	/* check to see if data type is seismic */
	seismic = ISSEISMIC(tr.trid);

	if (!seismic)
		warn("input is not seismic data, trid=%d", tr.trid);

	/* Get mode; note that imode is initialized to AMP */
	if (!getparstring("mode", &mode))	mode = "amp";

	if      (STREQ(mode, "phase"))  imode = ARG;
	else if (STREQ(mode, "freq"))	imode = FREQ;
	else if (STREQ(mode, "bandwidth")) imode = BANDWIDTH;
	else if (STREQ(mode, "normamp")) imode = NORMAMP;
	else if (STREQ(mode, "freqw")) imode = FREQW;
	else if (STREQ(mode, "thin")) imode = THIN;
	else if (STREQ(mode, "fdenv")) imode = FENV;
	else if (STREQ(mode, "sdenv")) imode = SENV;
	else if (STREQ(mode, "q")) imode = Q;
	else if (!STREQ(mode, "amp"))
		err("unknown mode=\"%s\", see self-doc", mode);

	/* getpar value of unwrap */
	switch(imode) {
	case FREQ:
		if (!getparfloat("unwrap", &unwrap))	unwrap=1;
	break;
	case Q:
		if (!getparfloat("unwrap", &unwrap))	unwrap=1;
	break;
	case FREQW:
		if (!getparfloat("unwrap", &unwrap))	unwrap=1;
		if (!getparint("wint", &wint))	wint=3; 
	break;
	case THIN:
		if (!getparfloat("unwrap", &unwrap))	unwrap=1;
		if (!getparint("wint", &wint))	wint=3;
	break;
	case ARG:
		if (!getparfloat("unwrap", &unwrap))	unwrap=0;
	break;
	}

	/* allocate space for data and hilbert transformed data, cmplx trace */
	data = ealloc1float(nt);
	hdata = ealloc1float(nt);
	ct = ealloc1complex(nt);


	/* Loop over traces */
	do {
		register int i;

		/* Get data from trace */
		for (i = 0; i < nt; ++i)  data[i] = tr.data[i];

		
		/* construct quadrature trace with hilbert transform */
		hilbert(nt, data, hdata);

		/* build the complex trace */
		for (i = 0; i < nt; ++i)  ct[i] = cmplx(data[i],hdata[i]);

		/* Form absolute value, phase, or frequency */
		switch(imode) {
		case AMP:
			for (i = 0; i < nt; ++i) {
				float re = ct[i].r;
				float im = ct[i].i;
				tr.data[i] = sqrt(re*re + im*im);
			}
			
			/* set trace id */
			tr.trid = ENVELOPE;
		break;
		case ARG:
		{
			float *phase = ealloc1float(nt);

			for (i = 0; i < nt; ++i) {
				float re = ct[i].r;
				float im = ct[i].i;
				if (re*re+im*im)  phase[i] = atan2(im, re);
				else              phase[i] = 0.0;
			}

			/* phase unwrapping */
			/* default unwrap=0 for this mode */
			if (unwrap!=0) unwrap_phase(nt, unwrap, phase);
			
			/* write phase values to tr.data */
			for (i = 0; i < nt; ++i) tr.data[i] = phase[i];
			
			/* set trace id */
			tr.trid = INSTPHASE;
		}
		break;
		case FREQ:
		{
			float *phase = ealloc1float(nt);
			float	fnyq = 0.5 / dt;

			for (i = 0; i < nt; ++i) {
				float re = ct[i].r;
				float im = ct[i].i;
				if (re*re+im*im) {
					phase[i] = atan2(im, re);
				} else {
					phase[i] = 0.0;
				}
				
			}

			/* unwrap the phase */
			if (unwrap!=0) unwrap_phase(nt, unwrap, phase);

			/* compute freq(t)=dphase/dt */
			differentiate(nt, 2.0*PI*dt, phase);
			
			/* correct values greater nyquist frequency */
			for (i=0 ; i < nt; ++i)	{
				if (phase[i] > fnyq)
					phase[i] = 2 * fnyq - phase[i];
			}
                                        
			/* write freq(t) values to tr.data */
			for (i=0 ; i < nt; ++i) tr.data[i] = phase[i];

			/* set trace id */
			tr.trid = INSTFREQ;
		}
		break;
		case FREQW:
		{
			float	fnyq = 0.5 / dt;
			float *freqw = ealloc1float(nt);
			float *phase = ealloc1float(nt);
			float *envelop = ealloc1float(nt);
			float *envelop2 = ealloc1float(nt);
			for (i = 0; i < nt; ++i) {
				float re = ct[i].r;
				float im = ct[i].i;
				if (re*re+im*im) {
					phase[i] = atan2(im, re);
					} else {
						phase[i] = 0.0;
						}
				envelop[i] = sqrt(re*re + im*im);
			}

			/* unwrap the phase */
			if (unwrap!=0) unwrap_phase(nt, unwrap, phase);

			/* compute freq(t)=dphase/dt */
			differentiate(nt, 2.0*PI*dt, phase);
			
			/* correct values greater nyquist frequency */
			for (i=0 ; i < nt; ++i)	{
				if (phase[i] > fnyq)
					phase[i] = 2 * fnyq - phase[i];
			envelop2[i]=envelop[i]*phase[i];
			}
			twindow(nt, wint, envelop);
			twindow(nt, wint, envelop2);
			/* correct values greater nyquist frequency */
			for (i=0 ; i < nt; ++i) {
			freqw[i] = (envelop[i] == 0.0) ? 0.0 :envelop2[i]/envelop[i];
			}
			/* write freq(t) values to tr.data */
			for (i=0 ; i < nt; ++i) tr.data[i] = freqw[i];
			
			/* set trace id */
			tr.trid = INSTFREQ;
		}
		break;
		case THIN:
		{
			float	fnyq = 0.5 / dt;
			float *phase = ealloc1float(nt);
			float *freqw = ealloc1float(nt);
			float *phase2 = ealloc1float(nt);



			for (i = 0; i < nt; ++i) {
				float re = ct[i].r;
				float im = ct[i].i;

				if (re*re+im*im) {
					phase[i] = atan2(im, re);
				} else {
					phase[i] = 0.0;
				}
			}

			/* unwrap the phase */
			if (unwrap!=0) unwrap_phase(nt, unwrap, phase);

			/* compute freq(t)=dphase/dt */
			differentiate(nt, 2.0*PI*dt, phase);