改进软硬阈值折衷法的小波去噪程序？求改正

dsfire

仿照着《小波去噪中软硬阈值的一种改良折衷法》一文写的程序，对比了软、硬以及软硬阈值折衷法的信噪比及均方根误差。仿真出来之后，原信号的信噪比是负的，不知道为什么？均方根误差的值也比较大，不知何解？信噪比和均方根误差是参照《基于matlab的小波去噪仿真》：吴伟，《信息与电子工程》 一文写的。
不知道此程序有没有问题，望高手给予指正！！

close all;
clc;
clear;
format short g;
Fs=50*128;
N=1024;
n=[0:1:N-1];
A= [220 2   15  9];
Ph=[20  60  45  80 ];
ff=[50  100 150 250];
Num=length(A);
y=zeros(1,length(n));
for i=1:Num
    myt=A(i)*cos(2*pi*ff(i)*n./Fs+Ph(i)*pi/180);
    y=y+myt;%%构造原始信号
end
figure(1);
subplot(511);
plot(y);
ylabel('幅值 A');
title('原始信号');
s=awgn(y,20,'measured');%加入高斯白噪声
subplot(512);
plot(s);
ylabel('幅值 A');
title('加噪信号');
wname='db3';%选db3小波基
lev=5;%5层分解
[c,l]=wavedec(s,lev,wname);
a5=appcoef(c,l,wname,lev);
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);
a=8500;%%%
b=13;%%%%a和b两值是根据《小波去噪中软硬阈值的一种改良折衷法》一文给出的
cD=[d1,d2,d3,d4,d5];
sigma=median(abs(cD))/0.6475;

%-----------软阈值处理--------------------------
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,wname);
subplot(513)
plot(ys);
title('软阈值处理');

%---------硬阈值处理-------------------------
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;
yh=waverec(c,l,wname);
subplot(514)
plot(yh);
title('硬阈值处理');
%%%%根据《小波去噪中软硬阈值的一种改良折衷法》一文来写程序 作者为：郭晓霞，杨惠中《智能系统学报》
thr1=(sigma*sqrt(2*(log10(length(d1)))))/(log10(1+1));
for i=1:length(d1)
    if(abs(d1(i))>=thr1)
        cD1(i)=sign(d1(i))*(abs(d1(i))-b*thr1/(a^(abs(abs(d1(i))-thr1))+b-1));%估计第一层小波系数
    else
        cD1(i)=0;
    end
end
thr2=(sigma*sqrt(2*(log10(length(d2)))))/(log10(2+1));
for i=1:length(d2)
    if(abs(d2(i))>=thr2)
        cD2(i)=sign(d2(i))*(abs(d2(i))-b*thr2/(a^(abs(abs(d2(i))-thr2))+b-1));%估计第二层小波系数
    else
        cD2(i)=0;
    end
end
thr3=(sigma*sqrt(2*(log10(length(d3)))))/(log10(3+1));
for i=1:length(d3)
    if(abs(d3(i))>=thr3)
        cD3(i)=sign(d3(i))*(abs(d3(i))-b*thr3/(a^(abs(abs(d3(i))-thr3))+b-1));%估计第三层小波系数
    else
        cD3(i)=0;
    end
end
thr4=(sigma*sqrt(2*(log10(length(d4)))))/(log10(4+1));
for i=1:length(d4)
    if(abs(d4(i))>=thr4)
        cD4(i)=sign(d4(i))*(abs(d4(i))-b*thr4/(a^(abs(abs(d4(i))-thr4))+b-1));%估计第四层小波系数
    else
        cD4(i)=0;
    end
end
thr5=(sigma*sqrt(2*(log10(length(d5)))))/(log10(5+1));
for i=1:length(d5)
    if(abs(d5(i))>=thr5)
        cD5(i)=sign(d5(i))*(abs(d5(i))-b*thr5/(a^(abs(abs(d5(i))-thr5))+b-1));%估计第五层小波系数
    else
        cD5(i)=0;
    end
