head	1.18;
access;
symbols;
locks; strict;
comment	@ * @;


1.18
date	2011.11.12.00.09.00;	author john;	state Exp;
branches;
next	1.17;

1.17
date	2006.11.07.22.58.42;	author john;	state Exp;
branches;
next	1.16;

1.16
date	2004.12.23.23.32.36;	author john;	state Exp;
branches;
next	1.15;

1.15
date	2003.06.09.16.17.07;	author john;	state Exp;
branches;
next	1.14;

1.14
date	2002.03.13.18.53.53;	author john;	state Exp;
branches;
next	1.13;

1.13
date	97.07.28.22.36.46;	author john;	state Exp;
branches;
next	1.12;

1.12
date	96.09.13.21.49.05;	author john;	state Exp;
branches;
next	1.11;

1.11
date	96.02.07.15.56.04;	author jkc;	state Exp;
branches;
next	1.10;

1.10
date	96.01.29.17.07.19;	author jkc;	state Exp;
branches;
next	1.9;

1.9
date	95.08.04.15.51.03;	author jkc;	state Exp;
branches;
next	1.8;

1.8
date	94.08.25.15.01.31;	author jkc;	state Exp;
branches;
next	1.7;

1.7
date	94.03.30.13.53.46;	author john;	state Exp;
branches;
next	1.6;

1.6
date	93.12.14.13.41.45;	author jkc;	state Exp;
branches;
next	1.5;

1.5
date	93.11.29.14.23.35;	author john;	state Exp;
branches;
next	1.4;

1.4
date	93.11.24.16.26.00;	author john;	state Exp;
branches;
next	1.3;

1.3
date	92.10.23.15.11.37;	author john;	state Exp;
branches;
next	1.2;

1.2
date	92.03.24.12.54.29;	author jkc;	state Exp;
branches;
next	1.1;

1.1
date	91.03.03.08.56.26;	author jkc;	state Exp;
branches;
next	;


desc
@Time variant band-pass filtering
@


1.18
log
@copyright
@
text
@/* Copyright (c) Colorado School of Mines, 2011.*/
/* All rights reserved.                       */

/* SUTVBAND: $Revision: 1.17 $ ; $Date: 2006/11/07 22:58:42 $      */

#include "su.h"
#include "segy.h"

/*********************** self documentation **********************/
char *sdoc[] = {
"                                                               ",
" SUTVBAND - time-variant bandpass filter (sine-squared taper)  ",
"                                                               ",
" sutvband <stdin >stdout tf= f=			        ",
"                                                               ",
" Required parameters:                                          ",
"       dt = (from header)      time sampling interval (sec)    ",
"       tf=             times for which f-vector is specified   ",
"       f=f1,f2,f3,f4   Corner frequencies corresponding to the ",
"                       times in tf. Specify as many f= as      ",
"                       there are entries in tf.                ",
"                                                               ",
" The filters are applied in frequency domain.                  ",
"                                                               ",
" Example:                                                      ",
" sutvband <data tf=.2,1.5 f=10,12.5,40,50 f=10,12.5,30,40 | ...",
NULL};

/* Credits:
 *      CWP: Jack, Ken
 *
 * Trace header fields accessed:  ns, dt, delrt
 */
/**************** end self doc ***********************************/


#define PIBY2   1.57079632679490
#define LOOKFAC 2       /* Look ahead factor for npfaro   */
#define PFA_MAX 720720  /* Largest allowed nfft           */

/* Prototypes */
void makefilter(float *f, int nfft, int nfreq, float dt, float *filter);
void bandpass(float *data, int ntime, int nfft, int nfreq, 
		float *filterj, float *ftracej);

segy tr;

