zhejiyunsuan.txt

Visual C++语言编写

/*  
     A Program of synthetic seismogram generation. 
                    Author: Zhang Junhua     Date: Sep.,2006     */ 

#include "stdio.h"		// standard input/output head file.
#include "math.h"		// math. head file.
#include "malloc.h"		// memory allocation head file.

#define fm  35			// The main frequency of Ricker wavelet, you can change it.   
#define dt  0.002		// The sample interval, you can also set it to 0.001 or 0.004.
#define	XL  0.060		// The length of wavelet(s), you can change it.
#define	HL  0.300		// The length of h(t), namely reflection coefficient.
                        // Here the length is 0.3 second, you can change it.

#define PI 3.1415926

main()
{ 
  FILE *fp1,*fp2,*fp3;

  float *x,*h,*y;	// define the dynamic array. 
  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.

// allocate dynamic array.
  x=(float *)calloc(M,sizeof(float));
  h=(float *)calloc(N,sizeof(float));
  y=(float *)calloc(M+N-1,sizeof(float));

// Generate Ricker wavelet.

  fp1=fopen("Ricker.txt","w");
  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);

// Generate reflection coefficient.

  fp2=fopen("Ksai.txt","w");
// set initial value to zero.
  for(i=0;i<=N-1;i++)h[i]=0.0;

// Generate different ksai[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);

// Do convolution 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("Seis.txt","w");
  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);

}

// The function of convolution.

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];
    }
}