www.gusucode.com > wavelet工具箱matlab源码程序 > wavelet/wavelet/utthrset_cmd.m
function varargout = utthrset_cmd(coefs,longs,nb_IntVal) %UTTHRSET_CMD Define interval dependent thresholds. % For a given wavelet decomposition (C,L) (see WAVEDEC), % UTTHRSET_CMD computes the intervals and the thresholds % used for an interval dependent denoising. % % [THRPAR,INT_DepThr] = UTTHRSET_CMD(C,L) returns % THRPAR is a cell array which length is the level of the % wavelet decomposition. For the level k, THRPAR{k} is an L-by-3 % array such that THRPAR{k}(i,1) and THRPAR{k}i,2) are the end % points of the ith interval and THRPAR{k}(i,3) is the corresponding % threshold. % The variable int_DepThr contains the interval locations % and the threshold values for a number of intervals from 1 to 6. % THRPAR corresponds to the best choice of the number of intervals. % % [THRPAR,INT_DepThr] = UTTHRSET_CMD(C,L,NBint) let you choose % the number of intervals. % M. Misiti, Y. Misiti, G. Oppenheim, J.M. Poggi 01-Oct-2008. % Last Revision: 13-Oct-2008. % Copyright 1995-2008 The MathWorks, Inc. details = wrepcoef(coefs,longs); maxTHR = max((abs(details)),[],2)'; level = size(details,1); nb_Max_Int = 6; if nargin<3 , nb_IntVal = NaN; end % Computing Intervals. %===================== % Extract the detail of order 1. %------------------------------- det = details(1,:); xdata = 1:length(det); % Replacing 2% of biggest values of by the mean. %----------------------------------------------- x = sort(abs(det)); v2p100 = x(fix(length(x)*0.98)); det(abs(det)>v2p100) = mean(det); lenDet = length(det); % Finding breaking points. %------------------------- d = 10; if lenDet>1024 ratio = ceil(lenDet/1024); [Rupt_Pts,nb_Opt_Rupt,Xidx] = wvarchg(det(1:ratio:end),nb_Max_Int,d); Xidx = min(ratio*Xidx,lenDet); else [Rupt_Pts,nb_Opt_Rupt,Xidx] = wvarchg(det,nb_Max_Int,d); end nb_Opt_Int = nb_Opt_Rupt+1; if isnan(nb_IntVal) , nb_IntVal = nb_Opt_Int; end % Computing denoising structure. %------------------------------- % Ensure that there are no zeros on the lower diagonal entries. [row,col] = find(~Xidx); idx = find(~Xidx); lowertriang = row>=col; Xidx(idx(lowertriang>0)) = 1; Xidx = [zeros(size(Xidx,1),1) Xidx]; norma = sqrt(2)*thselect(det,'minimaxi'); % sqrt(2) comes from the fact that if x is a white noise % of variance 1 the reconstructed detail_1 of x is of % variance 1/sqrt(2) int_DepThr = cell(1,nb_Max_Int); for nbint = 1:nb_Max_Int for j = 1:nbint sig = median(abs(det(Xidx(nbint,j)+1:Xidx(nbint,j+1))))/0.6745; thr = norma*sig; int_DepThr{nbint}(j,:) = ... [Xidx(nbint,j) , Xidx(nbint,j+1), thr]; end int_DepThr{nbint}(1,1) = 1; int_DepThr{nbint}(:,[1 2]) = xdata(int_DepThr{nbint}(:,[1 2])); end int_DepThr_Cell = {int_DepThr,nb_Opt_Int}; thrParams = cell(1,level); intervals = int_DepThr{nb_IntVal}; for k=1:level thrPAR = intervals; TMP = min(thrPAR(:,3),maxTHR(k)); thrPAR(:,3) = TMP; thrParams{k} = thrPAR; end varargout = {thrParams,int_DepThr_Cell};