end
%%%%开始重构
cd=[a5,cD5,cD4,cD3,cD2,cD1];
c=cd;
yhs=waverec(cd,l,wname);
subplot(515);
plot(ys,'LineWidth',1);
ylabel('幅值 A')
title('软硬阈值折衷去噪');
y1=sum(s.^2);
yxn=sum((s-y).^2);%
yys=sum((s-ys).^2);
yyh=sum((s-yh).^2);
yyhs=sum((s-yhs).^2);
snrxn=10*log10((y1/yxn));%原加噪信号信噪比
fprintf('原加噪信号信噪比%4.4f\n',snrxn);
snrys=10*log10((y1/yys));%软阈值去噪信号信噪比
fprintf('软阈值去噪信号信噪比%4.4f\n',snrys);
snryh=10*log10((y1/yyh));%硬阈值去噪信号信噪比
fprintf('硬阈值去噪信号信噪比%4.4f\n',snryh);
snrys=10*log10((y1/yyhs));%软硬阈值处理后的信噪比
fprintf('软硬阈值处理后的信噪比%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(yyhs/N);%软硬阈值处理后均方根误差
fprintf('软硬阈值处理后均方根误差%4.4f\n',REMys);

[ 本帖最后由 dsfire 于 2010-3-10 21:12 编辑 ]

dsfire

运行的结果也贴上去。
原加噪信号信噪比20.0155
软阈值去噪信号信噪比17.8229
硬阈值去噪信号信噪比19.6739
软硬阈值处理后的信噪比20.1870
原加噪信号均方根误差15.6866
软阈值去噪信号均方根误差20.1909
硬阈值去噪信号均方根误差16.3157
软硬阈值处理后均方根误差15.3798
均方根误差比原加噪的均方根误差都大！

[ 本帖最后由 dsfire 于 2010-3-10 21:14 编辑 ]

aprilcat

在软硬阈值折中的方法中，sigma=median(abs(cD))/0.6475;为什么不是sigma=median(abs(di))/0.6475呢？

dsfire

首先，非常感谢您的回复！关于是di还是cD，这个我也不太清楚！
根据《小波分析及其工程应用》-杨建国一书中第110页的表述为：sigma应为2^1尺度上的细节小波系数的绝对值的中值。2^1尺度不知道该怎么理解？

aprilcat

我也是才开始学习小波的，不懂，共同学习，期待高手来指点下。

无锡老板

回复 2 # dsfire 的帖子

我觉得应该是di，每层都应该不同的呀。
我想请问一下，譬如说thr1=(sigma*sqrt(2*(log10(N))))/(log10(2));
为什么最后还要除一个log10(i+1)?
不应该是sigma*sqrt(2*log10(N))吗？
谢谢！

forward39

期待讨论

SMALLTIGER

从结果看，降噪效果不是很明显啊！之所以是log10(i+1)是因为信号的小波变换系数随尺度的增大而增大，而高斯白噪声的却随尺度增大而减小。欢迎继续讨论。

dsfire

回复 6 # 无锡老板 的帖子

这个程序有一年啦，你把上面提到的文章看看，说不定就能知道为什么了！

zhangshun5233

回复 2 # dsfire 的帖子

原加噪信号信噪比20.0155
软阈值去噪信号信噪比17.8229
硬阈值去噪信号信噪比19.6739
软硬阈值去噪的信噪比都小于原加噪信号的信噪比，这不正确吧？
如果是这样的话，那岂不是去噪以后的噪声更大了？
应该是去噪以后的信噪比越大于原加噪信号的信噪比，去噪效果越好；你觉得呢？

dsfire

回复 10 # zhangshun5233 的帖子

你说的应该是对的，但是根据算法和程序的计算就得出了上面的结果！求教一二啊

zhangshun5233

回复 11 # dsfire 的帖子

可能要等九月份再讨论这个问题了~~！这个问题对你已经过期了么？

dsfire

回复 12 # zhangshun5233 的帖子

不能说过期吧，也有很长一段时间没看了。

piscesedwin

在你的软阈值处理中，各层细节系数用的都是同一个阈值，这样不对吧，各尺度上的sigma应该是不同的

kongming

请教楼主，为什么在软、硬阈值法时用thr1=(sigma*sqrt(2*(log10(N))))/(log10(2))，而用新阈值法时用thr1=(sigma*sqrt(2*(log10(length(d1)))))/(log10(1+1))，就是为什么要将N变成length（d1）？