www.gusucode.com > wavelet工具箱matlab源码程序 > wavelet/wavelet/dwtleader.m

    function [dh,h,cp,zq,leaders,structfunc] = dwtleader(x,varargin)
%Multifractal 1-D Wavelet Leader Estimates
% [DH,H] = DWTLEADER(X) returns the singularity spectrum, DH, and the
% Holder exponents, H, for the 1-D real-valued data, X. The singularity
% spectrum is estimated using structure functions determined for the
% linearly-spaced moments -5:5. The structure functions are computed based
% on the wavelet leaders obtained using the biorthogonal spline wavelet
% filter with 1 vanishing moment in the synthesis wavelet and 5 vanishing
% moments in the analysis wavelet ('bior1.5'). By default, multifractal
% estimates are derived from wavelet leaders at a minimum level of 3 and
% maximum level where there are at least 6 wavelet leaders. For the
% default wavelet and minimum regression level, you need a time series with
% at least 248 samples.
%
% [DH,H,CP] = DWTLEADER(X) returns the first three cumulants, CP. The first
% cumulant characterizes the local maximum of the singularity spectrum. For
% monofractal processes, the first cumulant measures the linear fit for the
% scaling exponents, while the second cumulant characterizes the departure
% from linearity.
%
% [DH,H,CP,TAUQ] = DWTLEADER(X) returns the scaling exponents for the
% linearly spaced moments -5:5.
%
% [DH,H,CP,TAUQ,LEADERS] = DWTLEADER(...) returns the wavelet leaders in
% LEADERS. LEADERS is a cell array with the i-th element containing the
% wavelet leaders at level i+1, or scale 2^(i+1). Wavelet leaders are not
% defined at level 1.
%
% [DH,H,CP,TAUQ,LEADERS,STRUCTFUNC] = DWTLEADER(...) returns the
% multiresolution structure functions, STRUCTFUNC. STRUCTFUNC is a
% structure with 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 DWTLEADER, Tq is a Ns-by-36 matrix where
%   Ns is the number of scales used in multifractal estimates. The first 11
%   columns of Tq are the scaling exponent estimates by scale for each of
%   the q-th moments, -5:5. The next 11 columns contain the singularity
%   spectrum estimates, DH, for each of the q-th moments. Columns 23 to 33
%   contain the Holder exponent estimates, H. The final three columns
%   contain the estimates for the 1st, 2nd, and 3rd order cumulants
%   respectively.
%
%   weights: Ns-by-1 vector of weights used in the regression. The
%   weights are the number of wavelet leaders by scale.
%
%   logscales: Ns-by-1 vector containing the base-2 logarithm of the scales
%   used as predictors in the regression.
%
% [...] = DWTLEADER(X,WNAME) uses the orthogonal or biorthogonal wavelet
% denoted by WNAME in the computation of the wavelet leaders and the
% fractal estimates. WNAME is a wavelet family short name and filter number
% recognized by the wavelet manager. You can query valid wavelet family
% short names using wavemngr('read'). To determine whether a particular
% wavelet is orthogonal or biorthogonal, you can use WAVEINFO with the
% wavelet family short name, waveinfo('db'), or use WAVEMNGR with the
% 'type' option for a specific wavelet, wavemngr('type','fk4'). A 1
% indicates an orthogonal wavelet and a 2 denotes a biorthogonal wavelet.
% The minimum-required data length depends on the wavelet filter and the
% levels used in the regression model.
%
% [...] = DWTLEADER(...,'MinRegressionLevel',MINLEV) uses only levels
% greater than or equal to MINLEV in the multifractal estimates. MINLEV
% is a positive integer greater than or equal to 2. DWTLEADER requires at
% least 6 wavelet leaders at the maximum level and two levels to be used in
% the multifractal estimates. The scale in the discrete wavelet transform
% corresponding to MINLEV is 2^MINLEV. If unspecified, MinRegressionLevel
% defaults to 3.
%
% [...] = DWTLEADER(...,'MaxRegressionLevel',MAXLEV) uses only levels less
% than or equal to MAXLEV in the multifractal estimates. MAXLEV is a
% positive integer greater than or equal to MINLEV+1. MAXLEV defaults to
% the largest level where there are at least 6 wavelet leaders. The scale
% in the discrete wavelet transform corresponding to MAXLEV is 2^MAXLEV.
% The MaxRegressionLevel name-value pair is intended for situations where
% you want to restrict the levels used in the regression to a value less
% than the default level. You can use the optional output argument,
% LEADERS, or the weights field of the optional output argument,
% STRUCTFUNC, to determine the number of wavelet leaders by level.
%
%   %Example 1:
%   %   Compute the singularity spectrum and cumulants for
%   %   a Brownian noise process.
%   rng(100);
%   x = cumsum(randn(2^15,1));
%   [dh,h,cp] = dwtleader(x);
%   plot(h,dh,'o-','MarkerFaceColor','b'); grid on;
%   title({'Singularity Spectrum'; ['First Cumulant ' num2str(cp(1))]});
%
%   %Example 2:
%   %   Compute the cumulants for a multifractal random walk. The
%   %   multifractal random walk is realization of a random process with a
%   %   theoretical first cumulant of 0.75 and second cumulant of -0.05.
%   %   The second cumulant value of -0.05 captures the fact that the
%   %   scaling exponents deviate from a linear function with slope 0.75.
%   load mrw07505;
%   [~,~,cp,tauq] = dwtleader(mrw07505);
%   cp([1 2])
%   plot(-5:5,tauq,'bo--'); title('Estimated Scaling Exponents');
%   grid on;
%   xlabel('Q-th Moments'); ylabel('\tau(q)');
%
%   %Example 3:
%   %   Compare multifractal spectrum of heart-rate variability data
%   %   before and after application of a drug that reduces heart dynamics.
%   load hrvDrug;
%   predrug = hrvDrug(1:4642);
%   postdrug = hrvDrug(4643:end);
%   [dhpre,hpre] = dwtleader(predrug);
%   [dhpost,hpost] = dwtleader(postdrug);
%   plot(hpre,dhpre,hpost,dhpost);
%   xlabel('h'); ylabel('D(h)');
%   grid on;
%   legend('Predrug','Postdrug');
%
% See also WTMM, WFBM



