X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?p=CreaPhase.git;a=blobdiff_plain;f=octave_packages%2Ftsa-4.2.4%2Facovf.m;fp=octave_packages%2Ftsa-4.2.4%2Facovf.m;h=d89902da962cb7f59d734c4528f714a80bbea2b9;hp=0000000000000000000000000000000000000000;hb=c880e8788dfc484bf23ce13fa2787f2c6bca4863;hpb=1705066eceaaea976f010f669ce8e972f3734b05 diff --git a/octave_packages/tsa-4.2.4/acovf.m b/octave_packages/tsa-4.2.4/acovf.m new file mode 100644 index 0000000..d89902d --- /dev/null +++ b/octave_packages/tsa-4.2.4/acovf.m @@ -0,0 +1,147 @@ +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 [default] +% '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. + +% $Id: acovf.m 7687 2010-09-08 18:39:23Z schloegl $ +% Copyright (C) 1998-2003,2008,2010 by Alois Schloegl +% This is part of the TSA-toolbox. See also +% http://biosig-consulting.com/matlab/tsa/ +% +% This program is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% This program 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 General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with this program. If not, see . + + + +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;