最小平方反滤波程序代码C语言以及MATLAB编写Word文件下载.docx
《最小平方反滤波程序代码C语言以及MATLAB编写Word文件下载.docx》由会员分享，可在线阅读，更多相关《最小平方反滤波程序代码C语言以及MATLAB编写Word文件下载.docx（12页珍藏版）》请在冰豆网上搜索。
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"
//子波与反射序列褶积的地震记录
fp_xcorr=fopen（"
xcorr.csv"
//地震记录自相关
fp_at=fopen（"
at.csv"
//反滤波因子
fp_yt=fopen（"
yt.csv"
//因子与地震记录褶积
for（i=0;
i<
P;
i++）
ref[i]=0.;
ref[50]=0.2;
ref[80]=0.4;
ref[150]=-0.7;
//构造反射序列//
{
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）;
Q;
fprintf（fp_con1,"
con1[i]）;
//褶积生成地震记录//
zixiangguan（wave,xcorr,200）;
fprintf（fp_xcorr,"
xcorr[i]）;
//子波自相关//
b[i]=0;
b[0]=1;
//构造方程组右侧结果数组//
jiefangchengzu（xcorr,P,b,at）;
fprintf（fp_at,"
at[i]）;
//解方程组，求反滤波因子//
}
con（con1,at,yt,Q,P）;
P+Q-2;
fprintf（fp_yt,"
yt[i]）;
//因子与地震记录褶积//
}
//普通褶积//
voidcon（floata[],floatb[],floatc[],intM,intL）
{
inti,j,N;
floattp=0;
N=M+L-1;
N;
i++）
for（j=0;
j<
M;
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）
floats=0;
=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;
=k-1;
beta=beta+R[j+1]*y[j];
q=q+R[k-j]*x[j];
if（fabs（a）+1.0==1.0）
free（y）;
printf（"
return（-1）;
c=-beta/a;
s[0]=c*y[k-1];
y[k]=y[k-1];
if（k!
=1）
for（i=1;
{s[i]=y[i-1]+c*y[k-i-1];
a=a+c*beta;
h=（b[k]-q）/a;
{x[i]=x[i]+h*s[i];
y[i]=s[i];
x[k]=h*y[k];
free（s）;
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）
ref（i）=0;
end
ref（50）=0.7;
ref（80）=-0.4;
ref（150）=0.5;
%构造反射序列
con1=conv（wave,ref）;
%褶积生成地震记录
wav（i）=wave（201-i）;
%把wave反一下
f3=conv（wave,wav）;
%wave与wav褶积相当于wave自相关
t（i）=f3（i+199）;
%利用自相关后200个数据
T=toeplitz（t）;
%用t构造托普利兹矩阵，相当于线性方程组的系数
b=zeros（200,1）;
%构造一个矩阵，线性方程组右侧的结果
b（1,1）=1;
%构造一个矩阵，相当于线性方程组右侧的结果
at=T\b;
%解方程组，求出反滤波因子at
yt=conv（con1,at）;
%地震记录与反滤波因子褶积
子波
地震记录
反射序列
子波自相关
截取自相关后半段
反滤波因子
Yt：
滤波因子at与地震记录褶积结果
与最初反射系数一致，实现反滤波
升级会员 