sucwt.c

su 的源代码库

	/* Build filters from summed wavelet */
	for(i=0; i<nscales; ++i) {
		nconv[i] = 1 + (int)(scales[i] * xmax);

		for(j=0; j<nconv[i]; ++j)
			index[j] = 1 + j / (scales[i] * dx);

		for(j=0; j<nconv[i]; ++j)
			f[j] = waveletsum[index[j]-1];

		/* flip left right */
		for(j=0,k=nconv[i]-1; j<nconv[i]; ++j,--k)
			filt[i][j] = f[k];
	}

	/* Verbose warning */
	if(verbose) {
		warn("Convolution Lengths");
		for(i=0; i<nscales; ++i) warn("%d ",nconv[i]);
	}
	if(verbose) warn("%d scales will be used for transforms",nscales);

	/* Get information from first trace */
	if(!gettr(&tr))
		err("Cannot get first trace\n");
	ns = tr.ns;

	/* Allocate temporary storage space */
	rt = ealloc1float(ns);
	qt = ealloc1float(ns);
	tmpdata = ealloc2float(nscales,ns);

	/* Zero out rt and qt */
	memset((void *) rt, 0, ns*FSIZE);
	memset((void *) qt, 0, ns*FSIZE);

	/* Alloc sucwt_buffer for longest convolution */
	sucwt_buff = ealloc1float(ns+nconv[nscales-1]+1);
	
	do {  /* main loop over traces */

		outtr.d2 = waveletinc;
		outtr.f2 = minscale;

		memcpy((void *)&outtr,(const void *)&tr,HDRBYTES);

		/* Apply filters to produce wavelet transform */
		for(i=0; i<nscales; ++i) { /* loop over scales */

			for(j=0; j<ns+nconv[nscales-1]+1; ++j)
			sucwt_buff[j] = 0;

			/* convolve wavelet with data */
			conv(ns,0,tr.data,nconv[i],0,
					filt[i],ns,0,sucwt_buff);

			for(j=0; j<ns; ++j) 
				rt[j] = sucwt_buff[j+nconv[i]/2-1];

			for(j=ns-1; j>0; --j) 
				rt[j] = rt[j] - rt[j-1];

			for(j=0; j<ns; ++j)
				rt[j] = -sqrt(scales[i]) * rt[j];

				/* form the hilbert transform of rt */
				hilbert(ns,rt,qt);

			/* If not holder, then output envelope */
			if (!holder) {
				
				for (j=0 ; j<ns; ++j) {			
			  		outtr.data[j] = sqrt(rt[j]*rt[j] 
						               + qt[j]*qt[j]);
				}

				outtr.cdpt = i + 1;
				puttr(&outtr);
			} else {
				/* compute the modulus */
				for (j=0 ; j<ns; ++j) {			
			  		tmpdata[j][i] = sqrt(rt[j]*rt[j] + qt[j]*qt[j]);
				}

			}
		}


		if (holder) { /* compute the Holder regularity traces */
			float *x;
			float *y;
			float lrcoeff[4];

			x = ealloc1float(nscales);
			y = ealloc1float(nscales);
			
	                /* Compute an estimate of the Lipschitz (Holder)
			* regularity. Following Mallat (1992)	
                        *				
                        * ln | Wf(x,s)| <  ln C + alpha * ln|s|
                        *					
                        * alpha here is the Holder or Lipschitz exponent
                        * s is the length scale, and Wf(x,s) is f in the
                        * wavelet basis.			         
                        *					         
			* Here we merely fit a straight line		 
			* through the log-log graph of the of our wavelet
			* transformed data and return the slope as      
			* the regularity measure. 			
                        *					         
			*/

                	for ( j =0 ; j< ns ; ++j ) {
				int icount=0;
                        	x[0]=0;
                        	for ( i = 1 ; i < nscales ; ++i ) {
				
					
			          /* We stay away from values that will make */
				  /*  NANs in the output */
				  if ((i>1) && 
				      (tmpdata[j][i-1] - tmpdata[j][1] > 0.0)) {
                               		y[icount] = log(ABS(tmpdata[j][i]
                                     	              - tmpdata[j][1]));
                                			x[icount] = log(scales[i]-scales[1]);
						
						++icount;
					}

                        	}
				--icount;

				/* straight line fit, return slope */
				if ((icount> 10) && (divisor==1.0) ) {
                        	   linear_regression(y, x, icount, lrcoeff);
                        	   /* lrcoeff[0] is the slope of the line */
				   /* which is the Holder (Lipschitz) */
				   /* exponent */

                        	   outtr.data[j] = lrcoeff[0];

				} else  if ((icount> 10) && (divisor>1.0) ) {

				   float maxalpha=0.0;
				   float interval=icount/divisor;

				   for ( k = interval; k < icount; k+=interval){
                        	   	   linear_regression(y, x, k, lrcoeff);
					   maxalpha = MAX(lrcoeff[0],maxalpha);
					}
				   outtr.data[j] = maxalpha;		

				} else if ((icount < 10) && (divisor>=1.0)) {
				   outtr.data[j] = 0.0;		
				} else if ( divisor < 1.0 ) {
				   err("divisor = %f < 1.0!", divisor);	
				}
				



                	}

			puttr(&outtr); /* output holder regularity traces */
		}
	} while(gettr(&tr));

	return(CWP_Exit());
}

