dgwaveform.c

计算高斯各阶导函数的C程序 Computing Gaussian derivative waveforms of any order. Dgwaveform efficiently cre

/*
 * Copyright (c) 2007 by the Society of Exploration Geophysicists.
 * For more information, go to http://software.seg.org/2007/0004 .
 * You must read and accept usage terms at:
 * http://software.seg.org/disclaimer.txt before use.
 *
 * Revision history:
 * Original SEG version by Werner M. Heigl, Apache E&P Technology,
 * February 2007.
 */
 
#include "su.h"
#include "segy.h"

/*************************** self documentation *************************/
char *sdoc[] = {
" 	   								",
" DGWAVEFORM - make Gaussian derivative waveform			",
" 	   								",
" dgwaveform >stdout  [optional parameters]				",
" 									",
"									",
" Optional parameters:							",
" n=2    	order of derivative (n>=1)				",
" fpeak=35	peak frequency						",
" nfpeak=n*n	max. frequency = nfpeak * fpeak				",
" nt=128	length of waveform					",
" shift=0	additional time shift in s (used for plotting)		",
" sign=1	use =-1 to change sign					",
" verbose=0	=0 don't display diagnostic messages			",
"               =1 display diagnostic messages				",
" Notes:								",
" This code computes a waveform that is the n-th order derivative of a	",
" Gaussian. The variance of the Gaussian is specified through its peak	",
" frequency, i.e. the frequency at which the amplitude spectrum of the	",
" Gaussian has a maximum. nfpeak is used to compute maximum frequency,	",
" which in turn is used to compute the sampling interval. Increasing	",
" nfpeak gives smoother plots. In order to have a (pseudo-) causal	",
" pulse, the program computes a time shift equal to sqrt(n)/fpeak. An	",
" additional shift can be applied with the parameter shift. A positive	",
" value shifts the waveform to the right.				",
"									",
" Examples:								",
" 2-loop Ricker: dgwaveform n=1	>ricker2.su				",
" 3-loop Ricker: dgwaveform n=2 >ricker3.su				",
" Sonic transducer pulse: dgwaveform n=10 fpeak=300 >sonic.su		",
"									",
" To display use suxgraph. For example:					",
" dgwaveform n=10 fpeak=300 | suxgraph style=normal &			",
"									",
NULL};

/* Credits:
 *
 *	Werner M. Heigl, February 2007
 *
*/
/************************ end self doc **********************************/


void
deriv_gauss(double dt, int nt, double t0, float fpeak, int n, double *w,
	    int sign, int verbose);

void
hermite(double *h, double *h0, double *h1, double *t, int nt, int n,
	double sigma);

static segy tr;	/* structure of type segy that contains the waveform */

int
main(int argc, char **argv)
{
	int i;			/* loop variable			*/
	int n;			/* order of derivative			*/
	int nfpeak;		/* fmax = nfpeak*fpeak			*/
	int nt;			/* length of time vector t[]		*/
	int sign;		/* scalar for polarity			*/
	int verbose;		/* flag for diagnostic messages		*/
	float fmax;		/* max. frequency			*/
	float fpeak;		/* peak frequency			*/
	double dt;		/* sampling interval 			*/
	double shift;		/* additional time shift		*/
	double t0;		/* time delay				*/
	double *w = NULL;	/* waveform and Gaussian		*/	

	/* Hook up getpar */
	initargs(argc, argv);
	requestdoc(0);
	
	/* Get parameters and do setup */
	if (!getparint("n",&n))			n = 2;
	if (!getparfloat("fpeak",&fpeak))	fpeak = 35;
	if (!getparint("sign",&sign))		sign = 1;
	if (!getparint("nfpeak",&nfpeak))	nfpeak = n * n;
	if (!getpardouble("shift",&shift))	shift = 0;
	if (!getparint("verbose",&verbose))	verbose = 0;	
	fmax = nfpeak * fpeak;
	dt = 0.5 / fmax;	
	t0 = shift + sqrt(n) / fpeak;
	if (!getparint("nt",&nt))	nt = (int) (2 * t0 / dt + 1);
	warn("n=%d fpeak=%.0f fmax=%.0f t0=%f nt=%d dt=%.12f",n,fpeak,fmax,t0,nt,dt);
	if (dt < 1e-6)	err("single-precision exceeded: reduce nfpeak or fpeak");

	/* allocate & initialize memory */
	w = ealloc1double(nt);
	memset((void *) w, '0', DSIZE * nt);
	if (verbose)	warn("memory for waveform allocated and initialized");
	
	/* compute waveform */
	if (n >= 1) {
		
		/* compute n-th derivative of Gaussian */
		deriv_gauss(dt, nt, t0, fpeak, n, w, sign, verbose);
	
		/* write out waveform */	

			for (i = 0; i < nt; ++i)	tr.data[i] = (float) w[i];
			tr.tracl = 1;
			tr.ns    = nt;
			tr.trid  = 1;
			tr.dt    = NINT(dt*1000000.0);
			tr.ntr   = 1;
			puttr(&tr);
			warn("waveform written to stdout");
			
		
	} else	err("specified n not >=1 !!");

	/* free memory */
	free1double(w);
	if (verbose)	warn("memory freed");
	
	return(CWP_Exit());
}


