sufxdecon.c

su 的源代码库

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

/* SUFXDECON: $Revision: 1.8 $ ; $Date: 2006/11/07 22:58:42 $	*/	

#include "su.h"
#include "segy.h"
#include "header.h"
#include <signal.h>

/*********************** self documentation **********************/
char *sdoc[] = {
" 									",
" SUFXDECON - random noise attenuation by FX-DECONvolution              ",
"									",
" sufxdecon <stdin >stdout [...]	                                ",
"									",
" Required Parameters:							",
"									",
" Optional Parameters:							",
" taper=.1	length of taper                                         ",
" fmin=6.       minimum frequency to process in Hz  (accord to twlen)   ",
" fmax=.6/(2*dt)  maximum frequency to process in Hz                    ",
" twlen=entire trace  time window length (minimum .3 for lower freqs)   ",
" ntrw=10       number of traces in window                              ",
" ntrf=4        number of traces for filter (smaller than ntrw)         ",
" verbose=0	=1 for diagnostic print					",
" tmpdir=	if non-empty, use the value as a directory path	prefix	",
"		for storing temporary files; else, if the CWP_TMPDIR	",
"		environment variable is set, use its value for the path;",
"		else use tmpfile()					",
"									",
" Notes: Each trace is transformed to the frequency domain.             ",
"        For each frequency, Wiener filtering, with unity prediction in ",
"        space, is used to predict the next sample.                     ",
"        At the end of the process, data is mapped back to t-x domain.  ", 
"									",
"                                                                       ",
NULL};

/* Credits:			
 *
 *	CWP: Carlos E. Theodoro (10/07/97)
 *
 * References:      							
 *		Canales(1984):'Random noise reduction' 54th. SEGM	
 *		Gulunay(1986):'FXDECON and complex Wiener Predicition   
 *                             filter' 56th. SEGM	                
 *		Galbraith(1991):'Random noise attenuation by F-X        
 *                             prediction: a tutorial' 61th. SEGM	
 *
 * Algorithm:
 *	- read data
 *	- loop over time windows
 *		- select data
 *		- FFT (t -> f)
 *		- loop over space windows
 *			- select data
 *			- loop over frequencies
 *				- autocorelation
 *				- matrix problem
 *				- construct filter
 *				- filter data
 *			- loop along space window
 *				- FFT (f -> t)
 *				- reconstruct data
 * 	- output data
 *
 * Trace header fields accessed: ns, dt, d1
 * Trace header fields modified: 
 */

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

#define LOOKFAC	2	/* Look ahead factor for npfaro	  */
#define PFA_MAX	720720	/* Largest allowed nfft	          */

/* Function Prototypes */
void cxcor (int lx, int ifx, complex *x,
            int ly, int ify, complex *y, 
            int lz, int ifz, complex *z);
void cconv (int lx, int ifx, complex *x,
	    int ly, int ify, complex *y,
	    int lz, int ifz, complex *z);
static void closefiles(void);

/* Globals (so can trap signal) defining temporary disk files */
char tracefile[BUFSIZ];	/* filename for the file of traces	*/
char headerfile[BUFSIZ];/* filename for the file of headers	*/
FILE *tracefp;		/* fp for trace storage file		*/
FILE *headerfp;		/* fp for header storage file		*/

/* segy trace */
segy tr;

