function [ccoef,sig,acr] = lagcorr(data,refer,lags,level,degf);
% function [ccoef,sig] = lagcorr(data,refer,[lags],[level],[degf]);
%
% calculate lagged correlations and their significance
%
% input  :	data		- data array
%		refer		- reference series
%		[lags]	[0]	- vector of lags (reference leads for positive)
%		[level]	[95]	- significance level in % to be checked
%		[degf]		- degrees of freedom of reference series
%
% output :	ccoef		- correlation coefficient
%		sig		- significant (1,-1) or not (0)
%
% Correlation is calculated between refer and the dim of data which agrees
% in length. If there is more than one, the first is taken.
% Default is simple correlation without lag. If lags are given, the results
% for the different lags are listed in the lagged dimension.
% The degrees of freedom of the series can be prescribed as an integer>0
% or as a fraction of the length, otherwise the length of the lagged
% series is taken.
%
% uses :	anorm.m
%
% version 0.1.2		last change 17.03.2004

% G.Krahmann, LODYC Paris, May 1998


% references :	Bendat, Piersol:   Random Data: Analysis and measurement
%                                  procedures
%		Wiley-Interscience, 1971, p.125-128

% history:
%
% removed bug 		G.Krahmann, LODYC Paris 1998	0.1.0 --> 0.1.1
% changed shiftdim in 5.2 ?
% catch correl of 1	G.K. LDEO 2004			0.1.1 --> 0.1.2


%
% check arguments
%
if nargin<2
  error('Not enough input arguments.')
end
refer = squeeze(refer);
lr = length(refer);
if ~(size(refer,1)==1 | size(refer,2)==1)
  error('Reference series must be a vector.')
end
if sum( lr==size(data) )==0
  error('Data must contain a dimension with the length of the reference.')
end
if nargin<5
  degf = [];
end
if nargin<4
  level = 95;
end
if nargin<3
  lags = 0;
end
if sum(isnan(refer))>1
  error('Reference series must not contain NaN''s')
end
if degf<0 | ( degf>1 & degf<4 )
  warning('Skipping degrees of freedom.')
  degf = [];
end


%
% check for simple case
%
simple = 0;
if length(lags)==1
  if lags==0
    simple = 1;
  end
end


%
% prepare output arrays
%
ds = size(data);
di = find(ds==lr);
di = di(1);
ns = ds;
ns(di) = length(lags);
ccoef = repmat(0,ns);
sig = repmat(0,ns);
acr = repmat(0,ns);


% 
% subtract mean and enlarge reference
%
dr = ones(1,length(ns));
dr(di) = size(data,di);
nr = ds;
nr(di) = 1;
data = data - repmat(nmean(data,di),dr);
refer = refer - nmean(refer);
ref = repmat(0,size(data));
ref = shiftdim(ref,di-1);		% I was not able to get repamt
for n=1:length(refer)			% to work properly here
  ref(n,:) = ref(n,:) + refer(n);
end
ref = shiftdim(ref,length(ds)-di+1);


%
% correlation level
%
al = anorm(100-level);


%
% do correlations
%
if simple==1

  ref2 = refer.^2;
  mref2 = mean(ref2);
  ccoef = nmean(data.*ref,di)./sqrt( mref2 * nmean(data.^2,di) );
  if nargout>1
    if isempty(degf)
      ng = sum( isfinite(data), di );
    elseif degf>3
      ng = floor(degf(1));
    else
      ng = sum( isfinite(data), di )*degf;
    end
    [ccoef,ind1] = nans(ccoef,nan,1,'==');
    [ccoef,ind2] = nans(ccoef,nan,-1,'==');
    acr = 0.5*sqrt(ng-3).*log( (1+ccoef)./(1-ccoef) );
    sig = (abs(acr)>al).*sign(real(acr));
    if ~isempty(ind1)
      sig(ind1) = 1;
      ccoef(ind1) = 1;
    end
    if ~isempty(ind2)
      sig(ind2) = -1;
      ccoef(ind2) = -1;
    end
  end 

