suradon.c

su 的源代码库

	crt=(complex *)VNDemalloc(nmax,
		"inverse_p_transform:crt");
	rt=(float *)crt;
	ctemp=(complex *)VNDemalloc(np*sizeof(complex),
		"inverse_p_transform:ctemp");

	fac=1./ntfft;
	ntfftny=ntfft/2+1;
	for(ip=0;ip<np;ip++) {
		V2Dr0(vnda,ip,(char *)rt,301);
		for(it=0;it<ntfft;it++) rt[it]*=fac;
		pfarc(1,ntfft,rt,crt);
		V2Dw0(vnda,ip,(char *)crt,302);
	}
	VNDr2c(vnda);

	fac=1./np;
	for(iw=0;iw<ntfftny;iw++) {
		wa=freqweight(iw,df,f1,f2);
		if(wa>0.) {
			w=iw*dw;
			V2Dr1(vnda,iw,(char *)crt,303);
			if(wa<1.) {
				for(ip=0;ip<np;ip++) crt[ip]=crmul(crt[ip],wa);
			}
			for(ip=0;ip<np;ip++) ctemp[ip]=crt[ip];
			for(ix=0;ix<nx;ix++) {
			    rsum = isum = 0.;
			    for(ip=ip1;ip<np;ip++) {
				p = pmin + ip*dp;
				tr = cos(w*p*g[ix]);
				ti = sin(w*p*g[ix]);
				dr = ctemp[ip].r;
				di = ctemp[ip].i;
				rsum += tr*dr - ti*di;
				isum += tr*di + ti*dr;
			    }
			    crt[ix].r   = fac*rsum;
			    crt[ix].i   = fac*isum;
			}
		}else{
			for(ix=0;ix<nx;ix++) crt[ix]=czero;
		}
		V2Dw1(vnda,iw,(char *)crt,304);
	}
	for(ix=0;ix<nx;ix++) {
		V2Dr0(vnda,ix,(char *)crt,305);
		pfacr(-1,ntfft,crt,rt);
		V2Dw0(vnda,ix,(char *)rt,306);
	}
	VNDc2r(vnda);
	VNDfree(crt,"inverse_p_transform: crt");
	VNDfree(ctemp,"inverse_p_transform: ctemp");
	return;
}

static void compute_r( float w, int nx, float *g, int np, float dp, complex *r)
/*******************************************************************
Compute the top row of the Hermitian Toeplitz Matrix
			+
		  R = B B

		  i w p g(x)
where B = (1/np) e	    for equal increments in p as

     +           -i w p g(x)
and B = (1/nx) e 

as used for the Discrete Radon Transform computation for
linear or parabolic tau-p.


		 nx-1	i w j dp g(x )
r[j] = 1/(nx*np) Sum	e	    k
		 k=0
						  2
g(x ) is initialized to  x  for linear tau-p or x   for the parabolic transform
   k		          k		         k
prior to calling this routine.  The use of g is intended to emphasize that the
spatial locations do not have to be equally spaced for either method.
In general, this routine can be called for g specified as any function
of spatial position only.  For a more general function of x, dp will
not correspond to an increment in slowness or slowness squared but
rather to a more general parameter.

******************************************************************
Function parameters:

float w	input as angular frequency component of interest
int   nx      number of spatial positions stored in g
float g[]     spatial function for this Radon Transform
int   np      number of slowness (or slowness squared) components
float dp      increment in slownes (or slowness squared)
float r[]     output vector of { real r0, imaginary r0, real r1, 
	      imaginary r1, ...}
******************************************************************
Author: John Anderson (visitor to CSM from Mobil) Spring 1993
******************************************************************/
{
	int j,k;
	float rsum, isum, fac;
	fac= 1./(nx*np);
	for(j=0;j<np;j++) {
		rsum=0.;
		isum=0.;
		for(k=0;k<nx;k++) {
			rsum = rsum+cos( w*j*dp*g[k] );
			isum = isum+sin( w*j*dp*g[k] );	
		}
		r[j].r    = fac*rsum;
		r[j].i    = fac*isum;
	}
}

static void compute_rhs( float w, int nx, float *g, complex *data, int np, 
		float pmin, float dp, complex *rhs)
