|
2414
|
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 |
|
7501
|
10 % Mode 'biased' : normalizes with N [default] |
|
2414
|
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$ |
|
12580
|
27 % Copyright (C) 1998-2003,2008,2010 by Alois Schloegl <alois.schloegl@gmail.com> |
|
7501
|
28 % This is part of the TSA-toolbox. See also |
|
|
29 % http://biosig-consulting.com/matlab/tsa/ |
|
4907
|
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/>. |
|
2414
|
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; |
