synthetic.c

su 的源代码库

						w=(iw+1)*PI*2/tsec;
						refw[iw+1]=response[iw][ix][iz];
					} else {
						w=0.0;
						refw[iw+1]=cmplx(0.0,0.0);
					}

					/* compute complex frequencies */
					wpie=cmplx(w,decay);
					wsq=cmul(cmplx(0.0,1.0),wpie);
					wsq=cipow(wsq,wtype);

					/* compute reducing velocity factor */
					red_vel_factor (x, red_vel, tlag, wpie, wsq, &rvfac);

					/* compute Hanning window */
					compute_Hanning_window (win, iw, if1, if2, nw, &window);

					/* apply Hanning window and rvfac */
					refw[iw+1]=cmul(refw[iw+1],cmul(rvfac,window));
				}

				/* construct time image */
				construct_tx_trace (nw1, nt, ntc, ntfft,nfilters, fil_phase,
					tsec, unexp, sphrd, refw, fil_type, f1, f2, dbpo,
					reflectivity[ix]);
			}
		}
	}

	/* free working space */
	free1float(zr);
	free1complex(refw);
}
		
/******************************************************************************

			Subroutine to compute the reducing velocity factor

******************************************************************************/
void red_vel_factor (float x, float red_vel, float tlag, complex wpie, 
	complex wsq, complex *rvfac)
/******************************************************************************
Input:
x			range (offset) of current trace 
red_vel		reducing velocity
tlag		time lag to be applied to the seismogram
wpie
wsq			complex factor

Output
rvfac		computed reducing velocity factor
******************************************************************************/
{

	/* compute reducing velocity factor */
	if (red_vel==0.) {
		*rvfac=crmul(cmul(cmplx(0.0,1.0),wpie),tlag);
	} else {
		*rvfac=crmul(cmul(cmplx(0.0,1.0),wpie),tlag-x/red_vel);
	}

	/* multiply by complex factor */
	*rvfac=cmul(cexp1(*rvfac),wsq);
}


/******************************************************************************	

			Subroutine to construct the trace in t-x domain

******************************************************************************/
void construct_tx_trace (int nw, int nt, int ntc, int ntfft, int nfilters,
	int *filters_phase, float tsec, float unexp, float sphrd, complex *refw,
	int *fil_type, float *f1, float *f2, float *dbpo, float *reft)
/******************************************************************************
Input:
nw				number of frequencies
nt				number of time samples
ix				trace index
nfilt			number of filters to be applied
tsec
unexp
sphrd			flattening earth correction factor
refw			array[nw+1] of complex trace

Output:
reft			array[nt] current trace in t-x domain
******************************************************************************/
{
	int it;			/* loop counter */
	float a,b,c;		/* auxiliary variables */
	
	/* if requested, apply filters */
	if (nfilters !=0)  {
		apply_filters (filters_phase, nfilters, nw, tsec, f1, f2, dbpo, 
			fil_type, refw);
	}

	/* inverse fft to transfer data to t-x domain */
	pfacc (-1, ntfft, refw);
	

	/* construct the t-x image */
	for (it=0; it<nt; it++) {
		a=(float)it/ntc;
		b=unexp*a;
		c=exp(b);
		reft[it]=refw[it].r*c*sphrd;
	}

	/* reinitialize working arrays */
	for (it=0; it<ntfft; it++) {
		refw[it]=cmplx(0.0,0.0);
	}
}
					
					
/******************************************************************************

				Subroutine to compute a Hanning window

******************************************************************************/
void compute_Hanning_window (int iwin, int iw, int if1, int if2, int nf, 
	complex *win)
/******************************************************************************
Input:
iwin		=1 for Hanning window
if1			lower f=window frequency
if2			higher window frequency
nw			number of frequencies in output data
iw			current frequency

Output:
win			computed Hanning window for given frequency
******************************************************************************/
{

	/* compute Hanning window */
	if (iwin==1) {
		if (iw<if1-1) {
			*win=crmul(cmplx(1.0,0.0),0.5*(1.+cos(PI*(if1-iw+1)/if1)));
		} else if((iw>=if1-1)&&(iw<=if2-1)) {
			*win=cmplx(1.0,0.0);
		} else if((iw>if2-1)&&(iw<nf-1)) {
			*win=crmul(cmplx(1.0,0.0),0.5*(1.+cos(PI*(iw-if2+1)/(nf-if2))));
		} else {
			*win=cmplx(0.0,0.0);
		}
	} else {
		if (iw<nf-1) {
			*win=crmul(cmplx(1.0,0.0),0.5*(1.+cos(PI*(iw+1)/nf)));
		} else {
			*win=cmplx(0.0,0.0);
		}
	}
}
			
/******************************************************************************

		Subroutine to apply hi-cut, lo-cut or notch zero or minimum
							phase filters

******************************************************************************/
void apply_filters (int *min_phase, int nfilters, int nw, float tsec, float *f1,
	float *f2, float *dbpo, int *filter_type, complex *refw)
