地震资料数字处理程序

程序

1. 褶积滤波

# include "stdio.h" # include "math.h" # include "conv.c" # define pi 3.1415926 # define N 100 # define dt 0.002 main() {

float x[100], h[101], h1[101], y_low[200], y_band[200]; 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,*fp5; fp1=fopen("INPUT1.DAT","r"); for(i=0;i<=N;i++) { fscanf(fp1,"%f",&x[i]); } fp4=fopen("h_low.dat","w"); //低通滤波设计 for(i=-50;i<=50;i++) { if(i==0) h[i+50]=2*pi*f/pi; else h[i+50]=sin(2*pi*f*i*dt)/(pi*i*dt); fprintf(fp4,"%f ",h[i+50]); //output lowpass filter } fp2=fopen("synthesisdata_lowpass.DAT","w"); conv(x,m,h,n,y_low,l); for(i=50;i<l-50;i++) { fprintf(fp2,"%f\n",y_low[i]); } //带通滤波器 fp5=fopen("h_band.dat","w"); for(i=-50;i<=50;i++) { if(i==0) h1[i+50]=140;

else h1[i+50]=sin(2*pi*f2*i*dt)/(pi*i*dt)-sin(2*pi*f1*i*dt)/(pi*i*dt); fprintf(fp5,"%f\n",h1[i+50]); // output bandpass filter } fclose(fp5); conv(x,m,h1,n,y_band,l); fp3=fopen("synthesisdata_bandpass.DAT","w"); for(i=50;i<l-50;i++) { fprintf(fp3,"%f\n",y_band[i]); } fclose(fp1);fclose(fp2);fclose(fp3); }

2. 快变滤波

# include "stdio.h" # include "math.h" # include "stdlib.h" # include "fft.c"

# define pi 3.1415926

main() { double *xr,*xi; float *H; int i,np,nfft,k; float t,dt,df,f,z,fc1,fc2; FILE *fpar,*fp1,*fp2,*fp3; //从参数文件中获得截至频率 fpar=fopen("lowpassfilter.par","r"); fscanf(fpar,"%f%f",&fc1,&fc2); np=100; k=log(np)/log(2); if(np>pow(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<100;i++) {

fscanf(fp1,"%f",&z); xr[i]=z; } fclose(fp1); //补零至128位 for(i=100;i<nfft;i++) { xr[i]=0.0; } for(i=0;i<nfft;i++) { xi[i]=0.0; } //transfer to frequency doamin fft(xr,xi,k,1); //generate lowpass filter(zero-phase) for(i=0;i<nfft/2;i++) { f=i*df; if(f<=fc2 && f>=fc1)H[i]=1.0; else H[i]=0.0; } //滤波器对称 for(i=nfft/2;i<nfft;i++) {

f=i*df; H[i]=H[nfft-i]; } //filtering in frequency domain for(i=0;i<nfft;i++) { xr[i]=xr[i]*H[i]; xi[i]=xi[i]*H[i]; } fft(xr,xi,k,-1); fp2=fopen("lowpass2.dat","w"); for(i=0;i<nfft;i++) { fprintf(fp2,"%f\n",xr[i]); } fclose(fp2);

//获取高通截至频率

fpar=fopen("bandpass.par","r"); fscanf(fpar,"%f%f",&fc1,&fc2); fp1=fopen("INPUT1.DAT","r"); for(i=0;i<100;i++) {

fscanf(fp1,"%f",&z); xr[i]=z; } for(i=100;i<nfft;i++) { xr[i]=0.0; } for(i=0;i<nfft;i++) { xi[i]=0.0; }

//transfer to frequency doamin fft(xr,xi,k,1); //generate lowpass filter(zero-phase) for(i=0;i<=nfft/2;i++) { f=i*df; if(f<=fc2 && f>=fc1)H[i]=1.0; else H[i]=0.0; }

for(i=nfft/2+1;i<nfft;i++) { f=i*df; H[i]=H[nfft-i]; } //filtering in frequency domain for(i=0;i<nfft;i++) { xr[i]=xr[i]*H[i]; xi[i]=xi[i]*H[i]; } fft(xr,xi,k,-1); fp3=fopen("bandpass2.dat","w"); for(i=0;i<nfft;i++) { fprintf(fp3,"%f\n",xr[i]); }

}

3. 褶积滤波与递归滤波 褶积滤波

# include "stdio.h" # include "math.h" # include "stdlib.h" # include "conv.c" # include "fft.c"

# define PI 3.1415926 main() { void conv(); float x[50],h[20],y[69],hreverse[20],hzero[39],yreverse[69]; float dt=0.002; int i,m,n,l,p,q; FILE *fp1,*fp2,*fp3,*fp4,*fp5,*fp6; m=50;n=20; l=m+n-1; //read x(n) fp1=fopen("INPUT3.DAT","r"); for(i=0;i<50;i++) { fscanf(fp1,"%f",&x[i]); } fclose(fp1); //read filterfactor h(n) fp2=fopen("hn.dat","r"); fp5=fopen("h_reverse.dat","w"); for(i=0;i<20;i++) { fscanf(fp2,"%f",&h[i]); hreverse[i]=h[19-i]; fprintf(fp5,"%f\n",hreverse[i]); } fclose(fp2);fclose(fp5); conv(x,m,h,n,y,l);//非零相位褶积滤波 fp3=fopen("con_filter.dat","w"); for(i=0;i<l;i++) { fprintf(fp3,"%f\n",y[i]); } fclose(fp3);

p=n+n-1; q=m+p-1; //构造零相位滤波因子 conv(h,n,hreverse,n,hzero,p); fp6=fopen("zerophasefilterfactor.dat","w"); for(i=0;i<p;i++) { fprintf(fp6,"%f\n",hzero[i]); } fclose(fp6); //零相位滤波

conv(x,m,hzero,p,yreverse,q); fp4=fopen("convfilterreverse.dat","w"); for(i=0;i<q;i++) { fprintf(fp4,"%f\n",yreverse[i]); } fclose(fp4); }

递归滤波

#include<stdio.h> #include<math.h> #include<stdlib.h> #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)); fil1=(float*)malloc(np*sizeof(float)); fil2=(float*)malloc(np*sizeof(float)); //输入地震数据

fp1=fopen("INPUT3.DAT","r"); for(i=0;i<np;i++)fscanf(fp1,"%f",x+i); fclose(fp1); //输入a数组值 fp2=fopen("a(n).txt","r"); for(i=0;i<5;i++)fscanf(fp2,"%f",a+i);

fclose(fp2); for(i=5;i<np;i++)a[i]=0.0; //输入b数组值 fp3=fopen("b(n).txt","r"); for(i=0;i<5;i++)fscanf(fp3,& …… 此处隐藏：4040字，全部文档内容请下载后查看。喜欢就下载吧 ……