function iscale=integralscale(xblocked,deltat,errorcheck)
% function iscale=integralscale(xblocked,deltat,errorcheck)
%
% This is a really quick implmentation of an integral scale calculation
%
% inputs:
%
% xblocked  array (n x m) with n=signal length and m=# of signals
%           Each xblocked(:,k) is a separate signal to calculate the inegral scale of.
%
% deltat    sampling rate (secends)
%
% errorcheck errorcheck==1 will cause some error checking to occur.
%
% outputs:
%
% iscale   iscale(k) is the estimated integral scale for signal xblocked(:,k)
%
xsize=size(xblocked);

iscale=zeros([xsize(2) 1]);
for m=1:xsize(2)  % calculate the integral scale at each block
   gamma=zeros([xsize(1) 1]);  % gamma(i)=cov(X(m),X(m+i-1))
   x=xblocked(:,m);
   xmean=mean(x);
   x=x-xmean;
   gamma(1)=mean( x(1:(xsize(1)-1+1)) .* x(1:(xsize(1)))  );
   gammalength=1;
   %   Use gamma(k)<.3*gamma(1) as the cut off.
   while (gamma(gammalength)>.3*gamma(1))&(gammalength<xsize(1)-1)
      gammalength=gammalength+1;
      gamma(gammalength)=mean( x(1:(xsize(1)-gammalength+1)) .* x(gammalength:(xsize(1)))  );
   end % end of k loop to calculate gamma
   gamma=gamma/gamma(1);  % Take care of sigma^2
   iscaleold=xmean*sum((gamma(1:100)))*deltat;
   tau=deltat*(0:(gammalength-1));
   lambda=(tau*log(gamma(1:gammalength)))/(tau*tau');      
   iscale(m)=-xmean/lambda;   % <-- This is the correct normalization
      if gammalength>xsize(1)/10
         disp('Warning in integralscale:  covariance estimate may be inaccurate.')
         iscale(m)=0;
   end
end % end of m loop to calculate the integral scale for block m
