24小时热门版块排行榜

返回列表

【奖励】本帖被评价2次，作者baobiao007增加金币 1.2 个

baobiao007

木虫 (职业作家)

应助: 201 (大学生)
金币: 6482.7
帖子: 3050
在线: 1009.9小时
虫号: 505962

[资源] 【分享】频率波数域二维扇形低通滤波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    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));
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[ i ],xti[ i ],k1,tag);
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 (*h)[j]=(*h)[n-j];
//设置第三象限
for(i=m2+1; i for(j=n2+1; j (*h)[j]=(*h)[m-i][n-j];
//设置第四象限
for(i=m2+1; i for(j=0; j (*h)[j]=(*h)[m-i][j];
}
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 for(j=0; j {
         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 ]