cwplib.doc

seismic software,very useful

Prototype and parameters:

int npfao (int nmin, int nmax);
nmin        i lower bound on returned value
nmax        i desired (but not guaranteed) upper bound on returned value

-----------------------------------------
	c. npfar  - find next allowed fft size for pfarc

Notes:
	npfar returns the smallest valid n not less than nmin for 
	real-to-complex or complex-to-real prime factor ffts

	Current implemenations of real-to-complex and complex-to-real
	prime factor ffts require that the transform length n be even
	and that n/2 be a valid length for a complex-to-complex prime
	factor fft.  The value returned by npfar satisfies these
	conditions.

	Also, see notes for npfa.

PARAMETERS:

int npfar (int nmin);
nmin        i lower bound on returned value

-----------------------------------------
	d. npfaro - find optimal fft size for pfarc
Notes:

	npfaro returns the optimal n between nmin and nmax for
	real-to-complex or complex-to-real prime factor ffts.

	Current implemenations of real-to-complex and complex-to-real
	prime factor ffts require that the transform length n be even
	and that n/2 be a valid length for a complex-to-complex prime
	factor fft.  The value returned by npfaro satisfies these
	conditions.

	Also, see notes for npfao.


Prototype and parameters:

int npfaro (int nmin, int nmax);
nmin        i lower bound on returned value
nmax        i desired (but not guaranteed) upper bound on returned value

-----------------------------------------
	e. pfacc  - prime factor complex-to-complex fft

Notes:
	pfacc is a complex-to-complex prime factor fft, in place.
 
	The fft size, n, must be factorable into mutually prime factors
	taken from the set {2,3,4,5,7,8,9,11,13,16}.  In other words,
		n = 2**p * 3**q * 5**r * 7**s * 11**t * 13**u
	where
		0 <= p <= 4,  0 <= q <= 2,  0 <= r,s,t,u <= 1
	is required for pfa to yield meaningful results.  This
	restriction implies that n is restricted to the range
	1 <= n <= 720720 (= 5*7*9*11*13*16).

References:
	Temperton, C., 1985, Implementation of a self-sorting in-place
	prime factor fft algorithm:  Journal of Computational Physics,
	v. 58, p. 283-299.

	Temperton, C., 1988, A new set of minimum-add rotated
	rotated dft modules: Journal of Computational Physics,
	v. 75, p. 190-198.

Prototype and parameters:

void pfacc (int isign, int n, complex z[]);
isign	i sign of isign is the sign of exponent in Fourier kernel
n	i length of transform
z	b array of complex numbers to be transformed in place

-----------------------------------------
	f. pfarc  - prime factor real-to-complex fft

Notes:
	pfarc is a prime factor real-to-complex fft.
 
	Because pfarc uses pfacc to do most of the work, n must be even
	and n/2 must be a valid length for pfacc.  The simplest way to
	obtain a valid n is via n = npfar(nmin).  A more optimal n can
	be obtained with npfaro.

References:  
	Press et al, 1988, Numerical Recipes in C, p. 417.

	Also, see notes and references for function pfacc.

Prototype and parameters:

void pfarc (int isign, int n, float rz[], complex cz[]);
isign       i sign of isign is the sign of exponent in Fourier kernel
n           i length of transform; must be even
rz          i array of n real values (may be equivalenced to cz)
cz          o array of n/2+1 complex values (may be equivalenced to rz)

-----------------------------------------
	g. pfacr  - prime factor complex-to-real fft

NOTES:
	pfacr is a complex-to-real prime factor fft.

	Because pfacr uses pfacc to do most of the work, n must be even
	and n/2 must be a valid length for pfacc.  The simplest way to
	obtain a valid n is via n = npfar(nmin).  A more optimal n can
	be obtained with npfaro.

References:  
	Press et al, 1988, Numerical Recipes in C, p. 417.

	Also, see notes and references for function pfacc.

Prototype and parameters:
void pfacr (int isign, int n, complex cz[], float rz[]);
isign       i sign of isign is the sign of exponent in Fourier kernel
n           i length of transform
cz          i array of n/2+1 complex values (may be equivalenced to rz)
rz          o array of n real values (may be equivalenced to cz)
 
-----------------------------------------
	h. pfa2cc - 2d prime factor complex-to-complex fft

