开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 图形图像 > 波变换 > 查看源码
wthrmngr.m
资源名称：wavelet_toolbox.rar　[点击查看]
上传用户：haiyisale
上传日期：2013-01-09
资源大小：3246k
文件大小：25k
源码类别：

波变换

开发平台：

Matlab

wthrmngr.m：源码内容
  1. function varargout = wthrmngr(option,varargin)
  2. %WTHRMNGR Threshold settings manager.
  3. %   THR = WTHRMNGR(OPTION,METHOD,VARARGIN) returns a 
  4. %   global threshold or level dependent thresholds
  5. %   depending on OPTION. The inputs VARARGIN depend
  6. %   on OPTION and METHOD.
  7. %
  8. %   This function returns the thresholds used throughout
  9. %   the MATLAB Wavelet Toolbox for de-noising and 
  10. %   compression tools (command line M-files or GUI 
  11. %   tools).
  12. %
  13. %   The options for the METHOD parameter are:
  14. %
  15. %   - 'scarcehi'            (see WDCBM or WDCBM2 with high 
  16. %                                predefined value of parameter M).
  17. %   - 'scarceme'            (see WDCBM or WDCBM2 with medium  
  18. %                                predefined value of parameter M).
  19. %   - 'scarcelo'            (see WDCBM or WDCBM2 with low  
  20. %                                predefined value of parameter M).
  21. %
  22. %   - 'sqtwolog'            (see THSELECT option 'sqtwolog' 
  23. %                                and see also WDEN).
  24. %   - 'sqtwologuwn'         (see THSELECT option 'sqtwolog' 
  25. %                                and see also WDEN option 'sln').
  26. %   - 'sqtwologswn'         (see THSELECT option 'sqtwolog' 
  27. %                                and see also WDEN option 'mln').
  28. %   - 'rigrsure'            (see THSELECT option 'rigrsure' 
  29. %                                and see also WDEN).
  30. %   - 'heursure'            (see THSELECT option 'heursure' 
  31. %                                and see also WDEN).
  32. %   - 'minimaxi'            (see THSELECT option 'minimaxi' 
  33. %                                and see also WDEN).
  34. %
  35. %   - 'penalhi'             (see WBMPEN or WPBMPEN with high 
  36. %                                value of parameter ALPHA).
  37. %   - 'penalme'             (see WBMPEN or WPBMPEN with medium  
  38. %                                value of parameter ALPHA).
  39. %   - 'penallo'             (see WBMPEN or WPBMPEN with low  
  40. %                                value of parameter ALPHA).
  41. %
  42. %   - 'rem_n0'              this option returns a threshold
  43. %                           close to 0, a typical THR value is  
  44. %                           median(abs(coefficients)).
  45. %
  46. %   - 'bal_sn'              this option returns a threshold
  47. %                           such that the percentages of retained
  48. %                           energy and number of zeros are 
  49. %                           the same.
  50. %
  51. %   - 'sqrtbal_sn'          this option returns a threshold
  52. %                           equal to the square root of the value
  53. %                           such that the percentages of retained
  54. %                           energy and number of zeros are 
  55. %                           the same.
  56. %
  57. %   #############################  
  58. %   Discrete Wavelet 1-D options:
  59. %   #############################  
  60. %
  61. %    Compression using a global threshold:
  62. %    -------------------------------------
  63. %    X is the signal to be compressed and [C,L] is the wavelet 
  64. %    decomposition structure of the signal to be compressed. 
  65. %     THR = WTHRMNGR('dw1dcompGBL','rem_n0',X) 
  66. %     THR = WTHRMNGR('dw1dcompGBL','bal_sn',C,L)
  67. %
  68. %    Compression using level dependent thresholds:
  69. %    ---------------------------------------------
  70. %    X is the signal to be compressed and [C,L] is the wavelet 
  71. %    decomposition structure of the signal to be compressed.
  72. %    ALFA is a sparsity parameter (see WDCBM).
  73. %
  74. %     THR = WTHRMNGR('dw1dcompLVL','scarcehi',C,L,ALFA)
  75. %            ALFA must be such that 2.5 < ALFA < 10
  76. %     THR = WTHRMNGR('dw1dcompLVL','scarceme',C,L,ALFA)
  77. %            ALFA must be such that 1.5 < ALFA < 2.5
  78. %     THR = WTHRMNGR('dw1dcompLVL','scarcelo',C,L,ALFA)
  79. %            ALFA must be such that 1 < ALFA < 2
  80. %
  81. %    De-noising using level dependent thresholds:
  82. %    --------------------------------------------
  83. %    [C,L] is the wavelet decomposition structure of the
  84. %    signal to be de-noised, SCAL defines the 
  85. %    multiplicative threshold rescaling (see WDEN) and
  86. %    ALFA is a sparsity parameter (see WBMPEN).
  87. %
  88. %     THR = WTHRMNGR('dw1ddenoLVL','sqtwolog',C,L,SCAL)
  89. %     THR = WTHRMNGR('dw1ddenoLVL','rigrsure',C,L,SCAL)
  90. %     THR = WTHRMNGR('dw1ddenoLVL','heursure',C,L,SCAL)
  91. %     THR = WTHRMNGR('dw1ddenoLVL','minimaxi',C,L,SCAL)
  92. %
  93. %     THR = WTHRMNGR('dw1ddenoLVL','penalhi',C,L,ALFA)
  94. %            ALFA must be such that 2.5 < ALFA < 10
  95. %     THR = WTHRMNGR('dw1ddenoLVL','penalme',C,L,ALFA)
  96. %            ALFA must be such that 1.5 < ALFA < 2.5
  97. %     THR = WTHRMNGR('dw1ddenoLVL','penallo',C,L,ALFA)
  98. %            ALFA must be such that 1 < ALFA < 2
  99. %
  100. %   ########################################    
  101. %   Discrete Stationary Wavelet 1-D options:
  102. %   ########################################    
  103. %
  104. %    De-noising using level dependent thresholds:
  105. %    --------------------------------------------
  106. %    SWTDEC is the stationary wavelet decomposition structure 
  107. %    of the signal to be de-noised, SCAL defines the 
  108. %    multiplicative threshold rescaling (see WDEN) and
  109. %    ALFA is a sparsity parameter (see WBMPEN).
  110. %     THR = WTHRMNGR('sw1ddenoLVL',METHOD,SWTDEC,SCAL)
  111. %     THR = WTHRMNGR('sw1ddenoLVL',METHOD,SWTDEC,ALFA)
  112. %     The options for METHOD are the same as in the 'dw1ddenoLVL'
  113. %     case.
  114. %
  115. %   #############################  
  116. %   Discrete Wavelet 2-D options:
  117. %   #############################  
  118. %
  119. %    Compression using a global threshold:
  120. %    -------------------------------------
  121. %    X is the image to be compressed and [C,S] is the wavelet 
  122. %    decomposition structure of the image to be compressed.
  123. %     THR = WTHRMNGR('dw2dcompGBL','rem_n0',X)
  124. %     THR = WTHRMNGR('dw2dcompGBL','bal_sn',C,S)
  125. %     THR = WTHRMNGR('dw2dcompGBL','sqrtbal_sn',C,S)
  126. %
  127. %    Compression using level dependent thresholds:
  128. %    ---------------------------------------------
  129. %    X is the image to be compressed and [C,S] is the wavelet 
  130. %    decomposition structure of the image to be compressed.
  131. %    ALFA is a sparsity parameter (see WDCBM2).
  132. %
  133. %     THR = WTHRMNGR('dw2dcompLVL','scarcehi',C,S,ALFA)
  134. %            ALFA must be such that 2.5 < ALFA < 10
  135. %     THR = WTHRMNGR('dw2dcompLVL','scarceme',C,S,ALFA)
  136. %            ALFA must be such that 1.5 < ALFA < 2.5
  137. %     THR = WTHRMNGR('dw2dcompLVL','scarcelo',C,S,ALFA)
  138. %            ALFA must be such that 1 < ALFA < 2
  139. %
  140. %    De-noising using level dependent thresholds:
  141. %    --------------------------------------------
  142. %    [C,S] is the wavelet decomposition structure of the
  143. %    image to be de-noised, SCAL defines the 
  144. %    multiplicative threshold rescaling (see WDEN) and
  145. %    ALFA is a sparsity parameter (see WBMPEN).
  146. %
  147. %     THR = WTHRMNGR('dw2ddenoLVL','penalhi',C,S,ALFA)
  148. %            ALFA must be such that 2.5 < ALFA < 10
  149. %     THR = WTHRMNGR('dw2ddenoLVL','penalme',C,S,ALFA)
  150. %            ALFA must be such that 1.5 < ALFA < 2.5
  151. %     THR = WTHRMNGR('dw2ddenoLVL','penallo',C,S,ALFA)
  152. %            ALFA must be such that 1 < ALFA < 2
  153. %
  154. %     THR = WTHRMNGR('dw2ddenoLVL','sqtwolog',C,S,SCAL)
  155. %     THR = WTHRMNGR('dw2ddenoLVL','sqrtbal_sn',C,S)
  156. %
  157. %   ########################################  
  158. %   Discrete Stationary Wavelet 2-D options:
  159. %   ########################################  
  160. %
  161. %    De-noising using level dependent thresholds:
  162. %    --------------------------------------------
  163. %    SWTDEC is the stationary wavelet decomposition structure 
  164. %    of the image to be de-noised, SCAL defines the 
  165. %    multiplicative threshold rescaling (see WDEN) and
  166. %    ALFA is a sparsity parameter (see WBMPEN).
  167. %     THR = WTHRMNGR('sw2ddenoLVL',METHOD,SWTDEC,SCAL)
  168. %     THR = WTHRMNGR('sw2ddenoLVL',METHOD,SWTDEC,ALFA)
  169. %     The options for METHOD are the same as in the 'dw2ddenoLVL'
  170. %     case.
  171. %
  172. %   ####################################  
  173. %   Discrete Wavelet Packet 1-D options:
  174. %   #################################### 
  175. % 
  176. %    Compression using a global threshold:
  177. %    -------------------------------------
  178. %    X is the signal to be compressed and WPT is the wavelet 
  179. %    packet decomposition structure of the signal to be compressed.
  180. %     THR = WTHRMNGR('wp1dcompGBL','bal_sn',WPT)
  181. %     THR = WTHRMNGR('wp1dcompGBL','rem_n0',X)
  182. %
  183. %    De-noising using a global threshold:
  184. %    ------------------------------------
  185. %    WPT is the wavelet packet decomposition structure of the signal
  186. %    to be de-noised.
  187. %     THR = WTHRMNGR('wp1ddenoGBL','sqtwologuwn',WPT)
  188. %     THR = WTHRMNGR('wp1ddenoGBL','sqtwologswn',WPT)
  189. %     THR = WTHRMNGR('wp1ddenoGBL','bal_sn',WPT)
  190. %
  191. %     THR = WTHRMNGR('wp1ddenoGBL','penalhi',WPT)
  192. %            see WPBMPEN with ALFA = 6.25
  193. %     THR = WTHRMNGR('wp1ddenoGBL','penalme',WPT)
  194. %            see WPBMPEN with ALFA = 2
  195. %     THR = WTHRMNGR('wp1ddenoGBL','penallo',WPT)
  196. %            see WPBMPEN with ALFA = 1.5
  197. %
  198. %   ####################################  
  199. %   Discrete Wavelet Packet 2-D options:
  200. %   #################################### 
  201. % 
  202. %    Compression using a global threshold:
  203. %    -------------------------------------
  204. %    X is the image to be compressed and WPT is the wavelet 
  205. %    packet decomposition structure of the image to be compressed.
  206. %     THR = WTHRMNGR('wp2dcompGBL','bal_sn',WPT)
  207. %     THR = WTHRMNGR('wp2dcompGBL','rem_n0',X)
  208. %     THR = WTHRMNGR('wp2dcompGBL','sqrtbal_sn',WPT)
  209. %
  210. %    De-noising using a global threshold:
  211. %    ------------------------------------
  212. %    WPT is the wavelet packet decomposition structure of the image
  213. %    to be de-noised.
  214. %     THR = WTHRMNGR('wp2ddenoGBL','sqtwologuwn',WPT)
  215. %     THR = WTHRMNGR('wp2ddenoGBL','sqtwologswn',WPT)
  216. %     THR = WTHRMNGR('wp2ddenoGBL','sqrtbal_sn',WPT)
  217. %
  218. %     THR = WTHRMNGR('wp2ddenoGBL','penalhi',WPT)
  219. %            see WPBMPEN with ALFA = 6.25
  220. %     THR = WTHRMNGR('wp2ddenoGBL','penalme',WPT)
  221. %            see WPBMPEN with ALFA = 2
  222. %     THR = WTHRMNGR('wp2ddenoGBL','penallo',WPT)
  223. %            see WPBMPEN with ALFA = 1.5
  224. %
  225. %   See also THSELECT, WBMPEN, WDCBM, WDCBM2, WDEN, WDENCMP,  
  226. %            WNOISEST, WPBMPEN, WPDENCMP.
  228. %   M. Misiti, Y. Misiti, G. Oppenheim, J.M. Poggi 21-Oct-98.
  229. %   Last Revision: 26-Oct-1999.
  230. %   Copyright 1995-2002 The MathWorks, Inc.
  231. %   $Revision: 1.7 $  $Date: 2002/04/14 19:49:35 $
  233. meth = varargin{1};
  234. switch option
  235.    case {'dw1dcompGBL','dw1dcompLVL','dw1ddenoLVL','dw1ddenoDEN'} ,
  236.    case {'dw2dcompGBL','dw2dcompLVL','dw2ddenoLVL'} ,
  237.    case {'wp1dcompGBL','wp2dcompGBL'} ,
  238.    case {'wp1ddenoGBL','wp2ddenoGBL'} ,
  239.      switch meth
  240.        case 'sqtwologswn' , meth = 'sqtwolog'; scal = 'sln';
  241.        case 'sqtwologuwn' , meth = 'sqtwolog'; scal = 'one';
  242.      end
  243.    case 'sw1ddenoLVL' ,
  244.    case 'sw2ddenoLVL' ,
  245.    otherwise ,  errargt(mfilename,'invalid option','msg'); error('*');
  246. end
  247. flgTYPE = option(1:2);
  248. flgDIM  = str2num(option(3));
  249. flgTOOL = option(5:8);
  250. flgMODE = option(9:11);
  252. if ~isequal(meth,'rem_n0')
  253.     switch option
  254.        case {'dw1dcompGBL','dw1dcompLVL','dw1ddenoLVL','dw1ddenoDEN'}
  255.          level = length(varargin{3})-2;
  257.        case {'dw2dcompGBL','dw2dcompLVL','dw2ddenoLVL'}
  258.          level = size(varargin{3},1)-2;
  260.        case {'wp1dcompGBL','wp1ddenoGBL','wp2dcompGBL','wp2ddenoGBL'}
  261.          level = treedpth(varargin{2});
  263.        case 'sw1ddenoLVL' , level = size(varargin{2},1)-1;
  264.        case 'sw2ddenoLVL' , level = (size(varargin{2},3)-1)/3;
  265.     end
  266. else
  267.    if length(varargin)>2 , level = varargin{3}; end
  268. end
  270. switch option
  271.    case 'sw1ddenoLVL'
  272.        tmp = varargin{2};
  273.        varargin{4} = varargin{3};
  274.        varargin{2} = [];
  275.        varargin{3} = size(tmp,2);
  276.        for k=1:level
  277.            cfs  = tmp(k,1:2^k:end);
  278.            varargin{2} = [cfs , varargin{2}];
  279.            varargin{3} = [length(cfs) , varargin{3}];
  280.        end
  281.        cfs = tmp(level+1,1:2^level:end);
  282.        varargin{2} = [cfs , varargin{2}];
  283.        varargin{3} = [length(cfs) , varargin{3}];
  285.    case 'sw2ddenoLVL'
  286.        tmp = varargin{2};
  287.        varargin{4} = varargin{3};
  288.        varargin{2} = [];
  289.        varargin{3} = size(tmp(:,:,1));
  290.        for k=1:level
  291.            cfs  = tmp(1:2^k:end,1:2^k:end,3*k);
  292.            varargin{2} = [cfs(:)' , varargin{2}];
  293.            cfs  = tmp(1:2^k:end,1:2^k:end,2*k);
  294.            varargin{2} = [cfs(:)' , varargin{2}];
  295.            cfs  = tmp(1:2^k:end,1:2^k:end,1*k);
  296.            varargin{2} = [cfs(:)' , varargin{2}];
  297.            varargin{3} = [size(cfs) ; varargin{3}];
  298.        end
  299.        cfs = tmp(1:2^level:end,1:2^level:end,end);
  300.        varargin{2} = [cfs(:)' , varargin{2}];
  301.        varargin{3} = [size(cfs) ; varargin{3}];
  302. end
  304. switch flgTOOL
  305.   %============================= COMPRESSION ==============================%
  306.   case 'comp'
  307.     switch flgMODE
  308.       case 'GBL'
  309.         switch meth
  310.            case 'rem_n0'
  311.              % sig = varargin{2};
  312.              %----------------------
  313.              varargout{1} = remNearZero(flgTOOL,flgTYPE,varargin{2});
  315.            case {'bal_sn','sqrtbal_sn'}
  316.              % coefs = varargin{2};
  317.              % sizes = varargin{3};
  318.              %------------------------------
  319.              % tree or wptree = varargin{2};
  320.              % data = varargin{3};
  321.              %-------------------------------
  322.              if length(varargin)<3 , varargin{3} = []; end 
  323.              [valTHR,maxTHR,thresVALUES,rl2SCR,n0SCR] = ...
  324.                    balanceSparsityNorm(meth,flgTYPE,varargin{2:3});
  325.              varargout = {valTHR,maxTHR,thresVALUES,rl2SCR,n0SCR};
  326.         end
  328.       case 'LVL'
  329.         switch meth
  330.            case {'scarcehi','scarceme','scarcelo'}
  331.              % coefs = varargin{2};
  332.              % sizes = varargin{3};
  333.              % alfa  = varargin{4};
  334.              %----------------------
  335.              varargout{1} = scarceStrategies(meth,flgDIM,varargin{2:4});
  337.            case 'rem_n0'
  338.              % sig = varargin{2};
  339.              % lev = varargin{3};
  340.              %-------------------
  341.              valTHR = remNearZero(flgTOOL,flgTYPE,varargin{2});
  342.              varargout{1} = expandTHR(valTHR,flgDIM,level);
  344.            case {'bal_sn','sqrtbal_sn'}
  345.              % coefs = varargin{2};
  346.              % sizes = varargin{3};
  347.              %----------------------
  348.              valTHR = balanceSparsityNorm(meth,flgTYPE,varargin{2:3});
  349.              varargout{1} = expandTHR(valTHR,flgDIM,level);
  350.        end
  351.     end
  352.   %========================================================================%
  353.  
  354.   %============================= DE-NOISING ===============================%
  355.   case 'deno'
  356.     switch flgMODE
  357.       case 'GBL'        % WP only.
  358.         % tree or wptree = varargin{2};
  359.         % data = varargin{3};
  360.         %-----------------------------
  361.         if length(varargin)==2 , varargin{3} = []; end
  362.         switch meth
  363.            case 'sqtwolog'
  364.              [valTHR,maxTHR,cfs] = fixedFormWP(flgDIM,varargin{2:3},scal);
  365.              varargout = {valTHR,maxTHR,cfs};
  367.            case {'bal_sn','sqrtbal_sn'}
  368.              [valTHR,maxTHR,thresVALUES] = ...
  369.                    balanceSparsityNorm(meth,flgTYPE,varargin{2:3});
  370.              varargout = {valTHR,maxTHR,thresVALUES};
  372.            case {'penalhi','penalme','penallo'}
  373.              [valTHR,maxTHR,cfs] = WPpenalStrategies(meth,flgDIM,varargin{2:3});
  374.              varargout = {valTHR,maxTHR,cfs};
  375.         end
  377.       case 'LVL'
  378.         switch meth
  379.            case 'sqtwolog', % DW & SW only.
  380.              switch flgDIM
  381.                 case 1 , varargout{1} = fixedForm1D(varargin{2:4});
  382.                 case 2 , varargout{1} = fixedForm2D(varargin{2:4},level);
  383.              end
  385.            case {'rigrsure','heursure','minimaxi'}  % DW1D & SW1D only.
  386.              % coefs = varargin{2};
  387.              % sizes = varargin{3};
  388.              % scal  = varargin{4};
  389.              %----------------------
  390.              coefs = detcoef(varargin{2:3},'all');
  391.              sigma = sigmaHAT(varargin{4},coefs);
  392.              varargout{1} = getTHR(meth,sigma,coefs);
  394.            case {'penalhi','penalme','penallo'}
  395.              % coefs = varargin{2};
  396.              % sizes = varargin{3};
  397.              % alfa  = varargin{4};
  398.              %----------------------
  399.              valTHR = penalStrategies(meth,flgDIM,varargin{2:4});
  400.              varargout{1} = expandTHR(valTHR,flgDIM,level);
  402.            case {'scarcehi','scarceme','scarcelo'}
  403.              % coefs = varargin{2};
  404.              % sizes = varargin{3};
  405.              % alfa  = varargin{4};
  406.              %----------------------
  407.              varargout{1} = scarceStrategies(meth,flgDIM,varargin{2:4});
  409.            case {'bal_sn','sqrtbal_sn'}
  410.              % coefs = varargin{2};
  411.              % sizes = varargin{3};
  412.              %----------------------
  413.              valTHR = balanceSparsityNorm(meth,flgTYPE,varargin{2:3});
  414.              varargout{1} = expandTHR(valTHR,flgDIM,level);
  415.         end
  418.       case 'DEN'  % estimation de densite
  419.         switch meth
  420.            case 'globalth', 
  421.              % coefs = varargin{2};
  422.              % sizes = varargin{3};
  423.              %----------------------
  424.              varargout{1} = GlobDens(varargin{2:3});
  426.            case {'bylevth1'}
  427.              % coefs = varargin{2};
  428.              % sizes = varargin{3};
  429.              %----------------------
  430.              varargout{1} = LvldDens(varargin{2:3},1);
  431.              
  432.            case {'bylevth2'}
  433.              % coefs = varargin{2};
  434.              % sizes = varargin{3};
  435.              %----------------------
  436.              varargout{1} = LvldDens(varargin{2:3},2);
  438.            case {'bylevsth'}
  439.              % coefs = varargin{2};
  440.              % sizes = varargin{3};
  441.              % alfa  = varargin{4};
  442.              %----------------------
  443.              varargout{1} = LvdsDens(varargin{2:4});
  444.         end
  445.     end
  446.   %========================================================================%
  448. end
  451. %=============================================================================%
  452. % INTERNAL FUNCTIONS
  453. %=============================================================================%
  454. %-----------------------------------------------------------------------------%
  455. function [valTHR,maxTHR,thresVALUES,rl2SCR,n0SCR] = ...
  456.                          balanceSparsityNorm(meth,flgTYPE,A,B)
  458. switch flgTYPE
  459.   case {'dw','sw'}
  460.     % coefs = A;
  461.     % sizes = B;
  462.     %-----------
  463.     [thresVALUES,rl2SCR,n0SCR,imin] = wcmpscr(A,B);
  465.  case 'wp'
  466.     % WP_Tree = varargin{2};
  467.     % WP_Data = varargin{3};
  468.     %-----------------------
  469.     if isa(A,'wptree')
  470.         [thresVALUES,rl2SCR,n0SCR,imin] = wpcmpscr(A,A);
  471.     else
  472.         [thresVALUES,rl2SCR,n0SCR,imin] = wpcmpscr(A,B);
  473.     end
  474. end
  475. valTHR = thresVALUES(imin);
  476. maxTHR = thresVALUES(end);
  477. if isequal(meth,'sqrtbal_sn') , valTHR = min(sqrt(valTHR),maxTHR); end
  478. %-----------------------------------------------------------------------------%
  479. function valTHR = remNearZero(flgTOOL,flgTYPE,X)
  481. switch flgTOOL
  482.   case 'comp' , argTOOL = 'cmp';
  483.   case 'deno' , argTOOL = 'den';
  484. end
  485. switch flgTYPE
  486.   case 'dw' , argTYPE = 'wv';
  487.   case 'wp' , argTYPE = 'wp';
  488. end
  489. valTHR = ddencmp(argTOOL,argTYPE,X);
  490. %-----------------------------------------------------------------------------%
  491. function valTHR = scarceStrategies(meth,flgDIM,coefs,sizes,alfa)
  493. switch flgDIM
  494.   case 1 , M = [1 , 1.5 ,   2] * sizes(1);
  495.   case 2 , M = 4 * [1 , 4/3 , 8/3] * prod(sizes(1,:));
  496. end
  497. switch meth
  498.     case 'scarcehi' , M = M(1);
  499.     case 'scarceme' , M = M(2);
  500.     case 'scarcelo' , M = M(3);
  501. end
  502. if flgDIM==1
  503.     valTHR = wdcbm(coefs,sizes,alfa,M);
  504. else
  505.     valTHR = wdcbm2(coefs,sizes,alfa,M);
  506. end
  507. %-----------------------------------------------------------------------------%
  508. function valTHR = penalStrategies(meth,flgDIM,coefs,sizes,sliBMVal)
  510. switch flgDIM
  511.   case 1
  512.     sigma = wnoisest(coefs,sizes,1);
  513.   case 2 
  514.     det   = detcoef2('compact',coefs,sizes,1);
  515.     sigma = wnoisest(det);
  516. end
  517. switch meth
  518.   case 'penalhi' , alfa = 5*(3*sliBMVal+1)/8;
  519.   case 'penalme' , alfa = (sliBMVal+5)/4;
  520.   case 'penallo' , alfa = (sliBMVal+3)/4;
  521. end
  522. valTHR = wbmpen(coefs,sizes,sigma,alfa);
  523. %-----------------------------------------------------------------------------%
  524. function [valTHR,maxTHR,cfs] = WPpenalStrategies(meth,flgDIM,wpt,wpd)
  526. sliBMVal = 3;
  527. switch meth
  528.   case 'penalhi' , alfa = 5*(3*sliBMVal+1)/8;
  529.   case 'penalme' , alfa = (sliBMVal+5)/4;
  530.   case 'penallo' , alfa = (sliBMVal+3)/4;
  531. end
  533. % Old Version
  534. %-------------
  535. if ~isa(wpt,'wptree') , wpt = convv1v2(wpt,wpd); end
  537. % Compute sigma.
  538. %---------------
  539. depth = treedpth(wpt);
  540. if depth==0 , valTHR = 0; return; end
  541. switch flgDIM
  542.   case 1
  543.     cD1 = wpcoef(wpt,[1,1]);
  544.     sigma = wnoisest(cD1);
  545.   case 2
  546.     cH1 = wpcoef(wpt,[1,1]);
  547.     cV1 = wpcoef(wpt,[1,2]);
  548.     cD1 = wpcoef(wpt,[1,3]);
  549.     sigma = wnoisest([cH1(:)',cV1(:)',cD1(:)']);
  550. end
  551. if isa(wpt,'wptree')
  552.     cfs = read(wpt,'allcfs');
  553. end
  554. valTHR = wpbmpen(wpt,sigma,alfa);
  555. maxTHR = max(abs(cfs));
  556. valTHR = min(valTHR,maxTHR);
  557. %-----------------------------------------------------------------------------%
  558. function valTHR = expandTHR(valTHR,flgDIM,nbLEV)
  559. switch flgDIM
  560.    case 1 , nbDIR = 1;
  561.    case 2 , nbDIR = 3;
  562. end
  563. valTHR = valTHR*ones(nbDIR,nbLEV);
  564. %-----------------------------------------------------------------------------%
  565. function s = sigmaHAT(scal,coefs);
  567. level = length(coefs);
  568. switch scal
  569.   case 'one' , s = ones(1,level);
  570.   case 'sln' , s = ones(1,level)*wnoisest(coefs{1});
  571.   case 'mln' , s = wnoisest(coefs);
  572. end
  573. %-----------------------------------------------------------------------------%
  574. function thr = getTHR(meth,s,coefs)
  576. switch meth
  577.     case 'minimaxi'
  578.       nbcfs = 0;
  579.       for k=1:length(s) , nbcfs = nbcfs+length(coefs{k}); end
  580.       if nbcfs <= 32
  581.           thr = 0*s;
  582.       else
  583.           thr = (0.3936 + 0.1829*(log(nbcfs)/log(2)))*s;
  584.       end
  586.     case {'rigrsure','heursure'}
  587.       thr = zeros(size(s));
  588.       for k=1:length(s)
  589.           mk = max(coefs{k});
  590.           if (mk<sqrt(eps)) | (s(k)<sqrt(eps)*mk)
  591.               thr(k) = 0;
  592.           else
  593.               thr(k) = sureTHR(meth,coefs{k}/s(k));
  594.           end
  595.       end
  596.       thr = thr.*s;
  597. end
  598. %-----------------------------------------------------------------------------%
  599. function thr = sureTHR(meth,x)
  601. x = x(:)';
  602. n = length(x);
  603. switch meth
  604.   case 'rigrsure'
  605.     sx2 = sort(abs(x)).^2;
  606.     risks = (n-(2*(1:n))+(cumsum(sx2)+(n-1:-1:0).*sx2))/n;
  607.     [risk,best] = min(risks);
  608.     thr = sqrt(sx2(best));
  610.   case 'heursure'
  611.     hthr = sqrt(2*log(n));
  612.     eta = (norm(x).^2-n)/n;
  613.     crit = (log(n)/log(2))^(1.5)/sqrt(n);
  614.     if eta < crit
  615.         thr = hthr;
  616.     else
  617.         thr = min(sureTHR('rigrsure',x),hthr);
  618.     end
  619. end
  620. %-----------------------------------------------------------------------------%
  621. function valTHR = fixedForm1D(coefs,sizes,scal)
  623. coefs  = detcoef(coefs,sizes,'all');
  624. sigma  = sigmaHAT(scal,coefs);
  625. nbcfs = 0;
  626. for k=1:length(coefs)
  627.     nbcfs = nbcfs+length(coefs{k});
  628. end
  629. valTHR = sqrt(2*log(nbcfs))*sigma;
  630. %-----------------------------------------------------------------------------%
  631. function valTHR = fixedForm2D(coefs,sizes,scal,level)
  633. strDET = ['h','d','v'];
  634. s = ones(3,level);
  635. switch scal
  636.   case 'one'
  637.   case 'sln'
  638.     det  = detcoef2('compact',coefs,sizes,1);
  639.     s = wnoisest(det) * s;  
  640.   case 'mln'
  641.     for k = 1:level
  642.       det = detcoef2('compact',coefs,sizes,k);
  643.       s(:,k) = wnoisest(det) * ones(3,1);
  644.     end
  645. end
  646. for d = 1:3
  647.   for k = 1:level
  648.     det = detcoef2(strDET(d),coefs,sizes,k);
  649.     univTHR     = sqrt(2*log(prod(size(det))));
  650.     valTHR(d,k) = univTHR*s(d,k);
  651.   end
  652. end
  653. %-----------------------------------------------------------------------------%
  654. function [valTHR,maxTHR,cfs] = fixedFormWP(flgDIM,A,B,scal)
  656. order = treeord(A);
  657. nodes = [2:order]'; % nodes for details of level 1.
  658. det = [];
  659. if isa(A,'wptree')
  660.    for k =1:length(nodes)
  661.        tmp = wpcoef(A,nodes(k));
  662.        det = [det , tmp(:)'];
  663.    end
  664.    cfs = read(A,'allcfs');
  665. else
  666.    for k =1:length(nodes)
  667.        tmp = wpcoef(A,B,nodes(k));
  668.        det = [det , tmp(:)'];
  669.    end
  670.    cfs = wdatamgr('rallcfs',B);
  671. end
  672. univTHR = sqrt(2*log(length(det)));
  673. switch scal
  674.   case 'one' , s = 1;
  675.   case 'sln' , s = wnoisest(det);
  676. end
  677. valTHR = s*univTHR;
  678. maxTHR = max(abs(cfs));
  679. valTHR = min(valTHR,maxTHR);
  680. %-----------------------------------------------------------------------------%
  681. function valTHR = GlobDens(coefs,sizes)
  683. n = sizes(end);
  684. J = size(sizes,2)-2;
  685. coefs = coefs(sizes(1)+1:end);
  686. valTHR = max(abs(coefs))*log(n)/sqrt(n);
  687. valTHR = expandTHR(valTHR,1,J);
  688. %-----------------------------------------------------------------------------%
  689. function valTHR = LvldDens(coefs,sizes,flag)
  691. n = sizes(end);
  692. J = size(sizes,2)-2;
  693. valTHR = zeros(1,J);
  694. for j=1:J
  695.     d = detcoef(coefs,sizes,j);
  696.     switch flag
  697.         case 1 , valTHR(j) = 0.4*max(abs(d));
  698.         case 2 , valTHR(j) = 0.8*max(abs(d));
  699.     end
  700. end
  701. %-----------------------------------------------------------------------------%
  702. function valTHR = LvdsDens(coefs,sizes,alfa)
  704. n = sizes(end);
  705. J = size(sizes,2)-2;
  706. valTHR = zeros(1,J);
  707. for j=1:J
  708.     d = detcoef(coefs,sizes,j);
  709.     valTHR(j) = max(abs(d));
  710. end
  711. valTHR = valTHR * (alfa/(5-sqrt(eps)));
  712. %-----------------------------------------------------------------------------%
  713. function medad = med_ad(x)
  715. medad = median(abs(x-median(x)));
  716. %-----------------------------------------------------------------------------%
  717. %=============================================================================%