/*********************************************************************
				     +
Compute the right-hand-side vector  B  data(x)

	+	    -i w p g(x)
where B   = (1/nx) e	        for equal increments in p as
used for the Discrete Radon Transform computation for
linear or parabolic tau-p.

Function parameters:

float w	input angular frequency of interest
int   nx	number of spatial positions ( defines length of g and data )
float g[]      spatial function corresponding to spatial locations of data
complex data[] data as a function of spatial position for a single
		angular frequency w as complex values 
int   np	number of output slownesses p (may be slowness squared
		or a more general function)
float pmin     starting value of output p
float dp	increment in output p
complex rhs[]  np complex values for the result 
*********************************************************************
Author: John Anderson (visitor to CSM from Mobil) Spring 1993
*********************************************************************/

{
	int ip, ix;
	float p, rsum, isum, dr, di, tr, ti, fac;
	fac=1./nx;
	for(ip=0;ip<np;ip++) {
		p = pmin + ip*dp;
		rsum = isum = 0.;
		for(ix=0;ix<nx;ix++) {
			tr = cos(w*p*g[ix]);
			ti = -sin(w*p*g[ix]);
			dr = data[ix].r;
			di = data[ix].i;
			rsum += tr*dr - ti*di;
			isum += tr*di + ti*dr;
		}
		rhs[ip].r   = fac*rsum;
		rhs[ip].i   = fac*isum;
	}
}

static int ctoep( int n, complex *r, complex *a, complex *b,
		 complex *f, complex *g )
/***********************************************************************	
Complex Hermitian Toeplitz Solver for

N-1
Sum  R	     A  = B      for i=0,1,2,...,N-1
j=0   (i-j)   j    i

where R is Hermitian Toeplitz and A and B are complex.  For
an example 4 x 4 system,  A returns as the solution of


   R0  R1  R2  R3	A0	     B0

     *
   R1  R0  R1  R2	A1	     B1
				=    
     *   *
   R2  R1  R0  R1	A2	     B2

     *   *   *
   R3  R2  R1  R0	A3	     B3

and


   R0  R1  R2  R3	F0	     1

     *
   R1  R0  R1  R2	F1	     0
				=    
     *   *
   R2  R1  R0  R1	F2	     0

     *   *   *
   R3  R2  R1  R0	F3	     0


***********************************************************************	
where the function parameters are defined by

n     dimension of system
*r    provides the top row of the Hermitian Toeplitz matrix R 
*a    returns the complex solution vector A
*b    input as complex vector B (not changed during call)
*f    returns the complex spiking filter F
      (may be needed later for Simpson's sideways recursion
      if do search for optimum filter lag)
*g    work space of length n complex values

The function value returns as the number of successfully
computed complex filter coefficients (up to n) if successful or
0 if no coefficients could be computed.
***********************************************************************	
Author: John Anderson (visitor to CSM from Mobil) Spring 1993
***********************************************************************/

