suharlan.c

su 的源代码库

		traces, anenv_traces, reliability);

	/* check constraint of unit area for histograms */
	if (verbose==1) {
		spd=spn=sps=spspn=0.0;
		for (ix=0;ix<nintlh;ix++) {
			spd +=pd[ix];
			spn +=pn[ix];
			sps +=ps[ix];
			spspn +=pspn[ix];
		}
	
		if (verbose) {
			warn("Area under the curve for pd=%g\n",spd);
			warn("Area under the curve for pn=%g\n",spn);
			warn("Area under the curve for ps=%g\n",sps);
			warn("Area under the curve for ps*pn=%g\n",spspn);
		}
	}	

	/* if requested, output plots of histograms, Esd and reliabilities */
	count++;
	if (plot==1||plot==3) {
		if (count==1) plotname="pd1";
		else if (count==2) plotname="pd2";
		plot_one_d (nintlh, xamps, pd, plotname);
		if (count==1) plotname="pn1";
		else if (count==2) plotname="pn2";
		plot_one_d (nintlh, xamps, pn, plotname);
		if (count==1) plotname="ps1";
		else if (count==2) plotname="ps2";
		plot_one_d (nintlh, xamps, ps, plotname);
		if (count==1) plotname="pspn1";
		else if (count==2) plotname="pspn2";
		plot_one_d (nintlh, xamps, pspn, plotname);
		if (count==1) plotname="Esd1";
		else if (count==2) plotname="Esd2";
		plot_one_d (nintlh, xamps, Esd, plotname);
		if (count==1) plotname="rel1";
		else if (count==2) plotname="rel2";
		plot_one_d (nintlh, xamps, reliability, plotname);
	}

	/* free allocated space */
	free2float(anenv_traces);
	free1float(pd);
	free1float(ps);
	free1float(pn);
	free1float(pspn);
	free1float(amps);
	free1float(xamps);
	free1float(Esd);
	free1float(reliability);
}

/******************************************************************************

	Subroutine to compute an array of expected values E(s|d)

******************************************************************************/
void expected_value (int nintlh, int nx, int nt, int rindex, float *pspn,
	float *pn, float *ps, float *xamps, float **traces, float *Esd)
/******************************************************************************
Input:
anenv		=0 do not use analytic envelopes
		!=0 use analytic envelopes
clipf		=1 use the maximum amplitude as the clip value for division
		=2 use the sample amplitude as the clip value for division
nintlh	  number of intervals (bins) per histogram
nx		number of horizontal samples (traces) (ignored)
nt		number of vertical (time) samples (ignored)
rindex		histogram index of mean data sample
pspn		1-D array[nintlh] of ps convolved with pn
ps		1-D array[nintlh] of signal probability density function
pn		1-D array[nintlh] of noise probability density function
pd		1-D array[nintlh] of data probability density function
traces    	2-D array[nx][nt] of input traces
anenv_traces    2-D array[nx][nt] of analytic envelopes traces

Output Parameters:
Esd		1-D array[nintlh] of E(s|d) values for reference amplitudes
*******************************************************************************
expected_value uses Harlan's description of the process. It only computes
the Esd indicator for the reference amplitudes, (bin centers) and applies
linear interpolation for all other amplitudes.
******************************************************************************/
{
	int ix;			/* loop counters */
	float *xps;		/* array of x*ps(x) */
	float max;		/* maxumum amplitude in input data */

	ix = 0*nt*nx; traces += ix;
	max=xamps[nintlh-1]+(xamps[1]-xamps[0]);

	/* allocate working space */
	xps = alloc1float(nintlh);

	/* compute the product x*ps(x) */
	for (ix=0; ix<nintlh; ix++) xps[ix] = xamps[ix]*ps[ix];

	/* compute the convolution of xps(x) with pn(x) */
	conv (nintlh, -rindex, xps, nintlh, -rindex, pn, nintlh, -rindex, Esd);

	/* compute Esd for reference amplitudes (histogram bin centers) */
	for (ix=0; ix<nintlh; ix++) {
	
		/* do the division */
		Esd[ix]=divide(max,Esd[ix],pspn[ix]);
	}

	/* free allocated space */
	free1float(xps);
}