int
main(int argc, char **argv)
{
        float **filter;         /* filter arrays                        */
        float *tf=NULL;		/* times at which filters are centered  */
        int *itf;               /* ... as integers                      */
	int jmin;		/* index of first filter itf value	*/
	int jmax;		/* index of last filter itf value	*/
        int nfft;     	        /* fft sizes in each time gate          */
        int nfreq;     	        /* number of frequencies  	        */
        float **ftrace;         /* filtered sub-traces                  */
        int nfilter;            /* number of filters specified          */
      	float dt;               /* sample spacing                       */
        float tmin;             /* first time on traces                 */
        int nt;                 /* number of points on input trace      */
	cwp_Bool seismic;	/* is this seismic data?		*/
	float *data;

        
        /* Initialize */
        initargs(argc, argv);
        requestdoc(1);


        /* Get info from first trace */ 
        if (!gettr(&tr))  err("can't get first trace");
	seismic = ISSEISMIC(tr.trid); 
		 
	if (!seismic)
		warn("input is not seismic data, trid=%d", tr.trid);
        nt = tr.ns;
        if (!getparfloat("dt", &dt))    dt = ((double) tr.dt)/1000000.0;
        if (!dt) err("dt field is zero and not getparred");
	tmin = tr.delrt/1000.0;


        /* Get number of filters and center times */
        if (!(nfilter = countparval("tf")))  MUSTGETPARFLOAT("tf", tf);
	if (countparname("f") != nfilter)
		err("must give one f 4-tuple for each"
		    " (%d) tf value", nfilter);
		
	/* Leave room for possibly missing filters at endpoints */
	tf = ealloc1float(nfilter+4);  /* never use first 2 or last 2 */
	itf = ealloc1int(nfilter+4);
        getparfloat("tf", tf+2); jmin = 2; jmax = nfilter + 1;
	{ register int j;
          for (j = jmin; j <= jmax; ++j)  itf[j] = NINT((tf[j] - tmin)/dt);
        }
	

        /* Make filters with scale for inverse transform */
	nfft = npfaro(nt, LOOKFAC * nt);
	if (nfft >= SU_NFLTS || nfft >= PFA_MAX)
		err("Padded nt=%d -- too big", nfft);
	nfreq = nfft/2 + 1;
	filter = ealloc2float(nfreq, nfilter+4); /* never use 1st & last */
        { register int j;
          for (j = jmin; j <= jmax; ++j) {
                float *f = ealloc1float(4);

                if (getnparfloat(j-jmin+1, "f", f) != 4)
                        err("must give 4 corner frequencies in f=");
        	if (f[0] < 0.0 || f[0] > f[1] ||
				  f[1] >= f[2] || f[2] > f[3])
               		err("Filter #%d has bad frequencies", j - jmin + 1);
                makefilter(f, nfft, nfreq, dt, filter[j]);
          }
        }
	

	/* User may not have given a filter for tmin and/or tmax--	*/
	/* Extend array so can always assume these filters are present.	*/
	/* Note don't really use any of the extra storage in **filter!	*/
	if (itf[jmin] > 0) {
		filter[jmin-1] = filter[jmin]; 
		itf[jmin-1] = 0;
		--jmin;
	}
	if (itf[jmax] < nt - 1) {
		filter[jmax+1] = filter[jmax];
		itf[jmax+1] = nt - 1;
		++jmax;
	}
	
	
	/* Extend array so can always consider time points to be interior */
	itf[jmin-1] = 0;      /* now jmin - 1 is a valid index */
	itf[jmax+1] = nt - 1; /* now jmax + 1 is a valid index */


        /* Main loop over traces */
        ftrace = ealloc2float(nt, nfilter+4); /* never use 1st & last */
	data = ealloc1float(nt);
        do {
                register int i, j;
		
		/* Construct filtered sub-traces */
		for (j = jmin; j <= jmax; ++j) {			
			memset( (void *) data, 0, nt*FSIZE);
			for (i = itf[j-1]; i <= itf[j+1]; ++i)
				data[i] = tr.data[i];
                        bandpass(data,nt,nfft,nfreq,filter[j],ftrace[j]);
                }
		               
               /* Compose filtered trace from sub-traces */
               for (j = jmin; j < jmax; ++j) {
	       		float fitfj;
                        for (fitfj = i = itf[j]; i <= itf[j+1]; ++i) {
                                float a = (i - fitfj)/(itf[j+1] - fitfj);
                                tr.data[i] = (1-a)*ftrace[j][i] +
                                                 a*ftrace[j+1][i];
			}
                }
                
                puttr(&tr);
        } while (gettr(&tr));

        return(CWP_Exit());
}