/******************************************************************************
Input:
min_phase		=1 for minimum phase filters
				=0 for zero phase filters
nfilters		number of filters to apply
nw				number of frequencies
tsec
f1				array[nfilters] of low frequencies
f2				array[nfilters] of high frequencies
dbpo				array[nfilters] of filter slopes in db/oct
filtype			array[nfilters] of filter types:
				=1 low cut filter
				=2 high cut filter
				=3 notch filter
refw			array[nw+1] of complex samples to filter

Output:
refw			array[nw+1] of complex filtered samples 
*******************************************************************************
Note:
refw contains only the positive frequencies, the negatives are simply their
complex conjugates and can be computed when necessary
******************************************************************************/
{
	int i,j,iw;			/* loop counters */
	int nw1=nw+1;		/* number of frequencies (positive+zero) */
	int nwt=2*nw;		/* total number of frequencies */
	int index1,index2;	/* auxiliary indices for notch filter */
	float delfq=0.0;
	float poles;		/* filter poles */
	float r,factor;		/* auxiliary variables */
	float maxabs;
	complex *filter;	/* scratch array for filter coefficients */
	complex *filter1;	/* scratch array for filter coefficients */


	/* defensive programming */
	if (nfilters==0) return;

	/* allocate working space */
	filter=alloc1complex(2*nwt);
	filter1=alloc1complex(2*nwt);

	/* initialize filter array */
	for (iw=0; iw<nw1; iw++) filter[iw]=cmplx(1.0,0.0);

	/* main loop over number of requested filters */
	for (i=0; i<nfilters; i++) {
		
		/* smearing range */
		if (min_phase[i]==1) delfq=2*PI/tsec;

		if (filter_type[i]==1) { 	/* if low-pass Butterworth */
			poles=dbpo[i]/6.+1;
			r=1./((f1[i]-delfq)*tsec);

			/* compute filter coefficients */
			for (iw=0; iw<nw1; iw++) {
				factor=1.0/sqrt(1.+pow(iw*r,poles));
				filter[iw]=crmul(filter[iw],factor);
			}

		} else if (filter_type[i]==2) {	/* hi-pass Butterworth */
			poles=dbpo[i]/6.+1;
			r=1./((f1[i]+delfq)*tsec);

			/* compute filter coefficients */
			for (iw=0; iw<nw1; iw++) {
				factor=pow(iw*r,poles);	
				filter[iw]=crmul(filter[iw],sqrt(factor/(1+factor)));
			}

		} else if (filter_type[i]==3) { /* notch filter */
			index1=(f1[i]-delfq)*tsec+1;
			index2=(f2[i]+delfq)*tsec+2;

			if ((index2-index1)<16) {
				warn("in notch filter %d f1 too close to f2, "
				"this filter skiped\n",i);
				continue;
			}

			/* db reduction Neils design */
			r=1.0-pow(10.0,-dbpo[i]/20.);
			for (j=index1; j<=index2; j++) {
				factor=cos(2*PI*(j-index1)/(index2-index1))-1;
				filter[j-1]=crmul(filter[j-1],(1+r*factor)/2.);
			}
		} else {
				warn("unknown filter type for filter number %d,"
				" this filter skiped\n",i);
		}

		/* apply minimum phase correction */
		if (min_phase[i]==1) {
			for (iw=0; iw<nw1; iw++) {
				if (filter[iw].r==0.0) {
					filter[iw]=cmplx(-30.0,0.0);
				} else {

					/* set A[w]=ln(f[w]) */
					filter[iw]=cmplx(log(filter[iw].r),0.0);
				}
			}


			/* take care of negative frequencies needed for FFT */
			for (iw=nw1; iw<nwt; iw++) {
				filter[iw]=conjg(filter[nwt-iw+1]);
			}

			/* take inverse FFT to go to time domain */
			pfacc (-1, nwt, filter);

			for (iw=1; iw<nw1; iw++) {
				r=(1.+cos(PI*iw/nw))/2.;	
				filter[iw]=crmul(filter[iw],r);
				filter[nwt-iw+1]=crmul(filter[nwt-iw+1],r);
			}

			/* construct filter coefficients for all frequencies */
			for (iw=0; iw<nw1; iw++) {
				filter1[iw]=filter[iw];
			}
			for (iw=nw1; iw<nwt; iw++) {
				filter1[iw]=cneg(conjg(filter[iw]));
			}

			/* take forward FFT to get back to frequency domain */
			pfacc (1, nwt, filter);
			pfacc (1, nwt, filter1);
	
			maxabs=0.0;
			for (iw=0; iw<nw1; iw++) {
				filter[iw]=cexp1(cadd(filter[iw],filter1[iw]));
				maxabs=MAX(rcabs(filter[iw]),maxabs);
			}
	
			/* normalize */
			maxabs=1.0/maxabs;
			for (iw=0; iw<nw1; iw++) {
				filter[iw]=crmul(filter[iw],maxabs);
			}
		}
	}

	/* apply the filter */
	for (iw=0; iw<nw1; iw++) {
		refw[iw]=cmul(refw[iw],filter[iw]);
	}
	for (iw=nw1; iw<nwt; iw++) {
		refw[iw]=conjg(refw[nwt-iw]);
	}

	/* free allocated space */
	free1complex(filter);
	free1complex(filter1);
}