www.gusucode.com > 时间序列分析工具箱 - tsa源码程序 > tsa/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: 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;