void bandpass(float *data, int nt, int nfft, int nfreq,
		float *filterj, float *ftracej)
{
	float *rt = ealloc1float(nfft);
	complex *ct = ealloc1complex(nfreq);
	register int i;

        /* Load trace into rt (zero-padded) */
        memcpy( (void *) rt, (const void *) data, nt*FSIZE);
        memset( (void *) (rt + nt), 0, (nfft-nt)*FSIZE);

        /* FFT, filter, inverse FFT */
        pfarc(1, nfft, rt, ct);
        for (i = 0; i < nfreq; ++i)  ct[i] = crmul(ct[i], filterj[i]);
        pfacr(-1, nfft, ct, rt);

        /* Load traces back in, recall filter had nfft factor */
        for (i = 0; i < nt; ++i)  ftracej[i] = rt[i]; /* ftracej = rt ?? */

	/* free temp arrays */
	free1float( rt );
	free1complex( ct );
}


void makefilter(float *f, int nfft, int nfreq, float dt, float *filter)
{
        float onfft = 1.0 / nfft;
        float df = onfft / dt;
        int nfreqm1 = nfreq - 1;
        int if1 = NINT(f[0]/df);
        int if2 = NINT(f[1]/df);
        int if3 = MIN(NINT(f[2]/df), nfreqm1);
        int if4 = MIN(NINT(f[3]/df), nfreqm1);


        { register int i;
	  register float c = PIBY2 / (if2 - if1 + 2);
	  for (i = if1; i <= if2; ++i) {
		register float s = sin(c*(i - if1 + 1));
		filter[i] = s * s * onfft;
          }
	 }

        { register int i;
	  register float c = PIBY2 / (if4 - if3 + 2);
	  for (i = if3; i <= if4; ++i) {
		register float s = sin(c*(if4 - i + 1));
		filter[i] = s * s * onfft;
	  }
        }

        { register int i;
          for (i = if2 + 1; i < if3;   ++i)  filter[i] = onfft; 
          for (i = 0;       i < if1;   ++i)  filter[i] = 0.0; 
          for (i = if4 + 1; i < nfreq; ++i)  filter[i] = 0.0; 
        }
}
@


1.17
log
@changed    memset( .., (int) '\0',...) to memset( .., 0,..) as per suggestion of Marion Zhu and Reg Beardsley
@
text
@d1 4
a4 1
/* SUTVBAND: $Revision: 1.16 $ ; $Date: 2004/12/23 23:32:36 $      */
@


1.16
log
@changed err to warn for data nonseismic
@
text
@d1 1
a1 1
/* SUTVBAND: $Revision: 1.15 $ ; $Date: 2003/06/09 16:17:07 $      */
d144 1
a144 1
			memset( (void *) data, (int) '\0', nt*FSIZE);
d175 1
a175 1
        memset( (void *) (rt + nt), (int) '\0', (nfft-nt)*FSIZE);
@


1.15
log
@9 June 2003
@
text
@d1 1
a1 1
/* SUTVBAND: $Revision: 1.14 $ ; $Date: 2002/03/13 18:53:53 $      */
d74 1
a74 1
		err("input is not seismic data, trid=%d", tr.trid);
@


1.14
log
@added Song Yeung Suh's changes
@
text
@d1 1
a1 1
/* SUTVBAND: $Revision: 1.13 $ ; $Date: 1997/07/28 22:36:46 $      */
d163 1
a163 1
        return EXIT_SUCCESS;
@


1.13
log
@changes suggested by Stew Levin to silence SGI compiler warnings, etc
@
text
@d1 1
a1 1
/* SUTVBAND: $Revision: 1.12 $ ; $Date: 1996/09/13 21:49:05 $      */
d61 1
d138 1
a139 1
		float *data = ealloc1float(nt);
d169 1
a169 1
	float *rt = ealloc1float(nt);
d184 4
@


