www.gusucode.com > wavelet工具箱matlab源码程序 > wavelet/wavelet/wtmm.m
function varargout = wtmm(x,varargin) %Wavelet Transform Modulus Maxima % HEXP = WTMM(X) returns an estimate of the global Holder exponent, HEXP, % for the real-valued 1-D input signal, X. X must have at least 128 % samples. HEXP is estimated by regressing the mean log of the CWT maxima % by scale on the log of the wavelet scales. The CWT is obtained using the % 2nd derivative of Gaussian wavelet with 10 voices per octave. An % equivalent syntax is: HEXP = WTMM(X,'ScalingExponent','global') % % [HEXP,TAUQ] = WTMM(X) returns an estimate of the the partition function % scaling exponents, TAUQ, for the linearly spaced Q-moments, -2:0.1:2. % An equivalent syntax is [HEXP,TAUQ] = WTMM(X,'ScalingExponent','global') % % [HEXP,TAUQ] = WTMM(X,'MinRegressionScale',SCALE) uses only scales greater % than or equal to SCALE in the estimate of the Holder exponent. SCALE is a % positive scalar greater than or equal to 4. The MinRegressionScale % name-value pair only affects global estimates. If unspecified, % MinRegressionScale defaults to 4. There must be at least two scales to % use in the estimation of the global Holder exponent. % % [HEXP,TAUQ,STRUCTFUNC] = WTMM(...) returns the multiresolution structure % functions, STRUCTFUNC, for the global estimates. STRUCTFUNC is a % structure array containing the following fields: % % Tq: matrix of multiresolution quantities that depend jointly on time % and scale. Tq provides measurements of the input X at various scales. % Scaling phenomena in X imply a power-law relationship between the % moments of Tq and scale. For WTMM, Tq is a Ns-by-44 matrix where Ns is % the number of scales. The first 41 columns of Tq constitute the % scaling exponent estimates for each of the q-moments -2:0.1:2 by scale. % The last three columns correspond to the 1st, 2nd, and 3rd order % cumulants respectively by scale. % % weights: Ns-by-1 vector of weights used in the regression estimates. % The weights correspond to the number of wavelet maxima at each scale. % % logscales: Ns-by-1 vector with the base-2 logarithm of the scales used % as predictors in the regression. % % [LOCALHEXP,WT,WAVSCALES] = WTMM(X,'ScalingExponent','local') returns the % continuous wavelet transform WT, the local Holder exponent estimates, and % the scales used in the CWT, WAVSCALES. WT is a numel(WAVSCALES)-by-N % matrix where N is the length of the input signal X. LOCALHEXP is a M-by-2 % array containing the sample of the local singularity in the first column % and the estimated local Holder exponent in the second column. If there % are no maxima lines that converge to the finest scale in the wavelet % transform, LOCALEXP is an empty array. % % [...] = WTMM(...,'VoicesPerOctave',NV) uses NV voices per octave to % determine the CWT, maxima lines, and fractal estimates. NV is an even % integer between 8 and 32. If you do not specify the number of voices per % octave, NV defaults to 10. % % [...] = WTMM(...,'NumOctaves',NO) uses NO octaves to determine the CWT, % maxima lines, and fractal estimates. NO is an integer greater than or % equal to 4. WTMM limits the number of octaves to be less than or equal to % floor(log2(N/(3*sqrt(1.1666)))) where N is the length of the input data. % The factor sqrt(1.1666) is the standard deviation of the 2nd derivative % of Gaussian wavelet. NO defaults to the minimum of 7 and % floor(log2(N/(3*sqrt(1.1666)))). If you specify NO greater than the % maximum number of octaves, WTMM uses the maximum supported number. % % WTMM(...,'ScalingExponent','local') with no output arguments plots the % wavelet maxima lines in the current figure. Estimates of the local Holder % exponents are displayed in a table to the right of the plot. % % % Example 1: Obtain local Holder exponents for cusp signal % % Note at sample 241, this signal has Holder exponent of 0.5. At % % sample 803, the signal has a Holder exponent of 0.3. % load cusp; % wtmm(cusp,'ScalingExponent','local'); % % % Example 2: Obtain local Holder exponents for two delta functions % % located at sample 200 and sample 500. The local Holder exponents % % of the delta functions should be -1 at these sample values. % x = zeros(1e3,1); % x([200 500]) = 1; % wtmm(x,'ScalingExponent','local','NumOctaves',6); % % % Example 3: Estimate global Holder exponent for Brownian motion. % % In theory, the Holder exponent should be around 0.5 for this % % monofractal process. % rng(100); % x = cumsum(randn(2^15,1)); % hexp = wtmm(x); % % See also DWTLEADER, WFBM narginchk(1,7); nargoutchk(0,3); % Determine the number of outputs nOutputs = nargout; % Parse the input p = inputParser; defaultScaling = 'global'; expectedScaling = {'global','local'}; defaultnv = 10; defaultno = 7; defaultminRegScale = 4; addRequired(p,'x',@(x)validateattributes(x,{'numeric'},{'real','finite',... 'vector'})); addParameter(p,'ScalingExponent',defaultScaling,... @(x)any(validatestring(x,expectedScaling))); addParameter(p,'VoicesPerOctave',defaultnv,@(x)validateattributes(x,... {'numeric'},{'scalar','integer','even','>=',8,'<=',32})); addParameter(p,'NumOctaves',defaultno,@(x)validateattributes(x,... {'numeric'},{'scalar','integer','>=',4,'finite'})); addParameter(p,'MinRegressionScale',defaultminRegScale,... @(x)validateattributes(x,{'numeric'},{'scalar','positive','>=',2})); % Parse the inputs to the function. Only required is the time series parse(p,x,varargin{:}); scalexp = p.Results.ScalingExponent; % Pad for global estimates if strcmpi(scalexp,'global') pad = true; else pad = false; end numoct = p.Results.NumOctaves; nv = p.Results.VoicesPerOctave; x = p.Results.x; j1 = p.Results.MinRegressionScale; % Detrend signal x = detrend(x,0); x = x(:)'; %Only supported 2nd derivative of Gaussian wavname = 'dergauss2'; % Record original signal length Norig = numel(x); if Norig < 128 error(message('Wavelet:mfa:WTMMLength')); end % 2 sigma_t for the 2nd derivative of Gaussian sigma2 = 2*sqrt(1.1666); % Determine maximum scale. For WTMM, keep the maximum scale from going % too large numoct = min(numoct,floor(log2(Norig/(3/2*sigma2)))); n = Norig; if pad padvalue = floor(Norig/2); x =[fliplr(x(1:padvalue)) x x(end:-1:end-padvalue+1)]; % Length of data plus any extension n = length(x); end ds = 1/nv; a0 = 2^ds; % Create scales wavscales = a0.^(1:numoct*nv); % Find the minimum scale greater than or equal to the MinRegressionScale % determine the cone of influence at finest scale % we will only trust local maxima that terminate inside the cone of % influence at the finest scale conofinfb = round(sigma2*wavscales(1)); % beginning and end of cone of influence conofinfb = [conofinfb Norig-conofinfb]; NbSc = numel(wavscales); % Create frequencies for real wavelet transform omega = (1:fix(n/2)); omega = omega.*(2*pi)/n; omega = [0, omega, -omega(fix((n-1)/2):-1:1)]; % Obtain the Fourier transform of the data and the wavelet filter bank f = fft(x); psift = wavelet.internal.waveft(wavname,omega,wavscales); % Obtain the CWT coefficients -- real-valued because the wavelet is % real-valued and the signal is real-valued cwtcfs = ifft(repmat(f,NbSc,1).*psift,[],2,'symmetric'); % Remove padding if necessary if pad cwtcfs = cwtcfs(:,padvalue+1:padvalue+Norig); end % Determine the maximum map. This is the input to the structure function % calculation and is used in the formation of the wavelet skeleton % We need this map for both the local and global exponents. maxmap = wtmaxmap(cwtcfs',pad); cfsmask = fliplr(maxmap)'; ncount = sum(cfsmask,2); % Mask the CWT coefficients using the local maxima wtmask = cwtcfs.*cfsmask; % If the scaling exponent calculation is global, we need % to compute the global estimates if strcmpi(scalexp,'global') % We keep the index of the scale, not the actual scale. j1 = find(wavscales>=j1,1,'first'); % Find terminal scale for regression J = find(ncount>=6,1,'last'); % Check that there are at least two scales for the regression estimates if J < j1+1 error(message('Wavelet:mfa:RegressionLevels')); end if ~any(wtmask(:)) error(message('Wavelet:mfa:AllZeroCWT')); end Nest = size(wtmask,1); param.q = -2:0.1:2; Nq = numel(param.q); param.cumulant = 3; zetaq = zeros(Nq,Nest); Cp = zeros(param.cumulant,Nest); for jj = 1:Nest [zetaq(:,jj),~, ~, Cp(:,jj)] = ... wavelet.internal.mfstructfunctions(abs(wtmask(jj,:)),param); end Y = [zetaq; Cp*log2(exp(1))]; xj = log2(wavscales); J = numel(wavscales); % Return scaling exponent results and structfunc.Tq = Y(:,j1:J)'; structfunc.weights = ncount(j1:J); structfunc.logscales = xj(j1:J)'; % Create design matrix X = ones(length(structfunc.logscales),2); X(:,2) = structfunc.logscales; % Weighted linear regression with multiple response variables betahat = lscov(X,structfunc.Tq,structfunc.weights); betahat = betahat(2,:); zq = betahat(1:Nq); cp = betahat(Nq+1:end); % If scaling exponents are local elseif strcmpi(scalexp,'local') maxchains = wtskeleton(maxmap,nv,conofinfb); if ~isempty(maxchains) hlocal = zeros(numel(maxchains),1); tsamp = zeros(numel(maxchains),1); for kk = 1:numel(maxchains) [hexp,samp] = wavelet.internal.localholderexp(kk,maxchains,cwtcfs,wavscales); hlocal(kk) = hexp; tsamp(kk) = samp; end localHdata = [tsamp hlocal]; else localHdata = []; end end if nOutputs == 0 && strcmpi(scalexp,'local') plotwavskeleton(cwtcfs,maxchains,wavscales,localHdata) end if strcmpi(scalexp,'global') varargout{1} = cp(1); varargout{2} = zq; varargout{3} = structfunc; elseif strcmpi(scalexp,'local') && nOutputs > 0 varargout{1} = localHdata; varargout{2} = cwtcfs; varargout{3} = wavscales; end function maxmap = wtmaxmap(cwtcfs,pad) % Get time and scale dimension of CWT % here n is the number of samples and nscale is the number of scales [n,nscale] = size(cwtcfs); maxmap = zeros(n,nscale); t = 1:n; tplus = circshift(t,1,2); tminus = circshift(t,-1,2); wtmag = fliplr(abs(cwtcfs)); for k = 1:nscale localmax = ... wtmag(:,k) >= wtmag(tplus,k) & wtmag(:,k) >= wtmag(tminus,k); y = localmax.* wtmag(:,k); nnzy = y(y>0); if pad medy = median(nnzy); maxmap(:,k) = (y>medy); else maxmap(:,k) = (y>0); end end function chains = wtskeleton(maxmap,nv,conofinfb) % Chain together Ridges of Wavelet Transform % A chain is a list of maxima at essentially the same position % across a range of scales. % It is identified from the maxmap data structure output by WTMM % by finding a root at coarse scales and identifying the closest % maxima at the next finest scale. % n is the number of time points % nscale is the number of scales % first column of maxmap is the coarsest scale [n,nscale] = size(maxmap); nchain = 0; chains = zeros(size(maxmap)); count = 0; while any(any(maxmap)) [i,j] = find(maxmap); % beginning chain formation iscale = j(1); ipos = i(1); nchain = nchain+1; % at chain n and scale iscale, ipos gets a time position chains(nchain,iscale) = ipos; % zero out the maxmap maxmap(ipos,iscale) = 0; count = count+1; while(iscale < nscale) iscale = iscale+1; j = find(maxmap(:,iscale))'; circdist = min([ abs(j-ipos) ; abs(j-ipos+n); abs(j-ipos-n) ]); [~,pos] = min(circdist); if ~isempty(pos) ipos = j(pos); chains(nchain,iscale) = ipos; maxmap(ipos,iscale) = 0; count = count+1; else iscale = nscale; end end end chaincell = cell(size(chains,1)); for kk = 1:size(chains,1) [~,j,v] = find(chains(kk,:)); chaincell{kk} = [j' v']; end % retain only nonempty cells chains = chaincell(~cellfun(@isempty,chaincell)); %chains = chaincell; chainLength = cellfun(@(x)numel(x(:,1)),chains); chainEnd = cellfun(@(x)x(end,1),chains); idxEnd = chainEnd == nscale; idxLen = chainLength>=nv; % what chains can we estimate Holder exponents for chains = chains(idxEnd & idxLen); rangechains = cellfun(@(x)max(x(:,2))-min(x(:,2)),chains); idx = rangechains<=floor(n/4); chains = chains(idx); %Find chains that terminate in the cone of influence and remove those cET = cellfun(@(x)(x(end,2)),chains); insidecf = cET>conofinfb(1) & cET <conofinfb(2); chains = chains(insidecf); % Sort by chain end time chainEndTime = cellfun(@(x)(x(end,2)),chains); [~,idxsort] = sort(chainEndTime); chains = chains(idxsort); function plotwavskeleton(wt,chains,wavscales,hexpdata) % Plot local maxima lines for zero output case for kk = 1:numel(chains) chains{kk}(:,1) = flipud(chains{kk}(:,1)); end f = gcf; clf(f,'reset'); ax1 = axes(f); f.Visible = 'off'; f.Units = 'normalized'; f.Position = [0.1 0.2 0.4 0.5]; movegui(f,'center'); f.Units = 'Pixels'; % Layout figure hl = wavelet.internal.layout.GridBagLayout(f); hl.add(ax1, 1, 1,'Fill','both'); imagesc(ax1,1:size(wt,2),log2(wavscales),abs(wt)); hold(ax1,'on'); for kk = 1:numel(chains) scales = chains{kk}(:,1); plot(ax1,chains{kk}(:,2),log2(wavscales(scales)),'w','linewidth',0.1); end ylim(ax1,[min(log2(wavscales)) max(log2(wavscales))]); xlim(ax1,[1 size(wt,2)]); set(ax1,'ydir','reverse'); ylbl = getString(message('Wavelet:mfa:WTMMYlabel')); ylabel(ax1,ylbl); xlbl = getString(message('Wavelet:mfa:WTMMXlabel')); xlabel(ax1,xlbl); titl = getString(message('Wavelet:mfa:WTMMTitle')); title(ax1,titl); grid(ax1,'on'); hold(ax1,'off'); colnames = {getString(message('Wavelet:mfa:Table1stCol')); ... getString(message('Wavelet:mfa:Table2ndCol')) }; tbl = uitable('Parent',f, 'Data', hexpdata, 'ColumnName', colnames,... 'units','pixels','visible','off'); tbl.BackgroundColor = [0.94 0.94 0.94]; hl.add(tbl, 1, 2,'Fill', 'Both', 'MinimumWidth', tbl.Extent(3)); hl.HorizontalWeights = [1 0]; hl.setConstraints(1, 2, 'TopInset', 5, 'RightInset', 7, 'BottomInset', 5); hl.setConstraints(1,2, 'MinimumHeight', tbl.Extent(4), 'Fill', 'None',... 'Anchor','North'); tbl.Visible = 'on'; f.Visible = 'on'; f.NextPlot = 'replace';