int
main(int argc, char **argv)
{
        int nsp;		/* number of samples in input trace	*/
	int nspw;		/* number of samples in time window     */
	int nspws;	       	/* number of samples in window (wo taper)*/
	int nspwf;	       	/* number of samples in first window    */
	int nspwi;	       	/* number of samples in intermed. window */
	int nstaper;	       	/* number of samples in taper		*/
	int ntrc;	       	/* number of input traces		*/
	int ntrw;	       	/* number of traces in window		*/
	int ntrwu;	       	/* updated number of traces in window	*/
	int ntw;	       	/* number of time windows		*/
	int nwx;	       	/* number of spacial windows		*/
	int nfft;		/* transform length			*/
	int nf;			/* number of frequencies		*/
	int ntrf;		/* number of traces for filter		*/
	int jx,jw,jr,itt;	/* summation indexes			*/
	int it,iw,ifq,ir;	/* summation indexes			*/
	int ig,ifv;		/* summation indexes			*/
	int verbose;		/* flag to get advisory messages	*/
	int *ipvt;		/* indices of pivot permutations	*/

	float *info;		/* index of last zero pivot		*/
	float dt;		/* sampling interval in secs		*/
	float d1;		/* sample interval in Hz		*/
	float fmin;		/* minimum frequency to process in Hz	*/
	float fmax;		/* maximum frequency to process in Hz	*/
	float taper;		/* length of taper			*/
	float twlen;       	/* time window length 			*/
	float **tdatain;      	/* real trace			   	*/
	float **tdataout;     	/* real trace			   	*/
	float *tidataw;       	/* real trace in time window - input   	*/
	float *ttidataw;      	/* real trace in time window - input   	*/
	float *todataw;       	/* real trace in time window - output   */
	float *ttodataw;      	/* real trace in time window - output   */
	float *rautoc;       	/* real part of autocorrelation   	*/
	float *iautoc;       	/* imaginary part of autocorretation 	*/
	float **rmatrix;      	/* complex autocorrelation matrix	*/
	float *gvector;       	/* advanced autocorrelation vector	*/

	complex *ct;		/* complex transformed trace (window)  	*/
	complex **fdata;	/* data - freq. domain          	*/
	complex **fdataw;	/* data - freq. domain - in window     	*/
	complex *sfreq;		/* single frequency vector 	  	*/
	complex **ffreq;	/* filtered frequency vector 	  	*/
	complex *freqv;		/* frequency vector		      	*/
	complex *sfreqout;	/* single frequency output vector      	*/
	complex *sfreqcj;	/* complex conjugate of sfreq	  	*/
	complex *autocorr;	/* autocorrelation output	  	*/
	complex *fvector;	/* filter               	  	*/

	char *tmpdir;		/* directory path for tmp files		*/
	cwp_Bool istmpdir=cwp_false;/* true for user-given path	*/
 
	/* Initialize */
	initargs(argc, argv);
	requestdoc(1);

	/* Get info from first trace */ 
	if (!gettr(&tr))  err("can't get first trace");
	nsp = tr.ns;

	if (!getparint("verbose", &verbose))	verbose=0;

	if (verbose) warn("nsp = %i",nsp);

	/* dt is used only to set output header value d1 */
	if (!getparfloat("dt", &dt))	dt = ((double) tr.dt)/1000000.0;
	if (!dt) {
	        dt = .004;
		if (verbose) warn("dt not set, assumed to be .004");
	}

	if (verbose) warn("dt = %f",dt);

	if (!getparfloat("taper", &taper)) taper=.1;
	if (taper==0.0) taper=.004; 

	if (!getparfloat("fmin", &fmin)) fmin= 6.;
	if (fmin==0.0) 
	        if (verbose) warn("using fmin=6 Hz");

	if (!getparfloat("fmax", &fmax)) fmax=.6/(2*dt);
	if (fmax==0.0) 
	        if (verbose) warn("using fmax=.6/(2*dt) Hz");

	if (!getparfloat("twlen", &twlen)) twlen=(float)(nsp-1)*dt;
	if (twlen<.3) {
	        twlen=.3;
	        if (verbose) warn("twlen cannot be less than .3s, using .3s");
	}

	/* setting taper and spatial and temporal windows */
	nstaper=NINT(taper/dt);
	nstaper = (nstaper%2 ? nstaper+1 : nstaper); 

	ntw = NINT((nsp-1)*dt/twlen);

	if (ntw==1) taper=0.0;

	nspws = NINT(twlen/dt) + 1;
	nspwf = nspws + nstaper/2;
	nspwi = nspws + nstaper;

	if (!getparint("ntrw", &ntrw)) ntrw=10;
	if (!ntrw) {
	        ntrw = 10;
		if (verbose) warn("ntrw cannot be zero, using 10 traces");
	}

	if (verbose) warn("ntrw = %i",ntrw);

	ntrwu = ntrw;

	if (!getparint("ntrf", &ntrf)) ntrf=4;
	if (!ntrf) 
	        if (verbose) warn("using ntrf=4");

	/* Look for user-supplied tmpdir */
	if (!getparstring("tmpdir",&tmpdir) &&
	        !(tmpdir = getenv("CWP_TMPDIR"))) tmpdir="";
	if (!STREQ(tmpdir, "") && access(tmpdir, WRITE_OK))
	        err("you can't write in %s (or it doesn't exist)", tmpdir);

	/* store traces and headers in tempfiles while getting a count */
	if (STREQ(tmpdir,"")) {
	        tracefp = etmpfile();
		headerfp = etmpfile();
		if (verbose) warn("using tmpfile() call");
	} else { /* user-supplied tmpdir */
	        char directory[BUFSIZ];
		strcpy(directory, tmpdir);
		strcpy(tracefile, temporary_filename(directory));
		strcpy(headerfile, temporary_filename(directory));
		/* Trap signals so can remove temp files */
		signal(SIGINT,  (void (*) (int)) closefiles);
		signal(SIGQUIT, (void (*) (int)) closefiles);
		signal(SIGHUP,  (void (*) (int)) closefiles);
		signal(SIGTERM, (void (*) (int)) closefiles);
		tracefp = efopen(tracefile, "w+");
		headerfp = efopen(headerfile, "w+");
		istmpdir=cwp_true;		
		if (verbose) warn("putting temporary files in %s", directory);
	}

	ntrc = 0;
	do {
	        ++ntrc;
		efwrite(&tr,HDRBYTES,1,headerfp);
		efwrite(tr.data, FSIZE, nsp, tracefp);
	} while (gettr(&tr));
	erewind(tracefp);
	erewind(headerfp);
	if (verbose) warn("%i traces input",ntrc);

	nwx = ntrc/ntrw;

	/* Set up pfa fft */
	nfft = npfaro(2*nspwf, LOOKFAC * 2*nspwf);

/*
	if (nfft >= SU_NFLTS || nfft >= PFA_MAX)  
	        err("Padded nsp=%d--too big", nfft);
*/
	if (nfft >= PFA_MAX)  
	        err("Padded nsp=%d--too big", nfft);

	nf = nfft/2 + 1;
	d1 = 1.0/(nfft*dt);

	/* space allocation */
	ct = ealloc1complex(nfft);
	tidataw=ealloc1float(nfft);
	ttidataw=ealloc1float(nfft);
	fdata = ealloc2complex(nf,ntrc);
	fdataw = ealloc2complex(nf,2*ntrw+2*ntrf);
	tdatain = ealloc2float(nsp,ntrc);
	tdataout = ealloc2float(nsp,ntrc);
	sfreq = ealloc1complex(2*ntrw+2*ntrf);
	sfreqcj = ealloc1complex(2*ntrw+2*ntrf);
	autocorr = ealloc1complex(ntrf+1);
	rautoc = ealloc1float(ntrf+1);
	iautoc = ealloc1float(ntrf+1);
	rmatrix = ealloc2float(2*ntrf,2*ntrf);
	gvector = ealloc1float(2*ntrf);
	fvector = ealloc1complex(2*ntrf+1);
	ffreq = ealloc2complex(nf,2*ntrw);
	sfreqout = ealloc1complex(2*ntrw);
	freqv = ealloc1complex(nfft);
	todataw = ealloc1float(nfft);
	ttodataw = ealloc1float(nfft);
	ipvt = ealloc1int(4*ntrw);
	info = ealloc1float(4*ntrw);

	/* zero output file */
	memset((void *) tdataout[0], 0, ntrc*nsp*FSIZE);

	/* load traces into the zero-offset array and close tmpfile */
	efread(*tdatain, FSIZE, nsp*ntrc, tracefp);
	efclose(tracefp);

	/* If dt not set, issue advisory on frequency step d1 */
	if (dt && verbose)  warn("d1=%f", 1.0/(nfft*dt));

	if (verbose)
	   warn("nf=%i, d1=%f, nfft=%i, nstaper=%i",
			nf,d1,nfft,nstaper);