应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 专业学科区 » 地学 » 地球物理 » 【分享】频率波数域FK偏移c程序【原创】
31/1返回列表
查看: 2418  |  回复: 2
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看
【奖励】 本帖被评价2次，作者baobiao007增加金币 1.4 个

baobiao007

木虫 (职业作家)

[资源] 【分享】频率波数域FK偏移c程序【原创】

#include
#include
#include
#include"FFT.h"

//将文件filename中的内容读入二维数组a
void MyReadFile(char *filename,float **a,int m,int n)
{
     int i,j;
         FILE *fp;
         fp=fopen(filename,"rb" );
         for(i=0; i                  for(j=0; j                          fread(&a[ i ][j],sizeof(float),1,fp);
         fclose(fp);
}
//将二维数组a中的数据写入文件filename中
void MyWriteFile(char *filename,float **a,int m,int n)
{
     int i,j;
         FILE *fp;
     fp=fopen(filename,"wb" );
         for(i=0; i                  for(j=0; j                          fwrite(&a[ i ][j],sizeof(float),1,fp);
         fclose(fp);
}
//申请二维动态数组的函数
float **MySpace(int m, int n)
{
        int i;
        float **p;
        p=(float **)calloc(m,sizeof(float *));
        if(p==NULL)
        {  printf("申请空间失败\n";
           exit(1);
        }
        for(i=0; i                 p=(float *)calloc(n,sizeof(float));
        return p;
}
//释放申请的二维动态数组
void FreeMySpace(float ***p, int m)
{
        int i;
        for(i=0; i                 free((*p)[ i ]);
        free(*p);
}
//矩阵转置a[m][n]->b[n][m]
void Zhuan(float ***a,float ***b,int m,int n)
{
        int i,j;
        for(i=0; i                 for(j=0; j                         *(*(*b+i)+j)=*(*(*a+j)+i);
}
//2D-FFT tag=1,正变换;tag=-1,反变换
void FFT2D(float **xr,float **xi,int m,int n,int tag)
{
        int i,k1,k2;
        float **xtr, **xti;
        xtr=MySpace(n,m);
        xti=MySpace(n,m);
        k1=log(m)/log(2.0);
        k2=log(n)/log(2.0);
        if(n>pow(2,k1)) k1=k1+1;
        if(m>pow(2,k2)) k2=k2+1;
    for(i=0; i         fft(xr[ i ],xi[ i ],k2,tag);
    Zhuan(&xr,&xtr,m,n);
    Zhuan(&xi,&xti,m,n);
    for(i=0; i         fft(xtr,xti,k1,tag);
        Zhuan(&xtr,&xr,n,m);
        Zhuan(&xti,&xi,n,m);
        FreeMySpace(&xtr,n);
        FreeMySpace(&xti,n);       
}
//sinc插值
void SincIns(float **pur,float **pui,float ***midr,float ***midi,
                           float v,int m,int n,float dkx,float dkz,float df)
{
        int i,j,l;
        float pi=3.1415926;
        float kx, kz;//圆波数
        float dw;//圆频率采样间隔
        float lm,k,a,b;
        dw=2.0*pi*df;
        for(i=0; i<=m/2; i++)
            for(j=0; j<=n/2; j++)               
                {
                        kx=i*2.0*pi*dkx;
                        kz=j*2.0*pi*dkz;
                        k=0.5*v*sqrt(kx*kx+kz*kz)/dw;//速度要减半(爆炸反射界面原理)
            l=(int)k;
            lm=k-l;
                        a=pi*lm;
                        b=1.0-lm;
                        if(fabs(a-0)<1e-6 || fabs(b-0)<1e-6) continue;
            (*midr)[ i ][j]=pur[ i ][ l ]*cos(a)*sin(a)/a
                                +sin(a)*sin(a)/a*pui[ i ][ l ]
                                +pur[ i ][l+1]*cos(b*pi)*sin(b*pi)/(b*pi)
                                -pui[ i ][l+1]*sin(b*pi)*sin(b*pi)/(b*pi);
                        (*midi)[j]=pui[ i ][ l ]*cos(a)*sin(a)/a
                                -sin(a)*sin(a)/a*pur[ i ][ l ]+pui[ i ][l+1]*cos(b*pi)*sin(b*pi)/(b*pi)
                                +pur[ i ][l+1]*sin(b*pi)*sin(b*pi)/(b*pi);

                }
        for(i=m/2+1; i                 for(j=0; j<=n/2; j++)
                {
                        (*midr)[ i ][j]= (*midr)[m-i][j];
                        (*midi)[ i ][j]=(*midi)[m-i][j];
                }
}
void main()
{
///初始参数设置区///////////////////////////////////
        const int DNUM=64;//零偏移距剖面接收道数       
        const int TNUM=256;//时间采样点数
    const int ZNUM=TNUM;//深度采样点数
        const float V=2000.0;//上覆介质速度
        const float DT=0.002;//时间采样间隔2ms
        const float DX=10.0;//道间距
        const float DZ=DT*V*0.5;//深度步长(速度要减半)
        const float pi=3.1415926;
    char fname1[]="before.dat";//原始地震记录文件名
        char fname2[]="after.dat";//偏移后的剖面文件名
        char fname3[]="mid.dat";//频谱图
////其他参数计算区//////////////////////////////////
        float **recordr, **recordi;
        float **midr, **midi;
        float df;//频率采样间隔
        float dkx;//x方向波数采样间隔
        float dkz;//z方向波数采样间隔
        float dex;
        int i,j;
    df=1.0/(DT*TNUM); dkx=1.0/(DX*DNUM);
        dkz=1.0/(DZ*ZNUM);
        //申请空间
        recordr=MySpace(DNUM,TNUM);
        recordi=MySpace(DNUM,TNUM);
        midr=MySpace(DNUM,TNUM);
        midi=MySpace(DNUM,TNUM);
        //读入原始剖面
        MyReadFile(fname1,recordr,DNUM,TNUM);
        //对原始剖面进行二维傅里叶变换
    FFT2D(recordr,recordi,DNUM,TNUM,1);
        MyWriteFile(fname3,recordr,DNUM,TNUM);
        //对变换结果进行插值
    SincIns(recordr,recordi,&midr,&midi,V,DNUM,TNUM,dkx,dkz,df);
        //乘因子
    for(i=0; i             for(j=0; j                 {   dex=0.5*V*j*dkz/sqrt(pow(i*dkx,2)+pow(j*dkz,2));
                        if(fabs(dex-0)<1e-6) continue;
                    midr[j]=midr[ i ][j]*dex;
            midi[ i ][j]=midi[ i ][j]*dex;                       
                }
        //傅里叶反变换
        FFT2D(midr,midi,DNUM,TNUM,-1);
        MyWriteFile(fname2,midr,DNUM,TNUM);
        //释放空间
        FreeMySpace(&recordr,DNUM);
        FreeMySpace(&recordi,DNUM);
        FreeMySpace(&midr,DNUM);
        FreeMySpace(&midi,DNUM);
}
五点脉冲的偏移结果：


[ Last edited by baobiao007 on 2011-2-18 at 03:48 ]
回复此楼

» 收录本帖的淘帖专辑推荐

科技写作与绘图

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
1楼 2011-02-18 03:33:54
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zmc

金虫 (正式写手)

★★★ 三星级,支持鼓励

沙发
回复此楼
2楼2011-03-21 10:13:04
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

biyadi

铁虫 (初入文坛)

★★★★★ 五星级,优秀推荐

c:\documents and settings\administrator\桌面\新建文件夹\c.c(4) : fatal error C1083: Cannot open include file: 'FFT.h': No such file or directory
执行 cl.exe 时出错.
我运行时总会出现以上问题 怎么解决呢
赞一下 回复此楼
3楼2011-12-14 15:56:15
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 baobiao007 的主题更新
31/1返回列表
☆ 无星级 ★ 一星级 ★★★ 三星级 ★★★★★ 五星级
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
信息提示
关闭
请填处理意见
关闭