narginchk(1,6)
nargoutchk(0,6)
validateattributes(x,{'numeric'},{'finite','vector','nonempty'});
params = parseinputs(varargin{:});
wname = params.wname;
% Mininum regression level
j1 = params.j1;

wtype = wavemngr('type',wname);
if wtype == 1 || wtype == 2
    [Lo,Hi] = wfilters(wname);
else
    error(message('Wavelet:mfa:OrthorBiorthWavelet'));
end


x = x(:)';
[leaders,scales,ncount] = wavelet.internal.dwtleaders(x,Lo,Hi);
% Determine default maximum level
J = find(ncount>=6,1,'last');
% Check if user has specified a maximum regression level
if isfield(params,'J') 
    if params.J <= J
        J = params.J;
    else 
        error(message('Wavelet:mfa:MaxRegressionLevel',J));
    end
end

if J < j1+1
    error(message('Wavelet:mfa:RegressionLevels'));
end

Nq = numel(params.q);
Nest = length(ncount);
Dq = zeros(Nq,Nest);
Hq = zeros(Nq,Nest);
Cp = zeros(params.cumulant,Nest);
zetaq = zeros(Nq,Nest);
% Compute wavelet multiresolution structure functions
% Structure function computation includes finest-scale "leaders"
% Regression does not
for jj = 1:Nest
    [zetaq(:,jj),Dq(:,jj), Hq(:,jj), Cp(:,jj)] = ...
        wavelet.internal.mfstructfunctions(abs(leaders{jj}),params);
end


Cp = Cp*log2(exp(1));
Y = [zetaq; Dq; Hq; Cp];
xj = log2(scales);

Y = Y(:,j1:J);
xj = xj(j1:J);
ncount = ncount(j1:J);

%Leaders are not defined at level 1
leaders = leaders(2:end);
% Forming multiresolution structure functions
structfunc.Tq = Y';
structfunc.weights = ncount;
structfunc.logscales = xj;
% Create design matrix
X = ones(length(structfunc.logscales),2);
X(:,2) = structfunc.logscales;
% Least-squares regression with weigths
betahat = lscov(X,structfunc.Tq,structfunc.weights);
% Ignore intercept terms -- use only slopes
betahat = betahat(2,:);
zq = betahat(1:Nq);
dh = betahat(Nq+1:2*Nq)+1;
h = betahat(2*Nq+1:3*Nq);
cp = betahat(3*Nq+1:end);




function params = parseinputs(varargin)

params.wname = 'bior1.5';
params.j1 = 3;
params.q = -5:5;
params.cumulant = 3;
if isempty(varargin)
    return;
end

tf = find(strncmpi(varargin,'MinRegressionLevel',3));
if nnz(tf) == 1
    params.j1 = varargin{tf+1};
    validateattributes(params.j1,{'numeric'},{'scalar','integer', '>=',2},...
        'dwtleader','MinRegressionLevel');
    varargin(tf:tf+1) = [];
    if isempty(varargin)
        return;
    end
end

tf = find(strncmpi(varargin,'MaxRegressionLevel',3));
if nnz(tf) == 1
    params.J = varargin{tf+1};
    validateattributes(params.J,{'numeric'},{'scalar','integer', ...
        '>=',params.j1+1},'dwtleader','MaxRegressionLevel');
    varargin(tf:tf+1) = [];
    if isempty(varargin)
        return;
    end
end

% Only remaining varargin argument must be wavelet
tf = cellfun(@(x)ischar(x),varargin);
if nnz(tf) == 1
    params.wname = varargin{tf>0};
    varargin(tf>0) = [];
else
    error('Invalid');
end

if ~isempty(varargin)
    error(message('Wavelet:mfa:UnsupportedSyntax'));
end