[已答复] 请教:matlab小波阈值去噪,阈值为什么取系数最大绝对值?

[复制链接]
Adel 发表于 2020-3-25 12:14:06
本帖最后由 Adel 于 2020-4-22 17:18 编辑

信号见附件。使用matlab自带的小波阈值去噪函数去噪,阈值计算函数thselect和阈值处理函数wthresh,结果阈值计算都很大,刚好是每层系数(细节)的最大绝对值,这样在下步阈值处理时,也就将所有系数处理成了0.想请教为什么会这样,看网上别人有些用同样方法并不会这样,而是得到很好的去噪效果。
而且,我换了信号,也能得到去噪效果,难得是信号不同么???

1、我的程序:

clear all; clc; close all;
fs=16000;                  %语音信号采样频率为16000
x1 = load('x1.mat');
x1 = x1.x1;
t = 40e-3;  % 40ms
x1 = x1(1:fs*t);
t_tick = 0:1/fs:t-1/fs;
y1=fft(x1,1024);           %对信号做1024点FFT变换
f=fs*(0:511)/1024;
figure(1)
plot(t_tick,x1)                   %做原始语音信号的时域图形
y=awgn(x1',10,'measured');   %加10db的高斯白噪声
[snr,mse]=snrmse(x1,y');%求得信噪比 均方误差
figure(2)
plot(t_tick,y)                   %做加噪语音信号的时域图形
[c,l]=wavedec(y,3,'db1');%多尺度一维分解
%用db1小波对信号进行3层分解并提取系数
a3=appcoef(c,l,'db1',3);
d3=detcoef(c,l,3);
d2=detcoef(c,l,2);
d1=detcoef(c,l,1);
%阈值获取,使用Stein的无偏风险估计原理
% 问题在这,为什么得到的阈值都是各层系数的最大绝对值呢?于是下步阈值处理时,就将所有系数都处理成0了!!!
thr1=thselect(d1,'rigrsure');
thr2=thselect(d2,'rigrsure');
thr3=thselect(d3,'rigrsure');
% 阈值处理,这里将系数全部处理成0了
gd1=wthresh(d1,'h',thr1);
gd2=wthresh(d2,'h',thr2);
gd3=wthresh(d3,'h',thr3);

c1=[a3 gd3 gd2 gd1];  % 利用新的阈值进行处理
y1=waverec(c1,l,'db1');%多尺度重构
[snr,mse]=snrmse(x1,y1');%求得信噪比 均方误差
figure(3)
plot(t_tick,y1);

下面是上述代码中的一个子函数代码,计算信噪比和MSE的:
function [snr,mse]=snrmse(I,In)
% 计算信噪比函数
% I :原始信号
% In:去噪后信号
snr=0;
Ps=sum(sum((I-mean(mean(I))).^2));%signal power
Pn=sum(sum((I-In).^2));           %noise power
snr=10*log10(Ps/Pn);
mse=Pn/length(I);
end

2、换信号,matlab中构造:


close all
clear,clc
snr=4;
init=2055615866;
%产生原始信号sref和高斯白噪声污染的信号s
[sref,s]=wnoise(1,11,snr,init);
[c,l]=wavedec(s,3,'db1');%多尺度一维分解
%用db1小波对信号进行3层分解并提取系数
a3=appcoef(c,l,'db1',3);
d3=detcoef(c,l,3);
d2=detcoef(c,l,2);
d1=detcoef(c,l,1);
thr1=thselect(d1,'rigrsure');%阈值获取,使用Stein的无偏风险估计原理
thr2=thselect(d2,'rigrsure');
thr3=thselect(d3,'rigrsure');
%利用改进阈值函数进行去噪处理。也相当于一种获取阈值,软阈值的一种
gd1=wthresh(d1,'h',thr1);
gd2=wthresh(d2,'h',thr2);
gd3=wthresh(d3,'h',thr3);
c1=[a3 gd3 gd2 gd1];  % 利用新的阈值进行处理
y1=waverec(c1,l,'db1');%多尺度重构
figure
plot(y1);
zhh =1;






x1.mat

483.23 KB, 下载次数: 26

17 条回复


Jasmine_cmOlJ 发表于 2020-4-19 10:05:28
你好,我最近也在做小波阈值去噪这部分,可以讨论下么?

Adel 发表于 2020-4-20 10:40:27
Jasmine_cmOlJ 发表于 2020-4-19 10:05
你好,我最近也在做小波阈值去噪这部分,可以讨论下么?

你好,可以啊 非常欢迎 我对基础理论也还是有很多疑问  
现在先在做C语言实现
你的问题是什么 可以发帖通知我,也可以直接在这里回

Jasmine_cmOlJ 发表于 2020-4-20 16:21:29
Adel 发表于 2020-4-20 10:40
你好,可以啊 非常欢迎 我对基础理论也还是有很多疑问  
现在先在做C语言实现
你的问题是什么 可以发帖通 ...

我找到了用matlab做小波阈值的一段代码,也运行出来了,但有些看不懂图,可以帮我看看吗,qq:3404676885

Adel 发表于 2020-4-21 11:36:59
Jasmine_cmOlJ 发表于 2020-4-20 16:21
我找到了用matlab做小波阈值的一段代码,也运行出来了,但有些看不懂图,可以帮我看看吗,qq:3404676885 ...

贴在这里撒 qq很少用 又麻烦
稍微耐心点 把自己的问题说清楚 能把问题描述的清清楚楚让别人最方便的帮助你,也是很有意义的
论坛这么好的平台 说不定有更多人帮你解决问题呢  

wx_N89FLl7S 发表于 2020-4-21 12:14:29
请问这个问题解决了吗,我的代码和您的差不多,也是每一层阈值选的是每一侧层系数的最大值。

Adel 发表于 2020-4-21 17:15:26
wx_N89FLl7S 发表于 2020-4-21 12:14
请问这个问题解决了吗,我的代码和您的差不多,也是每一层阈值选的是每一侧层系数的最大值。 ...

没呢 不知道matlab为啥算成这个样子,可能确实按照那几种规则的算法去算我们的系数,结果就是那样的
给你个建议,你去看下那几种规则是怎么算的,然后把你的系数带进去算,看下是不是就是最大绝对值?印象中还是有一些地方讲到这些准则的算法的  哈哈哈 当然我也可以去试试 但是我懒 或者暂时没去做哈:lol

Jasmine_cmOlJ 发表于 2020-4-21 17:45:08
Adel 发表于 2020-4-21 11:36
贴在这里撒 qq很少用 又麻烦
稍微耐心点 把自己的问题说清楚 能把问题描述的清清楚楚让别人最方便的帮助 ...

根据信噪比和均方根误差可以看出固定阈值方法比较好,但是看图的话用固定阈值的方法好像和原始语音信号一样,所以想问问还有没有其他的方法,是不是这个不对
clear;
t1=clock;
[YSJ,sr]= wavread('Voice.wav');
sr=8000;
%数据预处理,数据可能是存储在矩阵或者是EXCEL中的二维数据,衔接为一维的,如果数据是一维数据,此步骤也不会影响数据
[c,l]=size(YSJ);
Y=[];
for i=1:c
    Y=[Y,YSJ(i,:)];
end
[c1,l1]=size(Y);
X=[1:l1];
%绘制噪声信号图像
figure();
subplot(221);
plot(X,Y);
xlabel('');
ylabel('幅值');
title('原始信号');
%硬阈值处理
lev=3;
xd=wden(Y,'heursure','h','one',lev,'db4');%硬阈值去噪处理后的信号序列
%figure(2)
subplot(222);
plot(X,xd)
xlabel('');
ylabel('幅值');
title('硬阈值去噪处理')
set(gcf,'Color',[1 1 1])
%软阈值处理
lev=3;
xs=wden(Y,'heursure','s','one',lev,'db4');%软阈值去噪处理后的信号序列
%figure(3)
subplot(223);
plot(X,xs)
xlabel('');
ylabel('幅值');
title('软阈值去噪处理')
set(gcf,'Color',[1 1 1])
%固定阈值后的去噪处理
lev=3;
xz=wden(Y,'sqtwolog','s','sln',lev,'db4');%固定阈值去噪处理后的信号序列
%figure(4)
subplot(224);
plot(X,xz);
xlabel('');
ylabel('幅值');
title('固定阈值后的去噪处理')
set(gcf,'Color',[1 1 1])
%计算信噪比SNR
Psig=sum(Y*Y')/l1;
Pnoi1=sum((Y-xd)*(Y-xd)')/l1;
Pnoi2=sum((Y-xs)*(Y-xs)')/l1;
Pnoi3=sum((Y-xz)*(Y-xz)')/l1;
SNR1=10*log10(Psig/Pnoi1);
SNR2=10*log10(Psig/Pnoi2);
SNR3=10*log10(Psig/Pnoi3);
%计算均方根误差RMSE
RMSE1=sqrt(Pnoi1);
RMSE2=sqrt(Pnoi2);
RMSE3=sqrt(Pnoi3);
%输出结果
disp('三种阈值设定方式的降噪处理结果');
disp(['硬阈值去噪处理的SNR=',num2str(SNR1),',RMSE=',num2str(RMSE1)]);
disp(['软阈值去噪处理的SNR=',num2str(SNR2),',RMSE=',num2str(RMSE2)]);
disp(['固定阈值后的去噪处理SNR=',num2str(SNR3),',RMSE=',num2str(RMSE3)]);
t2=clock;
tim=etime(t2,t1);
disp(['运行耗时',num2str(tim),'秒'])

Adel 发表于 2020-4-22 11:30:04
Jasmine_cmOlJ 发表于 2020-4-21 17:45
根据信噪比和均方根误差可以看出固定阈值方法比较好,但是看图的话用固定阈值的方法好像和原始语音信号一 ...

你参数选sln不是固定阈值吧,这也是根据小波系数算来的 具体怎么算就要去查算法了 和楼上一样的建议 我有机会也去研究下 先把wdwn说明贴这  免得每次都去找help

固定阈值 是用wthresh函数处理 你试下效果  但是对于不固定的信号来说 固定阈值基本也没意义

wden

wden

Jasmine_cmOlJ 发表于 2020-4-22 14:21:50
Adel 发表于 2020-4-22 11:30
你参数选sln不是固定阈值吧,这也是根据小波系数算来的 具体怎么算就要去查算法了 和楼上一样的建议 我有 ...

好的,谢谢你哈,我再研究研究。

wx_N89FLl7S 发表于 2020-4-22 15:29:07
Adel 发表于 2020-4-21 17:15
没呢 不知道matlab为啥算成这个样子,可能确实按照那几种规则的算法去算我们的系数,结果就是那样的
给你 ...

可以 open thselect内有具体的计算公式,另外还存在一个问题是,这个代码对手动加入的白噪声(当前每一层细节小波系数都置零)有很好的去噪效果,但是对实际的含噪信号无能为力(不知是否是因为当前我们的错误)。

Adel 发表于 2020-4-22 16:54:11
Jasmine_cmOlJ 发表于 2020-4-22 14:21
好的,谢谢你哈,我再研究研究。

不客气 共同进步吧 有发现可以分享哈 哈哈 难得找到共同研究的伙伴

Adel 发表于 2020-4-22 17:03:16
wx_N89FLl7S 发表于 2020-4-22 15:29
可以 open thselect内有具体的计算公式,另外还存在一个问题是,这个代码对手动加入的白噪声(当前每一层 ...

嗯嗯 谢谢提醒 后面我看下按照公式算下  这种将所有细节系数置零 肯定是不行的
另一种思路,你可以去查找一些论文,有各种改进小波阈值的,可以按照其公式确定阈值,也是可以的,我这搜到下面两个,还可以:
心电信号智能分析关键技术研究_姚成-吉林-博士2012
远程心电监护诊断系统心电信号处理方法研究-苏丽-哈工程-博士2006
建议参考下,当然只是我的意见,如果不好下载 我可以传给你

Adel 发表于 2020-4-22 18:10:00
wx_N89FLl7S 发表于 2020-4-22 15:29
可以 open thselect内有具体的计算公式,另外还存在一个问题是,这个代码对手动加入的白噪声(当前每一层 ...

我重新调试了下程序,就是阈值出问题那个。看thselect函数,确实有计算公式,但是前面两种基本不知道意思,也就是算法原理。后面两种,公式很简单,算的是和系数长度的对数直接相关。
试了四种方式,计算出的阈值都是比系数的绝对值最大值要大,而且越是后面几种越大,所以最终结果一样,都是把所有细节系数置零了。
这里我也贴下matlab的thselect函数代码,免得每次再去找吧
  1. function thr = thselect(x,tptr)
  2. %THSELECT Threshold selection for de-noising.
  3. %   THR = THSELECT(X,TPTR) returns threshold X-adapted value
  4. %   using selection rule defined by string TPTR.
  5. %   
  6. %   Available selection rules are:
  7. %   TPTR = 'rigrsure', adaptive threshold selection using
  8. %       principle of Stein's Unbiased Risk Estimate.
  9. %   TPTR = 'heursure', heuristic variant of the first option.
  10. %   TPTR = 'sqtwolog', threshold is sqrt(2*log(length(X))).
  11. %   TPTR = 'minimaxi', minimax thresholding.
  12. %
  13. %   Threshold selection rules are based on the underlying
  14. %   model y = f(t) + e where e is a white noise N(0,1).
  15. %   Dealing with unscaled or nonwhite noise can be handled
  16. %   using rescaling output threshold THR (see SCAL parameter
  17. %   in WDEN).
  18. %
  19. %   See also WDEN.

  20. %   M. Misiti, Y. Misiti, G. Oppenheim, J.M. Poggi 12-Mar-96.
  21. %   Last Revision: 20-Dec-2010.
  22. %   Copyright 1995-2010 The MathWorks, Inc.

  23. x = x(:)';
  24. n = length(x);
  25. switch tptr
  26.     case 'rigrsure'
  27.         sx2 = sort(abs(x)).^2;
  28.         risks = (n-(2*(1:n))+(cumsum(sx2)+(n-1:-1:0).*sx2))/n;
  29.         [~,best] = min(risks);
  30.         thr = sqrt(sx2(best));

  31.     case 'heursure'
  32.         hthr = sqrt(2*log(n));
  33.         eta = (norm(x).^2-n)/n;
  34.         crit = (log(n)/log(2))^(1.5)/sqrt(n);
  35.         if eta < crit
  36.             thr = hthr;
  37.         else
  38.             thr = min(thselect(x,'rigrsure'),hthr);
  39.         end

  40.     case 'sqtwolog'
  41.         thr = sqrt(2*log(n));

  42.     case 'minimaxi'
  43.         if n <= 32
  44.             thr = 0;
  45.         else
  46.             thr = 0.3936 + 0.1829*(log(n)/log(2));
  47.         end

  48.     otherwise
  49.         error(message('Wavelet:FunctionArgVal:Invalid_ArgVal'))
  50. end
复制代码



Miss清河洛 发表于 2020-5-7 16:58:42
Adel 发表于 2020-4-22 18:10
我重新调试了下程序,就是阈值出问题那个。看thselect函数,确实有计算公式,但是前面两种基本不知道意思 ...

我在进行小波阈值去噪对信号进行处理的时候也产生了相同的问题,去噪之后的高频系数都置零了,请问你找到解决的办法了么

Adel 发表于 2020-5-22 18:14:23
Miss清河洛 发表于 2020-5-7 16:58
我在进行小波阈值去噪对信号进行处理的时候也产生了相同的问题,去噪之后的高频系数都置零了,请问你找到 ...

没有 后面两种计算方式 明显就是计算出很大的阈值 所以阈值处理就全部得0了

best-time 发表于 2021-8-20 19:48:36
请问系数置0的问题解决了吗:loveliness:

AlaricYZB 发表于 2022-5-12 15:05:10
2022年了,题主明白为啥都置0了吗,或者应该怎么做呢,我现在也是这种情况
您需要登录后才可以回帖 登录 | 注册

本版积分规则

相关帖子
相关文章
热门教程
站长推荐
快速回复 返回顶部 返回列表