kxwfilter.c

#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <genfft.h>

#define ISODD(n) ((n) & 01)
#define MAX(x,y) ((x) > (y) ? (x) : (y))
#define MIN(x,y) ((x) < (y) ? (x) : (y))

void kxwfilt(complex *data, float k, float dx, int nkx, float a1, float perc);

void kxwfilter(float *data, int nt, int nx, float dt, float dx, float fmin, float fmax, float angle, float cp, float perc)
{
	int ntfft, nfreq, nkx, ix, it, iomin, iomax, xorig, iom, ikx;
	float df, dkx, deltom, omin, omax, kp, om;
	float *pdata;
	complex *cdata;

    ntfft = optncr(nt);
    nfreq = ntfft/2+1;
    nkx = optncc(2*nx);

	df     = 1.0/((float)ntfft*dt);
	dkx    = 2.0*M_PI/(nkx*dx);
	deltom = 2.*M_PI*df;
	omin   = 2.*M_PI*fmin;
	omax   = 2.*M_PI*fmax;

	iomin  = (int)MIN((omin/deltom), (nfreq-1));
	iomin  = MAX(iomin, 1);
	iomax  = MIN((int)(omax/deltom), (nfreq-1));

    pdata = (float *)calloc(ntfft*nkx,sizeof(float));
    cdata = (complex *)malloc(nfreq*nkx*sizeof(complex));

    /* copy input data into extended arrays with padded zeroes */
    for (ix=0; ix<nx; ix++) {
        memcpy(&pdata[ix*ntfft],&data[ix*nt],nt*sizeof(float));
    }

    /* transform from t-x to kx-w */
    xorig = nkx/2;
    xt2wkx(pdata, cdata, ntfft, nkx, ntfft, nkx, xorig);


    for (iom = iomin; iom <= iomax; iom++) {
		om  = iom*deltom;
		kp  = om/cp;

		kxwfilt(&cdata[iom*nkx], kp, dx, nkx, angle, perc);
	}

    /* transform back to t-x */
    wkx2xt(cdata, pdata, ntfft, nkx, nkx, ntfft, xorig);

    /* reduce array to nt samples nx traces */
    for (ix=0; ix<nx; ix++) {
        memcpy(&data[ix*nt],&pdata[ix*ntfft],nt*sizeof(float));
    }

	free(pdata);
	free(cdata);

	return;
}


void kxwfilt(complex *data, float k, float dx, int nkx, float a1, float perc)
{
	int 	ikx, ik1, ik2, ntap;
	float 	kxnyq, dkx, kxfmax, kfilt;
	float 	kpos, band, li, *filter;

	kxnyq = M_PI/dx;
	dkx   = 2.0*M_PI/(nkx*dx);
	if (a1 > 90.0) kpos = kxnyq;
	else kpos = k*sin(M_PI*a1/180.0);

	filter = (float *)malloc(nkx*sizeof(float));

	band  = fabs(2*kpos);
	ntap  = (int)abs((int)(perc*band/dkx));
	kfilt = fabs(dkx*ntap);

	if (perc > 0) {
		if (kpos+kfilt < kxnyq) {
			kxfmax = kpos+kfilt;
		}
		else {
			kxfmax = kxnyq;
			ntap = (int)(0.15*nkx/2);
		}
	}
	else {
		kxfmax = MIN(kpos, kxnyq);
	}

	ik1 = (int)(kxfmax/dkx);
	ik2 = ik1 - ntap;

	if (perc < -0.5 || perc > 1.0) {
		if (kpos > 0.85*kxnyq) {
			kpos = 0.85*kxnyq;
		}
		ik1 = nkx/2-1;
		ik2 = (int)(kpos/dkx);
	}

	li = 1.0/(ik1-ik2);
	for (ikx = 0; ikx < ik2; ikx++) 
        filter[ikx] = 1.0;
	for (ikx = ik2; ikx < ik1; ikx++)
		filter[ikx] = 0.5*(cos(M_PI*(ikx-ik2)*li)+1);
	for (ikx = ik1; ikx <= nkx/2; ikx++)
		filter[ikx] = 0.0;

	for (ikx = 0; ikx <= (nkx/2); ikx++) {
		data[ikx].r *= filter[ikx];
		data[ikx].i *= filter[ikx];
	}
	for (ikx = (nkx/2+1); ikx < nkx; ikx++) {
		data[ikx].r *= filter[nkx-ikx];
		data[ikx].i *= filter[nkx-ikx];
	}

	free(filter);
	return;
}