else

  % shift dimensions to get the series into the first
  ref = shiftdim(ref,di-1);
  data = shiftdim(data,di-1);
  ccoef = shiftdim(ccoef,di-1);
  sig = shiftdim(sig,di-1);
  acr = shiftdim(acr,di-1);

  % reshape to get only two dimensions
  sr = size(ref);
  ref = reshape(ref,[sr(1),prod(sr(2:length(sr)))]);
  data = reshape(data,[sr(1),prod(sr(2:length(sr)))]);
  sc = size(ccoef);
  ccoef = reshape(ccoef,[sc(1),prod(sc(2:length(sc)))]);
  sig = reshape(sig,[sc(1),prod(sc(2:length(sc)))]);
  acr = reshape(acr,[sc(1),prod(sc(2:length(sc)))]);

  data2 = data.^2;
  for n=1:length(lags)
    if lags(n)<0
      ref2 = refer(1-lags(n):length(refer)).^2;
      mref2 = mean(ref2);
      d = data(1:size(data,1)+lags(n),:);
      r = ref(1-lags(n):size(ref,1),:);
      ccoef(n,:) = nmean(d.*r,1)./sqrt( mref2 * nmean(d.^2,1) );
      if isempty(degf)
        ng = sum( isfinite(d), 1 );
      elseif degf>3
        ng = floor(degf(1));
      else
        ng = sum( isfinite(d), 1 )*degf;
      end
      [ccoef,ind1] = nans(ccoef,nan,1,'==');
      [ccoef,ind2] = nans(ccoef,nan,-1,'==');
      acr(n,:) = 0.5*sqrt(ng-3).*log( (1+ccoef(n,:))./(1-ccoef(n,:)) );
      sig(n,:) = (abs(acr(n,:))>al).*sign(real(acr(n,:)));
      if ~isempty(ind1)
        sig(ind1) = 1;
        ccoef(ind1) = 1;
      end
      if ~isempty(ind2)
        sig(ind2) = -1;
        ccoef(ind2) = -1;
      end
    elseif lags(n)>0
      ref2 = refer(1:length(refer)-lags(n)).^2;
      mref2 = mean(ref2);
      d = data(1+lags(n):size(data,1),:);
      r = ref(1:size(ref,1)-lags(n),:);
      ccoef(n,:) = nmean(d.*r,1)./sqrt( mref2 * nmean(d.^2,1) );
      if isempty(degf)
        ng = sum( isfinite(d), 1 );
      elseif degf>3
        ng = floor(degf(1));
      else
        ng = sum( isfinite(d), 1 )*degf;
      end
      [ccoef,ind1] = nans(ccoef,nan,1,'==');
      [ccoef,ind2] = nans(ccoef,nan,-1,'==');
      acr(n,:) = 0.5*sqrt(ng-3).*log( (1+ccoef(n,:))./(1-ccoef(n,:)) );
      sig(n,:) = (abs(acr(n,:))>al).*sign(real(acr(n,:)));
      if ~isempty(ind1)
        sig(ind1) = 1;
        ccoef(ind1) = 1;
      end
      if ~isempty(ind2)
        sig(ind2) = -1;
        ccoef(ind2) = -1;
      end
    elseif lags(n)==0
      ref2 = refer.^2;
      mref2 = mean(ref2);
      d = data;
      r = ref;
      ccoef(n,:) = nmean(d.*r,1)./sqrt( mref2 * nmean(d.^2,1) );
      if isempty(degf)
        ng = sum( isfinite(d), 1 );
      elseif degf>3
        ng = floor(degf(1));
      else
        ng = sum( isfinite(d), 1 )*degf;
      end
      [ccoef,ind1] = nans(ccoef,nan,1,'==');
      [ccoef,ind2] = nans(ccoef,nan,-1,'==');
      acr(n,:) = 0.5*sqrt(ng-3).*log( (1+ccoef(n,:))./(1-ccoef(n,:)) );
      sig(n,:) = (abs(acr(n,:))>al).*sign(real(acr(n,:)));
      if ~isempty(ind1)
        sig(ind1) = 1;
        ccoef(ind1) = 1;
      end
      if ~isempty(ind2)
        sig(ind2) = -1;
        ccoef(ind2) = -1;
      end
    end
  end 

  % reshape back
  ref = reshape(ref,sr);
  data = reshape(data,sr);
  ccoef = reshape(ccoef,sc);
  sig = reshape(sig,sc);
  acr = reshape(acr,sc);

  % shift dimensions back
  ref = shiftdim(ref,length(ds)-di+1);
  data = shiftdim(data,length(ds)-di+1);
  ccoef = shiftdim(ccoef,length(ds)-di+1);
  sig = shiftdim(sig,length(ds)-di+1);
  acr = shiftdim(acr,length(ds)-di+1);

end