{
	float er, ei, vr, vi, cr, ci, vsq;
	int j;	  	/*  for the jth iteration, j=0,n-1 	*/
	int k;			/*  for the kth component, k=0,j-1 	*/
	int jmk;		/*  j-k 				*/
	if (r[0].r==0.) return 0;

	f[0].r = 1.0/r[0].r;
	f[0].i = 0.;
	a[0].r = b[0].r/r[0].r;
	a[0].i = b[0].i/r[0].r;
	vr=1.;
	vi=0.;
	vsq=1.;

	for(j=1;j<n;j++) {     	/* iteration loop for iteration j	*/
	/*  	Compute spiking filter that outputs {v,0,0,...} 
		for this iteration step j			*/
		f[j].r=0.;
		f[j].i=0.;
		er=ei=0.;
		for(k=0;k<j;k++) {
			jmk=j-k;
			er+=r[jmk].r*f[k].r+r[jmk].i*f[k].i;
			ei+=r[jmk].r*f[k].i-r[jmk].i*f[k].r;
		}
		cr  = (er*vr - ei*vi)/vsq;
		ci  = (er*vi + ei*vr)/vsq;
		vr  = vr - (cr*er+ci*ei);
		vi  = vi + (cr*ei-ci*er);
		vsq =  vr*vr + vi*vi;
		if (vsq <= 0.) break;
		for(k=0;k<=j;k++) {
			jmk=j-k;
			g[k].r = f[k].r - cr*f[jmk].r - ci*f[jmk].i;
			g[k].i = f[k].i + cr*f[jmk].i - ci*f[jmk].r;		
		}
		for(k=0;k<=j;k++) {
			f[k]=g[k];
		}

		/*  Compute shaping filter for this iteration */
		a[j].r=0.;
		a[j].i=0.;
		er=ei=0.;
		for(k=0;k<j;k++) {
			jmk=j-k;
			er+=r[jmk].r*a[k].r+r[jmk].i*a[k].i;
			ei+=r[jmk].r*a[k].i-r[jmk].i*a[k].r;
		}
		er  = er-b[j].r;
		ei  = ei-b[j].i;
		cr  = (er*vr - ei*vi)/vsq;
		ci  = (er*vi + ei*vr)/vsq;
		for(k=0;k<=j;k++) {
			jmk=j-k;
			a[k].r += - cr*f[jmk].r - ci*f[jmk].i;
			a[k].i += + cr*f[jmk].i - ci*f[jmk].r;		
		}	
	}

	/* Properly normalize the spiking filter so that R F = {1,0,0,...} */
	/* instead of {v,0,0,...}.  To be accurate, recompute vr,vi,vsq */ 
	vr=vi=0.;
	for(k=0;k<j;k++) {
		vr+=r[k].r*f[k].r-r[k].i*f[k].i;
		vi+=r[k].r*f[k].i+r[k].i*f[k].r;
	}

	vsq = vr*vr + vi*vi;

	/*  Compute (er,ei) = 1./(vr,vi)   */
	er = vr/vsq;
	ei = -vi/vsq;
	for(k=0;k<j;k++) {
		f[k].r = er*f[k].r - ei*f[k].i;
		f[k].i = er*f[k].i + ei*f[k].r;	
	}
	return (j);
}

static int ctoephcg( int niter, int n, complex *a, complex *x, complex *y, 
	complex *s, complex *ss, complex *g, complex *rr)

/*********************************************************************

Hestenes and Stiefel conjugate gradient algorithm 
specialized for solving Hermitian Toeplitz
system.  a[] is input as a vector defining the only the
top row of A.  x[] is the solution vector returned.
y[] is input.  niter is the maximum number of conjugate 
gradient steps to compute.  The function returns as
the number of steps actually computed.  The other 
vectors provide workspace.

Complex Hermitian Toeplitz Solver for

N-1
Sum  A	     x  = y      for i=0,1,2,...,N-1
j=0   (i-j)   j    i

where A is Hermitian Toeplitz and x and y are complex.  For
an example 4 x 4 system,  x returns as the solution of


   A0  A1  A2  A3	x0	     y0

     *
   A1  A0  A1  A2	x1	     y1
				=    
     *   *
   A2  A1  A0  A1	x2	     y2

     *   *   *
   A3  A2  A1  A0	x3	     y3

********************************************************************
Author: John Anderson (visitor to CSM from Mobil) Spring 1993
*********************************************************************/
{
	int j, iter;
	complex czero;
	float alpha, beta, gamma, gammam, rsq, rp, test;
	float eps=1.0e-6;
	float rcdot(int n, complex *a, complex *b);
	rp   = rcdot(n,y,y);
	test = n*eps*eps*rp;
	czero.r=czero.i=0.;

	for(j=0;j<n;j++) {
		x[j]=czero;
		rr[j]=y[j];
	}
	htmul(n,a,rr,g);	   /*  adjoint matrix multiply */

	for(j=0;j<n;j++) s[j]=g[j];
	gammam=rcdot(n,g,g);

	for(iter=0;iter<niter;iter++) { /* forward matrix multiply  */
		htmul(n,a,s,ss);  
		alpha  = gammam/rcdot(n,ss,ss);
		for(j=0;j<n;j++) {
			x[j] =cadd(x[j],crmul(s[j],alpha));
			rr[j]=csub(rr[j],crmul(ss[j],alpha));
		}
		rsq = rcdot(n,rr,rr);
		if ( iter>0 && ( rsq==rp || rsq<test ) ) return(iter-1);
		rp = rsq;

		htmul(n,a,rr,g);   /*  adjoint matrix multiply  */
		gamma  = rcdot(n,g,g);
		if (gamma<eps) break;
		beta   = gamma/gammam;
		gammam = gamma;