%======================================================================

%                    B E S T T L A G . M 

%                    doc: Wed Jul 19 15:34:53 2006

%                    dlm: Fri Mar  5 15:50:57 2010

%                    (c) 2006 M. Visbeck

%                    uE-Info: 51 16 NIL 0 0 72 0 2 4 NIL ofnI

%======================================================================

function [lag,co]=besttlag(a1,a2,nlag,npoint,nsect)

% function [lag,co]=besttlag(a1,a2,nlag,npoint,nsect)

%

% function to find the best shift for two vectors

% uses median difference

% 

%  input:  a1: first vector (time,variable) should be the high resolution one

%          a2: second vector(time,variable)

%          nlag = number of point to try shifiting 

%          npoint = use how many points

%          nsect = number of chunks of data to use

%

% M. Visbeck LDEO August-2002

% MODIFICATIONS BY ANT:

%	Jul 19, 2006: - removed NaNs before interp1() to avoid Matlab 7.2 warning

l1=size(a1,1);

l2=size(a2,1);

if nargin<4, npoint=min(length(a1),3*nlag); end 

if nargin<5, nsect=ceil(l1/(npoint)); end 

if nargin<3, nlag=fix(npoint/8); end

% loop over chunks of data to be used avoid begining and end

ismed=round(linspace(1,l1,nsect+4));

ismed([1,2, end-1, end])=[];

% add largest gradient to list

iok=fix(l1*0.25):(l1*.75);

[dum,ii]=find(maxnan(abs(diff(a1(iok,2)))));

ismed(end+1)=iok(ii);

for i=1:length(ismed)

% select data for section

 isect=[-fix(npoint/2):(npoint/2)]+ismed(i);

 iok=find(isect>0 & isect<l1);

 isect=isect(iok);

% interpolate a2 on a1 with regards to time

 igood = find(isfinite(a2(:,2)));

 a22=interp1(a2(igood,1),a2(igood,2),a1(isect,1),'nearest');

 a12=a1(isect,2);

 [lagv(i),i1,i2,cov(i)]=bestlag(a12,a22,nlag);

 if lagv(i)==0 & cov(i)==1, cov(i)=nan; end

 disp([' lag: ',int2str(lagv(i)),'  correlation: ',num2str(cov(i))])

end

if maxnan(cov)<0.97

 % try acceleration near large w-difference

 i=length(lagv)+1;

 [lagv(i),i1,i2,cov(i)]=bestlag(diff(a12),diff(a22),nlag);

 disp([' acceleration lag: ',int2str(lagv(i)),'  correlation: ',num2str(cov(i))])

end

if maxnan(cov)<0.97

 % one scan over whole length of max correlation was not great

 i=length(lagv)+1;

 isect=nlag:(length(a1)-nlag);

 a22=interp1(a2(:,1),a2(:,2),a1(isect,1),'nearest');

 a12=a1(isect,2);

 [lagv(i),i1,i2,cov(i)]=bestlag(a12,a22,nlag);

 disp([' all data lag: ',int2str(lagv(i)),'  correlation: ',num2str(cov(i))])

end

if maxnan(cov)<0.97

 % one scan over whole length of max correlation was not great

 i=length(lagv)+1;

 ii=find(~isfinite(a12)); a12(ii)=0; 

 ii=find(~isfinite(a22)); a22(ii)=0; 

 [lagv(i),i1,i2,cov(i)]=bestlag(cumsum(a12),cumsum(a22),nlag);

 disp([' all data integral lag: ',int2str(lagv(i)),'  correlation: ',num2str(cov(i))])

end

% choose the most likely lag

ii=find(~isfinite(cov));

cov(ii)=[];

lagv(ii)=[];

if length(lagv)>2

 lag0=median(lagv);

else

 lag0=nan;

end

disp([' median lag ',int2str(lag0)])

nnlag=-nlag:nlag;

lagh=hist(lagv,nnlag);

% prefer 0 lag slightly

[nlag1,i1]=max(lagh-abs(nnlag)/(2*max(nnlag)));

lag1=nnlag(i1);

iok=find(lag1==lagv);

co=meannan(cov(iok));

disp([' most popular lag ',int2str(lag1)])

[cos,is]=sort(cov);

lags=lagv(is);

lag2=lags(end);

disp([' best correlated lag ',int2str(lag2)])

% decide which one to use

if co*sqrt(length(iok)) > cos(end)

 lag=lag1;

else

 lag=lag2;

end

iok=find(lag==lagv);

co=maxnan(cov(iok));

disp([' BESTTLAG:  lag is: ',int2str(lag),'  for which ',...

      int2str(length(iok)/nsect*100),'% of ',int2str(nsect),' lags agree'])

return

%===============================

function y = median(x,dim)

%MEDIAN Median value.

%   For vectors, MEDIAN(X) is the median value of the elements in X.

%   For matrices, MEDIAN(X) is a row vector containing the median

%   value of each column.  For N-D arrays, MEDIAN(X) is the median

%   value of the elements along the first non-singleton dimension

%   of X.

%

%   MEDIAN(X,DIM) takes the median along the dimension DIM of X.

%

%   Example: If X = [0 1 2

%                    3 4 5]

%

%   then median(X,1) is [1.5 2.5 3.5] and median(X,2) is [1

%                                                         4]

%

%   See also MEAN, STD, MIN, MAX, COV.

%   Copyright 1984-2002 The MathWorks, Inc. 

%   $Revision: 5.15 $  $Date: 2002/06/05 17:06:39 $

if nargin==1, 

  dim = min(find(size(x)~=1)); 

  if isempty(dim), dim = 1; end

end

if isempty(x), y = []; return, end

siz = [size(x) ones(1,dim-ndims(x))];

n = size(x,dim);

% Permute and reshape so that DIM becomes the row dimension of a 2-D array

perm = [dim:max(length(size(x)),dim) 1:dim-1];

x = reshape(permute(x,perm),n,prod(siz)/n);

% Sort along first dimension

x = sort(x,1);

if rem(n,2) % Odd number of elements along DIM

  y = x((n+1)/2,:);

else % Even number of elements along DIM

  % y = (x(n/2,:) + x((n/2)+1,:))/2;

  y =  x(fix((n/2)+1),:);

end

% Check for NaNs

y(isnan(x(1,:)) | isnan(x(n,:))) = NaN;

% Permute and reshape back

siz(dim) = 1;

y = ipermute(reshape(y,siz(perm)),perm);