1.12
log
@silenced -Wall messages
@
text
@d1 1
a1 1
/* SUTVBAND: $Revision: 1.11 $ ; $Date: 1996/02/07 15:56:04 $      */
d97 1
a97 1
	if (nfft >= MIN(SU_NFLTS, PFA_MAX))
@


1.11
log
@Deal with header fields accessed/modified down to sudoc.
@
text
@d1 1
a1 1
/* SUTVBAND: $Revision: 1.10 $ ; $Date: 1996/01/29 17:07:19 $      */
d45 1
d49 1
a49 1
        float *tf;              /* times at which filters are centered  */
@


1.10
log
@Say sampling INTERVAL instead of rate to make Ken happy.
@
text
@d1 1
a1 1
/* SUTVBAND: $Revision: 1.9 $ ; $Date: 1995/08/04 15:51:03 $      */
a23 2
"                                                               ",
" Trace header fields accessed:  ns, dt, delrt.			",
d28 2
d86 1
a86 1
	tf = ealloc1float(nfilter+4);  /* never use ist2 or last 2 */
@


1.9
log
@Changed cast of tr.dt to double for consistency.
@
text
@d1 1
a1 1
/* SUTVBAND: $Revision: 1.8 $ ; $Date: 1994/08/25 15:01:31 $      */
d14 1
a14 1
"       dt = (from header)      time sampling rate (sec)        ",
@


1.8
log
@Removed copyright statement from codes that had it.
Moved self-doc line for sudoc.
@
text
@d1 1
a1 1
/* SUTVBAND: $Revision: 1.7 $ ; $Date: 94/03/30 13:53:46 $      */
d74 1
a74 1
        if (!getparfloat("dt", &dt))    dt = (float) tr.dt/1000000.0;
@


1.7
log
@added ISSEISMIC macro
@
text
@d1 1
a1 1
/* SUTVBAND: $Revision: 1.6 $ ; $Date: 93/12/14 13:41:45 $      */
a2 13
/*----------------------------------------------------------------------
 * Copyright (c) Colorado School of Mines, 1990.
 * All rights reserved.
 *
 * This code is part of SU.  SU stands for Seismic Unix, a processing line
 * developed at the Colorado School of Mines, partially based on Stanford
 * Exploration Project (SEP) software.  Inquiries should be addressed to:
 *
 *  Jack K. Cohen, Center for Wave Phenomena, Colorado School of Mines,
 *  Golden, CO 80401  (jkc@@dix.mines.colorado.edu)
 *----------------------------------------------------------------------
 */

a26 1
/**************** end self doc ***********************************/
d31 1
@


1.6
log
@Enlarged notion of seismic trace to include DEAD, etc.
@
text
@d1 1
a1 1
/* SUTVBAND: $Revision: 1.5 $ ; $Date: 93/11/29 14:23:35 $      */
d82 2
a83 2
	seismic = (tr.trid == 0 || tr.trid == TREAL ||
		   tr.trid == TDEAD || tr.trid == TDUMMY);
@


1.5
log
@put in an explicit cast
@
text
@d1 1
a1 1
/* SUTVBAND: $Revision: 1.4 $ ; $Date: 93/11/24 16:26:00 $      */
d72 1
d82 4
a85 2
        if (tr.trid && tr.trid != TREAL)
                err("input is not seismic data, trid=%d", tr.trid);
@


1.4
log
@changed bcopy to memset ; added explicit casts to memcpy
@
text
@d1 1
a1 1
/* SUTVBAND: $Revision: 1.3 $ ; $Date: 92/10/23 15:11:37 $      */
d84 1
a84 1
        if (!getparfloat("dt", &dt))    dt = tr.dt/1000000.0;
@


1.3
log
@new selfdoc installed
@
text
@d1 1
a1 1
/* SUTVBAND: $Revision: 1.2 $ ; $Date: 92/03/24 12:54:29 $      */
d152 1
a152 1
			bzero(data, nt*FSIZE);
d182 2
a183 2
        memcpy(rt, data, nt*FSIZE);
        bzero(rt + nt, (nfft-nt)*FSIZE);
@


