suharlan.c

su 的源代码库

analytic_envelopes	2-D array[nx][nt] of computed analytic envelops 	
******************************************************************************/
{
	int ix, it;			/* loop counters */
	float amp1,amp2,amp3;		/* auxiliary variables */
	float *hilbert_trace;		/* hilbert transform for 1 trace */

	/* allocate working space */
	hilbert_trace=alloc1float(nt);

	for (ix=0; ix<nx; ix++) {

		/* compute hilbert transform of current trace */
		hilbert (nt, traces[ix], hilbert_trace);

		for (it=0; it<nt; it++) {
			amp1=hilbert_trace[it];
			amp2=traces[ix][it];

			/* compute squared analytic envelopes */
			amp3=amp1*amp1+amp2*amp2;

			/* take care of square root */
			if (sgn==1||sgn==2||sgn==3) amp3=sqrt(amp3);

			/* take care of sign */
			if (sgn==3||sgn==6) analytic_envelopes[ix][it] =-amp3;
			else if (sgn==1||sgn==4) {
				if(amp2<0.0) analytic_envelopes[ix][it] =-amp3;
				else analytic_envelopes[ix][it] = amp3;
			} else {
				analytic_envelopes[ix][it]=amp3;
			}
		}
	}

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


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

	Subroutine to scale output traces to the same amplitude level
			of the input traces

******************************************************************************/
void scale_traces (int scl, int nx, int nt, float **in_traces, 
	float **out_traces)
/******************************************************************************
Input parameters:
scl		=1 for trace scaling only
		=2 for time sice scaling only
		=3 for both 
nx		number of traces
nt		number of samples per trace
in_traces	reference set of traces to scale out_traces
out_traces	set of traces to be scaled according to amplitudes of
		in_traces

Output parameters:
out_traces	set of scaled traces
******************************************************************************/
{
	int ix,it;

	/* if requested, apply trace by trace scaling */
	if (scl==1||scl==3) {

		/* loop over traces */
		for (ix=0; ix<nx; ix++) { 
			scale_one_trace (nt,in_traces[ix],out_traces[ix]);
		}
	}  

	/* if requested, apply scaling by time sices */
	if (scl==2||scl==3) {
		float **temp_in;		/* transpose of in_traces */
		float **temp_out;		/* transpose of out_traces */

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

		/* transpose input and output arrays */
		matrix_transpose (nx, nt, in_traces, temp_in);
		matrix_transpose (nx, nt, out_traces, temp_out);

		/* loop over time slices */
		for (it=0; it<nt; it++) {
			scale_one_trace (nx, temp_in[it], temp_out[it]);
		}

		/* transpose output back */
		matrix_transpose (nt, nx, temp_out, out_traces);

		/* free allocated space */
		free2float(temp_in);
		free2float(temp_out);
	}
}

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

	Subroutine to scale one trace to the same amplitude level of another

******************************************************************************/
void scale_one_trace (int ns, float *in_trace, float *out_trace)
/******************************************************************************
Input parameters:
ns		number of samples in trace to scale 
in_trace	reference trace to scale out_trace
out_trace 	trace to be scaled according to amplitudes of in_trace

Output parameters:
out_trace	scaled trace
******************************************************************************/	
{
	int is;
	float num,den;
	float scale;

	/* compute dot product of reference traces */
	num = dot_product (ns, in_trace, in_trace);

	/* compute dot product of output traces */	
	den = dot_product (ns, out_trace, out_trace);

	/* compute scale factor */
	if (den==0.0) {
		return;
	} else {
		scale = sqrt(num/den);	

		/* apply scale factor to each sample */
		for (is=0; is<ns; is++) {
			out_trace[is] *=scale;
		}
	}	
}

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

		Subroutine to transpose a matrix

******************************************************************************/
void matrix_transpose (int n1, int n2, float **matrix, float **tr_matrix)
/******************************************************************************
Input:
n1		number of number of columns in input matrix
n2		number of rows in input matrix
matrix		2-D array[n1][n2] to be transposed

Output:
tr_matrix	2-D array[n2][n1] of transposed matrix
******************************************************************************/
{
	int i1,i2;			/* loop counters */

	for (i1=0;i1<n1;i1++)
		for (i2=0;i2<n2;i2++)
			tr_matrix[i2][i1]=matrix[i1][i2];
}

/******************************************************************************
	
		compute the dot product of two time series

******************************************************************************/
float dot_product (int ns, float *vector1, float *vector2)
/******************************************************************************
Input:
ns		number of samples in vectors
vector1		array[ns] of first vector
vector2		array[ns] of second vector

Output:
		dot product of vector1 and vector2
******************************************************************************/
{
	int is;
	float sum=0.0;

	/* compute the dot product */
	for (is=0; is<ns; is++) sum +=vector1[is]*vector2[is]; 

	/* output result */
	return(sum);
}

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

	Subroutine to compute the maximum, minimum and sum of the
			samples in a 2-D array

******************************************************************************/
void compute_max_min_sum (int nx, int nt, float *min, float *max, float *sum,
	float **data)
/******************************************************************************
Input:
nx		number of horizontal samples (traces)
nt		number of vertical samples (samples per trace)
min		pointer to output minimum value
max		pointer to output maximum value
sum		pointer to output sum
data		2-D array[nx][nt] of traces

Output:
		max, min and sum values
******************************************************************************/
#define MAXVAL 9999999
{
	int ix,it;		/* loop counters */
	float d;		/* auxiliary variable for sample amplitude */

	/* initialize output variables */
	*min=MAXVAL;
	*max=-MAXVAL;
	*sum=0.0;
	
	for (ix=0; ix<nx; ix++)
		for (it=0; it<nt; it++) {

			/* get sample amplitude */
			d=data[ix][it];
	
			/* compute max, min and sum */
			if (d<*min) *min=d;
			if (d>*max) *max=d;
			*sum +=d;
		}
}

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

		Output data for a one dimensional plot

******************************************************************************/
void plot_one_d (int npoints, float *xamps, float *data, char *plotname)
/******************************************************************************
Input Parameters:
npoints		number of opints to plot
xamps		1-D array [npoints] of abscissa values
data		1-D array [npoints] of ordinate values
plotname	character array of plot name

Output:
		Binary file called plotname
******************************************************************************/
{
	FILE *out;		/* file pointer to output data */

	out=fopen(plotname,"w");
	fwrite (xamps, sizeof(float), npoints, out);
	fclose(out);
	out=fopen(plotname,"a");
	fwrite (data, sizeof(float), npoints, out);
	fclose(out);
}

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

		Subroutine to output a two-dimensional plot

******************************************************************************/
void plot_two_d (int npoints, float *data, char *plotname)
/******************************************************************************
Input Parameters:
npoints		number of points in output plot
data		1-D array [npoints] of data to plot
plotname	character array of plotname

Output:
		file called plotname ready to plot 
******************************************************************************/
{
	FILE *out;

	out=fopen(plotname,"w"); 
	fwrite (data, sizeof(float), npoints, out);
	fclose(out);
}

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

	Subroutine to compute convolution, correlation or bounded
		convolution of two input 1-D arrays

******************************************************************************/
void conv1 (int nx, int fx, float *x, int ny, int fy, float *y, int nz,
	int fz, float *z, int flag, float perc)
/******************************************************************************
Input:
nx		number of samples in first input array
fx		index of first sample of first input array
x[nx]		first input array
ny		number of samples in second input array
fy		index of first sample of second input array
y[ny]		second input array
nz		number of samples in output array
fz		index of first sample in output array
flag		=0 for convolution, =1 for correlation and =-1 for bounded
		convolution
perc		Desired percentage for bounded convolution.
		=1 for regular convolution and correlation

Output:
z[nz]		convolution or correlation of arrays x and y

Note: z cannot be equal to x or y
*******************************************************************************
This subroutine is a direct translation to C of a Fortran version written
by Dr. Bill Harlan, 1984.
******************************************************************************/
{
	int ix,iy,iz;
	int jx,jy,jz;
	float rl,rr,r,s;

	/* initialize output array */
	for (iz=0; iz<nz; iz++) z[iz]=0.0;

	/* compute convolution or correlation */	
	for (iz=1; iz<=nz; iz++) {
		jz=iz+fz-1;

		for (ix=1; ix<=nx; ix++) {
			
			jx=ix+fx-1;
			jy=jz-jx;

			/* correlation, change the sign */
			if (flag==1) jy=jz+jx;
		
			iy=jy-fy+1;

			/* exit inner loop if out of bounds */
			if (iy<1||iy>ny) continue;
	
			s=1;	/* for normal convolution or correlation */
	
			/* if bounded convolution is desired */
			if (flag==-1) {
				r=jz;
				r=ABS(r)*perc;
				rl=jy-0.5;
				rr=jy+0.5;
				if (rl>r) rl=r;
				if (rl<-r) rl=-r;
				if (rr>r) rr=r;
				if (rr<-r) rr=-r;
				s=rr-rl;
			}

			/* update sum */
			z[iz-1] += s*x[ix-1]*y[iy-1];
		}
	}
}
	
/******************************************************************************

	Subroutine to compute ps(x) as a deconvolution of pd(x)
		and pn(x) with constraints of positivity and
				unit area

******************************************************************************/
void deconvolve_histograms (int nsamples, int mean_index, int niter,
	float *pd, float *pn, float *ps)
/******************************************************************************
Input:
nsamples	Number of samples in histograms