地震资料数字处理程序(2)

fclose(fp2); printf("it is ok!\n"); }

5. 分析补零对振幅谱的影响 1、 不补零

# include "stdio.h" # include "math.h" # include "dft.c" # define N 60

# define PI 3.1415926 # define dt 0.004 main() { float x[N],xr[N],xi[N],w[N],wr[N],wi[N],z; int i,k; FILE *fp,*fp1,*fp2,*fp3;

fp3=fopen("sin.dat","w"); for(i=0;i<N;i++) { x[i]=sin(2.0*PI*(i+1)/30); fprintf(fp3,"%f\n",x[i]); } fclose(fp3); fp=fopen("WAVE.DAT","r"); for(i=0;i<N;i++) { fscanf(fp,"%f",&z); w[i]=z; } fclose(fp);

printf("it is ok !\n"); dft(x,xr,xi,N); dft(w,wr,wi,N); fp1=fopen("frequencydomain1_60.dat","w"); fp2=fopen("frequencydomain2_60.dat","w"); for(i=0;i<N;i++) {

fprintf(fp1,"%f\n",xr[i]); fprintf(fp2,"%f\n",wr[i]); } fclose(fp1);fclose(fp2); }

2、 补到64位 # include "stdio.h" # include "math.h" # include "stdlib.h" # include "fft.c" # define N 60

# define PI 3.1415926 # define dt 0.004 main() { void fft(); double *x1,*xi1,*x2,*xi2; float z; int i,k,nfft; FILE *fp,*fp1,*fp2,*fp3,*fp4; k=log(N)/log(2);

if(N>pow(2,k))k=k+1; nfft=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<N;i++) { x1[i]=sin(2*PI*(i+1)/30); } fp=fopen("WAVE.DAT","r"); for(i=0;i<N;i++) { fscanf(fp,"%f",&z); x2[i]=z; } for(i=0;i<N;i++) { xi1[i]=0.0; xi2[i]=0.0; } //补到64位 for(i=N;i<64;i++) { x1[i]=0.0; x2[i]=0.0; xi1[i]=0.0; xi2[i]=0.0; } fft(x1,xi1,6,1); fft(x2,xi2,6,1); fp1=fopen("frequencydomain1_64.dat","w"); fp2=fopen("frequencydomain2_64.dat","w"); for(i=0;i<nfft;i++) {

fprintf(fp1,"%f\n",sqrt(x1[i]*x1[i]+xi1[i]*xi1[i])*dt); fprintf(fp2,"%f\n",sqrt(x2[i]*x2[i]+xi2[i]*xi2[i])*dt); } }

3、 补到128位 # include "stdio.h" # include "math.h" # include "stdlib.h" # include "fft.c" # define N 60

# define PI 3.1415926 # define dt 0.004 main() { void fft(); double *x1,*xi1,*x2,*xi2; float z; int i,nfft; FILE *fp,*fp1,*fp2; nfft=128; 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<N;i++) { x1[i]=sin(2*PI*(i+1)/30); } fp=fopen("WAVE.DAT","r"); for(i=0;i<N;i++) { fscanf(fp,"%f",&z); x2[i]=z; } for(i=0;i<N;i++) { xi1[i]=0.0; xi2[i]=0.0; } //补到128位 for(i=N;i<128;i++)

{ x1[i]=0.0; x2[i]=0.0; xi1[i]=0.0; xi2[i]=0.0; } fft(x1,xi1,7,1); fft(x2,xi2,7,1); fp1=fopen("frequencydomain1_128.dat","w"); fp2=fopen("frequencydomain2_128.dat","w"); for(i=0;i<nfft;i++) {

fprintf(fp1,"%f\n",sqrt(x1[i]*x1[i]+xi1[i]*xi1[i])*dt); fprintf(fp2,"%f\n",sqrt(x2[i]*x2[i]+xi2[i]*xi2[i])*dt); } }

6. 线性褶积与循环褶积 # include "stdio.h" # include "math.h" # include "fft.c" # define dt 0.002

# define PI 3.1415926 # define L 101

main()

{ void conv(); void cir_conv(); float x[100],h[101],x1[L],h1[L],y[200],y1[100],df; int i,m,n,l,kc; FILE *fp,*fp1,*fp2; m=100;n=101;l=m+n-1; fp=fopen("INPUT1.DAT","r"); for(i=0;i<100;i++) { fscanf(fp,"%f",&x[i]); }

fclose(fp);

for(i=-50;i<=50;i++) { if(i==0)h[i+50]=140.0; else h[i+50]=sin(2*PI*70*i*dt)/(PI*i*dt); } //linear convolution

conv(x,m,h,n,y,l); fp1=fopen("linearconv.dat","w"); for(i=0;i<l;i++) { fprintf(fp1,"%f\n",y[i]); } //circus convolution fp=fopen("INPUT1.DAT","r"); for(i=0;i<100;i++) { fscanf(fp,"%f",&x1[i]); }

fclose(fp); df=1.0/(dt*100); kc=70.0/df; for(i=0;i<=50;i++) { if(i>=0 && i<=kc)h1[i]=1.0; else h1[i]=0.0; } for(i=50;i<=100;i++) { h1[i]=h1[100-i]; } if(L>=100) { for(i=100;i<L;i++){ x1[i]=0.0; h1[i]=0.0; } } cir_conv(x1,h1,y1,L); fp2=fopen("circusconv.dat","w"); for(i=0;i<L;i++) { fprintf(fp2,"%f\n",y1[i]); printf("%f\n",y1[i]); } fclose(fp2); }

//线性褶积子程序

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

//圆周褶积子程序

void cir_conv(float x[],float h[],float y[],int n) { int k,j; int temp; for(k=0;k<n/2;k++) { temp=h[k]; h[k]=h[n-k-1]; h[n-k-1]=temp; } for(k=0;k<n;k++) { temp=h[n-1]; for(j=n-2;j>=0;j--) h[j+1]=h[j]; h[0]=temp; y[k]=0.0; for(j=0;j<n;j++) y[k]=x[j]*h[j]*dt+y[k]; } return; }

7. 最小平方反滤波 # include "stdio.h" # include "math.h" # include "tlvs.c" # define N 200

//bn是地震子波序列 a是反射系数系列 main() {

void conv(); void autocorr(); float bn[60],a[200],x[211],rxx[211],a_cal[270]; double t[60],b[60],d[60]; int i,m,n,l; FILE *fp,*fp1,*fp2,*fp11,*fp12,*fp13; n=12;m=N;l=m+n-1; fp=fopen("INPUT8.DAT","r"); for(i=0;i<200;i++) { fscanf(fp,"%f",&a[i]); } fclose(fp); fp1=fopen("bn.dat","r"); for(i=0;i<12;i++) { fscanf(fp1,"%f",&bn[i]); } fclose(fp1); //synthetize s …… 此处隐藏：3925字，全部文档内容请下载后查看。喜欢就下载吧 ……