function [lambdaE,y] = powtr(x,seglen,fig)
%       COMPUTING ESTIMATE OF PARAMETER lambda FOR THE POWER TRANSFORM:
%                        y = x.^lambda   if lambda~=0
%                          = log(x)      if lambda =0  (natural logarithm)
% =============================================================================
% function [lambdaE,y] = powtr(x,seglen,fig)
% INPUT:
% x      ... is an n x 1 or 1 x n data vector of POSITIVE observed data
%            or vector of smoothed values (determnistic component) provided
%            that the residual error component r describing the decomposition
%            x = y + r is supplied instead of seglen
% seglen ... segment length used for computation of "std. deviation" growth
%            (recommended: 4<=seglen<=12)
%            or seglen=r is an n x 1 or 1 x n vector of the residual errors
%            provided that the deterministic component y was supplied instead
%            of x. In such a case the "std. deviation" growth will be evaluated
%            more precisely from the dependence of r on y.
%            optional: []=default segment length=8
% fig ... figure handle for the plot of estimated "std. deviation" versus mean
%         optional: []=default=no plot
% OUTPUT:
% lambdaE ... =[lambda,E] is and 1 x 2 vector, where
%              lambda is estimate of parameter lambda
%              [lambda-E,lambda+E] is its confidence interval of 95%
% y  ... =data transformed according to parameter lambda
%         vector of the same size as x

% (C) 1997, V. Vesely, Masaryk University of Brno (Czech Republic)
%
error(nargchk(1,3,nargin));
if nargin<3,        fig    = []; end
if nargin<2,        seglen = []; end
if isempty(seglen), seglen =  8; end

sz = size(x);
if min(sz)~=1
    error('x is not a vector');
end

n = length(x);

if length(seglen)==1,    % segment length supplied
   if seglen < 4
      disp('WARNING: segment length too short: 4 will be used');
      seglen = 4;
   end
   if seglen > 12
      disp('WARNING: segment length too long: 12 will be used');
      seglen = 12;
   end
   nseg = floor(n./seglen);
   if nseg < 5
      disp('WARNING: less than 5 segments: minimal segment length 4 will be used');
      seglen = 4;
      nseg = floor(n./seglen);
      if nseg < 5,
         error('Still less than 5 segments - estimate impossible');
      end
   end
   ns = nseg*seglen;
   n0 = n - ns;
   X  = reshape(x(1:ns),seglen,nseg); % Segmenting into matrix seglen x nseg
   r  = (max(X) - min(X)).';    % estimate of standard deviation in each segment
   m  = mean(X).';              % mean estimate in each segment
   if n0 > 4
      x0 = x(ns+1:n);           % append estimate for residual segment
      r  = [r;max(x0)-min(x0)];
      m  = [m;mean(x0)];
   end
else                            % residuals supplied
   r=seglen;
   if min(size(r))~=1 | length(r)~=n
      error('seglen should be a vector (of residuals) of the same length as x');
   end
   m=reshape(x,n,1);
   r=reshape(r,n,1);
   x=m+r;
   if sz(1)==1, x=x.'; end
   r=abs(r);
end

% Estimating theta for r(m)=sigma*m^theta using linear regression in
%                      log(r)=log(sigma)+theta*log(m)
D  = [log(m),ones(size(m))];
mm = linspace(min(m),max(m)).';
C  = [log(mm),ones(size(mm))];
[bE,rE,sE,PbE,CbE] = linregr(log(r),D,C);
                                % theta  = bE(1,1)   +- bE(1,2)
lambdaE = [1-bE(1,1),bE(1,2)];  % lambda = (1-theta) +- bE(1,2)

if ~isempty(fig)  % plot of r versus m + regression line + confidence interval
   rr  = exp(CbE(:,1));
   rr1 = exp(CbE(:,1)-CbE(:,2));
   rr2 = exp(CbE(:,1)+CbE(:,2));
   figure(fig); plot(m,r,'ow',mm,rr,'y-',mm,[rr1,rr2],'y--');
   title('Estimated standard deviation versus mean');
end

if nargout>1      % transform data according to lambda
   if abs(lambdaE(1)-1)<lambdaE(2)   % lambda=1 is in the confidence interval
       y = x;                        % then no transform
       disp('>>> CONSTANT VARIANCE (data not transformed): y=x <<<');
   elseif abs(lambdaE(1))<lambdaE(2) % lambda=0 in the confidence interval
       y=log(x);                     % then use logarithmic transform
       disp(['>>> LINEAR GROWTH OF VARIANCE (logarithmic transform): ',...
             'y=log(x) <<<']);
   else
       y=x.^lambdaE(1);              % else the power transform
       disp(['>>> POWER GROWTH OF VARIANCE (power transform): y=x.^',...
            num2str(lambdaE(1)),' <<<']);
   end
end
