dsfire 发表于 2010-3-10 20:47

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

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

close all;
clc;
clear;
format short g;
Fs=50*128;
N=1024;
n=;
A= ;
Ph=;
ff=;
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层分解
=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=;
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=;
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=;
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=;
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 发表于 2010-3-10 21:01

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

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

aprilcat 发表于 2010-3-31 22:27

回复 楼主 dsfire 的帖子

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

dsfire 发表于 2010-4-6 16:23

回复 板凳 aprilcat 的帖子

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

aprilcat 发表于 2010-4-8 16:53

回复 地板 dsfire 的帖子

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

无锡老板 发表于 2011-3-7 23:01

回复 2 # dsfire 的帖子

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

forward39 发表于 2011-3-10 23:14

期待讨论

SMALLTIGER 发表于 2011-3-12 15:39

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

dsfire 发表于 2011-3-13 19:47

回复 6 # 无锡老板 的帖子

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

zhangshun5233 发表于 2011-8-5 15:15

本帖最后由 zhangshun5233 于 2011-8-5 15:16 编辑

回复 2 # dsfire 的帖子

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

dsfire 发表于 2011-8-7 17:07

回复 10 # zhangshun5233 的帖子

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

zhangshun5233 发表于 2011-8-9 20:21

回复 11 # dsfire 的帖子

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

dsfire 发表于 2011-8-9 20:34

回复 12 # zhangshun5233 的帖子

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

piscesedwin 发表于 2011-11-14 21:54

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

kongming 发表于 2011-11-23 20:05

本帖最后由 kongming 于 2011-11-23 20:11 编辑

请教楼主,为什么在软、硬阈值法时用thr1=(sigma*sqrt(2*(log10(N))))/(log10(2)),而用新阈值法时用thr1=(sigma*sqrt(2*(log10(length(d1)))))/(log10(1+1)),就是为什么要将N变成length(d1)?
页: [1] 2
查看完整版本: 改进软硬阈值折衷法的小波去噪程序?求改正