最小平方反滤波程序代码C语言以及MATLAB编写

分C语言和MATLAB两部分第一部分：C语言#includestdio.h#includemath.h#includestdlib.h#includemalloc.h#includestring.h#definepi3.14#definefm25//主频#definedt0.001//采样间隔voidcon(floata[],floatb[],floatc[],intM,intL);//褶积voidzixiangguan(float*a,float*r,intP);//自相关intjiefangchengzu(floatR[],intn,floatb[],floatx[]);//解托普利茨方程组voidmain(){inti,j;intP=200;//子波长度intQ=399;//反射序列长度floata;//子波衰减因子floatwave[200],ref[200],con1[399],xcorr[200],b[200],at[200],yt[598];FILE*fp_wave,*fp_con1,*fp_xcorr,*fp_at,*fp_yt;fp_wave=fopen(wave.csv,w);//子波fp_con1=fopen(con1.csv,w);//子波与反射序列褶积的地震记录fp_xcorr=fopen(xcorr.csv,w);//地震记录自相关fp_at=fopen(at.csv,w);//反滤波因子fp_yt=fopen(yt.csv,w);//因子与地震记录褶积for(i=0;iP;i++)ref[i]=0.;ref[50]=0.2;ref[80]=0.4;ref[150]=-0.7;//构造反射序列//for(i=0;iP;i++){a=2*fm*fm*log(2)/3;wave[i]=cos(2*pi*fm*i*dt)*exp(-a*i*i*dt*dt);fprintf(fp_wave,%f\n,wave[i]);//构造子波//}con(wave,ref,con1,P,P);for(i=0;iQ;i++)fprintf(fp_con1,%f\n,con1[i]);//褶积生成地震记录//zixiangguan(wave,xcorr,200);for(i=0;iP;i++)fprintf(fp_xcorr,%f\n,xcorr[i]);//子波自相关//for(i=0;iP;i++)b[i]=0;b[0]=1;//构造方程组右侧结果数组//jiefangchengzu(xcorr,P,b,at);for(i=0;iP;i++){fprintf(fp_at,%f\n,at[i]);//解方程组，求反滤波因子//}con(con1,at,yt,Q,P);for(i=0;iP+Q-2;i++)fprintf(fp_yt,%f\n,yt[i]);//因子与地震记录褶积//}//普通褶积//voidcon(floata[],floatb[],floatc[],intM,intL){inti,j,N;floattp=0;N=M+L-1;for(i=0;iN;i++){for(j=0;jM;j++){if((i-j)=0&&(i-j)L)tp+=a[j]*b[(i-j)];}c[i]=tp;tp=0;}}//函数自相关//voidzixiangguan(float*a,float*r,intP){inti,j;floats=0;for(i=0;iP;i++){for(j=0;j=i;j++){s=s+a[j]*a[P-1-i+j];//从最左边开始，移到两者重合}r[P-1-i]=s;s=0;}}//莱文森递推解托普利茨方程组////R为矩阵的第一行或者第一列数据，b为方程组右侧结果，x为计算的反滤波因子intjiefangchengzu(floatR[],intn,floatb[],floatx[]){inti,j,k;floata,beta,q,c,h,*y,*s;s=(float*)calloc(n,sizeof(double));y=(float*)calloc(n,sizeof(double));a=R[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+R[j+1]*y[j];q=q+R[k-j]*x[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);}第二部分：MATLABf=25;%频率a=2/3*f*f*log(2);dt=0.001;%采样时间fori=1:200wave(i)=cos(2*pi*f*i*dt)*exp(-a*i*i*dt*dt);%生成子波endplot(wave)fori=1:200ref(i)=0;endref(50)=0.7;ref(80)=-0.4;ref(150)=0.5;%构造反射序列con1=conv(wave,ref);%褶积生成地震记录fori=1:200wav(i)=wave(201-i);%把wave反一下endf3=conv(wave,wav);%wave与wav褶积相当于wave自相关fori=1:200t(i)=f3(i+199);%利用自相关后200个数据endT=toeplitz(t);%用t构造托普利兹矩阵，相当于线性方程组的系数b=zeros(200,1);%构造一个矩阵，线性方程组右侧的结果b(1,1)=1;%构造一个矩阵，相当于线性方程组右侧的结果at=T\b;%解方程组，求出反滤波因子atyt=conv(con1,at);%地震记录与反滤波因子褶积子波地震记录反射序列子波自相关截取自相关后半段反滤波因子Yt：滤波因子at与地震记录褶积结果与最初反射系数一致，实现反滤波