基于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.1415926535898char *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);
}

基于SU的快速傅里叶变换（FFT）相关推荐

基于python的快速傅里叶变换FFT（二）
基于python的快速傅里叶变换FFT(二) 本文在上一篇博客的基础上进一步探究正弦函数及其FFT变换. 知识点 FFT变换,其实就是快速离散傅里叶变换,傅立叶变换是数字信号处理领域一种很重要的算 ...
基于python的快速傅里叶变换FFT（一）
基于python的快速傅里叶变换FFT(一) FFT可以将一个信号变换到频域.有些信号在时域上是很难看出什么特征的,但是如果变换到频域之后,就很容易看出特征了.这就是很多信号分析采用FFT变换的原因. ...
快速傅里叶变换python_基于python的快速傅里叶变换FFT（二）
基于python的快速傅里叶变换FFT(二) 本文在上一篇博客的基础上进一步探究正弦函数及其FFT变换. 知识点 FFT变换,其实就是快速离散傅里叶变换,傅立叶变换是数字信号处理领域一种很重要的算法. ...
基于MCU的快速傅里叶变换FFT调试记录
FFT主要用于将信号从时域变换到频域,时域变换到频域的主要目的是为了更方便分析信号,比如轴承振动异常的频率想通过时域信号去分析出来是哪个频率就很难分析,而通过把时域信号变化到频域的话就可以很清楚直接的 ...
MATLAB之傅里叶变换，快速傅里叶变换FFT
文章目录傅里叶变换及傅里叶逆变换定义窗函数/矩形脉冲信号的傅里叶变换基于MATLAB的快速傅里叶变换FFT 傅里叶变换及傅里叶逆变换定义能从时域的非周期连续信号转化到频域非周期连续信号. 窗函 ...
OpenCV快速傅里叶变换(FFT)用于图像和视讯流的模糊检测
OpenCV快速傅里叶变换(FFT)用于图像和视频流的模糊检测翻译自[OpenCV Fast Fourier Transform (FFT) for blur detection in images ...
基于FPGA的快速傅里叶变换加速（三）
基于FPGA的快速傅里叶变换加速(三) 硬件加速介绍及部分verilog代码实现 1. 硬件加速 1.1 FPGA 1.1.1 FPGA介绍概念: 基本结构: 工作原理: 1.1.2 开发板开发板 ...
MIT 线性代数 Linear Algebra 26:复矩阵，傅里叶矩阵, 快速傅里叶变换 FFT
这一讲我们来讲一下复矩阵.线性代数中,复矩阵是避免不了的话题,因为一个简单实矩阵都有可能有复数特征值. 复矩阵我们着重看一下复矩阵和实矩阵在运算上的区别. 距离首先,一个复数向量的的距离求法发生了 ...
Java中实现快速傅里叶变换FFT
Java中实现快速傅里叶变换FFT 一.概述 1.傅里叶变换(FT) 2.离散傅里叶变换(DFT) 3.快速傅里叶变换(FFT) 1)单位根 2)快速傅里叶变换的思想 3)蝶形图 4)快速傅里叶变换的 ...

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

基于SU的快速傅里叶变换（FFT）相关推荐

最新文章

热门文章