cwplib.doc

seismic software,very useful

	a. mkdiff - discrete Taylor series approximation to n'th derivative

Notes:
	The abscissae x of a smpled function f(x) can always be
	expressed as x = (j+a)*h, where j is an integer, a is a
	fraction, and h is the sampling interval.  To approximate the
	n'th order derivative fn(x) of the sampled function f(x) at x =
	(j+a)*h, use the m-l+1 coefficients in the output array d[] as
	follows:

	fn(x) = d[0]*f(j-l) + d[1]*f(j-l-1) +...+ d[m-l]*f(j-m)

	i.e., convolve the coefficients in d with the samples in f.

	m-l+1 (the number of coefficients) must not be greater than the
	NCMAX parameter specified in the code.

	For best approximations,
	when n is even, use a = 0.0, l = -m
	when n is odd, use a = 0.5, l = -m-1


Prototype and parameters:

void mkdiff (int n, float a, float h, int l, int m, float d[]);
n	i order of desired derivative (n>=0 && n<=m-l)
a	i fractional distance from integer sampling index
h	i sampling interval
l	i sampling index of first coefficient
m	i sampling index of last coefficient
d	o array of coefficients for n'th order differentiator

-----------------------------------------
	b. mkhdiff - filter approximating the bandlimited half-differentiator

Notes:
	The half-differentiator is defined by

					  pi
	    d[l+j] = sqrt(1/h)/(2pi) * integral dw sqrt(-iw)*exp(-iwj)
		  			 -pi

                                  pi
           = sqrt(2/h)/(2pi) * integral dw sqrt(w)*(cos(wj)-sin(wj))
                                  0 

	for j = -l, -l+1, ... , l.

	An alternative definition is that f'(j) = d(j)*d(j)*f(j), where
	f'(j) denotes the derivative of a sampled function f(j) and *
	denotes a convolution sum.


	The half-derivative g(j) of f(j) may be computed by the
	following sum:

	g(j) = d[0]*f(j+l) + d[1]*f(j+l-1) + ... + d[2*l]*f(j-l)

	The integral over frequency is evaluated numerically using
	Simpson's method.  Although the Filon method of numerical
	integration is more appropriate for this integral, the
	truncation of d[l+j] for |j| > l is probably the greatest
	source of error.  In any case, d[l+j] is cosine-tapered to
	reduce these truncation errors.


Prototype and parameters:

void mkhdiff (float h, int l, float d[]);
h	i sampling interval
l	i half-length of half-differentiator (length = 1+2*l is odd)
d	o array of 1+2*l coefficients for half-differentiator

------------------------------------------------------------------------
 8. General Signal Processing

	a. conv - compute z = x convolved with y

Notes:

           ifx+lx-1
    z[i] =   sum    x[j]*y[i-j]  ;  i = ifz,...,ifz+lz-1
            j=ifx


Prototype and parameters:

void conv (int lx, int ifx, float x[], int ly, int ify, float y[],
	int lz, int ifz, float z[]);
lx	i length of x array
ifx	i sample index of first x
x	i array to be convolved with y
ly	i length of y array
ify	i sample index of first y
y	i array with which x is to be convolved
lz	i length of z array
ifz	i sample index of first z
z	o array containing x convolved with y

-----------------------------------------
	b. xcor - compute z = x cross-correlated with y

Notes:

           ifx+lx-1
    z[i] =   sum    x[j]*y[i+j]  ;  i = ifz,...,ifz+lz-1
            j=ifx


Prototype and parameters:

void xcor (int lx, int ifx, float x[], int ly, int ify, float y[],
	int lz, int ifz, float z[]);
lx	i length of x array
ifx	i sample index of first x
x	i array to be cross-correlated with y
ly	i length of y array
ify	i sample index of first y
y	i array with which x is to be cross-correlated
lz	i length of z array
ifz	i sample index of first z
z	o array containing x cross-correlated with y

-----------------------------------------
	c. hilbert - compute Hilbert transform y of x

Notes:
	The Hilbert transform is computed by convolving x with a
	windowed (approximate) version of the ideal Hilbert
	transformer.


Prototype and parameters:

void hilbert (int n, float x[], float y[]);
n			i length of x and y
x			i function to be Hilbert transformed
y			o Hilbert transform of x

-----------------------------------------
	d. antialias - Anti-alias filter for use before sub-sampling

Notes:
	The anti-alias filter is a recursive (Butterworth) filter.  For
	zero-phase anti-alias filtering, the recursive filter is
	applied forwards and backwards.


Prototype and parameters:

void antialias (float frac, int phase, int n, float p[], float q[]);
frac	i current sampling interval / future interval (should be <= 1)
phase	i =0 for zero-phase filter; =1 for minimum-phase filter
n	i number of samples
p	i array[n] of input samples
q	o array[n] of output (anti-alias filtered) samples		

------------------------------------------------------------------------
 9. Discrete Abel transform package

Contents:

	a. abelalloc - allocate and return a pointer to an Abel transformer
	b. abelfree  - free an Abel transformer
	c. abel	     - compute the Abel transform

Notes:
	The Abel transform is defined by:

			 Infinity
		g(y) = 2 Integral dx f(x)/sqrt(1-(y/x)^2)
			   |y|

	Linear interpolation is used to define the continuous function
	f(x) corresponding to the samples in f[].  The first sample
	f[0] corresponds to f(x=0) and the sampling interval is assumed
	to be 1.  Therefore, the input samples correspond to 0 <= x <=
	n-1.  Samples of f(x) for x > n-1 are assumed to be zero.
	These conventions imply that

		g[0] = f[0] + 2*f[1] + 2*f[2] + ... + 2*f[n-1]

