sufxdecon.c

su 的源代码库

	/* loop over time windows */
	for (iw=0;iw<ntw;iw++) {

	        if (iw>0 && iw<ntw-1) nspw=nspwi; 
		else if (iw==0) 
		        if (ntw>1) nspw = nspwf;
			else        nspw = nsp;
		else
		        nspw = nsp - nspws*iw + nstaper/2;
 
		if (verbose)
		        warn("iw=%i, ntw=%i, nspw=%i, twlen=%f",
					iw,ntw,nspw,twlen); 

		/* zero fdata */
		memset((void *) fdata[0], 0, nf*ntrc*sizeof(complex));
     
		/* select data */
		for (jw=0;jw<ntrc;jw++) {
 
		        if (iw>0)
			   for (it=0;it<nspw;it++)
			     tidataw[it]=tdatain[jw][it + iw*nspws - nstaper/2];
			else
			        for (it=0;it<nspw;it++)  
				         tidataw[it]=tdatain[jw][it];	  

			memset((void *) (tidataw + nspw), 0, (nfft-nspw)*FSIZE);
			memset((void *) ct, 0, nfft*sizeof(complex));

			/* FFT from t to f */
			for (it=0;it<nfft;it++)
			        ttidataw[it]=(it%2 ? -tidataw[it] : tidataw[it]);
			pfarc(1, nfft, ttidataw, ct);

			/* Store values */    
			for (ifq = 0; ifq < nf; ifq++) { 
			        fdata[jw][ifq] = ct[nf-1-ifq];
			}
		}

		/* Loop over space windows */
		for (jx=0;jx<nwx;jx++){

		        /* to take care of a possible incomplete last window */
		        if (ntrc<jx*ntrw+2*ntrw) 
			        ntrwu = ntrc - jx*ntrw;
			else
			        ntrwu = ntrw;

			if (verbose) {
			        warn("jx=%i, ntrc=%i, ntrw=%i",
					jx,ntrc,ntrw);
				warn("ntrwu=%i, nwx=%i, ntrf=%i",
					ntrwu,nwx,ntrf);
			}

			/* zero fdataw */
			for (jw=0;jw<ntrw+2*ntrf;jw++)
			        memset((void *) fdataw[jw], 0, nf*sizeof(complex));

			/* select data */
			for (jw=0;jw<ntrwu+2*ntrf;jw++) 
			        for (ifq = 0; ifq < nf; ifq++) {

				   if (jx>0 && jx<nwx-1)  
				       fdataw[jw][ifq] = fdata[jw + jx*ntrw - ntrf][ifq];
					else if (jx==0)
					        if (jw>=ntrf && jw<ntrw+ntrf) 
						        fdataw[jw][ifq] = fdata[jw - ntrf][ifq];
						else if (jw<ntrf) 
						        fdataw[jw][ifq] = fdata[0][ifq];
						else 
						        if (nwx>1) 
							        fdataw[jw][ifq] = fdata[jw - ntrf][ifq];
							else 
							        fdataw[jw][ifq] = fdata[ntrc-1][ifq];
					else
					        if (jw<ntrwu+ntrf)
						        fdataw[jw][ifq] = fdata[jw + jx*ntrw - ntrf][ifq];
						else 
						        fdataw[jw][ifq] = fdata[ntrc-1][ifq];
	}

			/* loop over frequencies */
			for (ifq=0;ifq<nf;ifq++) {

			        if ((float)ifq*d1>=fmin && (float)ifq*d1<=fmax) {

				       /* Loop over space window */
				       memset((void *) sfreq, 0, (ntrwu+2*ntrf)*sizeof(complex));
				       memset((void *) sfreqcj, 0, (ntrwu+2*ntrf)*sizeof(complex));

				       for (jw=0;jw<ntrwu+2*ntrf;jw++) {
	  
					       sfreq[jw]=fdataw[jw][ifq];
					       sfreqcj[jw]=conjg(fdataw[jw][ifq]);
				       }

				       memset((void *) autocorr, 0, (ntrf+1)*sizeof(complex));

				       /* complex autocorrelation */
				       cxcor(ntrwu,0,sfreq,ntrwu,0,sfreq,ntrf+1,0,autocorr);

				       /* zeroing files */
				       memset((void *) rautoc, 0, (ntrf+1)*FSIZE);
				       memset((void *) iautoc, 0, (ntrf+1)*FSIZE);

				       /* taking real and imaginary parts */
				       for (jw=0;jw<ntrf+1;jw++) {
					       rautoc[jw]=autocorr[jw].r;
					       iautoc[jw]=autocorr[jw].i;
				       }

				       /* zeroing files */
				       memset((void *) gvector, 0, 2*ntrf*FSIZE);
				       memset((void *) fvector, 0, (2*ntrf+1)*sizeof(complex));
				       for (ir=0;ir<2*ntrf;ir++) 
					       memset((void *) rmatrix[ir], 0, 2*ntrf*FSIZE);

				       /* matrix problem */
				       for (ir=0;ir<ntrf;ir++) 
					       for (jr=0;jr<ntrf;jr++) { 
						       if (ir>=jr) rmatrix[ir][jr]=autocorr[ir-jr].r;
						       else        rmatrix[ir][jr]=autocorr[jr-ir].r;
					       }

				       for (ir=ntrf;ir<2*ntrf;ir++)
					       for (jr=0;jr<ntrf;jr++) {
						       if (ir-ntrf<jr) rmatrix[ir][jr]=-autocorr[jr-ir+ntrf].i;
						       else            rmatrix[ir][jr]= autocorr[ir-jr-ntrf].i;
					       }

				       for (ir=ntrf;ir<2*ntrf;ir++)
					       for (jr=ntrf;jr<2*ntrf;jr++)
						       rmatrix[ir][jr]=rmatrix[ir-ntrf][jr-ntrf];

				       for (ir=0;ir<ntrf;ir++)
					       for (jr=ntrf;jr<2*ntrf;jr++)
						       rmatrix[ir][jr]=-rmatrix[ir+ntrf][jr-ntrf];

				       for (ig=0;ig<2*ntrf;ig++) {
					       if (ig<ntrf) gvector[ig]=autocorr[ig+1].r;
				  else gvector[ig]=autocorr[ig-ntrf+1].i;
				       }

				       LU_decomposition(2*ntrf,rmatrix,ipvt,info);
				       backward_substitution(2*ntrf,rmatrix,ipvt,gvector);
      
				       /* construct filter */
				       for (ifv=0,ig=ntrf-1;ifv<ntrf;ifv++,ig--) 
					       fvector[ifv]=conjg(cmplx(gvector[ig]/2.,gvector[ig+ntrf]/2.));

				       for (ifv=ntrf+1,ig=0;ifv<2*ntrf+1;ifv++,ig++) 
					       fvector[ifv]=cmplx(gvector[ig]/2.,gvector[ig+ntrf]/2.);
	 
				       memset((void *) sfreqout, 0, ntrwu*sizeof(complex));

				       /* convolution of data with filter */
				       /* output is one sample ahead */
				       cconv(ntrwu+2*ntrf,-ntrf,sfreq,2*ntrf+1,-ntrf,fvector,ntrwu,0,sfreqout); 

				       /* store filtered values */
				       for (jw=0;jw<ntrwu;jw++) ffreq[jw][ifq]=sfreqout[jw];

				}
			} /* end of frequencies loop */

			/* loop along space windows */
			for (jw=0;jw<ntrwu;jw++) {
    
			        /* select data */
			        for (ifq=0,itt=nf-1;ifq<nf;ifq++,itt--)
				        freqv[ifq] = ffreq[jw][itt]; 

				memset((void *) (freqv+nf), 0, (nfft-nf)*sizeof(complex));
				memset((void *) todataw, 0, nfft*FSIZE);

				/* FFT back from f to t and scaling */
				pfacr(-1, nfft, freqv, todataw);
				for (it=0;it<MIN(nsp,nfft);it++)
					todataw[it]/=nfft; 
				for (it=0;it<MIN(nsp,nfft);it++)
				        ttodataw[it]=(it%2 ? -todataw[it] : todataw[it]); 
      
				/*loop along time */
				if (ntw>1) {
				        /* first portion of time window */
				        if (iw>0) 
					        for (it=0;it<nstaper;it++)
						        tdataout[jx*ntrw+jw][it+iw*nspws-nstaper/2]+=
							        ttodataw[it]*((float)(it)*dt/taper);
					else 
					        for (it=0;it<nstaper;it++)
						        tdataout[jx*ntrw+jw][it]=ttodataw[it];

					/* intermediate portion of time window */
					if (iw>0) 
					        for (it=nstaper;it<nspw-nstaper;it++)
						        tdataout[jx*ntrw+jw][it+iw*nspws-nstaper/2]=ttodataw[it];
					else 
					        for (it=nstaper;it<nspw-nstaper;it++)
						        tdataout[jx*ntrw+jw][it]=ttodataw[it];

					/* last portion of time window */
					if (iw>0 && iw<ntw-1) 
					        for (it=nspw-nstaper;it<nspw;it++)
						        tdataout[jx*ntrw+jw][it+iw*nspws-nstaper/2]+=
							        ttodataw[it]*(1.-((float)(it-nspw+nstaper))*dt/taper);
					else if (iw==ntw-1)
					        for (it=nspw-nstaper;it<nspw;it++)
						        tdataout[jx*ntrw+jw][it+iw*nspws-nstaper/2]=ttodataw[it];
					else 
					        for (it=nspw-nstaper;it<nspw;it++)
						        tdataout[jx*ntrw+jw][it]+=ttodataw[it]*(1.-((float)(it-nspw+nstaper))*dt/taper);
				}
				else {
				        for (it=0;it<nsp;it++) 
					        tdataout[jx*ntrw+jw][it]=ttodataw[it];
				}

			} /* end loop over space windows */

		} /* end loop over space windows */

	} /* end of time windows loop */

	for (jx=0;jx<ntrc;jx++) {
	        efread(&tr,HDRBYTES,1,headerfp);
		memcpy( (void *) tr.data, (const void *) tdataout[jx], nsp*FSIZE);
		puttr(&tr);
	}

	/* Clean up */
	free1float(tidataw);
	free1float(ttidataw);
	free2complex(fdataw);
	free1complex(sfreq);
	free1complex(sfreqcj);
	free1complex(autocorr);
	free1float(rautoc);
	free1float(iautoc);
	free2float(rmatrix);
	free1float(gvector);
	free1complex(fvector);
	free1complex(sfreqout);
	free1complex(freqv);
	free1float(todataw);
	free1float(ttodataw);
	free2complex(ffreq);
	free2float(tdatain);
 
	efclose(headerfp);
	if (istmpdir) eremove(headerfile);
	if (istmpdir) eremove(tracefile);

	return(CWP_Exit());

}

/* complex correlation */
void cxcor (int lx, int ifx, complex *x,
            int ly, int ify, complex *y, 
            int lz, int ifz, complex *z)
{
        int i,j;
	complex *xr;
	xr = alloc1complex(lx);
	for (i=0,j=lx-1; i<lx; ++i,--j)
	        xr[i] = conjg (x[j]);
	cconv(lx,1-ifx-lx,xr,ly,ify,y,lz,ifz,z);
	free1complex(xr);
}

/* complex convolution */
void cconv (int lx, int ifx, complex *x,
       	    int ly, int ify, complex *y,
	    int lz, int ifz, complex *z)
{
        int ilx=ifx+lx-1,ily=ify+ly-1,ilz=ifz+lz-1,i,j,jlow,jhigh;
	complex sum;
  
	x -= ifx;  y -= ify;  z -= ifz; 
	for (i=ifz; i<=ilz; ++i) {
	        jlow = i-ily;  if (jlow<ifx) jlow = ifx;
		jhigh = i-ify;  if (jhigh>ilx) jhigh = ilx;
		for (j=jlow,sum=cmplx(0.,0.); j<=jhigh; ++j){
		        sum = cadd(sum,cmul(x[j],y[i-j]));
		}
		z[i] = sum;
	}
}

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