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

A2.1 卡尔曼滤波程序： load xinshuju.mat %调入数据 h=; 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); end t=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=; 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)) figure plot(A0,r);title(the recontrust signal); %%%%%%%%%%%%%%%%%% [C,L] = wavedec(s,3,db8) % 三次小波分解 cA3 = appcoef(C,L,db8,3);[cD1,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); figure subplot(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)) figure plot(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); figure subplot(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;%给