References:

	Hansen, E. W., 1985, Fast Hankel transform algorithm:  IEEE
	Trans. on Acoustics, Speech and Signal Processing, v. ASSP-33,
	n. 3, p. 666-671.  (Beware of several errors in the equations
	in this paper!)


Prototypes:

void *abelalloc (int n);
void abelfree (void *at);
void abel (void *at, float f[], float g[]);

------------------------------------------------------------------------
10. Discrete Hankel transform package

Contents:

	a. hankelalloc - allocate and return a pointer to a Hankel transformer
	b. hankelfree  - free a Hankel transformer
	c. hankel0     - compute the zeroth-order Hankel transform
	d. hankel1     - compute the first-order Hankel transform

Notes:
	The zeroth-order Hankel transform is defined by:

			Infinity
		h0(k) = Integral dr r j0(k*r) f(r)
			   0

	where j0 denotes the zeroth-order Bessel function.

	The first-order Hankel transform is defined by:

			Infinity
		h1(k) = Integral dr r j1(k*r) f(r)
			   0

	where j1 denotes the first-order Bessel function.

	The Hankel transform is its own inverse.

	The Hankel transform is computed by an Abel transform followed
	by a Fourier transform.

References:
	Hansen, E. W., 1985, Fast Hankel transform algorithm:  IEEE
	Trans. on Acoustics, Speech and Signal Processing, v. ASSP-33,
	n. 3, p. 666-671.  (Beware of several errors in the equations
	in this paper!)

Prototypes:
	void *hankelalloc (int nfft);
	void hankelfree (void *ht);
	void hankel0 (void *ht, float f[], float h[]);
	void hankel1 (void *ht, float f[], float h[]);

------------------------------------------------------------------------
11. Sorting and Searching

	a. hpsort - sort an array so that a[0] <= a[1] <= ... <= a[n-1]

Notes:
	Adapted from Standish, T. A., Data Structure Techniques, p. 91.


Prototype and parameters:

void hpsort (int n, float a[]);
n	i number of elements in a
a	i array[n] to be sorted
a	o array[n] sorted

-----------------------------------------
	b. qksort - sort functions based on Hoare's quicksort algorithm

Contents:
	1. qksort - sort an array such that a[0] <= a[1] <= ... <= a[n-1]

	2. qkfind - partially sort an array so that the element a[m] has
		the value it would have if the entire array were sorted
		such that a[0] <= a[1] <= ... <= a[n-1]

Notes:
	This function is adapted from procedure partition by Hoare,
	C.A.R., 1961, Communications of the ACM, v. 4, p. 321.

Prototypes:
void qksort (int n, float a[]);
void qkfind (int m, int n, float a[]);

-----------------------------------------
	c. qkisort - index sort functions based on Hoare's quicksort algorithm

Contents:

	1. qkisort - sort an array of indices i[] so that 
		a[i[0]] <= a[i[1]] <= ... <= a[i[n-1]]

	2. qkifind - partially sort an array of indices i[] so that the 
		     index i[m] has the value it would have if the entire
		     array of indices were sorted such that 
		     a[i[0]] <= a[i[1]] <= ... <= a[i[n-1]]

Notes:
	This function is adapted from procedure partition by Hoare,
	C.A.R., 1961, Communications of the ACM, v. 4, p. 321.


Prototypes:

void qkisort (int n, float a[], int i[]);
void qkifind (int m, int n, float a[], int i[]);

------------------------------------------------------------------------
12. Statistics

	a. quest - return an estimate of a specified quantile 

Notes:
	The estimate should approach the sample quantile in the limit
	of large n.

	The estimate is most accurate for cumulative distribution
	functions that are smooth in the neighborhood of the quantile
	specified by p.

	This function is an implementation of the algorithm published
	by Jain, R. and Chlamtac, I., 1985, The PP algorithm for
	dynamic calculation of quantiles and histograms without storing
	observations:  Comm. ACM, v. 28, n. 10.


Prototype and parameters:

float quest (float p, int n, float x[]);
p	i quantile to be estimated (0.0<=p<=1.0 is required)
n	i number of samples in array x (n>=5 is required)
x	i array of floats

-----------------------------------------
	b. questinit - alloc,init,return pointer to state of quantile estimator 

Notes:
	This function must be called before calling function questUpdate.

	See also notes in questUpdate.

Prototype and parameters:

QuestState *questinit (float p, int n, float x[]);
p	i quantile to be estimated (0.0<=p<=1.0 is required)
n	i number of samples in array x (n>=5 is required)
x	i array of floats

-----------------------------------------
	b. questupdate - update and return a quantile estimate 

Notes:
	The estimate should approach the sample quantile in the limit
	of large n.

	The estimate is most accurate for cumulative distribution
	functions that are smooth in the neighborhood of the quantile
	specified by p.

	This function is an implementation of the algorithm published
	by Jain, R. and Chlamtac, I., 1985, The PP algorithm for
	dynamic calculation of quantiles and histograms without storing
	observations:  Comm. ACM, v. 28, n. 10.


Prototype and parameters:

float questupdate (QuestState *state, int n, float x[]);
s	i pointer to state of quantile estimator
n	i number of samples in array x
x	i array of floats

------------------------------------------------------------------------
13. Single precision BLAS (basic linear algebra subroutines)

Contents:
	a. isamax - return index of element with maximum absolute value
	b. sasum - return sum of absolute values
	c. saxpy - compute y[i] = a*x[i]+y[i]
	d. scopy - copy x[i] to y[i] (i.e., set y[i] = x[i])
	e. sdot - return sum of x[i]*y[i] (i.e., dot product of x and y)
	f. snrm2 - return square root of sum of squares of x[i]
	g. sscal - compute x[i] = a*x[i]
	h. sswap - swap x[i] and y[i]

Notes:
	Adapted from LINPACK FORTRAN.