supef.c

su 的源代码库

		
	/* Get parameters  */
	if (!getparint("showwiener",  &showwiener))	showwiener = 0;
	if (!getparfloat("pnoise",  &pnoise))	pnoise = PNOISE;
	
	/* .. mincorr and maxcorr */
	if (getparfloat("mincorr", &mincorr))	imincorr = NINT(mincorr/dt);
	else					imincorr = 0;
	if (imincorr < 0) err("mincorr=%g too small", mincorr);
	
	if (getparfloat("maxcorr", &maxcorr))	imaxcorr = NINT(maxcorr/dt);
	else					imaxcorr = nt;
	if (imaxcorr > nt) err("maxcorr=%g too large", maxcorr);

	if (imincorr >= imaxcorr)
		err("mincorr=%g, maxcorr=%g", mincorr, maxcorr);
	
	/*... Get mix weighting values values */
	if ((nmix = countparval("mix"))!=0) {
		mix = ealloc1float(nmix);
		getparfloat("mix",mix);
		
	} else { /* else use default values */
		nmix = 1;
		mix = ealloc1float(nmix);

		mix[0] = VAL0;
	}

	/* Divide mixing weight by number of traces to mix over */
	for (imix = 0; imix < nmix; ++imix)
		mix[imix]=mix[imix]/((float) nmix);

	/* Getpar interpolation method */
	getparstring("method",&method);

	

	/* compute filter sizes and correlation number */
	nlag  = imaxlag[0] - iminlag[0] + 1;
	ncorr = imaxcorr - imincorr + 1;
	lcorr = imaxlag[0] + 1;

	/* Compute byte sizes in wiener/spiker and autocorr */
	lagbytes = FSIZE*nlag;
	maxlagbytes = FSIZE*lcorr;
	mixbytes = maxlagbytes*nmix;

	/* Allocate memory */
	wiener	 = ealloc1float(nlag);
	spiker	 = ealloc1float(nlag);
	autocorr = ealloc1float(lcorr);
	temp = ealloc1float(lcorr);
	mixacorr = ealloc2float(lcorr,nmix);

	/* Set pointer to "cross" correlation */
	crosscorr = autocorr + iminlag[0];

	/* Zero out mixing array */
	memset((void *) mixacorr[0], 0, mixbytes);

	/* set old cdp  first trace */
	oldcdp = intrace.cdp;

	jcdp = 0;
	/* Main loop over traces */
	do {
		static int itr = 0;
		++itr;
		
		/* if neccessary, compute new filter parameters */
		if (intrace.cdp!=oldcdp && ncdp>1 && icdp<ncdp-1) {
			 int cdptotal = cdp[ncdp - 1] - intrace.cdp;

			++jcdp;

			free1float(cdpinterp);
			free1float(ominlag);
			free1float(omaxlag);

			

			 /* Compute uniformly sampled */
		         cdpinterp = ealloc1float(cdptotal);
			 ominlag = ealloc1float(cdptotal);
			 omaxlag = ealloc1float(cdptotal);
			
			 memset( (void *) omaxlag, 0, cdptotal * FSIZE);
                 	 memset( (void *) ominlag, 0, cdptotal * FSIZE);

        		 for(icdp=0; icdp<cdptotal; ++icdp)
                	 	cdpinterp[icdp] = intrace.cdp + icdp;

			 for(icdp=0; icdp<cdptotal; ++icdp)
                                ominlag[icdp] = minlag[icdp] + icdp;

			 for(icdp=0; icdp<cdptotal; ++icdp)
                                omaxlag[icdp] = maxlag[icdp] + icdp;
			

				
			 /* if linear interpolation or only one input sample */
   			if (method[0]=='l' || ncdp==1) {
                        	intlin(cdptotal,minlag,cdp,cdp[0],cdp[ncdp-1],ncdp,cdpinterp,ominlag);
				intlin(cdptotal,maxlag,cdp,cdp[0],cdp[ncdp-1],ncdp,cdpinterp,omaxlag);
				for (icdp = 0; icdp < ncdp; ++icdp) {
                        		iminlag[icdp] = NINT(ominlag[icdp]/dt);
                        		imaxlag[icdp] = NINT(omaxlag[icdp]/dt);
                		}

        		/* else, if monotonic interpolation */
        		} else if (method[0]=='m') {
                       		zind = (float (*)[4])ealloc1float(ncdp*4);
                        	cmonot(ncdp,minlag,cdp,zind);
                        	intcub(0,ncdp,minlag,zind,ncdp,cdpinterp,ominlag);
				cmonot(ncdp,maxlag,cdp,zind);
                                intcub(0,ncdp,maxlag,zind,ncdp,cdpinterp,omaxlag);			
				for (icdp = 0; icdp < ncdp; ++icdp) {
                                        iminlag[icdp] = NINT(ominlag[icdp]/dt);
                                        imaxlag[icdp] = NINT(omaxlag[icdp]/dt);
                                }

        		/* else, if Akima interpolation */
        		} else if (method[0]=='a') {
                        	zind = (float (*)[4])ealloc1float(ncdp*4);
                        	cakima(ncdp,minlag,cdp,zind);
                        	intcub(0,ncdp,minlag,zind,ncdp,cdpinterp,ominlag);
				cakima(ncdp,maxlag,cdp,zind);
                                intcub(0,ncdp,maxlag,zind,ncdp,cdpinterp,omaxlag);
				for (icdp = 0; icdp < ncdp; ++icdp) {
                                        iminlag[icdp] = NINT(ominlag[icdp]/dt);
                                        imaxlag[icdp] = NINT(omaxlag[icdp]/dt);
                                }

       			/* else, if cubic spline interpolation */
       			} else if (method[0]=='s') {
                       		zind = (float (*)[4])ealloc1float(ncdp*4);
                       		csplin(ncdp,minlag,cdp,zind);
                       		intcub(0,ncdp,minlag,zind,ncdp,cdpinterp,ominlag);
				csplin(ncdp,maxlag,cdp,zind);
                                intcub(0,ncdp,maxlag,zind,ncdp,cdpinterp,omaxlag);
				for (icdp = 0; icdp < ncdp; ++icdp) {
                                        iminlag[icdp] = NINT(ominlag[icdp]/dt);
                                        imaxlag[icdp] = NINT(omaxlag[icdp]/dt);
                                }

       			/* else, if unknown method specified */
       			} else {
               			err("%s is an unknown interpolation method!\n",method);
       			}


			/* compute filter sizes and correlation number */
			nlag  = imaxlag[jcdp] - iminlag[jcdp] + 1;
			ncorr = imaxcorr - imincorr + 1;
			lcorr = imaxlag[jcdp] + 1;

			/* Compute byte sizes in wiener/spiker and autocorr */
			/*lagbytes = FSIZE*nlag;*/
			maxlagbytes = FSIZE*lcorr;

			/* free memory */
			free1float(wiener);
			free1float(spiker);
			free1float(autocorr);
			free1float(temp);
			free2float(mixacorr);
			
			/* Allocate memory */
			wiener	 = ealloc1float(nlag);
			spiker	 = ealloc1float(nlag);
			autocorr = ealloc1float(lcorr);
			temp = ealloc1float(lcorr);
			mixacorr = ealloc2float(lcorr,nmix);

			/* Set pointer to "cross" correlation */
			crosscorr = autocorr + iminlag[jcdp];

			/* Zero out mixing array */
			memset((void *) mixacorr[0], 0, mixbytes);

		}

		/* zero out filter vectors */
		memset((void *) wiener, 0, lagbytes);
		memset((void *) spiker, 0, lagbytes);
		memset((void *) autocorr, 0, maxlagbytes);
		memset((void *) temp, 0, maxlagbytes);

		/* Form autocorrelation vector */
		xcor(ncorr, imincorr, intrace.data,
		     ncorr, imincorr, intrace.data,
		     lcorr, 0, autocorr);


		/* Leave trace alone if autocorr[0] vanishes */
		if (autocorr[0] == 0.0) {
			puttr(&intrace);
			if (showwiener)
				warn("NO Wiener filter, trace: %d", itr);

			continue;
		}


		/* Whiten */
		autocorr[0] *= 1.0 + pnoise;

		/* Read autocorr into first column of mixacorr[][] */
		memcpy( (void *) mixacorr[0], 
				(const void *) autocorr, maxlagbytes);

		/* Loop over values of the autocorrelation array */
		for (ilag = 0; ilag < lcorr; ++ilag) {

			/* Weighted moving average (mix) */
			for(imix=0; imix<nmix; ++imix)
				temp[ilag]+=mixacorr[imix][ilag]*mix[imix];

			/* put mixed data back in seismic trace */
			autocorr[ilag] = temp[ilag]; 

		}

		/* Bump columns of mixacorr[][] over by 1 */
		/* to make space for autocorr from next trace */
		for (imix=nmix-1; 0<imix; --imix)
			for (ilag=0; ilag<lcorr; ++ilag) 
				mixacorr[imix][ilag] = mixacorr[imix-1][ilag];


		/* Get inverse filter by Wiener-Levinson */
		stoepf(nlag, autocorr, crosscorr, wiener, spiker);
		

		/* Convolve pefilter with trace - don't do zero multiplies */
		for (i = 0; i < nt; ++i) {
			register int j;
			register int n = MIN(i, imaxlag[jcdp]); 
			register float sum = intrace.data[i];

			for (j = iminlag[jcdp]; j <= n; ++j)
				sum -= wiener[j-iminlag[jcdp]] * intrace.data[i-j];

			outtrace.data[i] = sum;
		}


		/* Show pefilter on request */
		if (showwiener && autocorr[0] != 0.0) {
			register int i;
			warn("Wiener filter, trace: %d", itr);


                        for (i = 0; i < imaxlag[icdp]; ++i) {
                                fprintf(outparfp, "%10g%c ", wiener[i],' ');
			}

                      	fprintf(outparfp, "\n");

			memcpy( (void *) &outtrace, 
				(const void *) &intrace, HDRBYTES);
		
			memset( (void *) outtrace.data, 0, nt*FSIZE);

			outtrace.data[0] = 1.0;

			for (i=0; i< imaxlag[icdp]; ++i)
				outtrace.data[iminlag[icdp]+i] = -wiener[i];

			puttr(&outtrace);

		} else {
			/* Output filtered trace */
			memcpy( (void *) &outtrace, 
				(const void *) &intrace, HDRBYTES);

			puttr(&outtrace);
		}

		/* update value of oldcdp */
		oldcdp = intrace.cdp;
		
	} while (gettr(&intrace));


	return(CWP_Exit());
}