supef.c

su 的源代码库

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

/* SUPEF: $Revision: 1.41 $ ; $Date: 2006/11/07 22:58:42 $		*/

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

/*********************** self documentation ******************************/
char *sdoc[] = {
" 									",
" SUPEF - Wiener predictive error filtering				",
" 									",
" supef <stdin >stdout  [optional parameters]				",
" 									",
" Required parameters:							",
" dt is mandatory if not set in header			 		",
" 									",
" Optional parameters:							",
" cdp= 			CDPs for which minlag, maxlag, pnoise, mincorr, ",
"			maxcorr are set	(see Notes)			",
" minlag=dt		first lag of prediction filter (sec)		",
" maxlag=last		lag default is (tmax-tmin)/20			",
" pnoise=0.001		relative additive noise level			",
" mincorr=tmin		start of autocorrelation window (sec)		",
" maxcorr=tmax		end of autocorrelation window (sec)		",
" showwiener=0		=1 to show Wiener filter on each trace		",
" mix=1,...	 	array of weights (floats) for moving		",
"				average of the autocorrelations		",
" outpar=/dev/null	output parameter file, contains the Wiener filter",
" 			if showwiener=1 is set				",
" method=linear	 for linear interpolation of cdp values			",
"		       =mono for monotonic cubic interpolation of cdps	",
"		       =akima for Akima's cubic interpolation of cdps	",
"		       =spline for cubic spline interpolation of cdps	",
" 									",
" Trace header fields accessed: ns, dt					",
" Trace header fields modified: none					",
" 									",
" Notes:								",
" 									",
" 1) To apply spiking decon (Wiener filtering with no gap):		",
"									",
" Run the following command						",
"									",
"    suacor < data.su | suximage perc=95				",
"									",
" You will see horizontal strip running across the center of your plot.	",
" This is the autocorrelation wavelet for each trace. The idea of spiking",
" decon is to apply a Wiener filter with no gap to the data to collapse	",
" the waveform into a spike. The idea is to pick the width of the	",
" autocorrelation waveform _from beginning to end_ (not trough to trough)",
" and use this time for MAXLAG_SPIKING:					",
"									",
"  supef < data.su maxlag=MAXLAG_SPIKING  > dataspiked.su		",
"									",
" 2) Prediction Error Filter (i.e. gapped Wiener filtering)		",
" The purpose of gapped decon is to suppress repetitions in the data	",
" such as those caused by water bottom multiples.			",
"									",
" To look for the period of the repetitions				",
" 									",
"    suacor ntout=1000 < dataspiked.su | suximage perc=95		",
" 									",
" The value of ntout must be larger than the default 100. The idea is	",
" to look for repetitions in the autocorrelation. These repetitions will",
" appear as a family of parallel stripes above and below the main	",
" autocorrelation waveform. Set MAXLAG_PEF to the period of the repetitions",
" Set MINLAG_PEF to be slightly larger than the value of MAXLAG_SPIKING	",
" that you used to spike the data.					",
" 									",
"  supef < dataspiked.su minlag=MINLAG_PEF maxlag=MAXLAG_PEF > datapef.su",
" 									",
" Some experimentation may be required to get a satisfactory result.	",
" 									",
" 3) It may be effective to sort your data into cdp gathers with susort,",
" and perform sunmo correction to the water speed with sunmo, prior to 	",
" attempts to suppress water bottom multiples. After applying supef, the",
" user should apply inverse nmo to undo the nmo to water speed prior to	",
" further processing.							",
"  									",
" For a filter expressed as a function of cdp, specify the array	",
"     cdp=cdp1,cdp2,...							",
" and for each cdp specified, specify the minlag and maxlag arrays as	",
"      minlag=min1,min2,...     maxlag=max1,max2,...   			",
"  									",
" It is required that the number of minlag and maxlag values be equal to",
" the number of cdp's specified.  If the number of			",
" values in these arrays does not equal the number of cdp's, only the first",
" value will be used.							",
"									",
" Caveat:								",
" The showweiner=1 option writes out the wiener filter to outpar, and   ",
" the prediction error filter to stdout, which is			", 
"     1,0,0,...,-wiener[0],...,-wiener[imaxlag-1] 			",
" where the sample value of -wiener[0], is  iminlag in the pe-filter.	",
" The pe-filter is output as a SU format datafile, one pe-filter for each",
" trace input.								",
"									",
 NULL};

/* Credits:
 *	CWP: Shuki Ronen, Jack K. Cohen, Ken Larner
 *      CWP: John Stockwell, added mixing feature (April 1998)
 *      CSM: Tanya Slota (September 2005) added cdp feature
 *
 *      Technical Reference:
 *	A. Ziolkowski, "Deconvolution", for value of maxlag default:
 *		page 91: imaxlag < nt/10.  I took nt/20.
 *
 * Notes:
 *	The prediction error filter is 1,0,0...,0,-wiener[0], ...,
 *	so no point in explicitly forming it.
 *
 *	If imaxlag < 2*iminlag - 1, then we don't need to compute the
 *	autocorrelation for lags:
 *		imaxlag-iminlag+1, ..., iminlag-1
 *	It doesn't seem worth the duplicated code to implement this.
 *
 * Trace header fields accessed: ns
 */
/**************** end self doc *******************************************/


/* External definitions */
#define PNOISE	0.001			/* default pnoise value		*/
#define VAL0    1.0			/* default weighting value	*/
#define OUTPAR_DEFAULT "/dev/null"	/* default output filename	*/




segy intrace, outtrace;

