succwt.c

su 的源代码库

		ia = 0;
		for (ioct=0; ioct<=noct-1; ioct++) {
			for (iv=0; iv<=nv-1; iv++) {
				e = ioct + (float) iv/nv;
				a[ia] = pow(2.0,e);
				if (verbose==1) fprintf(stderr,"aval %d %f %f \n",
							ia,a[ia],e);
				ia += 1;
			}
		}
	}
        
        /* echo extra info if verbose is requested */
	if (verbose==1) warn("fmax = %g, fmin = %g\n", fmax, fmax/a[na]);
	if (verbose==1) warn("scale_index / scale_value / frequency");
	if (verbose==1) {
		for (ia=0; ia<na; ++ia) {
			warn("  %i / %g / %g", ia, a[ia],fmax/a[ia]);
		}
	}

	/*construct basis wavelets*/
	for (ia=0; ia<na; ++ia) {
	if (k == 2) {
			cfourier(a[ia],fmax,c,p,dt,npad,ctmp);
	} else {
			cmorlet(a[ia],fmax,c,p,dt,npad,ctmp);
	}
	ishift=npad/2 ;
	for(i=0; i < ishift; i++){
		basis[ia][i]=ctmp[i+ishift];
		basis[ia][i+ishift]=ctmp[i];
	}
		pfacc(1,  npad, basis[ia]);
	}
        
        /* initialize counters for output traces and cwt spectra */
        itracl = tr.tracl;  /* preserve tracl */
        icdp = 1;
        
        /* loop over traces */
	do {
            
            /* load input trace into real part of complex trace */
            for (it=0,rms=0;it<nt;it++) {
                ctr[it].r = tr.data[it];
                ctr[it].i = 0.0;
                rms+=sqrt(tr.data[it]*tr.data[it]);
            }
            if (verbose==1) fprintf(stderr,"RMS IN %f \n",sqrt((double)rms/((double)nt)));
            for (it=nt;it<npad;it++)  { ctr[it].r=0.0; ctr[it].i=0.0;}
            pfacc(1,  npad, ctr);

            /* loop over scales */
            for (ia=0; ia<na; ia++) {
		for(i=0; i < npad;i++) {
			ctmp[i]= cmul(basis[ia][i],ctr[i]);
		}
		pfacc(-1,  npad, ctmp);
                        
               /* done with this scale, copy header from input trace */
                memcpy( (void *) &cwt, (const void *) &tr, 240);
                
                /* set output header values */
                cwt.tracl = itracl;
                cwt.cdp = icdp;
                cwt.cdpt = ia;
                cwt.unscale = a[ia];
                cwt.sut = (short) noct;
                cwt.gut = (short) noct*nv;
                cwt.trid = FUNPACKNYQ;
		cwt.ns=nt*2;

		memcpy((void*)cwt.data,(const void*)ctmp,nt*sizeof(float)*2);
                
                /* send output trace on its way */
                puttr(&cwt);
                
                /* bump output trace counter */
                itracl += 1;
            }
            
            /* bump output cwt spectrum counter */
            icdp += 1;
                
        } while (gettr(&tr));
        
        
        return(CWP_Exit());
}


static void cmorlet (float a, float fmax, float c, float p, float dt, int nt, complex *w)
/*****************************************************************************
Compute complex morlet wavelet trace for given translation (b) and scale (a)

    w[t,fmax,a,c] = a^p Exp( - i 2 pi (fmax*t/a) ) Exp( - ( t/(a*c) )^2 )
    
Ref:  Matlab wavelet toolbox www info
http://www.math.mcgill.ca/sysdocs/matlabr12/help/toolbox/wavelet/ch06_a37.html
   
******************************************************************************
Input:
a               wavelet scale
fmax            center frequency of wavelet (default = nyquist)
c               t-domain gaussian taper parameter (default = dt)
p               power of scale in CWT 
dt              time sample rate
nt              number of time samples
w[nt]           complex wavelet trace array
******************************************************************************
Notes: 
******************************************************************************
Author:  Chris Liner, UTulsa, 11/18/2003
******************************************************************************/
{
        float t;    /* time variable */
        float cc;       /* scale factor to normalize */
        float tc;   /* center time of wavelet */
        float fc;   /* wavelet center frequency at this scale */
        int it;     /* time index */
        float arg;        /* real argument of exponential */
        complex exparg;   /* complex exponential E^z */
        float rexp;       /* real exponential E^x */
         
        /* center at time index nt/2 - 1 */
        tc = (nt/2 - 1)*dt;
        
        /* constant scale factor 
           ... normalized to give ImpResp=1.0 at all scales */
	cc = pow(a,p);

	/* wavelet center frequency at this scale */
	fc = fmax / a;
                
        /* loop over time samples to build wavelet */
        for (it=1;it<=nt;++it) {

            	t = (it-1)*dt;

            	/* Set up args for complex exponential */
            	arg = - 2.0 * PI * fc * (t-tc);
            	exparg = cexp(crmul(I, arg));
	    
            	/* real exponential */
		arg = (t-tc)/(a*c);
		arg = - arg*arg;
            	rexp = cc*exp( arg );

		/* product of comples and real exponentials */
            	w[it] = crmul(exparg,rexp); 
        }
}

static void cfourier (float a, float fmax, float c, float p, float dt, int nt, complex *w)
/*****************************************************************************
Compute complex morlet wavelet trace for given translation (b) and scale (a)

    w[t,fmax,a] = Exp( - i 2 pi (fmax*t/a) )
    
******************************************************************************
Input:
a               wavelet scale
fmax            center frequency of wavelet (default = nyquist)
c               t-domain gaussian taper parameter (unused)
p               power of scale in CWT 
dt              time sample rate
nt              number of time samples
w[nt]           complex wavelet trace array
******************************************************************************
Notes: 
******************************************************************************
Author:  Chris Liner, UTulsa, 11/18/2003
******************************************************************************/
{
        float t;    /* time variable */
        float cc;       /* scale factor to normalize */
        float tc;   /* center time of wavelet */
        float fc;   /* wavelet center frequency at this scale */
        int it;     /* time index */
        float arg;        /* real argument of exponential */
        complex exparg;   /* complex exponential E^z */
         
        /* center at time index nt/2 - 1 */
        tc = (nt/2 - 1)*dt;
        
        /* constant scale factor 
           ... normalized to give ImpResp=1.0 at all scales */
	cc = pow(a,p);
        
	/* wavelet center frequency at this scale */
	fc = fmax / a;
                
        /* loop over time samples to build wavelet */
        for (it=1;it<=nt;++it) {

            	t = (it-1)*dt;

            	/* complex exponential */
            	arg = - 2.0 * PI * fc * (t-tc);
            	exparg = cexp(crmul(I, arg));
		w[it] = crmul(exparg,cc);

        }
}