合成地震记录.zip

/* A Program of synthetic seismogram generation. */ #include "stdio.h" // #include "math.h" // #include "malloc.h" // #define fm 35 // 雷克子波的主频 #define dt 0.002 // 采样间隔 0.001 or 0.004. #define XL 0.060 // 子波长度 #define HL 0.300 // 反射系数长度 // Here the length is 0.3 second, you can change it. #define PI 3.1415926 main() { FILE *fp1,*fp2,*fp3; float *x,*h,*y; // 定义动态数组 float t; int i,M,N; void conv(); M=XL/dt+1; // The number of x(t). For this example, M=31. N=HL/dt+1; // The number of h(t). For this example, N=151. // 分配动态数组 x=(float *)calloc(M,sizeof(float)); h=(float *)calloc(N,sizeof(float)); y=(float *)calloc(M+N-1,sizeof(float)); // 生成雷克子波 fp1=fopen("d:\\Ricker(2).txt","at+"); for(i=-M/2;i<=M/2;i++) { t=i*dt; x[i+M/2]=(1.0-2.0*pow(PI*fm*t,2.0))*exp(-pow(PI*fm*t,2.0)); fprintf(fp1,"%8.3f%10.4f\n",t*1000,x[i+M/2]); } fclose(fp1); // 生成反射系数 fp2=fopen("d:\\Ref.txt","at+"); // 将初始值设置为零 for(i=0;i<=N-1;i++)h[i]=0.0; // Generate different Ref[i]. // You can design different reflection coefficient sequence. h[20]=0.4; h[80]=-0.3; h[90]=-0.5; h[130]=0.5; for(i=0;i<=N-1;i++) { t=i*dt; fprintf(fp2,"%8.3f%10.4f\n",t*1000,h[i]); } fclose(fp2); // 做卷积 x(n)*h(n). conv(x,M,h,N,y,M+N-1); // Output the results. Let the length of y(n) be equal to h(t), cut y(n). fp3=fopen("d:\\Seismic.txt","at+"); for(i=M/2;i<=N+M/2-1;i++) { t=(i-M/2)*dt; fprintf(fp3,"%8.3f %10.4f\n",t*1000,y[i]); } fclose(fp3); free(x); free(h); free(y); } // 卷积 void conv(float x[],int m,float h[],int n,float y[],int l) { int k,i; for(k=0;k<l;k++) { y[k]=0.0; for(i=0;i<m;i++) if(k-i>=0&&k-i<=n-1)y[k]=y[k]+x[i]*h[k-i]; } }