/******************************************************************************

	Compute reliability via bounded convolution

******************************************************************************/
void compute_reliability (int norm, int nx, int nt, int nintlh, int rindex,
	float c, float *ps, float *pn, float *pspn, float **traces,
	float *reliability)
/******************************************************************************
Input:
anenv		=0 not to use analytic envelopes
		!=0 use analytic envelopes
nx		number of hotizontal samples (traces) (ignored)
nt		number of vertical samples (samples per trace)
nintlh		number of bins per histogram
rindex		index of reference amplitude for histograms
c		maximum allowed error in percent
ps		1-D array[nintlh] of signal probability density histogram
pn		1-D array[nintlh] of noise probability density histogram
pspn		1-D array[nintlh] of ps convolved with pn
traces		2-D array[nx][nt] of input traces (ignored)
anenv_races	2-D array[nx][nt] of input analytic envelopes traces


Output:
reliability	1-D array[nintlh] of reliability values for bin centers
*******************************************************************************
Note:
This subroutine computes the reliability as a bounded convolution without 
explicitly computing the expected value of the signal given the data (Esd
indicator), assuming that for high amplitude samples the signal and the data
amplitudes are about the same, so in the limits of integration for the 
reliability indicator the s (signal) is implicitly replaced by d (data).
The reliability of a given sample is assumed to be that of the histogram bin 
to which it belongs)
******************************************************************************/
{
	int i;				/* loop counters */
	float max=-999999999;
	float *numerator;		/* numnerator of reliability */

	i = nt*nx*0;
	traces += i;
	/* allocate working space */
	numerator=alloc1float(nintlh);

	/* compute numerator of reliability as a bounded convolution */
		conv1 (nintlh, -rindex, ps, nintlh, -rindex, pn, nintlh, -rindex,			numerator, -1, c);

	/* compute reliability for bin centers */
	for (i=0; i<nintlh; i++) 
	  	reliability[i]=divide (1.0, numerator[i], pspn[i]);
	
	/* if requested normalize reliability */
	if (norm==1) {	
		for (i=0; i<nintlh; i++) {
			if (max<reliability[i]) max=reliability[i];
		}

		for (i=0; i<nintlh; i++) {
			reliability[i] /=max;
		}
	}

	/* free allocated space */
	free1float(numerator);
}

/******************************************************************************

	Subroutine to zero out noisy samples based on their Esd value

******************************************************************************/
void zero_noisy_samples (int anenv, int sts, int nx, int nt, float r1,
	float r2, float r, int nintlh, float *amps, float **traces, 
	float **anenv_traces, float *rel)
/******************************************************************************
Input parameters:
anenv		=1 use analytic envelopes
		=0 don't use them
sts		flag for smoothing
nx		number of horizontal samples (traces)
nt		number of vertical samples 
r		samples with reliability >r are kept 
nintlh		number of intervals in histograms	
Esd		1-D array[nintlh] of expected values for reference amplitudes
traces		2-D array[nx][nt] of input traces
anenv_traces	2-D array[nx][nt] of analytic envelope traces
rel	 	1-D array[nintlh] of reliability values for reference amplitudes

traces		2-D array[nx][nt] of extracted (focused) signal
******************************************************************************/
{
	int it,ix;		/* loop counters */
	int index=0;		/* sample index */
	float d;		/* sample amplitude */	
	float **oz;		/* array of ones and zeros */

	/* allocate working space */
	oz=alloc2float(nt,nx);

	/* create array of ones and zeros depending on reliability value */
	for (ix=0; ix<nx; ix++)
		for (it=0; it<nt; it++) {

			/* get sample amplitude */
			if (anenv==0) d=traces[ix][it];	 
			else d=anenv_traces[ix][it];

			/* find to which bin sample d belongs */
			xindex (nintlh, amps, d, &index);
			
			/* flag sample with zero or one */
			if (rel[index]<r) oz[ix][it]=0.0;
			else oz[ix][it]=1.0;
		}

	/* if requested, smooth array of zeros and ones horizontally and vert*/
	if (sts !=0) {
		dlsq_smoothing (nt, nx, 0, nt, 0, nx, r1, r2, 0, oz);
	}

	/* multiply data by smoothed array of ones and zeros */
	for (ix=0; ix<nx; ix++)
		for (it=0; it<nt; it++) 
			traces[ix][it] *=oz[ix][it];

	/* free allocated space */
	free2float(oz);
}

