suattributes.c

su 的源代码库

			/* correct values greater nyquist frequency */
			for (i=0 ; i < nt; ++i)	{
				if (phase[i] > fnyq)
					phase[i] = 2 * fnyq - phase[i];
					phase2[i]= 2 * fnyq - phase[i];
			}
			/* Do windowing for Average Ins . Freq over wint*/
			twindow(nt, wint, phase2);

			for (i=0 ; i < nt; ++i)	{
				freqw[i] = phase[i] - phase2[i];
			/*	if (abs(freqw[i]) > fnyq)
				freqw[i] = 2 * fnyq - freqw[i];
			*/
			/* write Thin-Bed(t) values to tr.data */
				tr.data[i] = freqw[i];
				}
			/* set trace id */
			tr.trid = INSTFREQ;
		}
		break;
		case BANDWIDTH:
		{
			float *envelop = ealloc1float(nt);
			float *envelop2 = ealloc1float(nt);

		/* Bandwidth (Barnes 1992)

		          |d(envelope)/dt|
		band =abs |--------------|
		          |2 PI envelope |
	 	*/

			for (i = 0; i < nt; ++i) {
				float er = ct[i].r;
				float em = ct[i].i;
				envelop[i] = sqrt(er*er + em*em);
				envelop2[i]=sqrt(er*er + em*em);

			}
				differentiate(nt, dt, envelop);

				for (i = 0; i < ntout; ++i) {
				   if (2.0*PI*envelop2[i]!=0.0) {
					tr.data[i] = ABS(envelop[i]/(2.0*PI*envelop2[i]));
				   } else {
				        tr.data[i]=0.0;
				   }
				}
				tr.trid = ENVELOPE;
		}
		break;
		case NORMAMP:
		{
			float phase;
			float *na = ealloc1float(nt);
			for (i = 0; i < nt; ++i) {
				float re = ct[i].r;
				float im = ct[i].i;
				if (re*re+im*im)  phase = atan2(im, re);
				else              phase = 0.0;
				na[i] = cos(phase);
			}
			for (i=0 ; i < nt; ++i) tr.data[i] = na[i];
			
			/* set trace id */
			tr.trid = INSTPHASE;
			}
		break;
		case FENV:
		{
			float *amp = ealloc1float(nt);
			for (i = 0; i < nt; ++i) {
				float re = ct[i].r;
				float im = ct[i].i;
				amp[i] = sqrt(re*re + im*im);
			}
		/*conv(nt, 0, envelop, nt, 0, time, ntout, 0, ouput);*/

		differentiate(nt, 2.0*PI*dt, amp);
		for (i=0 ; i < nt; ++i) tr.data[i] = amp[i];
			/* set trace id */
			tr.trid = ENVELOPE;
		}
		break;
		case SENV:
		{
			float *amp = ealloc1float(nt);
			for (i = 0; i < nt; ++i) {
				float re = ct[i].r;
				float im = ct[i].i;
				amp[i] = sqrt(re*re + im*im);
			}

		differentiate(nt, 2.0*PI*dt, amp);
		differentiate(nt, 2.0*PI*dt, amp);
		for (i=0 ; i < nt; ++i) tr.data[i] = amp[i];
			/* set trace id */
			tr.trid = ENVELOPE;
		}
		break;

		case Q:
		{
			float *envelop = ealloc1float(nt);
			float *envelop2 = ealloc1float(nt);
			float *phase = ealloc1float(nt);
			float	fnyq = 0.5 / dt;

		/* Banswith (Barnes 1992)

		        -PI Freq(t) d(envelope)/dt
		band =  --------------------------
		                 envelope(t)
	 	*/

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

			}
			/* get envelope diff */
			differentiate(nt, dt, envelop);
			/* unwrap the phase */
			if (unwrap!=0) unwrap_phase(nt, unwrap, phase);
			/* compute freq(t)=dphase/dt */
			differentiate(nt, 2.0*PI*dt, phase);

			for (i=0 ; i < nt; ++i)	{
				if (phase[i] > fnyq)
					phase[i] = 2 * fnyq - phase[i];
			}

			for (i = 0; i < ntout; ++i) {
				if (envelop[i]!=0.0)
				tr.data[i] = -1*PI*phase[i]*envelop2[i]/envelop[i];
				else
				tr.data[i]=0.0;
				}
				tr.trid = INSTFREQ;
		}
		break;
		default:
			err("%s: mysterious mode=\"%s\"", __LINE__, mode);
		}


		tr.d1 = dt;   /* for graphics */
		puttr(&tr);

	} while (gettr(&tr));


	return(CWP_Exit());
}


