clc;
clear;
close all;
datafile='';
x=load(datafile);
subplot(4,2,1);
plot(x);
title('原始信号');
init=2055615866;
randn('seed',init);
nx=x+18*randn(size(x));
subplot(4,2,2);
plot(nx);
title('含噪声信号');
[c,l]=wavedec(nx,5,'db3');
a5=appcoef(c,l,'db3',5);
d5=detcoef(c,l,5);
d4=detcoef(c,l,4);
d3=detcoef(c,l,3);
d2=detcoef(c,l,2);
d1=detcoef(c,l,1);
d5=detcoef(c,l,5);
N=1024;
cD=[d1,d2,d3,d4,d5];
sigma=median(abs(cD))/0.6745;
thr1=(sigma*sqrt(2*(log10(N))))/(log10(2));
cD1=wthresh(d1,'h',thr1);
thr2=(sigma*sqrt(2*(log10(N))))/(log10(3));
cD2=wthresh(d2,'h',thr2);
thr3=(sigma*sqrt(2*(log10(N))))/(log10(4));
cD3=wthresh(d3,'h',thr3);
thr4=(sigma*sqrt(2*(log10(N))))/(log10(5));
cD4=wthresh(d4,'h',thr4);
thr5=(sigma*sqrt(2*(log10(N))))/(log10(6));
cD5=wthresh(d5,'h',thr5);
cd=[a5,cD5,cD4,cD3,cD2,cD1];
c=cd;
y=waverec(c,l,'db3');
subplot(4,2,3);
plot(y);
title('硬阈值去噪');
thr1=(sigma*sqrt(2*(log10(N))))/(log10(2));
cD1=wthresh(d1,'s',thr1);
thr2=(sigma*sqrt(2*(log10(N))))/(log10(3));
cD2=wthresh(d2,'s',thr2);
thr3=(sigma*sqrt(2*(log10(N))))/(log10(4));
cD3=wthresh(d3,'s',thr3);
thr4=(sigma*sqrt(2*(log10(N))))/(log10(5));
cD4=wthresh(d4,'s',thr4);
thr5=(sigma*sqrt(2*(log10(N))))/(log10(6));
cD5=wthresh(d5,'s',thr5);
cd=[a5,cD5,cD4,cD3,cD2,cD1];
c=cd;
ys=waverec(c,l,'db3');
subplot(4,2,4);
plot(ys);
title('软阈值去噪');
a1=0.7;
thr1=(sigma*sqrt(2*(log10(length(d1)))))/(log10(2));
for i=1:length(d1)
if(abs(d1(i))<=thr1)
cD1(i)=sign(d1(i))*(a1)*(abs(d1(i)));
else
cD1(i)=d1(i);
end
end
thr2=(sigma*sqrt(2*(log10(length(d2)))))/(log10(3));
for i=1:length(d2)
if(abs(d2(i))<=thr2)
cD2(i)=sign(d2(i))*(a1)*(abs(d2(i)));
else
cD2(i)=d2(i);
end
end
thr3=(sigma*sqrt(2*(log10(length(d3)))))/(log10(4));
for i=1:length(d3)
if(abs(d3(i))<=thr3)
cD3(i)=sign(d3(i))*(a1)*(abs(d3(i)));
else
cD3(i)=d3(i);
end
end
thr4=(sigma*sqrt(2*(log10(length(d4)))))/(log10(5));
for i=1:length(d4)
if(abs(d4(i))<=thr4)
cD4(i)=sign(d4(i))*(a1)*(abs(d4(i)));
else
cD4(i)=d4(i);
end
end
thr5=(sigma*sqrt(2*(log10(length(d5)))))/(log10(6));
for i=1:length(d5)
if(abs(d5(i))<=thr5)
cD5(i)=sign(d5(i))*(a1)*(abs(d5(i)));
else
cD5(i)=d5(i);
end
end
cd=[a5,cD5,cD4,cD3,cD2,cD1];
c=cd;
ysh1=waverec(cd,l,'db3');
subplot(4,2,5)
plot(ysh1);
title('改进的阈值去噪(a=0.7)');
a2=0.04;
thr1=(sigma*sqrt(2*(log10(length(d1)))))/(log10(2));
for i=1:length(d1)
if(abs(d1(i))<=thr1)
cD1(i)=sign(d1(i))*(a2)*(abs(d1(i)));
else
cD1(i)=d1(i);
end
end
thr2=(sigma*sqrt(2*(log10(length(d2)))))/(log10(3));
for i=1:length(d2)
if(abs(d2(i))<=thr2)
cD2(i)=sign(d2(i))*(a2)*(abs(d2(i)));
else
cD2(i)=d2(i);
end
end
thr3=(sigma*sqrt(2*(log10(length(d3)))))/(log10(4));
for i=1:length(d3)
if(abs(d3(i))<=thr3)
cD3(i)=sign(d3(i))*(a2)*(abs(d3(i)));
else
cD3(i)=d3(i);
end
end
thr4=(sigma*sqrt(2*(log10(length(d4)))))/(log10(5));
for i=1:length(d4)
if(abs(d4(i))<=thr4)
cD4(i)=sign(d4(i))*(a2)*(abs(d4(i)));
else
cD4(i)=d4(i);
end
end
thr5=(sigma*sqrt(2*(log10(length(d5)))))/(log10(6));
for i=1:length(d5)
if(abs(d5(i))<=thr5)
cD5(i)=sign(d5(i))*(a2)*(abs(d5(i)));
else
cD5(i)=d5(i);
end
end
cd=[a5,cD5,cD4,cD3,cD2,cD1];
c=cd;
ysh2=waverec(cd,l,'db3');
subplot(4,2,6)
plot(ysh2);
title('改进的阈值去噪(a=0.04)');
a3=0.007;
thr1=(sigma*sqrt(2*(log10(length(d1)))))/(log10(2));
for i=1:length(d1)
if(abs(d1(i))<=thr1)
cD1(i)=sign(d1(i))*(a3)*(abs(d1(i)));
else
cD1(i)=d1(i);
end
end
thr2=(sigma*sqrt(2*(log10(length(d2)))))/(log10(3));
for i=1:length(d2)
if(abs(d2(i))<=thr2)
cD2(i)=sign(d2(i))*(a3)*(abs(d2(i)));
else
cD2(i)=d2(i);
end
end
thr3=(sigma*sqrt(2*(log10(length(d3)))))/(log10(4));
for i=1:length(d3)
if(abs(d3(i))<=thr3)
cD3(i)=sign(d3(i))*(a3)*(abs(d3(i)));
else
cD3(i)=d3(i);
end
end
thr4=(sigma*sqrt(2*(log10(length(d4)))))/(log10(5));
for i=1:length(d4)
if(abs(d4(i))<=thr4)
cD4(i)=sign(d4(i))*(a3)*(abs(d4(i)));
else
cD4(i)=d4(i);
end
end
thr5=(sigma*sqrt(2*(log10(length(d5)))))/(log10(6));
for i=1:length(d5)
if(abs(d5(i))<=thr5)
cD5(i)=sign(d5(i))*(a3)*(abs(d5(i)));
else
cD5(i)=d5(i);
end
end
cd=[a5,cD5,cD4,cD3,cD2,cD1];
c=cd;
ysh3=waverec(cd,l,'db3');
subplot(4,2,7)
plot(ysh3);
title('改进的阈值去噪(a=0.007)');
a4=0.0001;
thr1=(sigma*sqrt(2*(log10(length(d1)))))/(log10(2));
for i=1:length(d1)
if(abs(d1(i))<=thr1)
cD1(i)=sign(d1(i))*(a4)*(abs(d1(i)));
else
cD1(i)=d1(i);
end
end
thr2=(sigma*sqrt(2*(log10(length(d2)))))/(log10(3));
for i=1:length(d2)
if(abs(d2(i))<=thr2)
cD2(i)=sign(d2(i))*(a4)*(abs(d2(i)));
else
cD2(i)=d2(i);
end
end
thr3=(sigma*sqrt(2*(log10(length(d3)))))/(log10(4));
for i=1:length(d3)
if(abs(d3(i))<=thr3)
cD3(i)=sign(d3(i))*(a4)*(abs(d3(i)));
else
cD3(i)=d3(i);
end
end
thr4=(sigma*sqrt(2*(log10(length(d4)))))/(log10(5));
for i=1:length(d4)
if(abs(d4(i))<=thr4)
cD4(i)=sign(d4(i))*(a4)*(abs(d4(i)));
else
cD4(i)=d4(i);
end
end
thr5=(sigma*sqrt(2*(log10(length(d5)))))/(log10(6));
for i=1:length(d5)
if(abs(d5(i))<=thr5)
cD5(i)=sign(d5(i))*(a4)*(abs(d5(i)));
else
cD5(i)=d5(i);
end
end
cd=[a5,cD5,cD4,cD3,cD2,cD1];
c=cd;
ysh4=waverec(cd,l,'db3');
subplot(4,2,8)
plot(ysh4);
title('改进的阈值去噪(a=0.0001)');
y1=sum(nx.^2);
yxn=sum((nx-x).^2);
y2=sum(ys.^2);
yys=sum((ys-x).^2);
y3=sum(y.^2);
yyh=sum((y-x).^2);
y41=sum(ysh1.^2);
yyhs1=sum((ysh1-x).^2);
y42=sum(ysh2.^2);
yyhs2=sum((ysh2-x).^2);
y43=sum(ysh3.^2);
yyhs3=sum((ysh3-x).^2);
y44=sum(ysh4.^2);
yyhs4=sum((ysh4-x).^2);
snrxn=10*log10((y1/yxn));
fprintf('原加噪信号信噪比%4.4f\n',snrxn);
snryh=10*log10((y2/yyh));
fprintf('硬阈值去噪信号信噪比%4.4f\n',snryh);
snrys=10*log10((y3/yys));
fprintf('软阈值去噪信号信噪比%4.4f\n',snrys);
snrys=10*log10((y41/yyhs1));
fprintf('(a=0.7)改进的阈值处理后的信噪比%4.4f\n',snrys);
snrys=10*log10((y42/yyhs2));
fprintf('(a=0.04)改进的阈值处理后的信噪比%4.4f\n',snrys);
snrys=10*log10((y43/yyhs3));
fprintf('(a=0.007)改进的阈值处理后的信噪比%4.4f\n',snrys);
snrys=10*log10((y44/yyhs4));
fprintf('(a=0.0001)改进的阈值处理后的信噪比%4.4f\n',snrys);
REMxn=sqrt(yxn/N);
fprintf('原加噪信号均方差%4.4f\n',REMxn);
REMys=sqrt(yys/N);
fprintf('软阈值去噪信号均方差%4.4f\n',REMys);
REMyh=sqrt(yyh/N);
fprintf('硬阈值去噪信号均方差%4.4f\n',REMyh);
REMys=sqrt(yyhs1/N);
fprintf('(a=0.7)改进软硬阈值处理后均方根误差%4.4f\n',REMys);
REMys=sqrt(yyhs2/N);
fprintf('(a=0.04)改进软硬阈值处理后均方根误差%4.4f\n',REMys);
REMys=sqrt(yyhs3/N);
fprintf('(a=0.007)改进软硬阈值处理后均方根�
