wfbmesti.m
上传用户:haiyisale
上传日期:2013-01-09
资源大小:3246k
文件大小:2k
源码类别:

波变换

开发平台:

Matlab

  1. function Hest = wfbmesti(x)
  2. %WFBMESTI Estimate fractal index.
  3. %   HEST = WFBMESTI(X) returns a row vector HEST which contains
  4. %   three estimates of the fractal index H of the signal X supposed
  5. %   to come from a fractional Brownian motion of parameter H.
  6. %
  7. %   The two first estimates are based on second order discrete 
  8. %   derivative, the second one is wavelet based.
  9. %   The third estimate is based on the linear regression in 
  10. %   loglog plot, of the variance of detail versus level.
  11. %   M. Misiti, Y. Misiti, G. Oppenheim, J.M. Poggi 22-May-2003.
  12. %   Last Revision: 11-Jul-2003.
  13. %   Copyright 1995-2004 The MathWorks, Inc.
  14. %   $Revision: 1.1.6.2 $  $Date: 2004/03/15 22:42:44 $ 
  15. % First method : second order discrete derivative.
  16. %-------------------------------------------------
  17. y = diff(x); % x is a Fractional Gaussian Noise
  18. y = cumsum(y(:)');  
  19. n = length(y);
  20. b1 = [1 -2 1];
  21. b2 = [1  0 -2 0 1];
  22. y1 = filter(b1,1,y);
  23. y1 = y1(length(b1):n);
  24. y2 = filter(b2,1,y);
  25. y2 = y2(length(b2):n);
  26. s1 = mean(y1.^2);
  27. s2 = mean(y2.^2);
  28. Hest(1) = 0.5*log2(s2/s1);
  29.         
  30. % Second method : second order discrete derivative (using wavelets).
  31. %-------------------------------------------------------------------
  32. [LO_D,c1,LO_R,HI_R] = wfilters('sym5'); 
  33. c2 = [c1;zeros(1,length(c1))];
  34. c2 = c2(:)';
  35. cy1 = filter(c1,1,y);
  36. cy1 = cy1(length(c1):n);
  37. cy2 = filter(c2,1,y);
  38. cy2 = cy2(length(c2):n);
  39. cs1 = mean(cy1.^2);
  40. cs2 = mean(cy2.^2);
  41. Hest(2) = 0.5*log2(cs2/cs1);
  42.         
  43. % Third method : variance versus level.
  44. %-------------------------------------------------
  45. % Wavelet decomposition of the fractal signal.
  46. levdec = min(wmaxlev(size(x),'haar'),6);
  47. [c,l]  = wavedec(x,levdec,'haar');
  48. % Robust estimates of the standard deviation of detail coeff.
  49. % Recall that x is supposed to be Gaussian.
  50. lvls  = [1:levdec];
  51. stdc  = wnoisest(c,l,lvls);
  52. % Perform regression and compute estimate of h.
  53. p = polyfit(lvls,log2(stdc.^2),1);
  54. Hest(3) = (p(1)-1)/2;