octave_packages/tsa-4.2.4/acovf.m

   1 function [ACF,NN] = acovf(Z,KMAX,Mode,Mode2);
   2 % ACOVF estimates autocovariance function (not normalized)
   3 % NaN's are interpreted as missing values.
   4 %
   5 % [ACF,NN] = acovf(Z,MAXLAG,Mode);
   6 %
   7 % Input:
   8 %  Z    Signal (one channel per row);
   9 %  MAXLAG  maximum lag
  10 %  Mode 'biased'  : normalizes with N [default]
  11 %       'unbiased': normalizes with N-lag
  12 %       'coeff'   : normalizes such that lag 0 is 1
  13 %        others   : no normalization
  14 %
  15 % Output:
  16 %  ACF autocovariance function
  17 %  NN  number of valid elements
  18 %
  19 % REFERENCES:
  20 %  A.V. Oppenheim and R.W. Schafer, Digital Signal Processing, Prentice-Hall, 1975.
  21 %  S. Haykin "Adaptive Filter Theory" 3ed. Prentice Hall, 1996.
  22 %  M.B. Priestley "Spectral Analysis and Time Series" Academic Press, 1981.
  23 %  W.S. Wei "Time Series Analysis" Addison Wesley, 1990.
  24 %  J.S. Bendat and A.G.Persol "Random Data: Analysis and Measurement procedures", Wiley, 1986.
  25
  26 %       $Id: acovf.m 7687 2010-09-08 18:39:23Z schloegl $
  27 %       Copyright (C) 1998-2003,2008,2010 by Alois Schloegl <a.schloegl@ieee.org>
  28 %       This is part of the TSA-toolbox. See also
  29 %       http://biosig-consulting.com/matlab/tsa/
  30 %
  31 %    This program is free software: you can redistribute it and/or modify
  32 %    it under the terms of the GNU General Public License as published by
  33 %    the Free Software Foundation, either version 3 of the License, or
  34 %    (at your option) any later version.
  35 %
  36 %    This program is distributed in the hope that it will be useful,
  37 %    but WITHOUT ANY WARRANTY; without even the implied warranty of
  38 %    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  39 %    GNU General Public License for more details.
  40 %
  41 %    You should have received a copy of the GNU General Public License
  42 %    along with this program.  If not, see <http://www.gnu.org/licenses/>.
  43
  44
  45
  46 if nargin<3, Mode='biased'; end;
  47
  48 [lr,lc] = size(Z);
  49
  50 MISSES = sum(isnan(Z)')';
  51 if any(MISSES); % missing values
  52         M = real(~isnan(Z));
  53         Z(isnan(Z))=0;
  54 end;
  55
  56 if (nargin == 1)
  57         KMAX = lc-1;
  58 elseif (KMAX >= lc-1)
  59         KMAX = lc-1;
  60 end;
  61
  62 ACF = zeros(lr,KMAX+1);
  63
  64 if nargin>3,            % for testing, use arg4 for comparing the methods,
  65
  66 elseif  (KMAX*KMAX > lc*log2(lc)) % & isempty(MISSES);
  67         Mode2 = 1;
  68 elseif  (10*KMAX > lc);
  69         Mode2 = 3;
  70 else
  71         Mode2 = 4;
  72 end;
  73
  74
  75 %%%%% ESTIMATION of non-normalized ACF %%%%%
  76
  77 % the following algorithms gve equivalent results, however, the computational effort is different,
  78 % depending on lr,lc and KMAX, a different algorithm is most efficient.
  79 if Mode2==1; % KMAX*KMAX > lc*log(lc);        % O(n.logn)+O(K²)
  80         tmp = fft(Z',2^nextpow2(size(Z,2))*2);
  81         tmp = ifft(tmp.*conj(tmp));
  82         ACF = tmp(1:KMAX+1,:)';
  83         if ~any(any(imag(Z))), ACF=real(ACF); end; % should not be neccessary, unfortunately it is.
  84 elseif Mode2==3; % (10*KMAX > lc)   % O(n*K)     % use fast Built-in filter function
  85         for L = 1:lr,
  86                 acf = filter(Z(L,lc:-1:1),1,Z(L,:));
  87                 ACF(L,:)= acf(lc:-1:lc-KMAX);
  88         end;
  89 else Mode2==4; % O(n*K)
  90         for L = 1:lr,
  91                 for K = 0:KMAX,
  92                         ACF(L,K+1) = Z(L,1:lc-K) * Z(L,1+K:lc)';
  93                 end;
  94         end;
  95 end;
  96
  97
  98 %%%%% GET number of elements used for estimating ACF - is needed for normalizing ACF %%%%%
  99
 100 if any(MISSES),
 101     % the following algorithms gve equivalent results, however, the computational effort is different,
 102     % depending on lr,lc and KMAX, a different algorithm is most efficient.
 103     if Mode2==1; % KMAX*KMAX > lc*log(lc);        % O(n.logn)+O(K²)
 104         tmp = fft(M',2^nextpow2(size(M,2))*2);
 105         tmp = ifft(tmp.*conj(tmp));
 106         NN = tmp(1:KMAX+1,:)';
 107         if ~any(any(imag(M))), NN=real(NN); end; % should not be neccessary, unfortunately it is.
 108     elseif Mode2==3; % (10*KMAX > lc)   % O(n*K)     % use fast Built-in filter function
 109         for L = 1:lr,
 110                 acf = filter(M(L,lc:-1:1),1,M(L,:));
 111                 NN(L,:)= acf(lc:-1:lc-KMAX);
 112         end;
 113     else Mode2==4; % O(n*K)
 114         for L = 1:lr,
 115                 for K = 0:KMAX,
 116                         NN(L,K+1) = M(L,1:lc-K) * M(L,1+K:lc)';
 117                 end;
 118         end;
 119     end;
 120 else
 121     NN = (ones(lr,1)*(lc:-1:lc-KMAX));
 122 end;
 123
 124
 125 if strcmp(Mode,'biased')
 126         if ~any(MISSES),
 127                 ACF=ACF/lc;
 128         else
 129                 %ACF=ACF./((lc-MISSES)*ones(1,KMAX+1));
 130                 ACF=ACF./max(NN + ones(lr,1)*(0:KMAX),0);
 131         end;
 132
 133 elseif strcmp(Mode,'unbiased')
 134         ACF=ACF./NN;
 135         %if ~any(MISSES),
 136         %       ACF=ACF./(ones(lr,1)*(lc:-1:lc-KMAX));
 137         %else
 138         %       ACF=ACF./((lc-MISSES)*ones(1,KMAX+1) - ones(lr,1)*(0:KMAX));
 139         %end;
 140
 141 elseif strcmp(Mode,'coeff')
 142         %ACF = ACF ./ ACF(:,ones(1,KMAX+1)) .* ((lc-MISSES)*ones(1,KMAX+1));
 143         ACF = ACF./NN;
 144         ACF = ACF./(ACF(:,1)*ones(1,size(ACF,2)));
 145 else
 146
 147 end;