请教：matlab小波阈值去噪，阈值为什么取系数最大绝对值？

信号见附件。使用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;

你好，我最近也在做小波阈值去噪这部分，可以讨论下么？

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

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

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

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

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

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

请问这个问题解决了吗，我的代码和您的差不多，也是每一层阈值选的是每一侧层系数的最大值。

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

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

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),'秒'])

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

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

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

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

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

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

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

Jasmine_cmOlJ 发表于 2020-4-22 14:21
好的，谢谢你哈，我再研究研究。

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

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

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

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. 
    
21. %   M. Misiti, Y. Misiti, G. Oppenheim, J.M. Poggi 12-Mar-96.
    
22. %   Last Revision: 20-Dec-2010.
    
23. %   Copyright 1995-2010 The MathWorks, Inc.
    
24. 
    
25. x = x(:)';
    
26. n = length(x);
    
27. switch tptr
    
28.     case 'rigrsure'
    
29.         sx2 = sort(abs(x)).^2;
    
30.         risks = (n-(2*(1:n))+(cumsum(sx2)+(n-1:-1:0).*sx2))/n;
    
31.         [~,best] = min(risks);
    
32.         thr = sqrt(sx2(best));
    
33. 
    
34.     case 'heursure'
    
35.         hthr = sqrt(2*log(n));
    
36.         eta = (norm(x).^2-n)/n;
    
37.         crit = (log(n)/log(2))^(1.5)/sqrt(n);
    
38.         if eta < crit
    
39.             thr = hthr;
    
40.         else
    
41.             thr = min(thselect(x,'rigrsure'),hthr);
    
42.         end
    
43. 
    
44.     case 'sqtwolog'
    
45.         thr = sqrt(2*log(n));
    
46. 
    
47.     case 'minimaxi'
    
48.         if n <= 32
    
49.             thr = 0;
    
50.         else
    
51.             thr = 0.3936 + 0.1829*(log(n)/log(2));
    
52.         end
    
53. 
    
54.     otherwise
    
55.         error(message('Wavelet:FunctionArgVal:Invalid_ArgVal'))
    
56. end
    

复制代码

：

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

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

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

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

请问系数置0的问题解决了吗:loveliness:

2022年了，题主明白为啥都置0了吗，或者应该怎么做呢，我现在也是这种情况