comparison util/Rice Wavelet Toolbox/denoise.m @ 78:f69ae88b8be5

added Rice Wavelet Toolbox with my modification, so it can be compiled on newer systems.

+function [xd,xn,option] = denoise(x,h,type,option)
+%    [xd,xn,option] = denoise(x,h,type,option); 
+%
+%    DENOISE is a generic program for wavelet based denoising.
+%    The program will denoise the signal x using the 2-band wavelet
+%    system described by the filter h using either the traditional 
+%    discrete wavelet transform (DWT) or the linear shift invariant 
+%    discrete wavelet transform (also known as the undecimated DWT
+%    (UDWT)). 
+%
+%    Input:  
+%       x         : 1D or 2D signal to be denoised
+%       h         : Scaling filter to be applied
+%       type      : Type of transform (Default: type = 0)
+%                   0 --> Discrete wavelet transform (DWT)
+%                   1 --> Undecimated DWT (UDWT)
+%       option    : Default settings is marked with '*':
+%                   *type = 0 --> option = [0 3.0 0 0 0 0]
+%                   type = 1 --> option = [0 3.6 0 1 0 0]
+%       option(1) : Whether to threshold low-pass part
+%                   0 --> Don't threshold low pass component 
+%                   1 --> Threshold low pass component
+%       option(2) : Threshold multiplier, c. The threshold is
+%                   computed as: 
+%                     thld = c*MAD(noise_estimate)). 
+%                   The default values are:
+%                     c = 3.0 for the DWT based denoising
+%                     c = 3.6 for the UDWT based denoising
+%       option(3) : Type of variance estimator
+%                   0 --> MAD (mean absolute deviation)
+%                   1 --> STD (classical numerical std estimate)
+%       option(4) : Type of thresholding
+%                   0 --> Soft thresholding
+%                   1 --> Hard thresholding
+%       option(5) : Number of levels, L, in wavelet decomposition. By
+%                   setting this to the default value '0' a maximal
+%                   decomposition is used.
+%       option(6) : Actual threshold to use (setting this to
+%                   anything but 0 will mean that option(3)
+%                   is ignored)
+%
+%    Output: 
+%       xd     : Estimate of noise free signal 
+%       xn     : The estimated noise signal (x-xd)
+%       option : A vector of actual parameters used by the
+%                program. The vector is configured the same way as
+%                the input option vector with one added element
+%                option(7) = type.
+%
+%  HERE'S AN EASY WAY TO RUN THE EXAMPLES:
+%  Cut-and-paste the example you want to run to a new file 
+%  called ex.m, for example. Delete out the % at the beginning 
+%  of each line in ex.m (Can use search-and-replace in your editor
+%  to replace it with a space). Type 'ex' in matlab and hit return.
+%
+%    Example 1: 
+%       h = daubcqf(6); [s,N] = makesig('Doppler'); n = randn(1,N);
+%       x = s + n/10;     % (approximately 10dB SNR)
+%       figure;plot(x);hold on;plot(s,'r');
+%
+%       %Denoise x with the default method based on the DWT
+%       [xd,xn,opt1] = denoise(x,h);
+%       figure;plot(xd);hold on;plot(s,'r');
+%
+%       %Denoise x using the undecimated (LSI) wavelet transform
+%       [yd,yn,opt2] = denoise(x,h,1);
+%       figure;plot(yd);hold on;plot(s,'r');
+%
+% Example 2: (on an image)  
+%      h = daubcqf(6);  load lena; 
+%      noisyLena = lena + 25 * randn(size(lena));
+%      figure; colormap(gray); imagesc(lena); title('Original Image');
+%       figure; colormap(gray); imagesc(noisyLena); title('Noisy Image'); 
+%       Denoise lena with the default method based on the DWT
+%      [denoisedLena,xn,opt1] = denoise(noisyLena,h);
+%      figure; colormap(gray); imagesc(denoisedLena); title('denoised Image');
+%       
+%
+%    See also: mdwt, midwt, mrdwt, mirdwt, SoftTh, HardTh, setopt
+%
+
+%File Name: denoise.m
+%Last Modification Date: 04/15/97	10:44:28
+%Current Version: denoise.m	2.4
+%File Creation Date: Mon Feb 20 08:33:15 1995
+%Author: Jan Erik Odegard  <[e-mail]>
+%
+%Copyright (c) 2000 RICE UNIVERSITY. All rights reserved.
+%Created by Jan Erik Odegard, Department of ECE, Rice University. 
+%
+%This software is distributed and licensed to you on a non-exclusive 
+%basis, free-of-charge. Redistribution and use in source and binary forms, 
+%with or without modification, are permitted provided that the following 
+%conditions are met:
+%
+%1. Redistribution of source code must retain the above copyright notice, 
+%   this list of conditions and the following disclaimer.
+%2. Redistribution in binary form must reproduce the above copyright notice, 
+%   this list of conditions and the following disclaimer in the 
+%   documentation and/or other materials provided with the distribution.
+%3. All advertising materials mentioning features or use of this software 
+%   must display the following acknowledgment: This product includes 
+%   software developed by Rice University, Houston, Texas and its contributors.
+%4. Neither the name of the University nor the names of its contributors 
+%   may be used to endorse or promote products derived from this software 
+%   without specific prior written permission.
+%
+%THIS SOFTWARE IS PROVIDED BY WILLIAM MARSH RICE UNIVERSITY, HOUSTON, TEXAS, 
+%AND CONTRIBUTORS AS IS AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, 
+%BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS 
+%FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL RICE UNIVERSITY 
+%OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 
+%EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 
+%PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; 
+%OR BUSINESS INTERRUPTIONS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, 
+%WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR 
+%OTHERWISE), PRODUCT LIABILITY, OR OTHERWISE ARISING IN ANY WAY OUT OF THE 
+%USE OF THIS SOFTWARE,  EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+%
+%For information on commercial licenses, contact Rice University's Office of 
+%Technology Transfer at [e-mail] or [phone]
+%
+%Change History: Fixed output of function and an error in the computation
+%                of the threshold for redundant denoising.
+%                <Jan Erik Odegard> <Mon Jul 31, 1995>
+%
+%                This code is composed of several of our old codes for
+%                wavelet based denoising. In an effort to make the mess
+%                more manageable we decided to create on code that would 
+%                handle all the various wavelet based denoising methods.
+%                However, only time will show (as we discover new and 
+%                improved forms of denoising) if we can succeed in our goals.
+%                <Jan Erik Odegard> <Thu May 11, 1995>
+%
+
+if(nargin < 2)
+  error('You need to provide at least 2 inputs: x and h');
+end;
+if(nargin < 3),
+  type = 0;
+  option = [];
+elseif(nargin < 4)
+  option = [];
+end;
+if(isempty(type)),
+  type = 0;
+end;
+if(type == 0),
+  default_opt = [0 3.0 0 0 0 0];
+elseif(type == 1),
+  default_opt = [0 3.6 0 1 0 0];
+else,
+  error(['Unknown denoising method',10,...
+	  'If it is any good we need to have a serious talk :-)']);
+end;
+option = setopt(option,default_opt);
+[mx,nx] = size(x);
+dim = min(mx,nx);
+if(dim == 1),
+  n = max(mx,nx);
+else,
+  n = dim;
+end;
+if(option(5) == 0),
+  L = floor(log2(n));
+else
+  L = option(5);
+end;
+if(type == 0), 			% Denoising by DWT
+  xd = mdwt(x,h,L);
+  if (option(6) == 0),
+    tmp = xd(floor(mx/2)+1:mx,floor(nx/2)+1:nx);
+    if(option(3) == 0),
+      thld = option(2)*median(abs(tmp(:)))/.67;
+    elseif(option(3) == 1),
+      thld = option(2)*std(tmp(:));
+    else
+      error('Unknown threshold estimator, Use either MAD or STD');
+    end;
+  else,
+    thld = option(6);
+  end;
+  if(dim == 1)
+    ix = 1:n/(2^L);
+    ykeep = xd(ix);
+  else
+    ix = 1:mx/(2^L);
+    jx = 1:nx/(2^L);
+    ykeep = xd(ix,jx);
+  end;
+  if(option(4) == 0),
+    xd = SoftTh(xd,thld);
+  elseif(option(4) == 1),
+    xd = HardTh(xd,thld);
+  else,
+    error('Unknown threshold rule. Use either Soft (0) or Hard (1)');
+  end;
+  if (option(1) == 0),
+    if(dim == 1),
+      xd(ix) = ykeep;
+    else,
+      xd(ix,jx) = ykeep;
+    end;
+  end;
+  xd = midwt(xd,h,L);
+elseif(type == 1), 			% Denoising by UDWT
+  [xl,xh] = mrdwt(x,h,L);
+  if(dim == 1),
+    c_offset = 1;
+  else,
+    c_offset = 2*nx + 1;
+  end;
+  if (option(6) == 0),
+    tmp = xh(:,c_offset:c_offset+nx-1);
+    if(option(3) == 0),
+      thld = option(2)*median(abs(tmp(:)))/.67;
+    elseif(option(3) == 1),
+      thld = option(2)*std(tmp(:));
+    else
+      error('Unknown threshold estimator, Use either MAD or STD');
+    end;
+  else,
+    thld = option(6);
+  end;
+  if(option(4) == 0),
+    xh = SoftTh(xh,thld);
+    if(option(1) == 1),
+      xl = SoftTh(xl,thld);
+    end;
+  elseif(option(4) == 1),
+    xh = HardTh(xh,thld);
+    if(option(1) == 1),
+      xl = HardTh(xl,thld);
+    end;
+  else,
+    error('Unknown threshold rule. Use either Soft (0) or Hard (1)');
+  end;
+  xd = mirdwt(xl,xh,h,L);
+else, 					% Denoising by unknown method
+  error(['Unknown denoising method',10,...
+         'If it is any good we need to have a serious talk :-)']);
+end;
+option(6) = thld;
+option(7) = type;
+xn = x - xd;