Logo Search packages:      
Sourcecode: octave-tsa version File versions  Download package

acovf.m

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: 4585 $
%     $Id: acovf.m 4585 2008-02-04 13:47:45Z adb014 $
%     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, see <http://www.gnu.org/licenses/>.


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;

Generated by  Doxygen 1.6.0   Back to index