《最小平方反滤波程序代码C语言以及MATLAB编写(共11页).doc》由会员分享，可在线阅读，更多相关《最小平方反滤波程序代码C语言以及MATLAB编写(共11页).doc（11页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、精选优质文档-倾情为你奉上分C语言和MATLAB两部分第一部分：C语言#include stdio.h#include math.h#include stdlib.h#include malloc.h#include string.h#define pi 3.14#define fm 25 /主频#define dt 0.001 /采样间隔void con(float a,float b,float c,int M,int L);/褶积void zixiangguan(float *a,float *r,int P);/自相关int jiefangchengzu(float R,int n,f
2、loat b,float x);/解托普利茨方程组void main() int i,j;int P=200;/子波长度int Q=399;/反射序列长度float a;/子波衰减因子float wave200,ref200,con1399,xcorr200,b200,at200,yt598;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);/地震记录自相关
3、fp_at=fopen(at.csv,w);/反滤波因子fp_yt=fopen(yt.csv,w);/因子与地震记录褶积for(i=0;iP;i+)refi=0.;ref50=0.2; ref80=0.4; ref150=-0.7; /构造反射序列/for(i=0;iP;i+)a=2*fm*fm*log(2)/3;wavei=cos(2*pi*fm*i*dt)*exp(-a*i*i*dt*dt);fprintf(fp_wave,%fn,wavei); / 构造子波 /con(wave,ref,con1,P,P);for(i=0;iQ;i+)fprintf(fp_con1,%fn,con1i);
4、 / 褶积生成地震记录 /zixiangguan(wave,xcorr,200);for(i=0;iP;i+)fprintf(fp_xcorr,%fn,xcorri); / 子波自相关 /for(i=0;iP;i+)bi=0;b0=1; / 构造方程组右侧结果数组 /jiefangchengzu(xcorr,P,b,at);for(i=0;iP;i+)fprintf(fp_at,%fn,ati); / 解方程组，求反滤波因子 /con(con1,at,yt,Q,P);for(i=0;iP+Q-2;i+)fprintf(fp_yt,%fn,yti); / 因子与地震记录褶积 / 普通褶积 /vo
5、id con(float a,float b,float c,int M,int L)int i,j,N;float tp=0;N=M+L-1;for(i=0;iN;i+) for(j=0;j=0&(i-j)L)tp+=aj*b(i-j); ci=tp;tp=0;/ 函数自相关 /void zixiangguan(float *a,float *r,int P) int i,j;float s=0;for(i=0;iP;i+)for(j=0;j=i;j+)s=s+aj*aP-1-i+j;/从最左边开始，移到两者重合rP-1-i=s;s=0;/ 莱文森递推解托普利茨方程组 / /R为矩阵的第一行
6、或者第一列数据，b为方程组右侧结果， x为计算的反滤波因子 int jiefangchengzu(float R,int n,float b,float x) int i,j,k; float a,beta,q,c,h,*y,*s; s=(float*)calloc(n,sizeof(double); y=(float*)calloc(n,sizeof(double); a=R0; if (fabs(a)+1.0=1.0) free(s); free(y); printf(failn); return(-1); y0=1.0;x0=b0/a; for(k=1;k=n-1;k+) beta=0.
7、0;q=0.0; for(j=0;j=k-1;j+) beta=beta+Rj+1*yj; q=q+Rk-j*xj; 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
8、=xi+h*si;yi=si; xk=h*yk; free(s); free(y); return(1); 第二部分：MATLABf=25;%频率a=2/3*f*f*log(2);dt=0.001;%采样时间for i=1:200wave(i)= cos(2*pi*f*i*dt)*exp(-a*i*i*dt*dt);%生成子波endplot(wave)for i=1:200ref(i)=0;endref(50)=0.7;ref(80)=-0.4;ref(150)=0.5;%构造反射序列con1=conv(wave,ref);%褶积生成地震记录for i=1:200wav(i)=wave(201
9、-i);%把wave反一下endf3=conv(wave,wav);%wave与wav褶积相当于wave自相关for i=1:200t(i)=f3(i+199);%利用自相关后200个数据endT=toeplitz(t);% 用t构造托普利兹矩阵，相当于线性方程组的系数b=zeros(200,1);%构造一个矩阵，线性方程组右侧的结果b(1,1)=1; %构造一个矩阵，相当于线性方程组右侧的结果at=Tb;%解方程组，求出反滤波因子atyt=conv(con1,at);%地震记录与反滤波因子褶积子波地震记录反射序列子波自相关截取自相关后半段反滤波因子Yt ：滤波因子at与地震记录褶积结果与最初反射系数一致，实现反滤波专心-专注-专业