一维信号小波去噪函数：wden
一维信号小波去噪函数：wden
function [xd,cxd,lxd,thrs] = wden(in1,in2,in3,in4,in5,in6,in7)
该函数实现了一维信号小波去噪，支持多种阈值选择规则，包括 MODWT（最大重叠离散小波变换）和 DWT（离散小波变换）。它还支持软阈值和硬阈值处理，并支持阈值的重新缩放。
函数参数：
in1：输入信号x或者 DWT 系数c，或者 MODWT 变换后的矩阵W。in2：可选参数，阈值选择规则，字符串类型，支持以下选项：- 'modwtsqtwolog'：使用 MODWT 进行去噪，并使用 Donoho 和 Johnstone 的通用阈值和层级依赖阈值。
- 'rigrsure'：使用 Stein 的无偏风险原理进行阈值选择。
- 'heursure'：Stein 的无偏风险原理的一种启发式变体。
- 'sqtwolog'：使用 Donoho 和 Johnstone 的通用阈值，但使用 DWT 进行去噪。
- 'minimaxi'：使用最小极大阈值。
in3：可选参数，阈值类型，字符串类型，'s' 表示软阈值，'h' 表示硬阈值。in4：可选参数，阈值重新缩放类型，字符串类型，支持以下选项：- 'one'：不进行重新缩放。
- 'sln'：使用基于第一层系数的噪声估计进行重新缩放。
- 'mln'：使用层级依赖的噪声估计进行重新缩放。'mln' 仅支持 MODWT 去噪。
in5：可选参数，小波变换层级，整数类型。in6：可选参数，小波名称，字符串类型，对于 MODWT 去噪，必须是正交小波。in7：可选参数，DWT 系数的层级信息，向量类型，由WAVEDEC函数返回。
函数返回值：
xd：去噪后的信号。cxd：去噪后的小波系数，对于 DWT 去噪，是一个向量；对于 MODWT 去噪，是一个矩阵。lxd：DWT 去噪时，返回每层系数的数量，仅支持 DWT 去噪。thrs：DWT 去噪时，返回每层的阈值；MODWT 去噪时，仅支持 'modwtsqtwolog' 阈值选择规则，返回每层的阈值。
代码实现：
function [xd,cxd,lxd,thrs] = wden(in1,in2,in3,in4,in5,in6,in7)
%WDEN Automatic 1-D denoising using wavelets.
% XD = WDEN(X,TPTR,SORH,SCAL,N,WNAME) returns a denoised version XD
% of the input signal X obtained by thresholding the wavelet
% coefficients.
% TPTR is the threshold selection rule specified as a string. Supported
% options for TPTR are:
% 'modwtsqtwolog' uses the maximal overlap discrete wavelet
% transform (MODWT) to denoise the signal with Donoho and Johnstone's
% universal threshold and level-dependent thresholding.
% The remaining TPTR options use the critically sampled DWT to denoise the
% signal:
% 'rigrsure' uses the principle of Stein's Unbiased Risk.
% 'heursure' is a heuristic variant of Stein's Unbiased Risk.
% 'sqtwolog' uses Donoho and Johnstone's universal threshold with the
% DWT.
% 'minimaxi' uses minimax thresholding.
%
% SORH specifies soft or hard thresholding with 's' or 'h'.
% SCAL defines the type of threshold rescaling:
% 'one' for no rescaling.
% 'sln' for rescaling using a noise estimate based on the first-level
% coefficients.
% 'mln' for rescaling using level-dependent estimates of the noise.
% 'mln' is the only option supported for MODWT denoising.
%
% N is the level of the wavelet transform and WNAME is the wavelet
% specified as a string. For MODWT denoising, WNAME must correspond to an
% orthogonal wavelet.
%
% XD = wden(C,L,TPTR,SORH,SCAL,N,WNAME) returns the denoised
% signal obtained by operating on the DWT coefficient vector C and number
% of DWT coefficients by level L. C and L are outputs of WAVEDEC. You
% must use the same wavelet in both WAVEDEC and WDEN.
%
% XD = wden(W,'modwtsqtwolog',SORH,'mln',N,WNAME) returns the
% denoised signal obtained by operating on the MODWT transform matrix W.
% W is the output of MODWT. You must use the same wavelet in both
% MODWT and WDEN.
%
% [XD,CXD] = WDEN(...) returns the denoised wavelet coefficients. For
% DWT denoising, CXD is a vector (see WAVEDEC). For MODWT denoising, CXD
% is a matrix with N+1 rows (see MODWT). The number of columns is equal
% to the length of the input signal X.
%
% [XD,CXD,LXD] = WDEN(...) returns the number of coefficients by level
% for DWT denoising. See the help for WAVEDEC for details. The LXD output
% is not supported for MODWT denoising.
%
% [XD,CXD,LXD,THR] = WDEN(...) returns the denoising thresholds by level
% for DWT denosing.
%
% [XD,CXD,THR] = WDEN(...) returns the denoising thresholds by level for
% MODWT denoising when you specify 'modwtsqtwolog'.
%
% %Example 1:
% % Denoise a signal consisting of a 2-Hz sine wave with transients
% % at 0.3 and 0.72 seconds. Use Donoho and Johnstone's universal
% % threshold with level-dependent estimation of the noise.
% % Obtain denoised versions using the DWT and MODWT. Compare the
% % results.
% N = 1000;
% t = linspace(0,1,N);
% x = 4*sin(4*pi*t);
% x = x - sign(t - .3) - sign(.72 - t);
% y = x+0.15*randn(size(t));
% xdDWT = wden(y,'sqtwolog','s','mln',3,'db2');
% xdMODWT = wden(y,'modwtsqtwolog','s','mln',3,'db2');
% subplot(2,1,1)
% plot(xdDWT), title('DWT Denoising'); axis tight;
% subplot(2,1,2)
% plot(xdMODWT), title('MODWT Denoising'); axis tight;
%
% %Example 2:
% % Denoise a blocky signal using the Haar wavelet with MODWT and DWT
% % denoising. Compare the L2 and L-infty norms of the difference
% % between the original signal and the denoised versions.
% [x,xn] = wnoise('blocks',10,3);
% xdMODWT = wden(xn,'modwtsqtwolog','s','mln',6,'haar');
% xd = wden(xn,'sqtwolog','s','mln',6,'haar');
% plot(x)
% hold on
% plot(xd,'r--')
% plot(xdMODWT,'k-.')
% legend('Original','DWT','MODWT')
% hold off
% norm(abs(x-xd),2), norm(abs(x-xd),Inf)
% norm(abs(x-xdMODWT),2), norm(abs(x-xdMODWT),Inf)
%
% See also THSELECT, MODWT, WAVEDEC, WDENCMP, WFILTERS, WTHRESH.
%
% Copyright 1995-2015 The MathWorks, Inc.
% Check arguments.
% nargoutchk(0,4);
nbIn = nargin;
switch nbIn
case {0,1,2,3,4,5}
error(message('Wavelet:FunctionInput:NotEnough_ArgNum'));
case 6
x = in1; tptr = in2; sorh = in3;
scal = in4; n = in5; w = in6;
case 7
c = in1; l = in2; tptr = in3;
sorh = in4; scal = in5; n = in6; w = in7;
end
if errargt(mfilename,tptr,'str')
error(message('Wavelet:FunctionArgVal:Invalid_ArgVal'));
end
if errargt(mfilename,sorh,'str')
error(message('Wavelet:FunctionArgVal:Invalid_ArgVal'));
end
if errargt(mfilename,scal,'str')
error(message('Wavelet:FunctionArgVal:Invalid_ArgVal'));
end
if errargt(mfilename,n,'int')
error(message('Wavelet:FunctionArgVal:Invalid_ArgVal'));
end
if errargt(mfilename,w,'str')
error(message('Wavelet:FunctionArgVal:Invalid_ArgVal'));
end
% Adding MODWT denoising
if strcmpi(tptr,'modwtsqtwolog')
if nargout>3
error(message('Wavelet:modwt:TooManyOutputs'));
end
[xd,cxd,lxd] = modwtdenoise1D(x,w,n,sorh,scal);
return;
end
if nbIn==6
% Wavelet decomposition of x.
[c,l] = wavedec(x,n,w);
end
% Threshold rescaling coefficients.
switch scal
case 'one' , s = ones(1,n);
case 'sln' , s = ones(1,n)*wnoisest(c,l,1);
case 'mln' , s = wnoisest(c,l,1:n);
otherwise
error(message('Wavelet:FunctionArgVal:Invalid_ArgVal'));
end
% Wavelet coefficients thresholding.
first = cumsum(l)+1;
first = first(end-2:-1:1);
ld = l(end-1:-1:2);
last = first+ld-1;
cxd = c;
lxd = l;
if iscolumn(c)
thrs = zeros(n,1);
else
thrs = zeros(1,n);
end
for k = 1:n
flk = first(k):last(k);
if strcmp(tptr,'sqtwolog') || strcmp(tptr,'minimaxi')
thr = thselect(c,tptr);
else
if s(k) < sqrt(eps) * max(c(flk))
thr = 0;
else
thr = thselect(c(flk)/s(k),tptr);
end
end % threshold.
thrs(k) = thr * s(k); % rescaled threshold.
cxd(flk) = wthresh(c(flk),sorh,thrs(k)); % thresholding or shrinking.
end
% Wavelet reconstruction of xd.
xd = waverec(cxd,lxd,w);
end
代码解释：
- 输入参数检查：检查输入参数的数量和类型，并根据输入参数确定是使用 DWT 还是 MODWT。
- 小波分解：如果输入的是原始信号
x，则使用wavedec函数进行 DWT 分解。 - 阈值重新缩放：根据选择的
scal参数，对小波系数进行重新缩放，以适应不同的噪声水平。 - 小波系数阈值处理：根据选择的阈值选择规则
tptr和阈值类型sorh，对小波系数进行阈值处理，得到去噪后的系数cxd。 - 小波重构：使用
waverec函数，将去噪后的系数cxd重构为去噪后的信号xd。
注意：
- 该函数使用
modwtdenoise1D函数进行 MODWT 去噪，需要单独定义modwtdenoise1D函数。 - 该函数使用
thselect函数选择阈值，需要单独定义thselect函数。 - 该函数使用
wnoisest函数估计噪声水平，需要单独定义wnoisest函数。 - 该函数使用
wthresh函数对小波系数进行阈值处理，需要单独定义wthresh函数。 - 该函数使用
waverec函数进行小波重构，需要单独定义waverec函数。
使用示例：
% 1. 使用 DWT 去噪一个包含 2Hz 正弦波和瞬态信号的信号
N = 1000;
t = linspace(0,1,N);
x = 4*sin(4*pi*t);
x = x - sign(t - 0.3) - sign(0.72 - t);
y = x + 0.15*randn(size(t));
xddwt = wden(y,'sqtwolog','s','mln',3,'db2');
% 2. 使用 MODWT 去噪一个包含方块信号的信号
[x,xn] = wnoise('blocks',10,3);
xdmodwt = wden(xn,'modwtsqtwolog','s','mln',6,'haar');
总结：
wden 函数提供了一种灵活且强大的方法，可以对一维信号进行小波去噪。它支持多种阈值选择规则、阈值类型和重新缩放方法，可以根据不同的信号特征和噪声水平选择最佳参数，以获得最佳去噪效果。
原文地址: https://www.cveoy.top/t/topic/n9Ov 著作权归作者所有。请勿转载和采集!