卡尔曼滤波和小波滤波程序.doc

吴国曾：石油钻井事故预报中数据的滤波方法研究A2.1 卡尔曼滤波程序：load xinshuju.mat %调入数据h=xinshuju;t1=14000;t2=15000;t=t1:t2;n=t2-t1+1;s=h(t);for i=1:n-1; b(i)=s(i+1)/s(i);end % a=mean(b);a=1;w(1)=0;w=randn(1,n);x(1)=s(1);for k=2:n x(k)=b(i)*x(k-1)+w(k-1);end v=randn(1,n);q1=std(v);rvv=q1.2;q2=std(w);rww=q2.2;c=1;z=c*s+v;p(1)=5;m(1)=0;for t=2:n %卡尔曼滤波过程 p1(t)=a.2*p(t-1)+rww; g(t)=c*p1(t)/(p1(t)+rvv); m(t)=a*m(t-1)+g(t)*(z(t)-a*c*m(t-1); p(t)= p1(t)-c*g(t)*p1(t);endt=1:n;figure subplot(211);plot(s);title(the orignal data); subplot(212);plot(m,g);title(the de-noising data) figure plot(s); hold on plot(m,r);A2.2 小波滤波：load xinshuju.mat %调入数据h=xinshuju;t1=1;t2=15000;t=t1:t2;n=t2-t1+1;s=h(t);l_s = length(s);cA1,cD1 = dwt(s,db8);%用db8小波对信号s进行单层分解a1=upcoef(a,cA1,db8,1,l_s); %重建小波系数至上一层a代表低通，d表高通，1表a作用在 d1 = upcoef(d,cD1,db8,1,l_s);% ca1上的次数subplot(311); plot(a1); title(Approximation A1)subplot(312); plot(d1); title(Detail D1) subplot(313);plot(s,g);A0 = idwt(cA1,cD1,db8,l_s);%用于单层小波重建err = max(abs(s-A0)figureplot(A0,r);title(the recontrust signal);%C,L = wavedec(s,3,db8) % 三次小波分解cA3 = appcoef(C,L,db8,3);%求得第三层的近似系数，长度60cD1,cD2,cD3 = detcoef(C,L,1,2,3);%求得1到3层的细节系数，长度依次减小A3 = wrcoef(a,C,L,db8,3);%重建小波系数，长度和源信号一样D1 = wrcoef(d,C,L,db8,1);D2 = wrcoef(d,C,L,db8,2);D3 = wrcoef(d,C,L,db8,3);figuresubplot(411); plot(A3); title(Approximation A3)subplot(412); plot(D1); title(Detail D1)subplot(413); plot(D2); title(Detail D2)subplot(414); plot(D3); title(Detail D3)%A0 = waverec(C,L,db8);%多层小波重建原始信号err1 = max(abs(s-A0)figureplot(A0,g);title(from the 3-dem-recontrust signal);% Remove noise by thresholding. thr,sorh,keepapp = ddencmp(den,wv,s);%den用于降噪，cmp用于压缩，wv使用小波变换，%wp使用小波包变换，thr求得的阈值，sorh，s%表软阈值，h表硬阈值，keepapp表保留的近似系数的层数。clean = wdencmp(gbl,s,db1,3,thr,sorh,keepapp);figuresubplot(211);plot(s);title(the orignal signal);subplot(212);plot(clean,g);title(Remove noise by thresholding);A2.3 MA 滤波：Max_m=59;Min_m=19;%给定启动新滤波段的阈值load EB_Da24.mat %调入数据h=EB_Da24;t1=20000;t2=20349;t=t1:t2;n=t2-t1+1;x=h(t);len=length(x);num1=10;for ki=1:num1 y(ki)=x(ki);endfor i=num1:length(x) f1=(abs(x(i)-x(i-1)Max_m); f2=(abs(x(i)-x(i-2)Max_m); if f1|f2 y(i)=x(i); y1(i)=x(i); else y(i)=0.10*x(i)+0.15*y(i-1)+0.75*y(i-2); y1(i)=mean(x(i-num1+1:i); end endt0=1:len; N=4; cA1,cD1=dwt(x,db8); A1=upcoef(a,cA1,db8,1,len); D1=upcoef(d,cD1,db8,1,len); C,L=wavedec(x,5,db8); %重构各层逼近系数 cA5=appcoef(C,L,db8,5);%提取一维近似系数 A5=wrcoef(a,C,L,db8,5);%由一维小波系数进行重构 for k=1:len dif1(k)=abs(x(k)-A5(k); dif2(k)=abs(x(k)-y(k); dif3(k)=abs(x(k)-y1(k); end figure(1);subplot(311);plot(t0,x(t0)title(the orignal signal)subplot(312);plot(t0,A5(t0)title(the processed signal)subplot(313);plot(t0,x(t0),r,t0,A5(t0),g);title(compare the two signal)figure(2);plot(t0,x(t0),t0,y(t0),r,