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

在地震勘探的数据处理过程中，我们经常遇到FFT的问题，我从SU中抽去了一部分代码，自己写了一个FFT的程序，可以做单道、多道、平均谱、和Log谱等，现在将代码贴出来，供大家参考：

/* 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 ( ir = 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);
}