denoise.m
上传用户:speoil
上传日期:2022-06-23
资源大小:224k
文件大小:9k
源码类别:

波变换

开发平台:

Matlab

  1. function [xd,xn,option] = denoise(x,h,type,option)
  2. %    [xd,xn,option] = denoise(x,h,type,option); 
  3. %
  4. %    DENOISE is a generic program for wavelet based denoising.
  5. %    The program will denoise the signal x using the 2-band wavelet
  6. %    system described by the filter h using either the traditional 
  7. %    discrete wavelet transform (DWT) or the linear shift invariant 
  8. %    discrete wavelet transform (also known as the undecimated DWT
  9. %    (UDWT)). 
  10. %
  11. %    Input:  
  12. %       x         : 1D or 2D signal to be denoised
  13. %       h         : Scaling filter to be applied
  14. %       type      : Type of transform (Default: type = 0)
  15. %                   0 --> Discrete wavelet transform (DWT)
  16. %                   1 --> Undecimated DWT (UDWT)
  17. %       option    : Default settings is marked with '*':
  18. %                   *type = 0 --> option = [0 3.0 0 0 0 0]
  19. %                   type = 1 --> option = [0 3.6 0 1 0 0]
  20. %       option(1) : Whether to threshold low-pass part
  21. %                   0 --> Don't threshold low pass component 
  22. %                   1 --> Threshold low pass component
  23. %       option(2) : Threshold multiplier, c. The threshold is
  24. %                   computed as: 
  25. %                     thld = c*MAD(noise_estimate)). 
  26. %                   The default values are:
  27. %                     c = 3.0 for the DWT based denoising
  28. %                     c = 3.6 for the UDWT based denoising
  29. %       option(3) : Type of variance estimator
  30. %                   0 --> MAD (mean absolute deviation)
  31. %                   1 --> STD (classical numerical std estimate)
  32. %       option(4) : Type of thresholding
  33. %                   0 --> Soft thresholding
  34. %                   1 --> Hard thresholding
  35. %       option(5) : Number of levels, L, in wavelet decomposition. By
  36. %                   setting this to the default value '0' a maximal
  37. %                   decomposition is used.
  38. %       option(6) : Actual threshold to use (setting this to
  39. %                   anything but 0 will mean that option(3)
  40. %                   is ignored)
  41. %
  42. %    Output: 
  43. %       xd     : Estimate of noise free signal 
  44. %       xn     : The estimated noise signal (x-xd)
  45. %       option : A vector of actual parameters used by the
  46. %                program. The vector is configured the same way as
  47. %                the input option vector with one added element
  48. %                option(7) = type.
  49. %
  50. %  HERE'S AN EASY WAY TO RUN THE EXAMPLES:
  51. %  Cut-and-paste the example you want to run to a new file 
  52. %  called ex.m, for example. Delete out the % at the beginning 
  53. %  of each line in ex.m (Can use search-and-replace in your editor
  54. %  to replace it with a space). Type 'ex' in matlab and hit return.
  55. %
  56. %    Example 1: 
  57. %       h = daubcqf(6); [s,N] = makesig('Doppler'); n = randn(1,N);
  58. %       x = s + n/10;     % (approximately 10dB SNR)
  59. %       figure;plot(x);hold on;plot(s,'r');
  60. %
  61. %       %Denoise x with the default method based on the DWT
  62. %       [xd,xn,opt1] = denoise(x,h);
  63. %       figure;plot(xd);hold on;plot(s,'r');
  64. %
  65. %       %Denoise x using the undecimated (LSI) wavelet transform
  66. %       [yd,yn,opt2] = denoise(x,h,1);
  67. %       figure;plot(yd);hold on;plot(s,'r');
  68. %
  69. % Example 2: (on an image)  
  70. %      h = daubcqf(6);  load lena; 
  71. %      noisyLena = lena + 25 * randn(size(lena));
  72. %      figure; colormap(gray); imagesc(lena); title('Original Image');
  73. %       figure; colormap(gray); imagesc(noisyLena); title('Noisy Image'); 
  74. %       Denoise lena with the default method based on the DWT
  75. %      [denoisedLena,xn,opt1] = denoise(noisyLena,h);
  76. %      figure; colormap(gray); imagesc(denoisedLena); title('denoised Image');
  77. %       
  78. %
  79. %    See also: mdwt, midwt, mrdwt, mirdwt, SoftTh, HardTh, setopt
  80. %
  81. %File Name: denoise.m
  82. %Last Modification Date: 04/15/97 10:44:28
  83. %Current Version: denoise.m 2.4
  84. %File Creation Date: Mon Feb 20 08:33:15 1995
  85. %Author: Jan Erik Odegard  <odegard@ece.rice.edu>
  86. %
  87. %Copyright (c) 2000 RICE UNIVERSITY. All rights reserved.
  88. %Created by Jan Erik Odegard, Department of ECE, Rice University. 
  89. %
  90. %This software is distributed and licensed to you on a non-exclusive 
  91. %basis, free-of-charge. Redistribution and use in source and binary forms, 
  92. %with or without modification, are permitted provided that the following 
  93. %conditions are met:
  94. %
  95. %1. Redistribution of source code must retain the above copyright notice, 
  96. %   this list of conditions and the following disclaimer.
  97. %2. Redistribution in binary form must reproduce the above copyright notice, 
  98. %   this list of conditions and the following disclaimer in the 
  99. %   documentation and/or other materials provided with the distribution.
  100. %3. All advertising materials mentioning features or use of this software 
  101. %   must display the following acknowledgment: This product includes 
  102. %   software developed by Rice University, Houston, Texas and its contributors.
  103. %4. Neither the name of the University nor the names of its contributors 
  104. %   may be used to endorse or promote products derived from this software 
  105. %   without specific prior written permission.
  106. %
  107. %THIS SOFTWARE IS PROVIDED BY WILLIAM MARSH RICE UNIVERSITY, HOUSTON, TEXAS, 
  108. %AND CONTRIBUTORS AS IS AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, 
  109. %BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS 
  110. %FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL RICE UNIVERSITY 
  111. %OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 
  112. %EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 
  113. %PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; 
  114. %OR BUSINESS INTERRUPTIONS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, 
  115. %WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR 
  116. %OTHERWISE), PRODUCT LIABILITY, OR OTHERWISE ARISING IN ANY WAY OUT OF THE 
  117. %USE OF THIS SOFTWARE,  EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  118. %
  119. %For information on commercial licenses, contact Rice University's Office of 
  120. %Technology Transfer at techtran@rice.edu or (713) 348-6173
  121. %
  122. %Change History: Fixed output of function and an error in the computation
  123. %                of the threshold for redundant denoising.
  124. %                <Jan Erik Odegard> <Mon Jul 31, 1995>
  125. %
  126. %                This code is composed of several of our old codes for
  127. %                wavelet based denoising. In an effort to make the mess
  128. %                more manageable we decided to create on code that would 
  129. %                handle all the various wavelet based denoising methods.
  130. %                However, only time will show (as we discover new and 
  131. %                improved forms of denoising) if we can succeed in our goals.
  132. %                <Jan Erik Odegard> <Thu May 11, 1995>
  133. %
  134. if(nargin < 2)
  135.   error('You need to provide at least 2 inputs: x and h');
  136. end;
  137. if(nargin < 3),
  138.   type = 0;
  139.   option = [];
  140. elseif(nargin < 4)
  141.   option = [];
  142. end;
  143. if(isempty(type)),
  144.   type = 0;
  145. end;
  146. if(type == 0),
  147.   default_opt = [0 3.0 0 0 0 0];
  148. elseif(type == 1),
  149.   default_opt = [0 3.6 0 1 0 0];
  150. else,
  151.   error(['Unknown denoising method',10,...
  152.   'If it is any good we need to have a serious talk :-)']);
  153. end;
  154. option = setopt(option,default_opt);
  155. [mx,nx] = size(x);
  156. dim = min(mx,nx);
  157. if(dim == 1),
  158.   n = max(mx,nx);
  159. else,
  160.   n = dim;
  161. end;
  162. if(option(5) == 0),
  163.   L = floor(log2(n));
  164. else
  165.   L = option(5);
  166. end;
  167. if(type == 0),  % Denoising by DWT
  168.   xd = mdwt(x,h,L);
  169.   if (option(6) == 0),
  170.     tmp = xd(floor(mx/2)+1:mx,floor(nx/2)+1:nx);
  171.     if(option(3) == 0),
  172.       thld = option(2)*median(abs(tmp(:)))/.67;
  173.     elseif(option(3) == 1),
  174.       thld = option(2)*std(tmp(:));
  175.     else
  176.       error('Unknown threshold estimator, Use either MAD or STD');
  177.     end;
  178.   else,
  179.     thld = option(6);
  180.   end;
  181.   if(dim == 1)
  182.     ix = 1:n/(2^L);
  183.     ykeep = xd(ix);
  184.   else
  185.     ix = 1:mx/(2^L);
  186.     jx = 1:nx/(2^L);
  187.     ykeep = xd(ix,jx);
  188.   end;
  189.   if(option(4) == 0),
  190.     xd = SoftTh(xd,thld);
  191.   elseif(option(4) == 1),
  192.     xd = HardTh(xd,thld);
  193.   else,
  194.     error('Unknown threshold rule. Use either Soft (0) or Hard (1)');
  195.   end;
  196.   if (option(1) == 0),
  197.     if(dim == 1),
  198.       xd(ix) = ykeep;
  199.     else,
  200.       xd(ix,jx) = ykeep;
  201.     end;
  202.   end;
  203.   xd = midwt(xd,h,L);
  204. elseif(type == 1),  % Denoising by UDWT
  205.   [xl,xh] = mrdwt(x,h,L);
  206.   if(dim == 1),
  207.     c_offset = 1;
  208.   else,
  209.     c_offset = 2*nx + 1;
  210.   end;
  211.   if (option(6) == 0),
  212.     tmp = xh(:,c_offset:c_offset+nx-1);
  213.     if(option(3) == 0),
  214.       thld = option(2)*median(abs(tmp(:)))/.67;
  215.     elseif(option(3) == 1),
  216.       thld = option(2)*std(tmp(:));
  217.     else
  218.       error('Unknown threshold estimator, Use either MAD or STD');
  219.     end;
  220.   else,
  221.     thld = option(6);
  222.   end;
  223.   if(option(4) == 0),
  224.     xh = SoftTh(xh,thld);
  225.     if(option(1) == 1),
  226.       xl = SoftTh(xl,thld);
  227.     end;
  228.   elseif(option(4) == 1),
  229.     xh = HardTh(xh,thld);
  230.     if(option(1) == 1),
  231.       xl = HardTh(xl,thld);
  232.     end;
  233.   else,
  234.     error('Unknown threshold rule. Use either Soft (0) or Hard (1)');
  235.   end;
  236.   xd = mirdwt(xl,xh,h,L);
  237. else,  % Denoising by unknown method
  238.   error(['Unknown denoising method',10,...
  239.          'If it is any good we need to have a serious talk :-)']);
  240. end;
  241. option(6) = thld;
  242. option(7) = type;
  243. xn = x - xd;