Wavelet Toolbox

wbmpen

Penalized threshold for wavelet 1-D or 2-D de-noising

Syntax

• THR = wbmpen(C,L,SIGMA,ALPHA)
  

Description

THR = wbmpen(C,L,SIGMA,ALPHA) returns global threshold THR for de-noising. THR is obtained by a wavelet coefficients selection rule using a penalization method provided by Birge-Massart.

[C,L] is the wavelet decomposition structure of the signal or image to be de-noised.

SIGMA is the standard deviation of the zero mean Gaussian white noise in de-noising model (see wnoisest for more information).

ALPHA is a tuning parameter for the penalty term. It must be a real number greater than 1. The sparsity of the wavelet representation of the de-noised signal or image grows with ALPHA. Typically ALPHA = 2.

THR minimizes the penalized criterion given by

let t^* be the minimizer of

• crit(t) = -sum(c(k)^2,k<=t) + 2*SIGMA^2*t*(ALPHA + log(n/t)) 
  

where c(k) are the wavelet coefficients sorted in decreasing order of their absolute value and n is the number of coefficients; then THR = c(t^*).

wbmpen(C,L,SIGMA,ALPHA,ARG) computes the global threshold and, in addition, plots three curves:

• 2*SIGMA^2*t*(ALPHA + log(n/t))
• sum(c(k)^2,k<=t)
• crit(t).

Examples

• % Example 1: Signal de-noising.
  % Load noisy bumps signal.
  load noisbump; x = noisbump;
  
  % Perform a wavelet decomposition of the signal
  % at level 5 using sym6.
  wname = 'sym6'; lev = 5;
  [c,l] = wavedec(x,lev,wname);
  % Estimate the noise standard deviation from the
  % detail coefficients at level 1, using wnoisest.
  sigma = wnoisest(c,l,1);
  
  % Use wbmpen for selecting global threshold  
  % for signal de-noising, using the tuning parameter.
  alpha = 2;
  thr = wbmpen(c,l,sigma,alpha)
  thr =
  
      2.7681
  
  % Use wdencmp for de-noising the signal using the above
  % threshold with soft thresholding and approximation kept.
  keepapp = 1;
  xd = wdencmp('gbl',c,l,wname,lev,thr,'s',keepapp);
  
  % Plot original and de-noised signals.
  figure(1)
  subplot(211), plot(x), title('Original signal')
  subplot(212), plot(xd), title('De-noised signal') 
  
  
  % Example 2: Image de-noising.
  % Load original image.
  load noiswom; 
  nbc = size(map,1);
  
  % Perform a wavelet decomposition of the image
  % at level 3 using coif2.
  wname = 'coif2'; lev = 3;
  [c,s] = wavedec2(X,lev,wname);
  
  % Estimate the noise standard deviation from the
  % detail coefficients at level 1.
  det1 = detcoef2('compact',c,s,1);
  sigma = median(abs(det1))/0.6745;
  
  % Use wbmpen for selecting global threshold  
  % for image de-noising.
  alpha = 1.2;
  thr = wbmpen(c,l,sigma,alpha)
  
  thr =
  
     36.0621
  
  % Use wdencmp for de-noising the image using the above
  % thresholds with soft thresholding and approximation kept.
  keepapp = 1;
  xd = wdencmp('gbl',c,s,wname,lev,thr,'s',keepapp);
  
  % Plot original and de-noised images.
  figure(2)
  colormap(pink(nbc));
  subplot(221), image(wcodemat(X,nbc))
  title('Original image')
  subplot(222), image(wcodemat(xd,nbc))
  title('De-noised image')
  
  
  

See Also
wden, wdencmp, wpbmpen, wpdencmp

waverec2

wcodemat