1.2
log
@First working (I think) version.
@
text
@d1 1
a1 1
/* SUTVBAND: $Revision: 1.10 $ ; $Date: 92/03/22 18:48:28 $      */
d20 20
a39 20
String sdoc =
"                                                               \n"
" SUTVBAND - time-variant bandpass filter (sine-squared taper)  \n"
"                                                               \n"
" sutvband <stdin >stdout tf= f=			        \n"
"                                                               \n"
" Required parameters:                                          \n"
"       dt = (from header)      time sampling rate (sec)        \n"
"       tf=             times for which f-vector is specified   \n"
"       f=f1,f2,f3,f4   Corner frequencies corresponding to the \n"
"                       times in tf. Specify as many f= as      \n"
"                       there are entries in tf.                \n"
"                                                               \n"
" The filters are applied in frequency domain.                  \n"
"                                                               \n"
" Example:                                                      \n"
" sutvband <data tf=.2,1.5 f=10,12.5,40,50 f=10,12.5,30,40 | ...\n"
"                                                               \n"
" Trace header fields accessed:  ns, dt, delrt.			\n"
;
d76 1
a76 1
        askdoc(1);
@


1.1
log
@Initial revision
@
text
@d1 1
a1 1
/* SUTVBAND: $Revision: 1.5 $ ; $Date: 90/12/17 17:50:27 $	*/
d20 19
a38 27
string sdoc =
" 								\n"
" SUTVBAND - time-variant bandpass filter (sine-squared taper)	\n"
" 								\n"
" sutvband <stdin >stdout [optional parameters]			\n"
" 							        \n"
" Required parameters:						\n"
" 	if dt is not set in header, then dt is mandatory	\n"
" 							        \n"
" Optional parameters: (nyquist calculated internally)		\n"
"	tfreq=0		times for which f is specified		\n"
"	f=f1,f2,f3,f4	Corner frequencies corresponding to the	\n"
"			times in tfreq. Specify as many f= as	\n"
"			there are entries in tfreq.  The 	\n"
"			default values of f1,f2,f3,f4 are:	\n"
" 							        \n"
" 	f1 = 0.10*(nyquist)	left  low  corner frequency (Hz)\n"
" 	f2 = 0.15*(nyquist)	left  high corner frequency (Hz)\n"
" 	f3 = 0.45*(nyquist)	right low  corner frequency (Hz)\n"
" 	f4 = 0.50*(nyquist)	right high corner frequency (Hz)\n"
" 	dt = (from header)	time sampling rate (sec)	\n"
" 							        \n"
" The filters are applied in frequency domain.			\n"
" 							        \n"
" Example:							\n"
" sutvband <data tfreq=.2,1.5 f=10,12.5,40,50 f=10,12.5,30,40 | ...	\n"
" 							        \n"
d43 1
a43 2
 *	CWP: Jack, Ken
 *
d47 3
a49 7
#define	piby2	1.57079632679490
#define FRAC1	0.10	/* Ratio of default f1 to Nyquist */
#define FRAC2	0.15	/* Ratio of default f2 to Nyquist */
#define FRAC3	0.45	/* Ratio of default f3 to Nyquist */
#define FRAC4	0.50	/* Ratio of default f4 to Nyquist */
#define LOOKFAC	2	/* Look ahead factor for npfao	  */
#define PFA_MAX	720720	/* Largest allowed nfft	          */
d52 3
a54 1
void makefilter(int if1, int if2, int if3, int if4, int nfft, float *filter);
d60 12
a71 16
	register float *rt;	/* real trace				*/
	register complex *ct;	/* complex transformed trace		*/
	float *filt;		/* filter array				*/
	float f1;		/* left lower corner frequency		*/
	float f2;		/* left upper corner frequency		*/
	float f4;		/* right lower corner frequency		*/
	float f3;		/* right upper corner frequency		*/
	int if1,if2,if3,if4;	/* integerizations of f1,f2,f3,f4	*/
	float dt;		/* sample spacing			*/
	float nyq;		/* nyquist frequency			*/
	int nt;			/* number of points on input trace	*/
	int nfft;		/* number of points for fft trace	*/
	int nf;			/* number of frequencies (incl Nyq)	*/
	int nfm1;		/* nf-1					*/
	float onfft;		/* reciprocal of nfft			*/
	float df;		/* frequency spacing (from dt)		*/
