aar.zip_AAR matlab code_adaptive kalman_autoregressive _matlab k

function [a,e,REV,TOC,CPUTIME,ESU] = aar(y, Mode, arg3, arg4, arg5, arg6, arg7, arg8, arg9); % Calculates adaptive autoregressive (AAR) and adaptive autoregressive moving average estimates (AARMA) % of real-valued data series using Kalman filter algorithm. % [a,e,REV] = aar(y, mode, MOP, UC, a0, A, W, V); % % The AAR process is described as following % y(k) - a(k,1)*y(t-1) -...- a(k,p)*y(t-p) = e(k); % The AARMA process is described as following % 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); % % Input: % y Signal (AR-Process) % Mode is a two-element vector [aMode, vMode], % aMode determines 1 (out of 12) methods for updating the co-variance matrix (see also [1]) % vMode determines 1 (out of 7) methods for estimating the innovation variance (see also [1]) % aMode=1, vmode=2 is the RLS algorithm as used in [2] % aMode=-1, LMS algorithm (signal normalized) % aMode=-2, LMS algorithm with adaptive normalization % % MOP model order, default [10,0] % MOP=[p] AAR(p) model. p AR parameters % MOP=[p,q] AARMA(p,q) model, p AR parameters and q MA coefficients % UC Update Coefficient, default 0 % a0 Initial AAR parameters [a(0,1), a(0,2), ..., a(0,p),b(0,1),b(0,2), ..., b(0,q)] % (row vector with p+q elements, default zeros(1,p) ) % A Initial Covariance matrix (positive definite pxp-matrix, default eye(p)) % W system noise (required for aMode==0) % V observation noise (required for vMode==0) % % Output: % a AAR(MA) estimates [a(k,1), a(k,2), ..., a(k,p),b(k,1),b(k,2), ..., b(k,q] % e error process (Adaptively filtered process) % REV relative error variance MSE/MSY % % % Hint: % The mean square (prediction) error of different variants is useful for determining the free parameters (Mode, MOP, UC) % % REFERENCE(S): % [1] A. Schloegl (2000), The electroencephalogram and the adaptive autoregressive model: theory and applications. % ISBN 3-8265-7640-3 Shaker Verlag, Aachen, Germany. % % More references can be found at % http://www.dpmi.tu-graz.ac.at/~schloegl/publications/ % % $Revision: 1.10 $ % $Id: aar.m,v 1.10 2005/05/25 02:52:38 pkienzle 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. % for compilation with the Matlab compiler mcc -V1.2 -ir %#realonly %#inbounds [nc,nr]=size(y); %if nc<nr y=y'; end; tmp=nr;nc=nr; nr=tmp;end; if nargin<2 Mode=0; end; % check Mode (argument2) if prod(size(Mode))==2 aMode=Mode(1); vMode=Mode(2); end; if any(aMode==(0:14)) & any(vMode==(0:7)), fprintf(1,['a' int2str(aMode) 'e' int2str(vMode) ' ']); else fprintf(2,'Error AAR.M: invalid Mode argument\n'); return; end; % check model order (argument3) if nargin<3 MOP=[10,0]; else MOP= arg3; end; if length(MOP)==0 p=10; q=0; MOP=p; elseif length(MOP)==1 p=MOP(1); q=0; MOP=p; elseif length(MOP)>=2 p=MOP(1); q=MOP(2); MOP=p+q; end; if nargin<4 UC=0; else UC= arg4; end; a0=zeros(1,MOP); A0=eye(MOP); if nargin>4, if all(size(arg5)==([1,1]*(MOP+1))); % extended covariance matrix of AAR parameters a0 = arg5(1,2:size(arg5,2)); A0 = arg5(2:size(arg5,1),2:size(arg5,2)) - a0'*a0; else a0 = arg5; if nargin>5 A0 = arg6; end; end; end; if nargin<7, W = []; else W = arg7; end; if all(size(W)==MOP), if aMode ~= 0, fprintf(1,'aMode should be 0, because W is given.\n'); end; elseif isempty(W), if aMode == 0, fprintf(1,'aMode must be non-zero, because W is not given.\n'); end; elseif any(size(W)~=MOP), fprintf(1,'size of W does not fit. It must be %i x %i.\n',MOP,MOP); return; end; if nargin<8, V0 = []; else V0 = arg8; end; if all(size(V0)==nr), if vMode ~= 0, fprintf(1,'vMode should be 0, because V is given.\n'); end; elseif isempty(V0), if aMode == 0, fprintf(1,'vMode must be non-zero, because V is not given.\n'); end; else fprintf(1,'size of V does not fit. It must be 1x1.\n'); return; end; % if nargin<7 TH=3; else TH = arg7; end; % TH=TH*var(y); % TH=TH*mean(detrend(y,0).^2); MSY=mean(detrend(y,0).^2); e=zeros(nc,1); Q=zeros(nc,1); V=zeros(nc,1); T=zeros(nc,1); %DET=zeros(nc,1); SPUR=zeros(nc,1); ESU=zeros(nc,1); a=a0(ones(nc,1),:); %a=zeros(nc,MOP); %b=zeros(nc,q); mu=1-UC; % Patomaeki 1995 lambda=(1-UC); % Schloegl 1996 arc=poly((1-UC*2)*[1;1]);b0=sum(arc); % Whale forgettting factor for Mode=258,(Bianci et al. 1997) dW=UC/MOP*eye(MOP); % Schloegl %------------------------------------------------ % First Iteration %------------------------------------------------ Y=zeros(MOP,1); C=zeros(MOP); %X=zeros(q,1); at=a0; A=A0; E=y(1); e(1)=E; if ~isempty(V0) V(1) = V0; else V(1) = (1-UC) + UC*E*E; end; ESU(1) = 1; %Y'*A*Y; A1=zeros(MOP);A2=A1; tic;CPUTIME=cputime; %------------------------------------------------ % Update Equations %------------------------------------------------ T0=2; for t=T0:nc, %Y=[y(t-1); Y(1:p-1); E ; Y(p+1:MOP-1)] if t<=p Y(1:t-1)=y(t-1:-1:1); % Autoregressive else Y(1:p)=y(t-1:-1:t-p); end; if t<=q Y(p+(1:t-1))=e(t-1:-1:1); % Moving Average else Y(p+1:MOP)=e(t-1:-1:t-q); end; % Prediction Error E = y(t) - a(t-1,:)*Y; e(t) = E; E2=E*E; AY=A*Y; esu=Y'*AY; ESU(t)=esu; if isnan(E), a(t,:)=a(t-1,:); else V(t) = V(t-1)*(1-UC)+UC*E2; if aMode == -1, % LMS % V(t) = V(t-1)*(1-UC)+UC*E2; a(t,:)=a(t-1,:) + (UC/MSY)*E*Y'; elseif aMode == -2, % LMS with adaptive estimation of the variance a(t,:)=a(t-1,:) + UC/V(t)*E*Y'; else % Kalman filtering (including RLS) if vMode==0, %eMode==4 Q(t) = (esu + V0); elseif vMode==1, %eMode==4 Q(t) = (esu + V(t)); elseif vMode==2, %eMode==2 Q(t) = (esu + 1); elseif vMode==3, %eMode==3 Q(t) = (esu + lambda); elseif vMode==4, %eMode==1 Q(t) = (esu + V(t-1)); elseif vMode==5, %eMode==6 if E2>esu V(t)=(1-UC)*V(t-1)+UC*(E2-esu); else V(t)=V(t-1); end;