1 function [a,e,REV,TOC,CPUTIME,ESU] = aar(y, Mode, arg3, arg4, arg5, arg6, arg7, arg8, arg9);
2 % Calculates adaptive autoregressive (AAR) and adaptive autoregressive moving average estimates (AARMA)
3 % of real-valued data series using Kalman filter algorithm.
4 % [a,e,REV] = aar(y, mode, MOP, UC, a0, A, W, V);
6 % The AAR process is described as following
7 % y(k) - a(k,1)*y(t-1) -...- a(k,p)*y(t-p) = e(k);
8 % The AARMA process is described as following
9 % y(k) - a(k,1)*y(t-1) -...- a(k,p)*y(t-p) = e(k) + b(k,1)*e(t-1) + ... + b(k,q)*e(t-q);
12 % y Signal (AR-Process)
13 % Mode is a two-element vector [aMode, vMode],
14 % aMode determines 1 (out of 12) methods for updating the co-variance matrix (see also [1])
15 % vMode determines 1 (out of 7) methods for estimating the innovation variance (see also [1])
16 % aMode=1, vmode=2 is the RLS algorithm as used in [2]
17 % aMode=-1, LMS algorithm (signal normalized)
18 % aMode=-2, LMS algorithm with adaptive normalization
20 % MOP model order, default [10,0]
21 % MOP=[p] AAR(p) model. p AR parameters
22 % MOP=[p,q] AARMA(p,q) model, p AR parameters and q MA coefficients
23 % UC Update Coefficient, default 0
24 % a0 Initial AAR parameters [a(0,1), a(0,2), ..., a(0,p),b(0,1),b(0,2), ..., b(0,q)]
25 % (row vector with p+q elements, default zeros(1,p) )
26 % A Initial Covariance matrix (positive definite pxp-matrix, default eye(p))
27 % W system noise (required for aMode==0)
28 % V observation noise (required for vMode==0)
31 % a AAR(MA) estimates [a(k,1), a(k,2), ..., a(k,p),b(k,1),b(k,2), ..., b(k,q]
32 % e error process (Adaptively filtered process)
33 % REV relative error variance MSE/MSY
37 % The mean square (prediction) error of different variants is useful for determining the free parameters (Mode, MOP, UC)
40 % [1] A. Schloegl (2000), The electroencephalogram and the adaptive autoregressive model: theory and applications.
41 % ISBN 3-8265-7640-3 Shaker Verlag, Aachen, Germany.
43 % More references can be found at
44 % http://www.dpmi.tu-graz.ac.at/~schloegl/publications/
47 % $Id: aar.m 8383 2011-07-16 20:06:59Z schloegl $
48 % Copyright (C) 1998-2003 by Alois Schloegl <a.schloegl@ieee.org>
50 % This program is free software: you can redistribute it and/or modify
51 % it under the terms of the GNU General Public License as published by
52 % the Free Software Foundation, either version 3 of the License, or
53 % (at your option) any later version.
55 % This program is distributed in the hope that it will be useful,
56 % but WITHOUT ANY WARRANTY; without even the implied warranty of
57 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
58 % GNU General Public License for more details.
60 % You should have received a copy of the GNU General Public License
61 % along with this program. If not, see <http://www.gnu.org/licenses/>.
65 %if nc<nr y=y'; end; tmp=nr;nc=nr; nr=tmp;end;
67 if nargin<2 Mode=0; end;
68 % check Mode (argument2)
69 if prod(size(Mode))==2
73 if any(aMode==(0:14)) && any(vMode==(0:7)),
74 fprintf(1,['a' int2str(aMode) 'e' int2str(vMode) ' ']);
76 fprintf(2,'Error AAR.M: invalid Mode argument\n');
80 % check model order (argument3)
81 if nargin<3 MOP=[10,0]; else MOP= arg3; end;
82 if length(MOP)==0 p=10; q=0; MOP=p;
83 elseif length(MOP)==1 p=MOP(1); q=0; MOP=p;
84 elseif length(MOP)>=2 p=MOP(1); q=MOP(2); MOP=p+q;
87 if nargin<4 UC=0; else UC= arg4; end;
92 if all(size(arg5)==([1,1]*(MOP+1))); % extended covariance matrix of AAR parameters
93 a0 = arg5(1,2:size(arg5,2));
94 A0 = arg5(2:size(arg5,1),2:size(arg5,2)) - a0'*a0;
103 if nargin<7, W = []; else W = arg7; end;
105 if all(size(W)==MOP),
107 fprintf(1,'aMode should be 0, because W is given.\n');
111 fprintf(1,'aMode must be non-zero, because W is not given.\n');
113 elseif any(size(W)~=MOP),
114 fprintf(1,'size of W does not fit. It must be %i x %i.\n',MOP,MOP);
118 if nargin<8, V0 = []; else V0 = arg8; end;
119 if all(size(V0)==nr),
121 fprintf(1,'vMode should be 0, because V is given.\n');
125 fprintf(1,'vMode must be non-zero, because V is not given.\n');
128 fprintf(1,'size of V does not fit. It must be 1x1.\n');
132 % if nargin<7 TH=3; else TH = arg7; end;
134 % TH=TH*mean(detrend(y,0).^2);
135 MSY=mean(detrend(y,0).^2);
148 mu=1-UC; % Patomaeki 1995
149 lambda=(1-UC); % Schloegl 1996
150 arc=poly((1-UC*2)*[1;1]);b0=sum(arc); % Whale forgettting factor for Mode=258,(Bianci et al. 1997)
152 dW=UC/MOP*eye(MOP); % Schloegl
155 %------------------------------------------------
157 %------------------------------------------------
168 V(1) = (1-UC) + UC*E*E;
174 %------------------------------------------------
176 %------------------------------------------------
181 %Y=[y(t-1); Y(1:p-1); E ; Y(p+1:MOP-1)]
183 if t<=p Y(1:t-1)=y(t-1:-1:1); % Autoregressive
184 else Y(1:p)=y(t-1:-1:t-p);
187 if t<=q Y(p+(1:t-1))=e(t-1:-1:1); % Moving Average
188 else Y(p+1:MOP)=e(t-1:-1:t-q);
192 E = y(t) - a(t-1,:)*Y;
203 V(t) = V(t-1)*(1-UC)+UC*E2;
204 if aMode == -1, % LMS
205 % V(t) = V(t-1)*(1-UC)+UC*E2;
206 a(t,:)=a(t-1,:) + (UC/MSY)*E*Y';
207 elseif aMode == -2, % LMS with adaptive estimation of the variance
208 a(t,:)=a(t-1,:) + UC/V(t)*E*Y';
210 else % Kalman filtering (including RLS)
211 if vMode==0, %eMode==4
213 elseif vMode==1, %eMode==4
215 elseif vMode==2, %eMode==2
217 elseif vMode==3, %eMode==3
218 Q(t) = (esu + lambda);
219 elseif vMode==4, %eMode==1
220 Q(t) = (esu + V(t-1));
221 elseif vMode==5, %eMode==6
223 V(t)=(1-UC)*V(t-1)+UC*(E2-esu);
228 elseif vMode==6, %eMode==7
230 V(t)=(1-UC)*V(t-1)+UC*(E2-esu);
234 Q(t) = (esu + V(t-1));
235 elseif vMode==7, %eMode==8
239 k = AY / Q(t); % Kalman Gain
240 a(t,:) = a(t-1,:) + k'*E;
242 if aMode==0, %AMode=0
243 A = A - k*AY' + W; % Schloegl et al. 2003
244 elseif aMode==1, %AMode=1
245 A = (1+UC)*(A - k*AY'); % Schloegl et al. 1997
246 elseif aMode==2, %AMode=11
248 A = A + sum(diag(A))*dW;
249 elseif aMode==3, %AMode=5
250 A = A - k*AY' + sum(diag(A))*dW;
251 elseif aMode==4, %AMode=6
252 A = A - k*AY' + UC*eye(MOP); % Schloegl 1998
253 elseif aMode==5, %AMode=2
254 A = A - k*AY' + UC*UC*eye(MOP);
255 elseif aMode==6, %AMode=2
256 T(t)=(1-UC)*T(t-1)+UC*(E2-Q(t))/(Y'*Y);
258 if T(t)>0 A=A+T(t)*eye(MOP); end;
259 elseif aMode==7, %AMode=6
260 T(t)=(1-UC)*T(t-1)+UC*(E2-Q(t))/(Y'*Y);
262 if T(t)>0 A=A+T(t)*eye(MOP); end;
263 elseif aMode==8, %AMode=5
264 Q_wo = (Y'*C*Y + V(t-1));
266 T(t)=(1-UC)*T(t-1)+UC*(E2-Q_wo)/(Y'*Y);
267 if T(t)>0 A=C+T(t)*eye(MOP); else A=C; end;
268 elseif aMode==9, %AMode=3
269 A = A - (1+UC)*k*AY';
270 A = A + sum(diag(A))*dW;
271 elseif aMode==10, %AMode=7
272 A = A - (1+UC)*k*AY' + sum(diag(A))*dW;
273 elseif aMode==11, %AMode=8
275 A = A - (1+UC)*k*AY' + UC*eye(MOP); % Schloegl 1998
276 elseif aMode==12, %AMode=4
277 A = A - (1+UC)*k*AY' + UC*UC*eye(MOP);
279 A = A - k*AY' + UC*diag(diag(A));
281 A = A - k*AY' + (UC*UC)*diag(diag(A));
289 CPUTIME = cputime - CPUTIME;
290 %REV = (e'*e)/(y'*y);
292 REV = mean(e.*e)./mean(y.*y);