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

    function Hest = wfbmesti(x)
%WFBMESTI Estimate fractal index.
%   HEST = WFBMESTI(X) returns a row vector HEST which contains
%   three estimates of the fractal index H of the signal X supposed
%   to come from a fractional Brownian motion of parameter H.
%
%   The two first estimates are based on second order discrete 
%   derivative, the second one is wavelet based.
%   The third estimate is based on the linear regression in 
%   loglog plot, of the variance of detail versus level.

%   M. Misiti, Y. Misiti, G. Oppenheim, J.M. Poggi 22-May-2003.
%   Last Revision: 11-Jul-2003.
%   Copyright 1995-2004 The MathWorks, Inc.

% First method : second order discrete derivative.
%-------------------------------------------------
y = diff(x); % x is a Fractional Gaussian Noise
y = cumsum(y(:)');  
n = length(y);
b1 = [1 -2 1];
b2 = [1  0 -2 0 1];
y1 = filter(b1,1,y);
y1 = y1(length(b1):n);
y2 = filter(b2,1,y);
y2 = y2(length(b2):n);
s1 = mean(y1.^2);
s2 = mean(y2.^2);
Hest(1) = 0.5*log2(s2/s1);
        
% Second method : second order discrete derivative (using wavelets).
%-------------------------------------------------------------------
[LO_D,c1,LO_R,HI_R] = wfilters('sym5'); 
c2 = [c1;zeros(1,length(c1))];
c2 = c2(:)';
cy1 = filter(c1,1,y);
cy1 = cy1(length(c1):n);
cy2 = filter(c2,1,y);
cy2 = cy2(length(c2):n);
cs1 = mean(cy1.^2);
cs2 = mean(cy2.^2);
Hest(2) = 0.5*log2(cs2/cs1);
        
% Third method : variance versus level.
%-------------------------------------------------

% Wavelet decomposition of the fractal signal.
levdec = min(wmaxlev(size(x),'haar'),6);
[c,l]  = wavedec(x,levdec,'haar');

% Robust estimates of the standard deviation of detail coeff.
% Recall that x is supposed to be Gaussian.
lvls  = [1:levdec];
stdc  = wnoisest(c,l,lvls);

% Perform regression and compute estimate of h.
p = polyfit(lvls,log2(stdc.^2),1);
Hest(3) = (p(1)-1)/2;