基于SU的快速傅里叶变换（FFT）

/* Copyright (c) China University of Petroleum , 2018.           */
/* All rights reserved.                                          */

/* CSIR: $Revision: 1.0 $ ; $Date: 2018/2/15  7:40:56 $        */

#include "par.h"
#include "su.h"
#include "segy.h"
#include "time.h"
#include "math.h"
#include "cwp.h"

#define LOOKFAC 2       /* Look ahead factor for npfaro   */
#define PFA_MAX 720720  /* Largest allowed nfft           */
#define N 222
#define pi 3.1415926535898

char *sdoc[] = {"   ",NULL};

/*
* Authors:
*
* References:
*
*/
/******************************** end self doc ********************************/

void FFT_tf(int nt,int nr,int nf,int nfft,float **dobs, complex **dxf);

segy instr;

int
main(int argc, char **argv)
{
time_t start,finish;
double duration;

int verbose;

int ir,it;
int nr,nt;
int nw;                 /* number of wavelet length          */
float dx,dt;

float **dobs;           /* arry of original seismic data     */
float **dobs_tw;        /* arry of original seismic data     */

/* file names */
char *amp_name ="";
char *dB_name  ="";

char fnm[BUFSIZ];
FILE *seisfp;           /* pointer to input obs seicmic data */
FILE *outpfp;           /* pointer to output obs seicmic data*/

/* Time window processing */
int itw;                /* counters                   */
int ntw;                /* The number of time windows */
int twlength;           /* The length of time window  */

int *tw;                /* arry of time windows       */

/* seismic data FFT T<->F */
int nfft;               /* number of points on output trace  */
int fnum;               /* number of frequencies             */
int nf;                 /* number of frequencies             */
float df;
float sum_amp;
float amp_max;
float dB_max;
complex **dxf;          /* arry of shot data in x-f domain   */
float **dxf_amp;        /* arry of shot data in x-f domain   */
float *spectrum_amp;
float *spectrum_dB;

/* hook up getpar to handle the parameters */
initargs(argc,argv);
requestdoc(0);
time(&start);

/* get requested parameters  */
if (!getparint("nr",&nr))           	err("must specify nr!");
if (!getparint("nt",&nt))           	err("must specify nt!");
if (!getparint("nw",&nw))           	err("must specify nw!");
if (!getparint("verbose",&verbose)) 	err("must specify verbose!");
if (!getparfloat("dx",&dx))         	err("must specify dx!");
if (!getparfloat("dt",&dt))         	err("must specify dt!");
if(!getparstring("amp_name",&_name)) err("must specify amp_name!");
if(!getparstring("dB_name",&dB_name)) 	err("must specify dB_name!");

warn("nr=%d,nt=%d;dx=%f,dt=%f",nr,nt,dx,dt);

/* read seismic data */
dobs  = alloc2float(nt,nr);
seisfp=stdin;
for ( i
4000
r = 0; ir < nr; ++ir)
{
fgettr(seisfp,&instr);
for (it = 0; it < nt; ++it)
{
dobs[ir][it] = instr.data[it];
}
}

if (verbose != 0 )
{
ntw = countparval("tw");
tw  = alloc1int(ntw);
if (!getparint("tw",tw))
{
err(" Using time window processing , need tw");
}
warn("The number of time windows is %d .",ntw/2);

for (itw = 0; itw < ntw/2; ++itw)
{
twlength = tw[2*itw+1]-tw[2*itw]+1;
warn("Time window %d : %d~%d",itw+1,tw[2*itw],tw[2*itw+1]);
warn("twlength=%d",twlength);

/* select data by time windows */
dobs_tw = alloc2float(twlength,nr);
for (ir = 0; ir < nr; ++ir)
{
for (it = tw[2*itw]-1; it < tw[2*itw+1]; ++it)
{
dobs_tw[ir][it-(tw[2*itw]-1)] = dobs[ir][it];
}
}

/* FFT the seismic data from t-x fo f-x */
/* Set up pfa fft T <-> F*/
nfft = npfaro(twlength, LOOKFAC * twlength);
if (nfft >= SU_NFLTS || nfft >= PFA_MAX)
err("Padded nt=%d--too big", nfft);
fnum= nfft/2 + 1;
nf = fnum;
df = 1.0/(twlength*dt);
warn("FFT:nf=%d; df=%f",nf,df);

/* Allocate fft arrays */
dxf = alloc2complex(nf,nr);
dxf_amp = alloc2float(nf,nr);
spectrum_amp = alloc1float(nf);
spectrum_dB  = alloc1float(nf);

/* FFT from t-x fo f-x */
memset((void *) dxf[0], 0, sizeof(complex)*nf*nr);
FFT_tf(twlength,nr,nf,nfft,dobs_tw,dxf);

/* calculate spectrum */
for (ir = 0; ir < nr ; ++ir)
{
for (it = 0; it < nf; ++it)
{
dxf_amp[ir][it] = (float)sqrt( pow(dxf[ir][it].r,2) +
pow(dxf[ir][it].i,2)) / 2.0;
}
}

/* calculate average amplitude spectrum */
for (it = 0; it < nf; ++it)
{
sum_amp = 0;
for (ir = 0; ir < nr ; ++ir)
{
sum_amp += dxf_amp[ir][it];
}
spectrum_amp[it] = sum_amp/nr;
spectrum_dB [it] =  20*log10(spectrum_amp[it]);
}

/* norm */
amp_max = 0;
dB_max = 0;
for (it = 0; it < nf; ++it)
{
amp_max = MAX(amp_max,spectrum_amp[it]);
dB_max  = MAX(dB_max,spectrum_dB[it]);
}

for (it = 0; it < nf; ++it)
{
spectrum_amp[it] = spectrum_amp[it]/amp_max;
spectrum_dB [it] = spectrum_dB [it]/dB_max;
}

/* Output data */
sprintf(fnm,"./Data/tmpt/spectrum_amp_%s_%d.bin",amp_name,itw+1);
outpfp=fopen(fnm,"wb");
fwrite(spectrum_amp,sizeof(float),nf,outpfp);
fclose(outpfp);

/* calculate dB spectrum */
sprintf(fnm,"./Data/tmpt/spectrum_dB_%s_%d.bin",dB_name,itw+1);
outpfp=fopen(fnm,"wb");
fwrite(spectrum_dB,sizeof(float),nf,outpfp);
fclose(outpfp);

free2float(dobs_tw);
free1float(spectrum_amp);
free1float(spectrum_dB);
free2float(dxf_amp);
free2complex(dxf);
}
}else{
/* FFT the seismic data from t-x fo f-x */
/* Set up pfa fft T <-> F*/
nfft = npfaro(nt, LOOKFAC * nt);
if (nfft >= SU_NFLTS || nfft >= PFA_MAX)
err("Padded nt=%d--too big", nfft);
fnum= nfft/2 + 1;
nf = fnum;
df = 1.0/(nt*dt);
warn("FFT:nf=%d; df=%f",nf,df);

/* Allocate fft arrays */
dxf = alloc2complex(nf,nr);
dxf_amp = alloc2float(nf,nr);
spectrum_amp = alloc1float(nf);
spectrum_dB  = alloc1float(nf);

/* FFT from t-x fo f-x */
memset((void *) dxf[0], 0, sizeof(complex)*nf*nr);
FFT_tf(nt,nr,nf,nfft,dobs,dxf);

/* calculate spectrum */
for (ir = 0; ir < nr ; ++ir)
{
for (it = 0; it < nf; ++it)
{
dxf_amp[ir][it] = (float)sqrt( pow(dxf[ir][it].r,2) +
pow(dxf[ir][it].i,2)) / 2.0;
}
}

/* calculate average amplitude spectrum */
for (it = 0; it < nf; ++it)
{
sum_amp = 0;
for (ir = 0; ir < nr ; ++ir)
{
sum_amp += dxf_amp[ir][it];
}
spectrum_amp[it] = sum_amp/nr;
spectrum_dB [it] =  20*log10(spectrum_amp[it]);
}

/* norm */
amp_max = 0;
dB_max = 0;
for (it = 0; it < nf; ++it)
{
amp_max = MAX(amp_max,spectrum_amp[it]);
dB_max  = MAX(dB_max,spectrum_dB[it]);
}

for (it = 0; it < nf; ++it)
{
spectrum_amp[it] = spectrum_amp[it]/amp_max;
spectrum_dB [it] = spectrum_dB [it]/dB_max;
}

/* Output data */
sprintf(fnm,"./Data/tmpt/spectrum_amp_%s_all.bin",amp_name);
outpfp=fopen(fnm,"wb");
fwrite(spectrum_amp,sizeof(float),nf,outpfp);
fclose(outpfp);

/* calculate dB spectrum */
sprintf(fnm,"./Data/tmpt/spectrum_dB_%s_all.bin",dB_name);
outpfp=fopen(fnm,"wb");
fwrite(spectrum_dB,sizeof(float),nf,outpfp);
fclose(outpfp);

free1float(spectrum_amp);
free1float(spectrum_dB);
free2float(dxf_amp);
free2complex(dxf);
}
free1int(tw);
free2float(dobs);

time(&finish);
duration=difftime(finish,start);
warn("Program Complete! It takes %4.2f MS.",1000.0*duration/60);
warn("***********CSIR End***************");
return(CWP_Exit());
}

void FFT_tf(int nt,int nr,int nf,int nfft,float **dobs, complex **dxf)
{
register float *rt1;    /* real trace               */
register complex *ct1;  /* complex transformed trace        */
int ir,iff,it;

/* Allocate fft arrays */
rt1 = ealloc1float(nfft);
ct1 = ealloc1complex(nf);

/* Main loop over traces */
for(ir=0;ir<nr;++ir)
{
/* Load trace into rt (zero-padded) */
for(it=0;it<nt;++it)
{
rt1[it]=dobs[ir][it];
}

memset((void *) (rt1 + nt), 0, (nfft-nt)*FSIZE);
/* FFT (T->F) */
pfarc(1, nfft, rt1, ct1);
/* output f-x domain shot  data */
for(iff=0;iff<nf;++iff)
{
dxf[ir][iff].r=ct1[iff].r;
dxf[ir][iff].i=ct1[iff].i;
}
}
free1float(rt1);
free1complex(ct1);
}