void differentiate(int n, float h, float *f)
/************************************************************************
differentiate - compute the 1st derivative of a function f[]
************************************************************************
Input:
n		number of samples
h		sample rate
f		array[n] of input values

Output:
f		array[n], the derivative of f
************************************************************************
Notes:
This is a simple 2 point centered-difference differentiator.
The derivatives at the endpoints are computed via 2 point leading and
lagging differences. 
************************************************************************
Author: John Stockwell, CWP, 1994
************************************************************************/
{
	int i;	
	float *temp;
	float h2=2*h;

	/* allocate space in temporary vector */
	temp = ealloc1float(n);

	/* do first as a leading difference */
	temp[0] = (f[1] - f[0])/h;

	/* do the middle values as a centered difference */
	for (i=1; i<n-1; ++i) temp[i] = (f[i+1] - f[i-1])/h2;

	/* do last value as a lagging difference */
	temp[n-1] = (f[n-1] - f[n-2])/h;

	for (i=0 ; i < n ; ++i) f[i] = temp[i];

	free1float(temp);
}

void unwrap_phase(int n, float w, float *phase)
/************************************************************************
unwrap_phase - unwrap the phase
*************************************************************************
Input:
n		number of samples
w		unwrapping flag; returns an error if w=0
phase		array[n] of input phase values

Output:
phase		array[n] of output phase values
*************************************************************************
Notes:
The phase is assumed to be continuously increasing. The strategy is
to look at the change in phase (dphase) with each time step. If it is larger
than PI/w, then use the previous value of dphase. No attempt is
made at smoothing the dphase curve.
*************************************************************************
Author: John Stockwell, CWP, 1994
************************************************************************/
{
	int i;
	float pibyw=0.0;
	float *dphase;
	float *temp;

	/* prevent division by zero in PI/w */
	if (w==0)  err("wrapping parameter is zero");
	else       pibyw = PI/w;

	/* allocate space */
	dphase = ealloc1float(n);
	temp = ealloc1float(n);

	/* initialize */
	temp[0]=phase[0];
	dphase[0]=0.0;

	/* compute unwrapped phase at each time step */
	for (i = 1; i < n; ++i) {

		/* compute jump in phase */
		dphase[i] = ABS(phase[i] - phase[i-1]);

		/* if dphase >= PI/w, use previous dphase value */
		if (ABS(dphase[i] - dphase[i-1]) >= pibyw )
			dphase[i] = dphase[i-1];

		/* sum up values in temporary vector */
		temp[i] = temp[i-1] + dphase[i];
	}

	/* assign values of temporary vector to phase[i] */
	for (i=0; i<n; ++i) phase[i] = temp[i];

	/* free space */
	free1float(temp);
	free1float(dphase);
}

void twindow(int nt, int wtime, float *data)
/************************************************************
twindow - simple time gating
*************************************************************
Input:
nt	number of time samples
wtime	= n*dt   where n are integer ex=1,2,3,4,5,...
          wtime=3 as default
 used for Frequency Weighted and Thin-bed attributes
*************************************************************
Author:	UGM (Geophysics Students): Agung Wiyono, 2005
************************************************************/
{
	float val;
	float *temp;
	int i;
	float sum;
	int nwin;
	
	nwin=2*wtime+1;
	temp = ealloc1float(nt);
	sum=0.0;
	for (i = 0; i< wtime+1; ++i) {
		val = data[i];
		sum +=val;
	}
	/* weighted */
	temp[0] = sum/nwin;
	
	/* dt<wtime */
	for (i = 1; i < wtime; ++i) {
		val = data[i+wtime];
		sum+=val;
		++nwin;
		temp[i] = sum/nwin;
		}
	/*wtime<dt<dt-wtime */
	for (i = wtime ; i < nt-wtime; ++i) {
		val = data[i+wtime];
		sum += val;
		val = data[i-wtime];
		sum -=val;
		temp[i] = sum/nwin;
	}

	/*dt-wtime<dt*/
	for (i = nt - wtime; i < nt; ++i) {
		val = data[i-wtime];
		sum -= val;
		--nwin;
		temp[i] = sum/nwin;
	}
	
	
	for (i=0; i<nt; ++i) data[i] = temp[i];

	/* Memori free */
	free1float(temp);
}