%GPS单历元信号变形监测数据去噪处理(220页例子)基于阈值化数据去噪处理
clc
load leleccum;
s=leleccum(1:3920);
s(500)=300
ls=length(s);
figure
% subplot(221);
plot(s);
title('原始形变数据序列');grid on;
xlabel('历元数')
ylabel('x方向形变量')
[c,l]=wavedec(s,3,'db1');
ca3=appcoef(c,l,'db1',3);
cd3=detcoef(c,l,3);
cd2=detcoef(c,l,2);
cd1=detcoef(c,l,1);
%下面对信号进行强制消噪处理
cdd3=zeros(1,length(cd3));
cdd2=zeros(1,length(cd2));
cdd1=zeros(1,length(cd1));
cl=[ca3,cdd3,cdd2,cdd1];
s1=waverec(cl,l,'db1');
figure
subplot(311);plot(s1);grid on;
title('强制消噪后形变数据序列');
xlabel('历元数')
ylabel('x方向形变量')
% %下面利用默认值消噪处理
% [thr,sorth,keepapp]=ddencmp('den','wv',s);
% s2=wdencmp('gbl',c,l,'db1',3,thr,sorth,keepapp);
% subplot(212);plot(s2);
% title('默认阈值消噪处理');grid
% %下面给定软阈值进行消噪处理
%确定阈值LAMTA
% N1=length(cd1)
% sigma1=std(cd1)
% lamta1=sigma1*sqrt(2*log(N1))
%
% N2=length(cd2)
% sigma2=std(cd2)
% lamta2=sigma2*sqrt(2*log(N2))
%
% N3=length(cd3)
% sigma3=std(cd3)
% lamta3=sigma3*sqrt(2*log(N3))
N1=length(cd1)
sigma1=median(cd1)/0.6754
lamta1=sigma1*sqrt(2*log(N1)/N1)
N2=length(cd2)
sigma2=median(cd2)/0.6754
lamta2=sigma2*sqrt(2*log(N2)/N2)
N3=length(cd3)
sigma3=median(cd3)/0.6754
lamta3=sigma3*sqrt(2*log(N3)/N3)
cd1soft=mydefine(cd1,'s',lamta1);
cd2soft=mydefine(cd2,'s',lamta2);
cd3soft=mydefine(cd3,'s',lamta3);
% cd1soft=mydefine(cd1,'s',1.456) ;
% cd2soft=mydefine(cd2,'s',1.832);
% cd3soft=mydefine(cd3,'s',2.786);
c2=[ca3,cd3soft,cd2soft,cd1soft];
s3=waverec(c2,l,'db1');
subplot(312);plot(s3);
title('软阈值进行消噪处理');grid on;
xlabel('历元数')
ylabel('x方向形变量')
%基于小波系数的中值滤波去噪
cd1m=medfilt1(cd1,9)
cd2m=medfilt1(cd2,9)
cd3m=medfilt1(cd3,9)
c2=[ca3,cd3m,cd2m,cd1m];
s3=waverec(c2,l,'db1');grid on
subplot(313);plot(s3);
xlabel('历元数')
ylabel('x方向形变量')
title('小波系数中值滤波消噪处理');grid;