/******************************************************************************

	Subroutine to compute minimum, maximum, interval width and 
			amplitudes for histograms

******************************************************************************/
void compute_histogram_stuff (int verbose, int nx, int nt, int nintlh,
	float **traces, float **rand_traces, int *rindex, float *amps,
	float *xamps)
/******************************************************************************
Input:
nx		number of traces in input array
nt		number of samples per trace
nintlh		number of intervals (bins) in a histogram
rindex		pointer to reference index
traces		2-D array[nx][nt] of input traces
rand_traces	2-D array[nx][nt] of randomly reversed traces

Output:
amps		1-D array[nintlh] of left end bin amplitudes in histograms
xamps		1-D array[nintlh] of bin centers in histograms
		computed min, max and del_int values 
******************************************************************************/
{
	int ix;			/* loop counter */
	float del_int;		/* histogram bin width */
	float min1,max1;	/* minimum and maximum amplitudes of data */
	float min2,max2;	/* minimum and maximum amplitudes of noise */
	float fint;		/* auxiliary variable */
	float sum1=0.0;
	float sum2=0.0;
	float min;
	float max;

	/* compute min, max and sum of input array */
	compute_max_min_sum (nx, nt, &min1, &max1, &sum1, traces);  

	/* compute min, max and sum of randomly reversed array */
	compute_max_min_sum (nx, nt, &min2, &max2, &sum2, rand_traces);
	
	/* compute minimum and maximummbin amplitudes and bin interval */
	min = min1-0.5*(max2-min2);
	max = max1+0.5*(max2-min2);
	del_int = (max-min)/(nintlh-1);
	fint = min;
	min -= del_int/2.0;
	max += del_int/2.0;
		
	/* compute array of amplitudes for histograms */
	for (ix=0; ix<nintlh; ix++) {
		amps[ix]=min+del_int*ix;
		xamps[ix]=fint+del_int*ix;
	}

	/* find interval for the mean value of data histogram */
	xindex (nintlh, amps, sum1/(nt*nx), rindex);

	/* if requested, print processin information */
	if (verbose==1) {
		warn("for histograms: min=%g max=%g",min,max);
		warn("first bin=%g binwidth=%g\n",fint,del_int);
	}
}


/******************************************************************************

	Subroutine to compute one-dimensional histograms

******************************************************************************/
void make_histogram (int nintlh, int nt, int nx, float **traces, float *amps, 
	float *pdf)
/******************************************************************************
Input:
nintlh	  number of intervals per local histogram
nt		number of time samples
nx		number of horizontal samples
amps		1-D array[nintlh] of histogram amplitudes
traces		2-D array[nx][nt] of input traces

Output:
pdf		1-D array[nintlh] of computed  probabilty density function
******************************************************************************/
{
	int i,it,ix;   			/* loop counters */
	int index=0;			/* search index */
	int ns=nt*nx;			/* total number of samples */

	/* initialize pdf array */
	for (i=0; i<nintlh; i++) pdf[i]=0.0;

	/* compute histogram of data samples */
	for (ix=0; ix<nx; ix++) 
		for (it=0; it<nt; it++) {

			/* find interval to which sample belongs */
			xindex (nintlh, amps, traces[ix][it], &index);

			/* update interval frequencies */
			pdf[index] +=1.0;
		}

	/* normalize frequencies to get probability density distribution */
	for (i=0; i<nintlh; i++) pdf[i] /=ns;
}

/******************************************************************************

	Subroutine to compute the analytic envelopeis of an input 2-D
				array of traces

******************************************************************************/
void compute_analytic_envelopes (int sgn, int nx, int nt, float **traces, 
	float **analytic_envelopes)
/******************************************************************************
Input:
sgn			=1 compute signed analytic envelopes
			=2 compute positive analytic envelopes
			=3 compute negative analytic envelopes
			=4 compute signed square of analytic envelope
			=5 compute positive square analytic envelopes
			=6 compute negative analytic envelopes
nx			number of traces
nt			number of samples per trace
traces			2-D array[nx][nt] of input traces

Output: