suk1k2filter.c

su 的源代码库

		/* default is a trapezoidal bandpass filter */
		for(iamps=1; iamps<namps1-1; iamps++) amps1[iamps]=1.;

		amps1[0]=0.;
		amps1[namps1-1]=0.;
        }
	if (!(namps1==npoly1)) 
		err("number of k1 values must = number of amps1 values");

        /* Check amps1 values */
        for(iamps = 0, icount=0; iamps < namps1 ; iamps++) {
		icount+=amps1[iamps];
                if( amps1[iamps] < 0.) err("amps1 values must be positive");
        }
        if (icount==0) err("All amps1 values are zero");
        for(iamps = 0, icount=0; iamps < namps1-1 ; ++iamps) {
			if(!(amps1[iamps]==amps1[iamps+1])) ++icount;
	}
        if (icount==0) warn("All amps1 values are the same");


        /* Get wavenumbers that define the k2 filter */
        if ((npoly2 = countparval("k2"))!=0) {
                k2 = ealloc1float(npoly2);
                getparfloat("k2",k2);
        } else {
                npoly2 = 4;
                k2 = ealloc1float(npoly2);

                k2[0] = FRAC0 * nyq2;
                k2[1] = FRAC1 * nyq2;
                k2[2] = FRAC2 * nyq2;
                k2[3] = FRAC3 * nyq2;
        }

	/* Check k2 values */
	if(npoly2 < 2) warn("Only %d value defining filter",npoly2);
        for(ik=0; ik < npoly2 - 1; ++ik)
		if(k2[ik] < 0.0 || k2[ik] > k2[ik+1])
                                err("Bad filter parameters");

	/* Get k2 filter amplitude values */
        if ((namps2 = countparval("amps2"))!=0) {
                amps2 = ealloc1float(namps2);
                getparfloat("amps2",amps2);
        } else {
                namps2 = npoly2;
                amps2 = ealloc1float(namps2);

		/* default is a trapezoidal bandpass filter */
		for(iamps=1; iamps<namps2-1; ++iamps) amps2[iamps]=1.;
		amps2[0]=0.;
		amps2[namps2-1]=0.;
        }
	if (!(namps2==npoly2)) 
		err("number of k2 values must = number of amps2 values");
        
	
        /* Check amps2 values */
        for(iamps = 0, icount=0; iamps < namps2 ; ++iamps) {
		icount+=amps2[iamps];
                if( amps2[iamps] < 0.) err("amps2 values must be positive");
        }
        if (icount==0) err("All amps2 values are zero");
        for(iamps = 0, icount=0; iamps < namps2-1 ; ++iamps) {
			if(!(amps2[iamps]==amps2[iamps+1])) ++icount;
	}
        if (icount==0) warn("All amps2 values are the same");

        if (!getparint("quad",&quad))		quad=0;

	/* Allocate space */
	rt = alloc2float(nx1fft, nx2fft);
	ct = alloc2complex(nK1,nx2fft);
	kfilter = alloc2float(nx1fft,nx2fft);
	k1filt = alloc1float(nK1);
	k2filt = alloc1float(nK2);

	/* Zero all arrays */
	memset((void *) rt[0], 0, nx1fft*nx2fft*FSIZE);
	memset((void *) kfilter[0], 0, nx1fft*nx2fft*FSIZE);
	memset((void *) ct[0], 0, nK1*nx2fft*sizeof(complex));
	memset((void *) k1filt, 0, nK1*FSIZE);
	memset((void *) k2filt, 0, nK2*FSIZE);

	/* Build Filters */
	polygonalFilter(k1,amps1,npoly1,nx1fft,dx1,k1filt);
	polygonalFilter(k2,amps2,npoly2,nx2fft,dx2,k2filt);

	/* There are only positive k1 values, but both pos and neg k2 values */	
	/* positive k1, positive k2 */
	if (quad==0 || quad==1)
		for (ik2=0; ik2<nK2; ++ik2) 
			for (ik1=0; ik1<nK1; ++ik1) 
				kfilter[ik2][ik1]=k1filt[ik1]*k2filt[ik2];

	/* positive k1, negative k2 */
	if (quad==0 || quad==2)
		for (ik2=nK2; ik2<nx2fft; ++ik2) 
			for (ik1=0; ik1<nK1; ++ik1) 
				kfilter[ik2][ik1]=k1filt[ik1]*k2filt[nx2fft-ik2];

	/* Load traces into fft arrays and close tmpfile */
	rewind(tracefp);
	for (ix2=0; ix2<nx2; ++ix2)
		efread(rt[ix2], FSIZE, nx1, tracefp);
	
	/* Fourier transform dimension 1 */
	pfa2rc(-1,1,nx1fft,nx2,rt[0],ct[0]);

	/* Fourier transform dimension 2 */
	pfa2cc(-1,2,nK1,nx2fft,ct[0]);

	/* Apply filter */
	for (ik2=0; ik2 < nx2fft; ++ik2)
		for (ik1=0; ik1 < nK1; ++ik1)
			ct[ik2][ik1] = crmul(ct[ik2][ik1], kfilter[ik2][ik1]);

	/* Inverse Fourier transform dimension 2 */
	pfa2cc(1,2,nK1,nx2fft,ct[0]);
	
	/* Inverse Fourier transform dimension 1 */
	pfa2cr(1,1,nx1fft,nx2,ct[0],rt[0]);

	erewind(headerfp);
	/* Output filtered traces */
	for (ix2=0; ix2 < nx2; ++ix2) { 

		efread(&tr, HDRBYTES, 1, headerfp);
		for (ix1=0; ix1 <nx1 ; ++ix1) 
			tr.data[ix1] = rt[ix2][ix1];

		puttr(&tr);
	}

	/* Clean up */
	efclose(headerfp);
	if (istmpdir) eremove(headerfile);
	efclose(tracefp);
	if (istmpdir) eremove(tracefile);

	return(CWP_Exit());
}

