地震记录合成代码

工程物探上有一个实验要求合成地震记录，还是比较简单的，就是反射系数序列和雷克子波卷积就可。卷积公式这儿就不写出来了。反射系数序列可以自己设置，200个点同时赋值为0，地下有几层，就把几个点设置为其他的值。对于雷克子波，有公式就可以做出来，对于工程物探，频率尽量大一些，这样地震记录的图要好看一些。

下边把代码奉献给大家。

#include<stdio.h>
#include<math.h>
#include<stdlib.h>
//#define dz1 "file1.txt"
//#define dz2 "file2.txt"

main ()
{ FILE *fp1,*fp2;
double dt=0.002,t[71];int fm1=70,i,j;double zibo[71],xishu[200],dzjl[300];
// fp1=fopen(dz1,"w+");
// fp2=fopen(dz2,"w+");
// fclose(fp1);
//fclose(fp2);

fp1=fopen("file1.txt","r+");
fp2=fopen("file2.txt","r+");
for(i=0;i<=35;i++)

{
t[i+35]=i*dt;

zibo[i+35]=exp(-pow(3.141592653*fm1*t[i+35],2))*(1-2*pow(3.141592653*fm1*t[i+35],2));

}
//for(i=0;i<71;i++)
// printf("%f\n",zibo[i]);
for(i=0;i<=34;i++)
{
t[i]=-t[70-i];
zibo[i]=zibo[70-i];

}
for(i=0;i<=70;i++)
{printf("%f\t%f\n",t[i],zibo[i]);
fprintf(fp1," %f\t%f\n",t[i],zibo[i]);

}

for(i=0;i<=199;i++)
{
xishu[i]=0;
}
xishu[10]=0.2;xishu[65]=0.4;xishu[110]=-0.5;xishu[155]=0.7;xishu[190]=-0.36;
for(i=0;i<200;i++)

{ dzjl[i]=0.0;

for(j=0;j<70;j++)

{ if((i-j)>0)

dzjl[i]=dzjl[i]+zibo[j]*xishu[i-j] ; }

}
for(i=0;i<200;i++)
{
// printf("%d\t%f\n",i,dzjl[i]);
fprintf(fp2,"%d\t%f\n",i,dzjl[i]);
}
fclose(fp1);
fclose(fp2);

}