clear all
data=load('e:\ceshi.txt');
size(data);
save data;
ceshidata=data';
save ceshidata;
load ceshidata;
s=ceshidata(1:563);
ls=length(s);
[C,L]=wavedec(s,5,'db1');
cD1=detcoef(C,L,1);
cA1=appcoef(C,L,'db1',1);
D1=wrcoef('d',C,L,'db1',1);
D2=wrcoef('d',C,L,'db1',2);
D3=wrcoef('d',C,L,'db1',3);
D4=wrcoef('d',C,L,'db1',4);
D5=wrcoef('d',C,L,'db1',5);
A=wrcoef('a',C,L,'db1',1);
D=0*D1;
D1(1,564)=0;
D2(1,564)=0;
D3(1,564)=0;
D4(1,564)=0;
D5(1,564)=0;
for i=1:563
if abs(D1(1,i+1))>abs(D1(1,i))
if abs(D1(1,i+1))>abs(D1(1,i-1))
D(1,i+1)=D1(1,i+1);
Dn1=D;
A1=max(Dn1(1,:))/2;
T1=wthresh(Dn1(1,:),'s',A1);
end
end
end
for i=1:563
if abs(D2(1,i+1))>abs(D2(1,i))
if abs(D2(1,i+1))>abs(D2(1,i-1))
D(1,i+1)=D2(1,i+1);
Dn2=D;
A2=max(Dn2(1,:))/2;
T2=wthresh(Dn2(1,:),'s',A2);
end
end
end
for i=1:563
if abs(D3(1,i+1))>abs(D3(1,i))
if abs(D3(1,i+1))>abs(D3(1,i-1))
D(1,i+1)=D3(1,i+1);
Dn3=D;
A3=max(Dn3(1,:))/2;
T3=wthresh(Dn3(1,:),'s',A3);
end
end
end
for i=1:563
if abs(D4(1,i+1))>abs(D4(1,i))
if abs(D4(1,i+1))>abs(D4(1,i-1))
D(1,i+1)=D4(1,i+1);
Dn4=D;
A4=max(Dn4(1,:))/2;
T4=wthresh(Dn4(1,:),'s',A4);
end
end
end
for i=1:563
if abs(D5(1,i+1))>abs(D5(1,i))
if abs(D5(1,i+1))>abs(D5(1,i-1))
D(1,i+1)=D5(1,i+1);
Dn5=D;
A5=max(Dn5(1,:))/2;
T5=wthresh(Dn5(1,:),'s',A5);
end
end
end
p=T1+T2+T3+T4+T5+A;
mdata=mean(p(1,1:68));
canchadata=p(1,1:68)-mdata;
sizedata=size(p);
d=canchadata.^2/67;
d1=sum(d,2);
deta=sqrt(d1);
for i=1:68
if abs(p(1,i+1)-mdata)>2*deta
display('可疑值为：')
p(1,i+1)
p(1,i+1)=p(1,i)
end
end
mdata=mean(p(1,68:85));
canchadata=p(1,68:85)-mdata;
sizedata=size(p);
d=canchadata.^2/(85-68-1);
d1=sum(d,2);
deta=sqrt(d1);
for i=68:85
if abs(p(1,i+1)-mdata)>2*deta
display('可疑值为：')
p(1,i+1)
p(1,i+1)=p(1,i)
end
end
mdata=mean(p(1,85:135));
canchadata=p(1,85:135)-mdata;
sizedata=size(p);
d=canchadata.^2/(135-85-1);
d1=sum(d,2);
deta=sqrt(d1);
for i=85:135
if abs(p(1,i+1)-mdata)>2*deta
display('可疑值为：')
p(1,i+1)
p(1,i+1)=p(1,i)
end
end
mdata=mean(p(1,135:563));
canchadata=p(1,135:563)-mdata;
sizedata=size(p);
d=canchadata.^2/(563-135-1);
d1=sum(d,2);
deta=sqrt(d1);
p(1,564)=0;
for i=135:563
if abs(p(1,i+1)-mdata)>2*deta
display('可疑值为：')
p(1,i+1)
p(1,i+1)=p(1,i)
end
end
subplot(2,1,1);plot(ceshidata);title('原始信号');
subplot(2,1,2);plot(p);title('去噪信号');