1function [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); 5% 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); 10% 11% Input: 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 19% 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) 29% 30% Output: 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 34% 35% 36% Hint: 37% The mean square (prediction) error of different variants is useful for determining the free parameters (Mode, MOP, UC) 38% 39% REFERENCE(S): 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. 42% 43% More references can be found at 44% http://pub.ist.ac.at/~schloegl/publications/ 45 46% 47% $Id$ 48% Copyright (C) 1998-2003 by Alois Schloegl <alois.schloegl@gmail.com> 49% 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. 54% 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. 59% 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/>. 62 63 64[nc,nr]=size(y); 65%if nc<nr y=y'; end; tmp=nr;nc=nr; nr=tmp;end; 66 67if nargin<2 Mode=0; end; 68% check Mode (argument2) 69if prod(size(Mode))==2 70 aMode=Mode(1); 71 vMode=Mode(2); 72end; 73if any(aMode==(0:14)) && any(vMode==(0:7)), 74 fprintf(1,['a' int2str(aMode) 'e' int2str(vMode) ' ']); 75else 76 fprintf(2,'Error AAR.M: invalid Mode argument\n'); 77 return; 78end; 79 80% check model order (argument3) 81if nargin<3 MOP=[10,0]; else MOP= arg3; end; 82if length(MOP)==0 p=10; q=0; MOP=p; 83elseif length(MOP)==1 p=MOP(1); q=0; MOP=p; 84elseif length(MOP)>=2 p=MOP(1); q=MOP(2); MOP=p+q; 85end; 86 87if nargin<4 UC=0; else UC= arg4; end; 88 89a0=zeros(1,MOP); 90A0=eye(MOP); 91if nargin>4, 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; 95 else 96 a0 = arg5; 97 if nargin>5 98 A0 = arg6; 99 end; 100 end; 101end; 102 103if nargin<7, W = []; else W = arg7; end; 104 105if all(size(W)==MOP), 106 if aMode ~= 0, 107 fprintf(1,'aMode should be 0, because W is given.\n'); 108 end; 109elseif isempty(W), 110 if aMode == 0, 111 fprintf(1,'aMode must be non-zero, because W is not given.\n'); 112 end; 113elseif any(size(W)~=MOP), 114 fprintf(1,'size of W does not fit. It must be %i x %i.\n',MOP,MOP); 115 return; 116end; 117 118if nargin<8, V0 = []; else V0 = arg8; end; 119if all(size(V0)==nr), 120 if vMode ~= 0, 121 fprintf(1,'vMode should be 0, because V is given.\n'); 122 end; 123elseif isempty(V0), 124 if aMode == 0, 125 fprintf(1,'vMode must be non-zero, because V is not given.\n'); 126 end; 127else 128 fprintf(1,'size of V does not fit. It must be 1x1.\n'); 129 return; 130end; 131 132% if nargin<7 TH=3; else TH = arg7; end; 133% TH=TH*var(y); 134% TH=TH*mean(detrend(y,0).^2); 135MSY=mean(detrend(y,0).^2); 136 137e=zeros(nc,1); 138Q=zeros(nc,1); 139V=zeros(nc,1); 140T=zeros(nc,1); 141%DET=zeros(nc,1); 142SPUR=zeros(nc,1); 143ESU=zeros(nc,1); 144a=a0(ones(nc,1),:); 145%a=zeros(nc,MOP); 146%b=zeros(nc,q); 147 148mu=1-UC; % Patomaeki 1995 149lambda=(1-UC); % Schloegl 1996 150arc=poly((1-UC*2)*[1;1]);b0=sum(arc); % Whale forgettting factor for Mode=258,(Bianci et al. 1997) 151 152dW=UC/MOP*eye(MOP); % Schloegl 153 154 155%------------------------------------------------ 156% First Iteration 157%------------------------------------------------ 158Y=zeros(MOP,1); 159C=zeros(MOP); 160%X=zeros(q,1); 161at=a0; 162A=A0; 163E=y(1); 164e(1)=E; 165if ~isempty(V0) 166 V(1) = V0; 167else 168 V(1) = (1-UC) + UC*E*E; 169end; 170ESU(1) = 1; %Y'*A*Y; 171 172A1=zeros(MOP);A2=A1; 173tic;CPUTIME=cputime; 174%------------------------------------------------ 175% Update Equations 176%------------------------------------------------ 177T0=2; 178 179for t=T0:nc, 180 181 %Y=[y(t-1); Y(1:p-1); E ; Y(p+1:MOP-1)] 182 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); 185 end; 186 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); 189 end; 190 191 % Prediction Error 192 E = y(t) - a(t-1,:)*Y; 193 e(t) = E; 194 E2=E*E; 195 196 AY=A*Y; 197 esu=Y'*AY; 198 ESU(t)=esu; 199 200 if isnan(E), 201 a(t,:)=a(t-1,:); 202 else 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'; 209 210 else % Kalman filtering (including RLS) 211 if vMode==0, %eMode==4 212 Q(t) = (esu + V0); 213 elseif vMode==1, %eMode==4 214 Q(t) = (esu + V(t)); 215 elseif vMode==2, %eMode==2 216 Q(t) = (esu + 1); 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 222 if E2>esu 223 V(t)=(1-UC)*V(t-1)+UC*(E2-esu); 224 else 225 V(t)=V(t-1); 226 end; 227 Q(t) = (esu + V(t)); 228 elseif vMode==6, %eMode==7 229 if E2>esu 230 V(t)=(1-UC)*V(t-1)+UC*(E2-esu); 231 else 232 V(t)=V(t-1); 233 end; 234 Q(t) = (esu + V(t-1)); 235 elseif vMode==7, %eMode==8 236 Q(t) = esu; 237 end; 238 239 k = AY / Q(t); % Kalman Gain 240 a(t,:) = a(t-1,:) + k'*E; 241 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 247 A = A - k*AY'; 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); 257 A=A*V(t-1)/Q(t); 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); 261 A=A*V(t)/Q(t); 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)); 265 C=A-k*AY'; 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 274 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); 278 elseif aMode==13 279 A = A - k*AY' + UC*diag(diag(A)); 280 elseif aMode==14 281 A = A - k*AY' + (UC*UC)*diag(diag(A)); 282 end; 283 end; 284 end; 285end; 286 287%a=a(end,:); 288TOC = toc; 289CPUTIME = cputime - CPUTIME; 290%REV = (e'*e)/(y'*y); 291 292REV = mean(e.*e)./mean(y.*y); 293