Notes:
	pfa2cc does multiple 2-D complex to complex ffts, in place.

	Only one (either the 1st or 2nd) dimension of the 2-D array is
	transformed.

	If idim equals 1, then n2 transforms of n1 complex elements are
	performed; else, if idim equals 2, then n1 transforms of n2
	complex elements are performed.

	Although z appears in the argument list as a one-dimensional
	array, z may be viewed as an n1 by n2 two-dimensional array:
	z[n2][n1].

	Let n denote the transform length, either n1 or n2, depending
	on idim.  Then, n must be factorable into mutually prime
	factors taken from the set {2,3,4,5,7,8,9,11,13,16}, thus:
	    n = 2**p * 3**q * 5**r * 7**s * 11**t * 13**u
	where
	    0 <= p <= 4,  0 <= q <= 2,  0 <= r,s,t,u <= 1
	is required for pfa2cc to yield meaningful results.  This
	restriction implies that n is restricted to the range
	    1 <= n <= 720720 (= 5*7*9*11*13*16)

	To perform a two-dimensional transform of an n1 by n2 complex
	array (assuming that both n1 and n2 are valid "n"), stored with
	n1 fast and n2 slow:
		pfa2cc(isign,1,n1,n2,z);
		pfa2cc(isign,2,n1,n2,z);

References:  
	Temperton, C., 1985, Implementation of a self-sorting
	in-place prime factor fft algorithm:  Journal of
	Computational Physics, v. 58, p. 283-299.

	Temperton, C., 1988, A new set of minimum-add rotated
	rotated dft modules: Journal of Computational Physics,
	v. 75, p. 190-198.

Prototypes and parameters:

void pfa2cc (int isign, int idim, int n1, int n2, complex z[]);
isign       i sign of isign is the sign of exponent in Fourier kernel
idim        i dimension to transform, which must be either 1 or 2 (see notes)
n1          i 1st (fast) dimension of array to be transformed (see notes)
n2          i 2nd (slow) dimension of array to be transformed (see notes)
z           b array of complex elements to be transformed in place
 
-----------------------------------------
	i. pfa2rc - 2d prime factor real-to-complex fft

Notes:
	pfa2rc does multiple 2-D real-to-complex transforms.

 
	If idim equals 1, then n2 transforms of n1 real elements to
	n1/2+1 complex elements are performed; else, if idim equals 2,
	then n1 transforms of n2 real elements to n2/2+1 complex
	elements are performed.

	Although rz appears in the argument list as a one-dimensional
	array, rz may be viewed as an n1 by n2 two-dimensional array:
	rz[n2][n1].  Likewise, depending on idim, cz may be viewed as
	either an n1/2+1 by n2 or an n1 by n2/2+1 two-dimensional array
	of complex elements.

	Let n denote the transform length, either n1 or n2, depending
	on idim.  Because pfa2rc uses pfa2cc to do most of the work, n
	must be even and n/2 must be a valid length for pfa2cc.  The
	simplest way to obtain a valid n is via n = npfar(nmin).  A
	more optimal n can be obtained with npfaro.

References:
	Press et al, 1988, Numerical Recipes in C, p. 417.

	Also, see notes and references for function pfa2cc.


Prototypes and parameters:

void pfa2rc (int isign, int idim, int n1, int n2, float rz[], complex cz[]);
isign       i sign of isign is the sign of exponent in Fourier kernel
idim        i dimension to transform, which must be either 1 or 2 (see notes)
n1          i 1st (fast) dimension of array to be transformed (see notes)
n2          i 2nd (slow) dimension of array to be transformed (see notes)
rz          i array of real values (may be equivalenced to cz)
cz          o array of complex values (may be equivalenced to rz)

-----------------------------------------
	j. pfa2cr - 2d prime factor complex-to-real fft

Notes:
	pfa2cr does multiple complex-to-real prime factor ffts.  

	Because pfacr uses pfacc to do most of the work, n must be even
	and n/2 must be a valid length for pfacc.  The simplest way to
	obtain a valid n is via n = npfar(nmin).  A more optimal n can
	be obtained with npfaro.

References:  
	Press et al, 1988, Numerical Recipes in C, p. 417.

	Also, see notes and references for function pfacc.

Prototype and parameters:
void pfa2cr (int isign, int idim, int n1, int n2, complex cz[], float rz[]);
isign       i sign of isign is the sign of exponent in Fourier kernel
n           i length of transform (see notes below)
cz          i array of n/2+1 complex values (may be equivalenced to rz)
rz          o array of n real values (may be equivalenced to cz)

-----------------------------------------
	k. pfamcc - multiple dimension prime factor complex-to-complex fft

