最小平方反濾波程序代碼C語(yǔ)言以及MATLAB編寫(xiě)

1、精選優(yōu)質(zhì)文檔-傾情為你奉上精選優(yōu)質(zhì)文檔-傾情為你奉上專心-專注-專業(yè)專心-專注-專業(yè)精選優(yōu)質(zhì)文檔-傾情為你奉上專心-專注-專業(yè)分C語(yǔ)言和MATLAB兩部分第一部分：C語(yǔ)言#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,floa

2、t *r,int P);/自相關(guān)int jiefangchengzu(float R,int n,float b,float x);/解托普利茨方程組void main() int i,j;int P=200;/子波長(zhǎng)度int Q=399;/反射序列長(zhǎng)度f(wàn)loat 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);/

3、子波與反射序列褶積的地震記錄fp_xcorr=fopen(xcorr.csv,w);/地震記錄自相關(guān)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; /構(gòu)造反射序列/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); / 構(gòu)造子波 /con(wave,ref,c

4、on1,P,P);for(i=0;iQ;i+)fprintf(fp_con1,%fn,con1i); / 褶積生成地震記錄 /zixiangguan(wave,xcorr,200);for(i=0;iP;i+)fprintf(fp_xcorr,%fn,xcorri); / 子波自相關(guān) /for(i=0;iP;i+)bi=0;b0=1; / 構(gòu)造方程組右側(cè)結(jié)果數(shù)組 /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

5、;i+)fprintf(fp_yt,%fn,yti); / 因子與地震記錄褶積 / 普通褶積 /void 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;/ 函數(shù)自相關(guān) /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;/從最左

6、邊開(kāi)始，移到兩者重合rP-1-i=s;s=0;/ 萊文森遞推解托普利茨方程組 / /R為矩陣的第一行或者第一列數(shù)據(jù)，b為方程組右側(cè)結(jié)果， x為計(jì)算的反濾波因子 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); ret

7、urn(-1); y0=1.0;x0=b0/a; for(k=1;k=n-1;k+) beta=0.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(f

8、ailn);return(-1); h=(bk-q)/a; for(i=0;i=k-1;i+) xi=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;%采樣時(shí)間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;%構(gòu)造反射序列con1=conv(wave,ref);%褶積生成地震記錄for i=1:200wav(i)=wave(201-i);%把wave反一下endf3=conv(wave,wav);%wave與wav褶積相當(dāng)于wave自相關(guān)for i=1:200t(i)=f3(i+199);%利用自相關(guān)后200個(gè)數(shù)據(jù)endT=toeplitz(t);% 用t構(gòu)造托普利茲矩陣，相當(dāng)于線性方程組的系數(shù)b=zeros(200,1);%構(gòu)造一個(gè)矩陣，線性方程組右側(cè)的結(jié)果b(1,1)