function [ACF,NN] = acovf(Z,KMAX,Mode,Mode2);
% ACOVF estimates autocovariance function (not normalized)
% NaN's are interpreted as missing values. 
%
% [ACF,NN] = acovf(Z,MAXLAG,Mode);
%
% Input:
%  Z    Signal (one channel per row);
%  MAXLAG  maximum lag
%  Mode	'biased'  : normalizes with N
%	'unbiased': normalizes with N-lag
%	'coeff'	  : normalizes such that lag 0 is 1	
%        others	  : no normalization
%
% Output:
%  ACF autocovariance function
%  NN  number of valid elements 
%
%
% REFERENCES:
%  A.V. Oppenheim and R.W. Schafer, Digital Signal Processing, Prentice-Hall, 1975.
%  S. Haykin "Adaptive Filter Theory" 3ed. Prentice Hall, 1996.
%  M.B. Priestley "Spectral Analysis and Time Series" Academic Press, 1981. 
%  W.S. Wei "Time Series Analysis" Addison Wesley, 1990.
%  J.S. Bendat and A.G.Persol "Random Data: Analysis and Measurement procedures", Wiley, 1986.

%	$Revision: 1.10 $
%	$Id: acovf.m,v 1.10 2005/05/31 14:30:57 qspencer Exp $
%	Copyright (C) 1998-2003 by Alois Schloegl <a.schloegl@ieee.org>		

% This library is free software; you can redistribute it and/or
% modify it under the terms of the GNU Library General Public
% License as published by the Free Software Foundation; either
% version 2 of the License, or (at your option) any later version.
% 
% This library is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% Library General Public License for more details.
%
% You should have received a copy of the GNU Library General Public
% License along with this library; if not, write to the
% Free Software Foundation, Inc., 59 Temple Place - Suite 330,
% Boston, MA  02111-1307, USA.


if nargin<3, Mode='biased'; end;

[lr,lc] = size(Z);

MISSES = sum(isnan(Z)')';
if any(MISSES); % missing values
	M = real(~isnan(Z));
	Z(isnan(Z))=0;
end;

if (nargin == 1) 
	KMAX = lc-1; 
elseif (KMAX >= lc-1) 
	KMAX = lc-1;
end;

ACF = zeros(lr,KMAX+1);

if nargin>3,		% for testing, use arg4 for comparing the methods,

elseif 	(KMAX*KMAX > lc*log2(lc)) % & isempty(MISSES);	
	Mode2 = 1;
elseif 	(10*KMAX > lc);
	Mode2 = 3;
else
	Mode2 = 4;
end;


%%%%% ESTIMATION of non-normalized ACF %%%%%

% the following algorithms gve equivalent results, however, the computational effort is different,
% depending on lr,lc and KMAX, a different algorithm is most efficient.
if Mode2==1; % KMAX*KMAX > lc*log(lc);        % O(n.logn)+O(KČ)
        tmp = fft(Z',2^nextpow2(size(Z,2))*2);
        tmp = ifft(tmp.*conj(tmp));
        ACF = tmp(1:KMAX+1,:)'; 
        if ~any(any(imag(Z))), ACF=real(ACF); end; % should not be neccessary, unfortunately it is.
elseif Mode2==3; % (10*KMAX > lc)   % O(n*K)     % use fast Built-in filter function
        for L = 1:lr,
                acf = filter(Z(L,lc:-1:1),1,Z(L,:));
                ACF(L,:)= acf(lc:-1:lc-KMAX);
        end;    
else Mode2==4; % O(n*K)
        for L = 1:lr,
                for K = 0:KMAX, 
                        ACF(L,K+1) = Z(L,1:lc-K) * Z(L,1+K:lc)';
                end;
        end;    
end;


%%%%% GET number of elements used for estimating ACF - is needed for normalizing ACF %%%%%

if any(MISSES),
    % the following algorithms gve equivalent results, however, the computational effort is different,
    % depending on lr,lc and KMAX, a different algorithm is most efficient.
    if Mode2==1; % KMAX*KMAX > lc*log(lc);        % O(n.logn)+O(KČ)
        tmp = fft(M',2^nextpow2(size(M,2))*2);
        tmp = ifft(tmp.*conj(tmp));
        NN = tmp(1:KMAX+1,:)'; 
        if ~any(any(imag(M))), NN=real(NN); end; % should not be neccessary, unfortunately it is.
    elseif Mode2==3; % (10*KMAX > lc)   % O(n*K)     % use fast Built-in filter function
        for L = 1:lr,
                acf = filter(M(L,lc:-1:1),1,M(L,:));
                NN(L,:)= acf(lc:-1:lc-KMAX);
        end;    
    else Mode2==4; % O(n*K)
        for L = 1:lr,
                for K = 0:KMAX, 
                        NN(L,K+1) = M(L,1:lc-K) * M(L,1+K:lc)';
                end;
        end;    
    end;
else
    NN = (ones(lr,1)*(lc:-1:lc-KMAX));
end;


if strcmp(Mode,'biased')
	if ~any(MISSES),
	        ACF=ACF/lc;
	else
	        %ACF=ACF./((lc-MISSES)*ones(1,KMAX+1));
	        ACF=ACF./max(NN + ones(lr,1)*(0:KMAX),0);
	end;

elseif strcmp(Mode,'unbiased')
        ACF=ACF./NN; 
	%if ~any(MISSES),
	%       ACF=ACF./(ones(lr,1)*(lc:-1:lc-KMAX));
	%else
	%	ACF=ACF./((lc-MISSES)*ones(1,KMAX+1) - ones(lr,1)*(0:KMAX));
	%end;

elseif strcmp(Mode,'coeff')
        %ACF = ACF ./ ACF(:,ones(1,KMAX+1)) .* ((lc-MISSES)*ones(1,KMAX+1));
        ACF = ACF./NN; 
	ACF = ACF./(ACF(:,1)*ones(1,size(ACF,2)));
else 

end;