d73 4
a76 4
	
	/* Initialize */
	initargs(argc, argv);
	askdoc(1);
d79 8
a86 8
	/* Get info from first trace */ 
	if (!gettr(&tr))  err("can't get first trace");
	if (tr.trid && tr.trid != TREAL)
		err("input is not seismic data, trid=%d", tr.trid);
	nt = tr.ns;
	if (!getparfloat("dt", &dt))	dt = tr.dt/1000000.0;
	if (!dt) err("dt field is zero and not getparred");
	nyq = 0.5/dt;
d89 16
a104 1
	/* Set up FFT parameters */
d108 5
d114 9
a122 3
	nf = nfft/2 + 1;
	nfm1 = nf - 1;
	onfft = 1.0 / (float) nfft;
d124 18
a142 7
	/* Get corner frequencies */
	if (!getparfloat("f1", &f1))	f1 = FRAC1 * nyq;
	if (!getparfloat("f2", &f2))	f2 = FRAC2 * nyq;
	if (!getparfloat("f3", &f3))	f3 = FRAC3 * nyq;
	if (!getparfloat("f4", &f4))	f4 = FRAC4 * nyq;
	if (f1 < 0.0 || f1 > f2 || f2 >= f3 || f3 > f4)
		err("Bad filter parameters");
d144 26
a169 5

	/* Allocate fft arrays */
	rt   = ealloc1float(nfft);
	ct   = ealloc1complex(nf);
	filt = ealloc1float(nf);
d171 2
d174 6
a179 8
	/* Compute integer frequencies */
	df = onfft / dt;
	if1 = NINT(f1/df);
	if2 = NINT(f2/df);
	if3 = NINT(f3/df);
	if (if3 > nfm1) if3 = nfm1;
	if4 = NINT(f4/df);
	if (if4 > nfm1) if4 = nfm1;
d181 3
d185 4
a188 2
	/* Make filter with scale for inverse transform */
	makefilter(if1,if2,if3,if3,nfft,filt);
d190 2
a191 23



	/* Main loop over traces */
	do {
		register int i;

		/* Load trace into rt (zero-padded) */
		memcpy(rt, tr.data, nt*FSIZE);
		bzero(rt + nt, (nfft-nt)*FSIZE);

		/* FFT, filter, inverse FFT */
		pfarc(1, nfft, rt, ct);
		for (i = 0; i < nf; ++i)  ct[i] = crmul(ct[i], filt[i]);
		pfacr(-1, nfft, ct, rt);

		/* Load traces back in, recall filter had nfft factor */
		for (i = 0; i < nt; ++i)  tr.data[i] = rt[i];

		puttr(&tr);
	} while (gettr(&tr));

	return EXIT_SUCCESS;
d195 1
a195 2

void makefilter(int if1, int if2, int if3, int if4, int nfft, float *filter)
d197 7
a203 4
	int nf = nfft/2 + 1;
	float onfft = 1.0 / (float) nfft;
	register float s;
	register int i;
a204 4
	if (if1 == if2) {
		filter[if1] = onfft;
	} else {
		register float c = piby2 / ((float)(if2 - if1));
d206 7
a212 5
		for (i = if1; i <= if2; ++i) {
			s = sin(c*(float)(i - if1));
			filter[i] = s * s * onfft;
		}
	}
d214 7
a220 4
	if (if3 == if4) {
		filter[if3] = onfft;
	} else {
		register float c = piby2 / ((float)(if4 - if3));
d222 5
a226 11
		for (i = if3; i <= if4; ++i) {
			s = sin(c * (float)(if4 - i));
			filter[i] = s * s * onfft;
		}
	}

	{ register int i;
	  for (i = if2 + 1; i < if3;     ++i)  filter[i] = onfft; 
	  for (i = 0;       i < if1;     ++i)  filter[i] = 0.0; 
	  for (i = if4 + 1; i < nf; ++i)       filter[i] = 0.0; 
	}
@