void
MexicanHatFunction(int nwavelet, float xmin, float xcenter, 
			float xmax, float sigma, float *wavelet) 
/***********************************************************************
MexicanHat - the classic Mexican hat function of length nwavelet 
************************************************************************
Input:
nwavelet	number of points total
xmin		minimum  x value
xcenter		central x value
xmax		maximum  x value
Output:
wavelet[]	array of floats of length nwavelet
************************************************************************
Notes:  

Gaussian:  g(x) = 1
		------------------- exp[ - ( x - xcenter)^2/2*sigma^2) ]
	       sigma*sqrt(2 * PI)

1st derivative of Gaussian:
g'(x) =  -(x - xcenter) 
	------------------- exp[ - (x - xcenter)^2/(2*sigma^2) ]
        (sigma^3)*sqrt(2 * PI)

Mexican Hat (2nd derivative of a Gaussian):
g''(x) = 1  
	--- 
      (sigma^3)* sqrt(2*PI)

           / ( x - xcenter)^2         \
	* |----------------   -  1  |
	   \ (sigma^2)              /
   
        * exp[ - (x - xcenter)^2/(2*sigma^2) ]

3rd derivative of Gaussian:
g'''(x) = 1  
	--- 
      (sigma^3)* sqrt(2*PI)

	* (x - xcenter)

           /      3 * ( x - xcenter)^3   \
	* |  1 -  ----------------        |
	   \         (sigma^2)           /
   
        * exp[ - (x - xcenter)^2/(2*sigma^2) ]

4th derivative of Gaussian (Witches' hat)
 (iv)
g  (x) =  1
	----
      (sigma^5)* sqrt(2*PI)
  
	   /      10 *( x - xcenter)^2         3 * ( x - xcenter)^4    \
        * | 1 -  -----------------------   +  ----------------------   |
           \           (sigma^2)                   (sigma^4)           /

        * exp[ - (x - xcenter)^2/(2*sigma^2) ]

************************************************************************
Author: CWP: John Stockwell (Nov 2004) 
************************************************************************/
{
	int i;			/* counter */
	double dxwavelet;	/* sampling interval in x */
	double mult;		/* multiplier on function */
	double twopi=2*PI;	/* 2 PI */

	/* check for error in nwavelet */
	if (nwavelet<=1) err("nwavelet must be greater than 1!");

	/* increment in x */
	dxwavelet = (xmax - xmin)/( nwavelet - 1);


	/* generate mexican hat function */
        /* ... multiplier ....*/
        mult = (1.0/(sigma*sigma*sigma * sqrt(twopi)));
	for(i=0; i<nwavelet; ++i) {
		float x = i*dxwavelet - xcenter; 

		wavelet[i] = mult * (x*x/(sigma*sigma) - 1.0) 
				* exp(- x*x/(2.0*sigma*sigma) );
	}
}