1 function [z,e,REV,ESU,V,Z,SPUR] = aarmam(y, Mode, MOP, UC, z0, Z0, V0, W);
2 % Estimating Adaptive AutoRegressive-Moving-Average-and-mean model (includes mean term)
4 % !! This function is obsolete and is replaced by AMARMA
6 % [z,E,REV,ESU,V,Z,SPUR] = aarmam(y, mode, MOP, UC, z0, Z0, V0, W);
7 % Estimates AAR parameters with Kalman filter algorithm
8 % y(t) = sum_i(a_i(t)*y(t-i)) + m(t) + e(t) + sum_i(b_i(t)*e(t-i))
11 % z(t) = G*z(t-1) + w(t) w(t)=N(0,W)
12 % y(t) = H*z(t) + v(t) v(t)=N(0,V)
15 % z = [m(t),a_1(t-1),..,a_p(t-p),b_1(t-1),...,b_q(t-q)];
16 % H = [1,y(t-1),..,y(t-p),e(t-1),...,e(t-q)];
17 % W = E{(z(t)-G*z(t-1))*(z(t)-G*z(t-1))'}
18 % V = E{(y(t)-H*z(t-1))*(y(t)-H*z(t-1))'}
22 % y Signal (AR-Process)
23 % Mode determines the type of algorithm
25 % MOP Model order [m,p,q], default [0,10,0]
26 % m=1 includes the mean term, m=0 does not.
27 % p and q must be positive integers
28 % it is recommended to set q=0.
29 % UC Update Coefficient, default 0
30 % z0 Initial state vector
31 % Z0 Initial Covariance matrix
35 % E error process (Adaptively filtered process)
36 % REV relative error variance MSE/MSY
39 % [1] A. Schloegl (2000), The electroencephalogram and the adaptive autoregressive model: theory and applications.
40 % ISBN 3-8265-7640-3 Shaker Verlag, Aachen, Germany.
42 % More references can be found at
43 % http://pub.ist.ac.at/~schloegl/publications/
45 % $Id: aarmam.m 9609 2012-02-10 10:18:00Z schloegl $
46 % Copyright (C) 1998-2002,2008,2012 by Alois Schloegl <alois.schloegl@ist.ac.at>
47 % This is part of the TSA-toolbox. See also
48 % http://pub.ist.ac.at/~schloegl/matlab/tsa/
49 % http://octave.sourceforge.net/
50 % http://biosig.sourceforge.net/
52 % This program is free software: you can redistribute it and/or modify
53 % it under the terms of the GNU General Public License as published by
54 % the Free Software Foundation, either version 3 of the License, or
55 % (at your option) any later version.
57 % This program is distributed in the hope that it will be useful,
58 % but WITHOUT ANY WARRANTY; without even the implied warranty of
59 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
60 % GNU General Public License for more details.
62 % You should have received a copy of the GNU General Public License
63 % along with this program. If not, see <http://www.gnu.org/licenses/>.
68 warning('AARMAM is obsolete. Use AMARMA instead!')
73 elseif ischar(Mode) Mode=bin2dec(Mode);
74 elseif isnan(Mode) return; end;
75 if nargin<3 MOP=[0,10,0]; end;
76 if length(MOP)==0, m=0;p=10; q=0; MOP=p;
77 elseif length(MOP)==1, m=0;p=MOP(1); q=0; MOP=p;
78 elseif length(MOP)==2, fprintf(1,'Error AMARMA: MOP is ambiguos\n');
79 elseif length(MOP)>2, m=MOP(1); p=MOP(2); q=MOP(3);MOP=m+p+q;
86 %fprintf(1,['a' int2str(aMode) 'e' int2str(eMode) ' ']);
90 V=zeros(nc,1);V(1)=V0;
96 arc=poly((1-UC*2)*[1;1]);b0=sum(arc); % Whale forgetting factor for Mode=258,(Bianci et al. 1997)
98 dW=UC/MOP*eye(MOP); % Schloegl
101 %------------------------------------------------
103 %------------------------------------------------
109 fprintf(2,'Warning AARMAM: m must be 0 or 1\n');
113 if (p<0) || (q<0) || (round(p)~=p) || (round(q)~=q),
114 fprintf(2,'Error AARMAM: p and q must be positive integers\n');
122 A1 = zeros(MOP); A2 = A1;
126 %------------------------------------------------
128 %------------------------------------------------
132 % make measurement matrix
137 H=[1; y_1; H(m+(1:p-1)'); E(1:min(1,q-1)) ; H(p+m+(1:q-1)')]; % shift y and e
139 else % this seem to be slightly faster
140 if t<=p, H(m+(1:t-1)) = y(t-1:-1:1); % Autoregressive
141 else H(m+(1:p)) = y(t-1:-1:t-p);
144 if t<=q, H(m+p+(1:t-1)) = e(t-1:-1:1); % Moving Average
145 else H(m+p+(1:q)) = e(t-1:-1:t-q);
161 V(t) = V0*(1-UC)+UC*E2;
164 V(t) = V0; %V(t-1)*(1-UC)+UC*E2;
167 V(t) = V0; %(t-1)*(1-UC)+UC*E2;
169 V0 = V0*(1-UC)+UC*E2;
176 V0=(1-UC)*V0+UC*(E2-ESU(t));
182 V(t) = (1-UC)*V0+UC*(E2-ESU(t));
188 V(t) = V0; % (t-1)*(1-UC)+UC*E2;
191 k = AY / (ESU(t) + V0); % Kalman Gain
196 T(t)=(1-UC)*T(t-1)+UC*(E2-Q(t))/(H'*H); % Roberts I 1998
198 if T(t)>0 W=T(t)*eye(MOP); else W=zeros(MOP);end;
200 Q_wo = (H'*C*H + V(t-1)); % Roberts II 1998
201 T(t)=(1-UC)*T(t-1)+UC*(E2-Q_wo)/(H'*H);
202 if T(t)>0 W=T(t)*eye(MOP); else W=zeros(MOP); end;
204 T(t)=(1-UC)*T(t-1)+UC*(E2-Q(t))/(H'*H);
206 if T(t)>0 W=T(t)*eye(MOP); else W=zeros(MOP); end;
213 W = UC*diag(diag(Z));
215 W = (UC*UC)*diag(diag(Z));
219 W = UC*eye(MOP); % Schloegl 1998
224 Z = Z - k*AY'; % Schloegl 1998
231 if any(any(isnan(W))), W=UC*Z;end;
234 Z = Z + W; % Schloegl 1998
238 REV = mean(e.*e)/mean(y.*y);
239 if any(~isfinite(Z(:))), REV=inf; end;