void
deriv_gauss(double dt, int nt, double t0, float fpeak, int n, double *w,
	    int sign, int verbose)
/****************************************************************************
  deriv_gauss:  Compute n-th order derivative of a Gaussian.
*****************************************************************************
Input:
dt		sampling interval
nt		length of waveform in samples
t0		time shift for (pseudo-) causality
fpeak		maximum frequency
n		order of derivative
sign		multiplier for polarity of waveform
verbose		flag for diagnostic messages
*****************************************************************************
Output:
w	array of size nt containing the waveform
*****************************************************************************
Return: none
*****************************************************************************
Author: Werner M. Heigl, Apache Corporation, E&P Technology, November 2006
*****************************************************************************/
{
	int i;			/* loop variable			*/
	double sigma;		/* temporal variance of Gaussian	*/
	double C;		/* normalization constant		*/
	double *h  = NULL;	/* Hermite polynomial			*/
	double *h0 = NULL;	/* temp array for H_{n-1}		*/
	double *h1 = NULL;	/* temp array for H_{n}			*/
	double *t  = NULL;	/* time vector				*/

	
	/* allocate & initialize memory */
	t  = ealloc1double(nt);
	h  = ealloc1double(nt);
	h0 = ealloc1double(nt);
	h1 = ealloc1double(nt);

	memset((void *) t , '0', DSIZE * nt);
	memset((void *) h , '0', DSIZE * nt);
	memset((void *) h0, '0', DSIZE * nt);
	memset((void *) h1, '0', DSIZE * nt);
	if (verbose)
		warn("memory allocated and initialized");

	/* initialize time vector */
	for (i = 0; i < nt; ++i)	t[i] = i * dt - t0;
	if (verbose)	warn("t[] initialized");
	
	/* compute Gaussian */
	sigma = n / ( 4 * PI * PI * fpeak * fpeak );
	if (verbose)	warn("sigma=%f",sigma);
	for (i = 0; i < nt; ++i)	w[i] = exp( - t[i] * t[i] / (2 * sigma) );
	if (verbose)	warn("Gaussian computed");
	
	/* compute Hermite polynomial */
	for (i = 0; i < nt; ++i) {
		h0[i] = 1.0;
		h1[i] = t[i] / sigma;
	}
	if (n==1)	memcpy((void *) h, (const void *) h1, DSIZE * nt);
	if (n>1)	hermite(h, h0, h1, t, nt, n, sigma);
	if (verbose)	warn("Hermite polynomial H_%d computed",n);
	
	/* compute waveform */
	for (i = 0; i < nt;++i)		w[i] = h[i] * w[i];
	if (verbose)	warn("waveform computed");
	
	/* find normalization constant */
	C = fabs(w[0]);
	for (i = 1; i < nt; ++i)
		if (fabs(w[i]) > C)	C = fabs(w[i]);
	if (ISODD(n))	C = -C;	/* to account for (-1)^n */
	if (verbose)	warn("C=%f",C);
	
	/* and finally normalize */
	for (i = 0; i < nt; ++i)	w[i] = sign * w[i] / C;
	if (verbose)	warn("waveform normalized");

	/* check amplitude a t=0 */
	if (verbose)	warn("w[o]=%.12f",w[0]);
	
	/* free memory */
	free1double(h); free1double(h0); free1double(h1);
	free1double(t);
	if (verbose)	warn("memory freed");
}	

void
hermite(double *h, double *h0, double *h1, double *t, int nt, int n, double sigma)
/*********************************************************************************
  hermite:  Compute n-th order generalized Hermite polynomial.
**********************************************************************************
Input:
h0		array of size nt holding H_{n-1}
h1		array of size nt holding H_{n}
t		array of size nt holding time vector
nt		size of arrays, no. of samples
n		order of polynomial
sigma		variance
**********************************************************************************
Output:
h		array of size nt holding H_{n+1}
**********************************************************************************
Note:	Note that n in the function call is the order of the derivative and
        j in the code below is the n in the recurrence relation
**********************************************************************************
Author: Werner M. Heigl, Apache Corporation, E&P Technology, December 2006
*********************************************************************************/
{
	int i;		/* loop variable */
	int j=1;	/* recurrence counter */

	/* as long as necessary use recurrence relation */
	do {
		/* current instance of recurrence relation */
		for (i = 0; i < nt; ++i) 		
			h[i] = ( t[i] * h1[i] - j * h0[i] ) / sigma;
		
		/* update inputs to recurrence relation */
		memcpy((void *) h0, (const void *) h1, DSIZE * nt);
		memcpy((void *) h1, (const void *) h , DSIZE * nt);
		
		/* update counter */
		++j;
		
	} while (j < n);
	
}