《地震资料数字处理程序》由会员分享，可在线阅读，更多相关《地震资料数字处理程序（26页珍藏版）》请在金锄头文库上搜索。
1、 第 1 页 程序程序 1. 褶积滤波褶积滤波 # include “stdio.h“ # include “math.h“ # include “conv.c“ # define pi 3.1415926 # define N 100 # define dt 0.002 main() float x100, h101, h1101, y_low200, y_band200; float df; int i,m=100,n=101,l=m+n- 1; float f=70.0; float f1=10.0; float f2=80.0; FILE *fp1,*fp2,*fp3,*fp4,*fp
2、5; fp1=fopen(“INPUT1.DAT“,“r“); for(i=0;ipow(2,k)k=k+1; nfft=pow(2,k); dt=0.002; df=1.0/(nfft*dt); xr=(double*)calloc(nfft,sizeof(double); xi=(double*)calloc(nfft,sizeof(double); H=(float*)calloc(nfft,sizeof(float); / read x(n) fp1=fopen(“INPUT1.DAT“,“r“); for(i=0;i=fc1)Hi=1.0; else Hi=0.0; /滤波器对称 f
3、or(i=nfft/2;i=fc1)Hi=1.0; else Hi=0.0; for(i=nfft/2+1;i #include #include #define np 50 void main() float *x,*a,*b,*fil1,*fil2; int i; void recur1(); void recur2(); FILE *fp1,*fp2,*fp3,*fp4,*fp5; x=(float*)malloc(np*sizeof(float); a=(float*)malloc(np*sizeof(float); b=(float*)malloc(np*sizeof(float);
4、 fil1=(float*)malloc(np*sizeof(float); fil2=(float*)malloc(np*sizeof(float); /输入地震数据 fp1=fopen(“INPUT3.DAT“,“r“); for(i=0;i=0;i- - ) y3i=0.0; y4i=0.0; for(j=0;jpow(2,k)k=k+1; nfft=pow(2,k); h=(double*)malloc(nfft*sizeof(double); hi=(double*)malloc(nfft*sizeof(double); for(i=- 50;ipow(2,k)k=k+1; nfft
5、=pow(2,k); x1=(double*)malloc(nfft*sizeof(double); xi1=(double*)malloc(nfft*sizeof(double); x2=(double*)malloc(nfft*sizeof(double) xi2=(double*)malloc(nfft*sizeof(double); for(i=0;i=0 i=0j- - ) hj+1=hj; h0=temp; yk=0.0; for(j=0;j=0 taoi=j; taoi=taoi- 99; printf(“%dn“,taoi); /静校正处理 for(i=0;i=0) for(j
6、=0;j=0;j- - ) k=j+taoi; if(k=0)xij=xik; fp2=fopen(“processeddata.DAt“,“w“); for(i=0;i=0linp;li+=4) j=li;k=j+2; t1=srj- srk;t2=sxj- sxk; srj+=srk; sxj+=sxk; srk=t1;sxk=t2; +j;+k; t1=srj- srk;t2=sxj- sxk; srj+=srk;sxj+=sxk; srk=inv*t2;sxk=- inv*t1; /* 2 points DFT */ for(li=0;linp;li+=2) j=li;k=j+1; t
7、1=srj- srk;t2=sxj- sxk; srj+=srk;sxj+=sxk; srk=t1;sxk=t2; /* sort according to bit reversal */ lmx=np/2;j=0; for(i=1;inp- 1;+i) k=lmx; while(k=j) 第 25 页 j- =k;k/=2; j+=k; if(ij) t1=srj;srj=sri;sri=t1; t1=sxj;sxj=sxi;sxi=t1; /* if Inverse FFT, multiply 1.0/np */ if(inv!=- 1) return; t1=1.0/np; for(i=
8、0;inp;+i) sri*=t1;sxi*=t1; 求解 tolpriz 方程 /函数名：tlvs.c # include “stdio.h“ # include “math.h“ # include “stdlib.h“ int tlvs(double t,int n,double b,double x) int i,j,k; double a,beta,q,c,h,*y,*s; s=malloc(n*sizeof(double); y=malloc(n*sizeof(double); a=t0; if (fabs(a)+1.0=1.0) free(s); free(y); printf(
9、“failn“); return(- 1); y0=1.0; x0=b0/a; for (k=1; k=n- 1; k+) beta=0.0; q=0.0; for (j=0; j=k- 1; j+) 第 26 页 beta=beta+yj*tj+1; q=q+xj*tk- j; if (fabs(a)+1.0=1.0) free(s); free(y); printf(“failn“); return(- 1); c=- beta/a; s0=c*yk- 1; yk=yk- 1; if (k!=1) for (i=1; i=k- 1; i+) si=yi- 1+c*yk- i- 1; a=a+c*beta; if (fabs(a)+1.0=1.0) free(s); free(y); printf(“failn“); return(- 1); h=(bk- q)/a; for (i=0; i=k- 1; i+) xi=xi+h*si; yi=si; xk=h*yk; free(s); free(y); return(1)
