| 查看: 1921 | 回复: 3 | |||
| 【奖励】 本帖被评价2次，作者baobiao007增加金币 1.2 个 | |||
[资源]
【分享】频率波数域二维扇形低通滤波c程序示例【原创】
|
|||
|
利用二维傅里叶变换将地震记录由x-t转换到f-k域在进行滤波，这种方法主要利用有效波与干扰波之间的视速度差异来分离两类信号。 对于学习地震数据处理的朋友，牢固地掌握好数字信号处理是很有好处的。尤其是傅里叶变换，用途非常广泛。这个程序进一步展示了离散傅里叶频谱的性质（共轭对称性），这种性质在编程实现时是尤其要注意的。理论可以参照任何一本地震数据处理教材。 #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); } //设置扇形滤波器频谱,一定要注意二维傅里叶变换的共轭对称性 void FilterFan(float v,float fn,float df,float dk,float ***h,int m,int n) { int i, j; int m2, n2; float b; m2=m/2; n2=n/2; //设置第一象限 for(i=0; i<=m2; i++) for(j=0; j<=n2; j++) { b=i*df/(v*dk); if(i*df<=fn && j<=b) (*h)[j]=1.0;//想得到另一条线,只需if else的1.0和0.0对调 else (*h)[j]=0.0; } //设置第二象限 for(i=0; i<=m2; i++) for(j=n2+1; j //设置第三象限 for(i=m2+1; i //设置第四象限 for(i=m2+1; i } void main() {//////////设置参数,根据情况自行修改//////////////// const int M=64;//行数 const int N=256;//列数 const float V=3700.0;//视速度，重要滤波参数 const float Fn=100.0;//频率边界,根据效果调节 const float dt=0.002;//2ms时间采样间隔，重要参数 const float dx=10.0;//接收道间距10米，重要参数 //////////////////////////////////////////////////// int i,j; float **record1;//要滤波的二维信号的实部 float **record2;//虚部 float **mid1;//中间结果 float **mid2;//中间结果 float **filter;//滤波器 float df;//频率采样间隔 float dk;//波数采样间隔 char fname[]="test-64-256.dat";//待分析信号文件名 char fresult[]="after64-256.dat";//滤波结果 char ffuliye[]="pu-64-256.dat";//信号的二位频谱 //离散傅里叶变换必须满足的条件 df=1.0/(N*dt); dk=1.0/(M*dx); printf("%f\n",df); //申请空间 record1=MySpace(M,N); record2=MySpace(M,N); mid1=MySpace(N,M); mid2=MySpace(N,M); filter=MySpace(N,M); //读入待滤波信号 MyReadFile(fname,record1,M,N); //对待滤波信号进行二维傅里叶变换 FFT2D(record1,record2,M,N,1); MyWriteFile(ffuliye,record1,M,N); //产生扇形滤波器 FilterFan(V,Fn,df,dk,&filter,N,M); //转置傅里叶变换后的信号频谱 Zhuan(&record1,&mid1,M,N); Zhuan(&record2,&mid2,M,N); //频谱相乘 for(i=0; i mid1[j]=mid1[j]*filter[j]; mid2[j]=mid2[j]*filter[j]; } //再一次矩阵转置 Zhuan(&mid1,&record1,N,M); Zhuan(&mid2,&record2,N,M); //傅里叶反变换 FFT2D(record1,record2,M,N,-1); //输出滤波结果 MyWriteFile(fresult,record1,M,N); //释放空间 FreeMySpace(&record1,M); FreeMySpace(&record2,M); FreeMySpace(&mid1,N); FreeMySpace(&mid2,N); FreeMySpace(&filter,N); } 原始地震剖面是两条交叉的同相轴：  它的二维频谱：  分离后的图像：   分离效果不太理想，原因是多方面的,既与参数的选择合理性有关，还可能因为没有考虑吉布斯效应等，有待进一步研究。 [ Last edited by baobiao007 on 2011-2-16 at 08:48 ] |
回复此楼» 收录本帖的淘帖专辑推荐
 科技写作与绘图 |
» 猜你喜欢
首次在隐伏矿体上方发现金、硫酸铅、硝酸铅、三氧化钨等含金属纳米微粒
已经有8人回复
-
求助：有没有大神有流体包裹体计算软件
已经有0人回复
-
地球化学论文润色/翻译怎么收费?
已经有83人回复
-
青C最后一年，散金祈福吧
已经有48人回复
-
没中，有点难受，科研也许就到这了，回归生活吧
已经有15人回复
-
282求调剂
已经有0人回复
-
地质学复试调剂
已经有1人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- c程序设计语言 非常适合初学者的
已经有175人回复
- 急求一个matlab编码，各位如果知道一定要帮我一把啊
已经有10人回复
- C++编程，关于循环结构的，大家看看我这程序哪里错了？
已经有21人回复
- 请问二维地震剖面与三维地震剖面有什么区别？！
已经有6人回复
- 求程序 模拟点蚀的元胞自动机程序
已经有8人回复
- Matlab求解偏微分方程组
已经有22人回复
- 安卓王牌软件三连发2（安卓系统一键ROOT，快速重启，二维码识别）
已经有72人回复
- 怎么判断晶体结构是一维？二维？三维啊？
已经有15人回复
- 波特图上频率发生了位移应该怎么解释啊？
已经有6人回复
- 请问，结晶过程是几维增长是怎样回事？例如一维，二维（细节请进）
已经有3人回复
- 求一道有具体解题步骤的二维核磁共振波普例题
已经有6人回复
- 请教一个关于滤波方面的问题
已经有4人回复
- 如何对方波进行FIR滤波生成此波形
已经有9人回复
- 求助C语言有限元一维、二维、矩形网格以及三角形网格程序
已经有4人回复
2楼2011-06-25 12:35:11
3楼2016-01-19 21:57:58
简单回复
2016-01-19 23:21
回复
yes
