www.gusucode.com > wavelet工具箱matlab源码程序 > wavelet/wavelet/wden.m

    function [xd,cxd,lxd] = 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. 
%
%   %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.
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>2
            error(message('Wavelet:modwt:TooManyOutputs'));
        end
        [xd,cxd] = 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;
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.
    thr      = thr * s(k);                  % rescaled threshold.
    cxd(flk) = wthresh(c(flk),sorh,thr);    % thresholding or shrinking.
end

% Wavelet reconstruction of xd.
xd = waverec(cxd,lxd,w);