int
main(int argc, char **argv)
{
	int nt;			/* number of points on trace		*/
	int i,ilag;		/* counters				*/
	int ncdp;		/* number of cdp's specified		*/
	int nminlag;		/* number of minlags specified		*/
	int nmaxlag;		/* number of maxlag specified		*/
	int icdp=0;		/* counter of cdp's			*/
	int jcdp=0;		/* counter of cdp's			*/

	float *cdp=NULL;	/* array[ncdp] of cdps specified 	*/
	float dt;		/* time sample interval (sec)		*/
	float *wiener=NULL;	/* Wiener error filter coefficients	*/
	float *spiker=NULL;	/* spiking decon filter			*/

	float pnoise;		/* pef additive noise level		*/

	float *minlag=NULL;	/* start of error filter (sec)		*/
	int *iminlag=NULL;	/* ... in samples			*/
	float *maxlag=NULL;	/* end of error filter (sec)		*/
	int *imaxlag=NULL;	/* ... in samples			*/
	int nlag;		/* length of error filter in samples	*/
	int ncorr;		/* length of corr window in samples	*/
	int lcorr;		/* length of autocorr in samples	*/

	long oldcdp;		/* cdp of previous trace */

	float *crosscorr=NULL;	/* right hand side of Wiener eqs	*/
	float *autocorr=NULL;	/* vector of autocorrelations		*/

	float mincorr;		/* start time of correlation window	*/
	int imincorr;		/* .. in samples			*/
	float maxcorr;		/* end time of correlation window	*/
	int imaxcorr;		/* .. in samples			*/
	int showwiener;		/* flag to display pred. error filter	*/
	
	size_t lagbytes;	/* bytes in wiener and spiker filters	*/
	size_t maxlagbytes;	/* bytes in autocorrelation		*/
	
	int imix;		/* mixing counter			*/
	int nmix;		/* number of traces to average over	*/
	size_t mixbytes;	/* number of bytes = maxlagbytes*nmix	*/ 
	float *mix=NULL;	/* array of averaging weights		*/
	float **mixacorr=NULL;	/* mixing array				*/
	float *temp=NULL;	/* temporary array			*/

	/* Interpolation */
	float (*zind)[4]=NULL;	/* array of interpolation coefficients  */
	char  *method="linear";	/* interpolation method			*/
	float *cdpinterp=NULL;	/* interpolated cdps			*/
	float *ominlag=NULL;	/* interpolated minlags			*/ 
	float *omaxlag=NULL;	/* interpolated maxlags			*/

        char *outpar=NULL;	/* name of file holding output parfile	*/
        FILE *outparfp=NULL;	/* ... its file pointer			*/


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



	/* Get info from first trace */ 
	if (!gettr(&intrace)) err("can't get first trace");
	nt = intrace.ns;
	dt = ((double) intrace.dt)/1000000.0;
	if (!dt) MUSTGETPARFLOAT ("dt", &dt); 

	/* Get parameters and do set up */
	if (!getparstring("outpar", &outpar))	outpar = OUTPAR_DEFAULT;
	outparfp = efopen(outpar, "w");


	/* Count parameters and allocate space */
	ncdp = countparval("cdp");
	if (ncdp>0) {
		if ((nminlag = countparval("minlag"))!=0) {
			if (!(nminlag==ncdp)) 
			   err("number of cdp values must equal number of minlag values");
		} else {
			err("cdp set, but minlag not getparred!");
		}

		if ((nmaxlag = countparval("maxlag"))!=0) {
			if (!(nmaxlag==ncdp)) 
			   err("number of cdp values must equal number of maxlag values");
		} else {
			err("cdp set, but maxlag not getparred!");
		}

		/* allocate space */
		cdp = ealloc1float(ncdp);
		minlag = ealloc1float(ncdp);
		maxlag = ealloc1float(ncdp);
		
		/* get parameters */
		getparfloat("minlag",minlag);
		getparfloat("maxlag",maxlag);
		getparfloat("cdp",cdp);

		/* check cdp values */
		for (icdp = 0; icdp < ncdp-1; ++icdp)  {
			if(cdp[icdp] > cdp[icdp+1])
				err("cdp values must increase monotonically!");
		}


		iminlag = ealloc1int(ncdp);
		imaxlag = ealloc1int(ncdp);

		for (icdp = 0; icdp < ncdp; ++icdp) {
			iminlag[icdp] = NINT(minlag[icdp]/dt);	
			imaxlag[icdp] = NINT(maxlag[icdp]/dt);
			if (iminlag[icdp] < 1) 
				err("minlag[%d]=%g too small",
					icdp, minlag[icdp]);
			if (imaxlag[icdp] >= nt)
				err("maxlag[%d]=%g too large",
					icdp, maxlag[icdp]);
			if (iminlag >= imaxlag)
				err("minlag[%d]=%g, maxlag[%d]=%g",
					icdp, minlag[icdp], icdp, maxlag[icdp]);
		}

	} else {
		ncdp = 1;
		cdp = ealloc1float(ncdp);
		minlag = ealloc1float(ncdp);
		maxlag = ealloc1float(ncdp);
		iminlag = ealloc1int(ncdp);
		imaxlag = ealloc1int(ncdp);

		/* zero out array values */
		memset( (void *) minlag , 0, ncdp * FSIZE);
		memset( (void *) maxlag , 0, ncdp * FSIZE);
		memset( (void *) imaxlag , 0, ncdp * ISIZE);
		memset( (void *) iminlag , 0, ncdp * ISIZE);


		cdp[0] = intrace.cdp;

		if (getparfloat("minlag", minlag)) {
			iminlag[0] = NINT(minlag[0]/dt);
		} else {
			iminlag[0] = 1;
		}

		if (getparfloat("maxlag", maxlag)) {
			imaxlag[0] = NINT(maxlag[0]/dt);
		} else {
			imaxlag[0] = NINT(0.05 * nt);
		}

	}