Notes:
	pfamcc does multiple complex to complex prime factor ffts, in place.

	n must be factorable into mutually prime factors taken 
	from the set {2,3,4,5,7,8,9,11,13,16}.  Thus:
	    n = 2**p * 3**q * 5**r * 7**s * 11**t * 13**u
	where
	    0 <= p <= 4,  0 <= q <= 2,  0 <= r,s,t,u <= 1
	is required for pfamcc to yield meaningful results.  This
	restriction implies that n is restricted to the range
	    1 <= n <= 720720 (= 5*7*9*11*13*16)

	To perform a two-dimensional transform of an n1 by n2 complex
	array (assuming that both n1 and n2 are valid "n"), stored with
	n1 fast and n2 slow:
	    pfamcc(isign,n1,n2,1,n1,z); (to transform 1st dimension)
	    pfamcc(isign,n2,n1,n1,1,z); (to transform 2nd dimension)
 
References:  
	Temperton, C., 1985, Implementation of a self-sorting
	in-place prime factor fft algorithm:  Journal of
	Computational Physics, v. 58, p. 283-299.

	Temperton, C., 1988, A new set of minimum-add rotated
	rotated dft modules: Journal of Computational Physics,
	v. 75, p. 190-198.

Prototype and parameters:
void pfamcc (int isign, int n, int nt, int k, int kt, complex z[]);
isign       i sign of isign is the sign of exponent in Fourier kernel
n           i number of complex elements per transform
nt          i number of transforms
k           i stride in complex elements within transforms
kt          i stride in complex elements between transforms
z           b array of complex elements to be transformed in place
 

------------------------------------------------------------------------
 5. Interpolation routines


  a. fsinc - single precision sinc function
  b. dsinc - double precision sinc function

float fsinc (float x);
double dsinc (double x);

-----------------------------------------
  c. mksinc - compute least-squares optimal sinc interpolation coefficients

Notes:

	The coefficients are a least-squares-best approximation to the
	ideal sinc function for frequencies from zero up to a computed
	maximum frequency.  For a given interpolator length, lsinc,
	mksinc computes the maximum frequency, fmax (expressed as a
	fraction of the nyquist frequency), using the following
	empirically derived relation (from a Western Geophysical
	Technical Memorandum by Ken Larner):

		fmax = min(0.066+0.265*log(lsinc),1.0)

	Note that fmax increases as lsinc increases, up to a maximum of
	1.0.  Use the coefficients to interpolate a uniformly-sampled
	function y(i) as follows:

		    lsinc-1
	    y(i+d) =  sum  sinc[j]*y(i+j+1-lsinc/2)
		      j=0

	Interpolation error is greatest for d=0.5, but for frequencies
	less than fmax, the error should be less than 1.0 percent.


Prototype and parameters:

void mksinc (float d, int lsinc, float sinc[]);
d		i fractional distance to interpolation point; 0.0<=d<=1.0
lsinc		i length of sinc approximation; lsinc%2==0 and lsinc<=20
sinc		o array containing lsinc interpolation coefficients

-----------------------------------------
  d. ints8r - interpolation of a uniformly-sampled real function y(x)
  e. ints8c - interpolation of a uniformly-sampled complex function y(x)


Notes:
	These routines use a table of 8-coefficient sinc approximations;
	maximum error for frequencies less than 0.6 nyquist is less than
	one percent.

	Because extrapolation of the input function y(x) is defined by
	the left and right values yinl and yinr, the xout values are not
	restricted to lie within the range of sample locations defined by
	nxin, dxin, and fxin.


Prototype and parameters:

void ints8r (int nxin, float dxin, float fxin, float yin[], 
	float yinl, float yinr, int nxout, float xout[], float yout[]);
void ints8c (int nxin, float dxin, float fxin, complex yin[], 
	complex yinl, complex yinr, int nxout, float xout[], complex yout[]);

nxin		i number of x values at which y(x) is input
dxin		i x sampling interval for input y(x)
fxin		i x value of first sample input
yin		i array of input y(x) values:  yin[0] = y(fxin), etc.
yinl		i value used to extrapolate yin values to left of yin[0]
yinr		i value used to extrapolate yin values to right of yin[nxin-1]
nxout		i number of x values a which y(x) is output
xout		i array of x values at which y(x) is output
yout		o array of output y(x) values:  yout[0] = y(xout[0]), etc.

-----------------------------------------
  f. intt8r - interpolation of a uniformly-sampled real function y(x)
  g. intt8c - interpolation of a uniformly-sampled complex function y(x)

Notes:

	These routines use a table of 8-coefficient interpolators.

	Because extrapolation of the input function y(x) is defined by
	the left and right values yinl and yinr, the xout values are
	not restricted to lie within the range of sample locations
	defined by nxin, dxin, and fxin.


Prototype and parameters:

void intt8r (int ntable, float table[][8],
	int nxin, float dxin, float fxin, float yin[], 
	float yinl, float yinr, int nxout, float xout[], float yout[]);
ntable		i number of tabulated interpolation operators; ntable>=2
table		i array of tabulated 8-point interpolation operators