最小平方反滤波程序代码C语言以及MATLAB编写.docx
- 文档编号:2824081
- 上传时间:2022-11-15
- 格式:DOCX
- 页数:12
- 大小:297.33KB
最小平方反滤波程序代码C语言以及MATLAB编写.docx
《最小平方反滤波程序代码C语言以及MATLAB编写.docx》由会员分享,可在线阅读,更多相关《最小平方反滤波程序代码C语言以及MATLAB编写.docx(12页珍藏版)》请在冰豆网上搜索。
最小平方反滤波程序代码C语言以及MATLAB编写
分C语言和MATLAB两部分
第一部分:
C语言
#include"stdio.h"
#include"math.h"
#include"stdlib.h"
#include"malloc.h"
#include"string.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;i
ref[i]=0.;
ref[50]=0.2;ref[80]=0.4;ref[150]=-0.7;//构造反射序列//
for(i=0;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;i fprintf(fp_con1,"%f\n",con1[i]);//褶积生成地震记录// zixiangguan(wave,xcorr,200); for(i=0;i fprintf(fp_xcorr,"%f\n",xcorr[i]);//子波自相关// for(i=0;i b[i]=0; b[0]=1;//构造方程组右侧结果数组// jiefangchengzu(xcorr,P,b,at); for(i=0;i { fprintf(fp_at,"%f\n",at[i]);//解方程组,求反滤波因子// } con(con1,at,yt,Q,P); for(i=0;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;i { for(j=0;j { if((i-j)>=0&&(i-j) tp+=a[j]*b[(i-j)]; } c[i]=tp; tp=0; } } //函数自相关// voidzixiangguan(float*a,float*r,intP) { inti,j; floats=0; for(i=0;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); } 第二部分: MATLAB f=25;%频率 a=2/3*f*f*log (2); dt=0.001;%采样时间 fori=1: 200 wave(i)=cos(2*pi*f*i*dt)*exp(-a*i*i*dt*dt);%生成子波 end plot(wave) fori=1: 200 ref(i)=0; end ref(50)=0.7; ref(80)=-0.4; ref(150)=0.5;%构造反射序列 con1=conv(wave,ref);%褶积生成地震记录 fori=1: 200 wav(i)=wave(201-i);%把wave反一下 end f3=conv(wave,wav);%wave与wav褶积相当于wave自相关 fori=1: 200 t(i)=f3(i+199);%利用自相关后200个数据 end T=toeplitz(t);%用t构造托普利兹矩阵,相当于线性方程组的系数 b=zeros(200,1);%构造一个矩阵,线性方程组右侧的结果 b(1,1)=1;%构造一个矩阵,相当于线性方程组右侧的结果 at=T\b;%解方程组,求出反滤波因子at yt=conv(con1,at);%地震记录与反滤波因子褶积 子波 地震记录 反射序列 子波自相关 截取自相关后半段 反滤波因子 Yt: 滤波因子at与地震记录褶积结果 与最初反射系数一致,实现反滤波
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 最小 平方 滤波 程序代码 语言 以及 MATLAB 编写