matlab中的wden函数（Automatic 1-D Denoising Using Wavelets）
matlab中的wden函数（Automatic 1-D Denoising Using Wavelets）
为了处理数学建模中的数据进行必要的学习:小波降噪
（本文参考的是matlab官方文档，仅供学习使用）
目录
wden函数
XD = wden(X,TPTR,SORH,SCAL,N,wname)
` returns a denoised version `XD` of the signal `X`. The function uses an `N`-level wavelet decomposition of `X` using the specified orthogonal or biorthogonal wavelet `wname` to obtain the wavelet coefficients. The thresholding selection rule `TPTR` is applied to the wavelet decomposition. `SORH` and `SCAL` define how the rule is applied.
[XD,CXD,LXD,THR] = wden(___)
returns the denoising thresholds by level for DWT denoising
对上面重要的参数进行解释
TPTR — Threshold selection rule character array
Threshold selection rule to apply to the wavelet decomposition structure of X:
'rigsure'— Use the principle of Stein's Unbiased Risk.'heursure'— Use a heuristic variant of Stein's Unbiased Risk.'sqtwolog— Use the universal threshold G2ln(length(x)).'minimaxi'— Use minimax thresholding. (Seethselectfor more information.)
SORH — Type of thresholding 's' | 'h'
Type of thresholding to perform:
's'— Soft thresholding'h'— Hard thresholding
SCAL — Multiplicative threshold rescaling 'one' | 'sln' | 'mln'
Multiplicative threshold rescaling:
'one'— No rescaling （网上有人说不随噪声水平变化）'sln'— Rescaling using a single estimation of level noise based on first-level coefficients'mln'— Rescaling using a level-dependent estimation of level noise
噪声处理样例
噪声生成
首先我们先随便制造个噪声：
% Create a signal consisting of a 2 Hz sine wave with transients at 0.3 and 0.72 seconds. Add randomly generated noise to the signal and plot the result.
N = 1000;
t = linspace(0,1,N);
x = 4*sin(4*pi*t);
x = x - sign(t-0.3) - sign(0.72-t);
sig = x + 0.5*randn(size(t));
plot(t,sig)
title('Signal')
grid on
细节说明：
这里可能很多是数字信号处理的内容，不过提前了解也无妨；下面对上面噪声生成的细节进行补充说明:
sign函数（sgn函数、符号函数）：
plot(t,sign(t-0.5))
axis([0,1,-2,2])
这个函数的值域是从-1到1，为了方便显示，向右将这个信号平移了0.5个单位。假设我们对另一个信号x(t)加入这个信号sign(t-0.5)，容易理解，作用后的函数x'(t)=x(t)+sign(t-0.5)在t=0.5处断开了，之前的函数图像向下平移一个单位，之后的函数图像向上平移一个单位
细节：这里不要和阶跃函数混淆，阶跃函数的范围是0到1
x = 4*sin(4*pi*t);
x = x - sign(t-0.3) - sign(0.72-t);
这里上面已经说了这里生成一个2Hz的噪声，我们要注意的是下面这行代码，实际上它给x(t)增加是这样的：
噪声
rng('default')
sig = x + 0.5*randn(size(t));
这里是用默认种子去生成伪随机数，噪声我们是每一个t值都对应了一个伪随机数，加了噪声之后我们得到了我们最后的结果：
降噪
硬阈值
% Using the sym8 wavelet, perform a level 5 wavelet decomposition of the signal and denoise it by applying three different threshold selection rules to the wavelet coefficients: SURE, minimax, and Donoho and Johnstone's universal threshold with level-dependent estimation of the noise. In each case, apply hard thresholding.
lev = 5;
wname = 'sym8';
[dnsig1,c1,l1,threshold_SURE] = wden(sig,'rigrsure','h','mln',lev,wname);
[dnsig2,c2,l2,threshold_Minimax] = wden(sig,'minimaxi','h','mln',lev,wname);
[dnsig3,c3,l3,threshold_DJ] = wden(sig,'sqtwolog','h','mln',lev,wname);
[dnsig4,c4,l4,threshold_heursure] = wden(sig,'heursure','h','mln',lev,wname);
这里原理先不做进一步分析，目前只需要知道这样用就够了，以后有时间再回来写分析
subplot(4,1,1)
plot(t,dnsig1)
title('Denoised Signal - SURE')
grid on
subplot(4,1,2)
plot(t,dnsig2)
title('Denoised Signal - Minimax')
grid on
subplot(4,1,3)
plot(t,dnsig3)
title('Denoised Signal - Donoho-Johnstone')
grid on
subplot(4,1,4)
plot(t,dnsig4)
title('Denoised Signal - Heursure')
grid on
最后我们返回四个阈值看一下
threshold_SURE %stein无偏估计
threshold_Minimax %极大极小值阈值
threshold_DJ %固定式阈值
threshold_heursure %启发式阈值
软阈值
噪声生成部分还是上面那种生成方式
% Using the db2 wavelet, perform a level 3 wavelet decomposition of the signal and denoise it using Donoho and Johnstone's universal threshold with level-dependent estimation of the noise. Obtain denoised versions using DWT and MODWT, both with soft thresholding.
wname = 'db2';
lev = 3;
[xdDWT,c1,l1,threshold_DWT] = wden(sig,'sqtwolog','s','mln',lev,wname);
[xdMODWT,c2,threshold_MODWT] = wden(sig,'modwtsqtwolog','s','mln',lev,wname);
画图
subplot(2,1,1)
plot(t,xdDWT)
grid on
title('DWT Denoising')
subplot(2,1,2)
plot(t,xdMODWT)
grid on
title('MODWT Denoising')
阈值
块状信号噪声处理
噪声生成
%Compare DWT and MODWT Denoising of a Blocky Signal
%This example denoises a blocky signal using the Haar wavelet with DWT and MODWT denoising. It compares the results with plots and metrics for the original and denoised versions.
%First, to ensure reproducibility of results, set a seed that will be used to generate random noise.
rng('default')
%Generate a signal and a noisy version with the square root of the signal-to-noise ratio equal to 3. Plot and compare each.
[osig,nsig] = wnoise('blocks',10,3);
plot(nsig,'r')
hold on
plot(osig,'b')
legend('Noisy Signal','Original Signal')
噪声处理（软阈值）
%Using the Haar wavelet, perform a level 6 wavelet decomposition of the noisy signal and denoise it using Donoho and Johnstone's universal threshold with level-dependent estimation of the noise. Obtain denoised versions using DWT and MODWT, both with soft thresholding.
wname = 'haar';
lev = 6 ;
[xdDWT,c1,l1] = wden(nsig,'sqtwolog','s','mln',lev,wname);
[xdMODWT,c2] = wden(nsig,'modwtsqtwolog','s','mln',lev,wname);
画图
figure
plot(osig,'b')
hold on
plot(xdDWT,'r--')
plot(xdMODWT,'k-.')
legend('Original','DWT','MODWT')
hold off
硬软阈值的区别
定义：
- 硬阈值
\[\begin{align}
\eta_H(x,th_H) & = x\cdot \varepsilon (|x|-th_H)\\
& =
\left \{
\begin{matrix}
x,\quad |x|>th_H\\
0,\quad|x|<th_H
\end{matrix}
\right.\\
(where\;\varepsilon(t)\;is\;&Step\;Function(阶跃函数) )
\end{align}
\]
- 软阈值
\[\begin{align}
\eta_S(x,th_S) & = sgn(x)\cdot (|x|-th_s)\\
& =
\left \{
\begin{matrix}
x+th_S,\quad x<-th_S\\
0,\quad|x|<th_S\\
x-th_S,\quad x>th_S
\end{matrix}
\right.\\
&(sgn(t)是符号函数)
\end{align}
\]
简单画一下图是这样的
如果觉得这个画不好的话，我们也可以用matlab来画
y=linspace(-1,1,100);
>> sigH=wthresh(y,'h',0.4);
>> sigS=wthresh(y,'s',0.4);
>> subplot(1,3,1)
>> plot(y,'b')
>> title('Origin Signal')
>> subplot(1,3,2)
>> plot(sigH,'r')
>> title('Hardthresh Signal')
>> subplot(1,3,3)
>> plot(sigS,'g')
>> title('Softthresh Signal')
更进一步的区别可以参考一下大佬的话：
