1 function [ACF,NN] = acovf(Z,KMAX,Mode,Mode2);
2 % ACOVF estimates autocovariance function (not normalized)
3 % NaN's are interpreted as missing values.
5 % [ACF,NN] = acovf(Z,MAXLAG,Mode);
8 % Z Signal (one channel per row);
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
16 % ACF autocovariance function
17 % NN number of valid elements
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.
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/
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.
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.
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/>.
46 if nargin<3, Mode='biased'; end;
50 MISSES = sum(isnan(Z)')';
51 if any(MISSES); % missing values
62 ACF = zeros(lr,KMAX+1);
64 if nargin>3, % for testing, use arg4 for comparing the methods,
66 elseif (KMAX*KMAX > lc*log2(lc)) % & isempty(MISSES);
68 elseif (10*KMAX > lc);
75 %%%%% ESTIMATION of non-normalized ACF %%%%%
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
86 acf = filter(Z(L,lc:-1:1),1,Z(L,:));
87 ACF(L,:)= acf(lc:-1:lc-KMAX);
89 else Mode2==4; % O(n*K)
92 ACF(L,K+1) = Z(L,1:lc-K) * Z(L,1+K:lc)';
98 %%%%% GET number of elements used for estimating ACF - is needed for normalizing ACF %%%%%
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
110 acf = filter(M(L,lc:-1:1),1,M(L,:));
111 NN(L,:)= acf(lc:-1:lc-KMAX);
113 else Mode2==4; % O(n*K)
116 NN(L,K+1) = M(L,1:lc-K) * M(L,1+K:lc)';
121 NN = (ones(lr,1)*(lc:-1:lc-KMAX));
125 if strcmp(Mode,'biased')
129 %ACF=ACF./((lc-MISSES)*ones(1,KMAX+1));
130 ACF=ACF./max(NN + ones(lr,1)*(0:KMAX),0);
133 elseif strcmp(Mode,'unbiased')
136 % ACF=ACF./(ones(lr,1)*(lc:-1:lc-KMAX));
138 % ACF=ACF./((lc-MISSES)*ones(1,KMAX+1) - ones(lr,1)*(0:KMAX));
141 elseif strcmp(Mode,'coeff')
142 %ACF = ACF ./ ACF(:,ones(1,KMAX+1)) .* ((lc-MISSES)*ones(1,KMAX+1));
144 ACF = ACF./(ACF(:,1)*ones(1,size(ACF,2)));