| 查看: 2465 | 回复: 2 | ||||
| 【奖励】 本帖被评价2次，作者baobiao007增加金币 1.4 个 | ||||
[资源]
【分享】频率波数域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 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 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 return p; } //释放申请的二维动态数组 void FreeMySpace(float ***p, int m) { int i; for(i=0; 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 } //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 Zhuan(&xr,&xtr,m,n); Zhuan(&xi,&xti,m,n); for(i=0; i 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 { (*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 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 ] |
回复此楼» 收录本帖的淘帖专辑推荐
 科技写作与绘图 |
» 猜你喜欢
地球物理学和空间物理学论文润色/翻译怎么收费?
已经有50人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- 在聚甲基丙烯酸酯类聚合物上的羟基引入C=C双键
已经有5人回复
- 一个根据线长、线径、绕匝等参数计算【电感】的小程序
已经有5人回复
- Matlab求解偏微分方程组
已经有22人回复
- 支持向量机训练好的参数如何用C代码写测试程序
已经有5人回复
- 求《环境工程学报》发表的有关程序~~~有编辑部电话的最好
已经有8人回复
- C语言变量问题，求助大侠！！
已经有17人回复
- 波特图上频率发生了位移应该怎么解释啊？
已经有6人回复
- 请教Windows系统下perl脚本的ARGV应用问题，谢谢
已经有5人回复
- 关于瞬时频率估计
已经有12人回复
- 同事说c语言之父去世了
已经有30人回复
- 【求助】求助关于XPS单质硅的峰偏移问题
已经有10人回复
- 【讨论】关于阻抗频率
已经有8人回复
- 【求助】帮我编个程啊~C语言，C++，Matlab，VB...都行~
已经有14人回复
- 求助C语言有限元一维、二维、矩形网格以及三角形网格程序
已经有4人回复
2楼2011-03-21 10:13:04
3楼2011-12-14 15:56:15