void polygonalFilter(float *f, float *amps, int npoly,
				int nfft, float dt, float *filter)
/*************************************************************************
polygonalFilter -- polygonal filter with sin^2 tapering
**************************************************************************
Input:
f		array[npoly] of frequencies defining the filter
amps		array[npoly] of amplitude values
npoly		size of input f and amps arrays
dt		time sampling interval
nfft		number of points in the fft

Output:
filter		array[nfft] filter values
**************************************************************************
Notes: Filter is to be applied in the frequency domain
**************************************************************************
Author:  CWP: John Stockwell   1992
*************************************************************************/
#define PIBY2   1.57079632679490
{
        int *intfr;             /* .... integerizations of f		*/
        int icount,ifs;		/* loop counting variables              */
	int taper=0;		/* flag counter				*/
        int nf;                 /* number of frequencies (incl Nyq)     */
        int nfm1;               /* nf-1                                 */
        float onfft;            /* reciprocal of nfft                   */
        float df;               /* frequency spacing (from dt)          */

        
	intfr=alloc1int(npoly);

        nf = nfft/2 + 1;
        nfm1 = nf - 1;
        onfft = 1.0 / nfft;

        /* Compute array of integerized frequencies that define the filter*/
        df = onfft / dt;
        for(ifs=0; ifs < npoly ; ++ifs) {
                intfr[ifs] = NINT(f[ifs]/df);
                if (intfr[ifs] > nfm1) intfr[ifs] = nfm1;
        }

	/* Build filter, with scale, and taper specified by amps[] values*/
	/* Do low frequency end first*/
	for(icount=0; icount < intfr[0] ; ++icount) 
		filter[icount] = amps[0] * onfft;

	/* now do the middle frequencies */
	for(ifs=0 ; ifs<npoly-1 ; ++ifs){
	   if(amps[ifs] < amps[ifs+1]) {	
		++taper;
		for(icount=intfr[ifs]; icount<=intfr[ifs+1]; ++icount) {
		    float c = PIBY2 / (intfr[ifs+1] - intfr[ifs] + 2);
		    float s = sin(c*(icount - intfr[ifs] + 1));
		    float adiff = amps[ifs+1] - amps[ifs];
		    filter[icount] = (amps[ifs] + adiff*s*s) * onfft;
		}
	   } else if (amps[ifs] > amps[ifs+1]) {	
		++taper;
		for(icount=intfr[ifs]; icount<=intfr[ifs+1]; ++icount) {
			   float c = PIBY2 / (intfr[ifs+1] - intfr[ifs] + 2);
                	   float s = sin(c*(intfr[ifs+1] - icount + 1));
			   float adiff = amps[ifs] - amps[ifs+1];
                	   filter[icount] = (amps[ifs+1] + adiff*s*s) * onfft;
		  }
	   } else 
		if(!(taper)){
		for(icount=intfr[ifs]; icount <= intfr[ifs+1]; ++icount)
		   	   filter[icount] = amps[ifs] * onfft;
		} else {
		for(icount=intfr[ifs]+1; icount <= intfr[ifs+1]; ++icount)
		   	   filter[icount] = amps[ifs] * onfft;
		}
	}

	/* finally do the high frequency end */
	for(icount=intfr[npoly-1]+1; icount<nf; ++icount){
		filter[icount] = amps[npoly-1] * onfft;
	}

}


/* for graceful interrupt termination */
static void closefiles(void)
{
	efclose(headerfp);
	efclose(tracefp);
	eremove(headerfile);
	eremove(tracefile);
	exit(EXIT_FAILURE);
}