地震资料处理实习报告

《地震资料数据处理》课程设计总结报告专业班级：地球物理学1302班姓名：学号：成绩：2016年12月31日目录一、设计内容………………………………………………………………（1）褶积滤波……………………………………………………………（2）快变滤波……………………………………………………………（3）褶积滤波与快变滤波的比较………………………………………（4）设计高通滤波因子…………………………………………………（5）频谱分析……………………………………………………………（6）分析补零对振幅谱的影响…………………………………………（7）线性褶积与循环褶积………………………………………………（8）最小平方反滤波……………………………………………………（9）零相位转换…………………………………………………………（11）静校正………………………………………………………………二、程序…………………………………………………………………………一、设计内容1、褶积滤波理想低通滤波因子理想带通滤波因子原始信号低通滤波结果带通滤波结果2、快变滤波原始数据快变低通滤波结果快变带通滤波结果3．褶积滤波与递归滤波的比较原始数据零相位褶积滤波结果非零相位褶积滤波结果零相位递归滤波非零相位递归滤波4．设计高通滤波因子原始数据高通滤波因子频率域高通滤波因子5．频谱分析正弦函数频谱尖脉冲频谱地震波频谱6．分析补零对振幅谱的影响正弦函数n=60正弦函数n=128\非周期波形n=60 非周期波形n=64非周期波形n=128 7．线性褶积与圆周褶积线性褶积模型圆周褶积结果长度与圆周褶积相等的线性褶积8．最小平方反滤波原始反射系数序列求出的反射系数序列9．零相位转换非零相位子波零相位子波11、地震记录原始数据选择第一道为参考到静校正结果二、程序1. 褶积滤波# include "stdio.h"# include "math.h"# include "conv.c"# define pi 3.1415926# define N 100# define dt 0.002main(){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;elseh[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;elseh1[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.1415926main(){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 doaminfft(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 domainfor(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 doaminfft(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;elseH[i]=0.0;}for(i=nfft/2+1;i<nfft;i++){f=i*df;H[i]=H[nfft-i];}//filtering in frequency domainfor(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.1415926main(){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 50void 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,"%f",b+i);fclose(fp3);for(i=5;i<np;i++)b[i]=0.0;//正向递归滤波recur1(x,a,b,fil1); fp4=fopen("正向递归结果.DAT","wb"); for(i=0;i<np;i++){fprintf(fp4,"%12.4f\n",fil1[i]);}fclose(fp4);for(i=0;i<np;i++){printf("%12.4f\n",fil1[i]);}printf("\n");//反向递归滤波recur2(fil1,a,b,fil2); fp5=fopen("反向递归结果.DAT","wb"); for(i=0;i<np;i++){fprintf(fp5,"%12.4f\n",fil2[i]);}fclose(fp5);for(i=0;i<np;i++){printf("%12.4f\n",fil2[i]);}}void recur1(float x[],float a[],float b[],float fil1[]) {int i,j,k;float y1[np],y2[np];for(i=0;i<np;i++){ y1[i]=0.0;y2[i]=0.0;for(j=0;j<=i;j++)y1[i]=y1[i]+a[j]*x[i-j];if(i==0)y2[i]=0.0;else{ for(k=1;k<=i;k++)y2[i]=y2[i]+b[k]*fil1[i-k];}fil1[i]=y1[i]-y2[i];}}void recur2(float fil1[],float a[],float b[],float fil2[]) {int i,j,k;float y3[np],y4[np];for(i=np-1;i>=0;i--){else{} y3[i]=0.0;y4[i]=0.0;for(j=0;j<=np-1-i;j++)y3[i]=y3[i]+a[j]*fil1[i+j];if(i==np-1)y4[i]=0.0;for(k=1;k<=np-1-i;k++)y4[i]=y4[i]+b[k]*fil2[i+k]; }fil2[i]=y3[i]-y4[i];}4. 设计高通滤波因子# include "stdio.h"# include "math.h"# include "stdlib.h"# include "fft.c"# define N 101# define dt 0.004# define PI 3.1415926main(){double *h,*hi;int i,k,nfft;FILE *fp1,*fp2;k=log(N)/log(2);if(N>pow(2,k))k=k+1;nfft=pow(2,k);h=(double*)malloc(nfft*sizeof(double));hi=(double*)malloc(nfft*sizeof(double));for(i=-50;i<=50;i++){if(i==0)h[i+50]=1.0/dt-60;elseh[i+50]=-sin(2*PI*30.0*i*dt)/(PI*i*dt);}for(i=100;i<128;i++){h[i]=0.0;}for(i=0;i<128;i++){hi[i]=0.0;}fp1=fopen("timedomain.dat","w");for(i=0;i<128;i++){fprintf(fp1,"%f\n",h[i]);}fclose(fp1);fft(h,hi,k,1);fp2=fopen("frequencydoamin.dat","w");for(i=0;i<128;i++){fprintf(fp2,"%f\n",sqrt(h[i]*h[i]+hi[i]*hi[i]));}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.004main(){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.004main(){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.004main(){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 101main(){ 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;elseh[i+50]=sin(2*PI*70*i*dt)/(PI*i*dt);}//linear convolutionconv(x,m,h,n,y,l);fp1=fopen("linearconv.dat","w");for(i=0;i<l;i++){fprintf(fp1,"%f\n",y[i]);}//circus convolutionfp=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;elseh1[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 seismic recordsconv(bn,n,a,m,x,l);fp11=fopen("synthseismicdata.dat","w");for(i=0;i<211;i++){fprintf(fp11,"%f\n",x[i]);}autocorr(x,rxx,l);fp12=fopen("autocorrdata.dat","w");for(i=0;i<l;i++){fprintf(fp12,"%f\n",rxx[i]);}for(i=0;i<60;i++){if(i==0)d[i]=1.0;elsed[i]=0.0;t[i]=rxx[i];}tlvs(t,60,d,b);//输出反子波fp13=fopen("waveletreverse.dat","w");for(i=0;i<60;i++){fprintf(fp13,"%f\n",b[i]);}conv(x,l,b,60,a_cal,l+60-1);fp2=fopen("sigma.dat","w");for(i=0;i<270;i++){fprintf(fp2,"%f\n",a_cal[i]);}}//褶积子程序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]*0.004;}}//自相关子程序void autocorr(float x[],float y[],int n){int i,j;for(i=0;i<n;i++){y[i]=0.0;for(j=i;j<n;j++)y[i]=x[j]*x[j-i]+y[i];}return;}8. 零相位转换# include "stdio.h"# include "math.h"# include "stdlib.h"# include "fft.c"# define PI 3.1415926main(){void fft();double *xr,*xi;float z,w[25];int i,nfft=32;FILE *fp,*fp1,*fp2;xr=(double*)malloc(nfft*sizeof(double));xi=(double*)malloc(nfft*sizeof(double));fp=fopen("wavelet.dat","r");for(i=0;i<25;i++){fscanf(fp,"%f",&z);xr[i]=z;}for(i=25;i<32;i++){xr[i]=0.0;}for(i=0;i<32;i++){xi[i]=0.0;}fft(xr,xi,5,1);fp1=fopen("amplitude.dat","w");for(i=0;i<32;i++){fprintf(fp1,"%f\n",sqrt(xr[i]*xr[i]+xi[i]*xi[i]));}for(i=0;i<32;i++){ xr[i]=sqrt(xr[i]*xr[i]+xi[i]*xi[i]); xi[i]=0.0;}fft(xr,xi,5,-1);fp2=fopen("changedwavelet.dat","w");for(i=0;i<32;i++){fprintf(fp2,"%f\n",xr[i]);}}10. 最小相位转换//最小相位转换# include "stdio.h"# include "math.h"# include "tlvs.c"main(){void conv();void autocorr();double temp,bo[73],ao[25],g[49],gr[49];float z;double b[25],rbb[25],d[25],y[25],c[49];int i;FILE *fp,*fp1,*fp2;fp=fopen("wavelet.dat","r");for(i=0;i<25;i++){fscanf(fp,"%f",&z);b[i]=z;//printf("%f\n",b[i]);}autocorr(b,25,b,25,rbb);for(i=0;i<25;i++){if(i==0)d[i]=1.0;elsed[i]=0.0;}//求y(t)tlvs(rbb,25,d,y);//printf("it is ok !\n");fp1=fopen("yt.dat","w");for(i=0;i<25;i++){fprintf(fp1,"%f\n",y[i]);}//求bo(0)conv(b,25,y,25,c,49);//c(t)=b(t)*g(t)temp=0.0;for(i=0;i<49;i++){temp=temp+c[i]*c[i];printf("%f\n",c[i]);}bo[0]=1.0/sqrt(temp);printf("bo(0)=%f\n",bo[0]);//求ao(t)for(i=0;i<25;i++){ao[i]=bo[0]*y[i];}conv(b,25,ao,25,g,49);//g(t)=b(t)*ao(t)//求g(-t)for(i=0;i<49;i++){gr[i]=g[48-i];}//求b o(t)=b(t)*g(-t)conv(b,25,gr,49,bo,73);fp2=fopen("minphasew.dat","w");for(i=0;i<73;i++){fprintf(fp2,"%f\n",bo[i]);}}//褶积子程序void conv(double x[],int m,double h[],int n,double 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]*0.004;}}//自相关子程序void autocorr(double x[],int m,double h[],int n,double y[] ) {int i,j;for(i=0;i<m;i++){y[i]=0.0;for(j=1;j<n;j++)if(i+j<m)y[i]=x[j]*h[j+i]+y[i];}return;}11. 静校正处理# include "stdio.h"# include "math.h"main(){void ccor();float x[10][100],y[100],rxy[10][199];float z,max;int i,j,k,tao[10];FILE *fp,*fp1,*fp2;//从文件中读取地震数据fp=fopen("seis.dat","r");for(i=0;i<10;i++){for(j=0;j<100;j++){fscanf(fp,"%f",&z);x[i][j]=z;}}fclose(fp);//求参考道for(j=0;j<100;j++){y[j]=0.0;for(i=0;i<10;i++){y[j]=y[j]+x[i][j];}y[j]=0.1*y[j];//printf("%f\n",y[j]);}//互相关求取静校正量值for(i=0;i<10;i++){ccor(x[i],100,y,100,rxy[i]);}fp1=fopen("rxy.dat","w");for(i=0;i<10;i++){for(j=0;j<199;j++){ fprintf(fp1,"%f\n",rxy[i][j]);}}fclose(fp1);//求取tao 值for(i=0;i<10;i++){max=rxy[i][0];tao[i]=0;for(j=0;j<199;j++){if(rxy[i][j]>max){max=rxy[i][j];tao[i]=j;}}tao[i]=tao[i]-99;printf("%d\n",tao[i]);}//静校正处理for(i=0;i<10;i++){if(tao[i]>=0){ for(j=0;j<100;j++){k=j+tao[i];if(k<100)x[i][j]=x[i][k];}}if(tao[i]<0){for(j=99;j>=0;j--){k=j+tao[i];if(k>=0)x[i][j]=x[i][k];}}}fp2=fopen("processeddata.DAt","w");for(i=0;i<10;i++){for(j=0;j<100;j++){ fprintf(fp2,"%f\n",x[i][j]);}}fclose(fp2);}//自相关程序void ccor(float x[],int m,float h[],int n,float y[]) {int i,j;for(i=-n+1;i<=m-1;i++){y[i+n-1]=0.0;for(j=0;j<=m-1;j++)if(j-i>=0&&j-i<=n-1)y[i+n-1]=y[i+n-1]+x[j]*h[j-i];}}部分子程序：//程序名：fft.c#include "stdio.h"#include "math.h"#include "stdlib.h"#define PI 3.1415926/* sr,sx：双精度型一维数组,输入（输出）信号的实部和虚部*/ /* m0: 2 的次方数, 2**m0=nfft *//* inv=1 forward transform; inv=-1 inverse transform */void fft(double sr[],double sx[],int m0,int inv){int i,j,lm,li,k,lmx,np,lix,mm2;double t1,t2,c,s,cv,st,ct;if(m0<0)return;lmx=1;for(i=1;i<=m0;++i)lmx+=lmx; //form 2**m0cv=2.0*PI/(double)lmx;ct=cos(cv); st=-inv*sin(cv);np=lmx;mm2=m0-2;/* fft butterfly numeration */for(i=1;i<=mm2;++i){lix=lmx;lmx/=2;c=ct;s=st;for(li=0;li<np;li+=lix){j=li;k=j+lmx;t1=sr[j]-sr[k];t2=sx[j]-sx[k];sr[j]+=sr[k];sx[j]+=sx[k];sr[k]=t1;sx[k]=t2;++j;++k;t1=sr[j]-sr[k];t2=sx[j]-sx[k];sr[j]+=sr[k];sx[j]+=sx[k];sr[k]=c*t1-s*t2;sx[k]=s*t1+c*t2;}for(lm=2;lm<lmx;++lm){cv=c;c=ct*c-st*s;s=st*cv+ct*s;for(li=0;li<np;li+=lix){j=li+lm;k=lmx+j;t1=sr[j]-sr[k];t2=sx[j]-sx[k];sr[j]+=sr[k];sx[j]+=sx[k];sr[k]=c*t1-s*t2;sx[k]=s*t1+c*t2;}}cv=ct;ct=2.0*ct*ct-1.0;st=2.0*st*cv;}/* 4 points DFT */if(m0>=2)for(li=0;li<np;li+=4){j=li;k=j+2;t1=sr[j]-sr[k];t2=sx[j]-sx[k];sr[j]+=sr[k];sx[j]+=sx[k];sr[k]=t1;sx[k]=t2;++j;++k;t1=sr[j]-sr[k];t2=sx[j]-sx[k];sr[j]+=sr[k];sx[j]+=sx[k];sr[k]=inv*t2;sx[k]=-inv*t1;}/* 2 points DFT */for(li=0;li<np;li+=2){j=li;k=j+1;t1=sr[j]-sr[k];t2=sx[j]-sx[k];sr[j]+=sr[k];sx[j]+=sx[k];sr[k]=t1;sx[k]=t2;}/* sort according to bit reversal */lmx=np/2;j=0;for(i=1;i<np-1;++i){k=lmx;while(k<=j){j-=k;k/=2;}j+=k;if(i<j){t1=sr[j];sr[j]=sr[i];sr[i]=t1;t1=sx[j];sx[j]=sx[i];sx[i]=t1;}}/* if Inverse FFT, multiply 1.0/np */if(inv!=-1)return;t1=1.0/np;for(i=0;i<np;++i){sr[i]*=t1;sx[i]*=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=t[0];if (fabs(a)+1.0==1.0){free(s);free(y);printf("fail\n");return(-1);}y[0]=1.0;x[0]=b[0]/a;for (k=1; k<=n-1; k++){beta=0.0; q=0.0;for (j=0; j<=k-1; j++){beta=beta+y[j]*t[j+1];q=q+x[j]*t[k-j];}if (fabs(a)+1.0==1.0){free(s);free(y);printf("fail\n");return(-1);}c=-beta/a; s[0]=c*y[k-1]; y[k]=y[k-1];if (k!=1)for (i=1; i<=k-1; i++)s[i]=y[i-1]+c*y[k-i-1];a=a+c*beta;if (fabs(a)+1.0==1.0){free(s);free(y);printf("fail\n");return(-1);}h=(b[k]-q)/a;for (i=0; i<=k-1; i++){x[i]=x[i]+h*s[i];y[i]=s[i];}x[k]=h*y[k];}free(s); free(y);return(1)}。