1 function [z,e,REV,ESU,V,Z,SPUR] = amarma(y, Mode, MOP, UC, z0, Z0, V0, W);
2 % Adaptive Mean-AutoRegressive-Moving-Average model estimation
3 % [z,E,ESU,REV,V,Z,SPUR] = amarma(y, mode, MOP, UC, z0, Z0, V0, W);
4 % Estimates AAR parameters with Kalman filter algorithm
5 % y(t) = sum_i(a(i,t)*y(t-i)) + mu(t) + E(t)
8 % z(t)=G*z(t-1) + w(t) w(t)=N(0,W)
9 % y(t)=H*z(t) + v(t) v(t)=N(0,V)
12 % z = [µ(t)/(1-sum_i(a(i,t))),a_1(t-1),..,a_p(t-p),b_1(t-1),...,b_q(t-q)];
13 % H = [1,y(t-1),..,y(t-p),e(t-1),...,e(t-q)];
14 % W = E{(z(t)-G*z(t-1))*(z(t)-G*z(t-1))'}
15 % V = E{(y(t)-H*z(t-1))*(y(t)-H*z(t-1))'}
18 % y Signal (AR-Process)
22 % MOP Model order [m,p,q], default [0,10,0]
23 % m=1 includes the mean term, m=0 does not.
24 % p and q must be positive integers
25 % it is recommended to set q=0.
26 % UC Update Coefficient, default 0
27 % z0 Initial state vector
28 % Z0 Initial Covariance matrix
32 % E error process (Adaptively filtered process)
33 % 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.
41 % [2] Schlögl A, Lee FY, Bischof H, Pfurtscheller G
42 % Characterization of Four-Class Motor Imagery EEG Data for the BCI-Competition 2005.
43 % Journal of neural engineering 2 (2005) 4, S. L14-L22
45 % More references can be found at
46 % http://www.dpmi.tu-graz.ac.at/~schloegl/publications/
48 % $Id: amarma.m 5376 2008-10-13 15:53:47Z schloegl $
49 % Copyright (C) 1998-2002,2005,2006,2007,2008 by Alois Schloegl <a.schloegl@ieee.org>
51 % This program is free software: you can redistribute it and/or modify
52 % it under the terms of the GNU General Public License as published by
53 % the Free Software Foundation, either version 3 of the License, or
54 % (at your option) any later version.
56 % This program is distributed in the hope that it will be useful,
57 % but WITHOUT ANY WARRANTY; without even the implied warranty of
58 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
59 % GNU General Public License for more details.
61 % You should have received a copy of the GNU General Public License
62 % along with this program. If not, see <http://www.gnu.org/licenses/>.
68 elseif isnan(Mode) return; end;
69 if nargin<3, MOP=[0,10,0]; end;
70 if nargin<8, W = nan ; end;
71 if length(MOP)==0, m=0;p=10; q=0; MOP=p;
72 elseif length(MOP)==1, m=0;p=MOP(1); q=0; MOP=p;
73 elseif length(MOP)==2, fprintf(1,'Error AMARMA: MOP is ambiguos\n');
74 elseif length(MOP)>2, m=MOP(1); p=MOP(2); q=MOP(3);MOP=m+p+q;
81 %fprintf(1,['a' int2str(aMode) 'e' int2str(eMode) ' ']);
85 V = zeros(nc,1);V(1)=V0;
87 ESU = zeros(nc,1)+nan;
88 SPUR = zeros(nc,1)+nan;
91 dW = UC/MOP*eye(MOP); % Schloegl
93 %------------------------------------------------
95 %------------------------------------------------
99 %M0 = z0(1)/(1-sum(z0(2:p+1))); %transformierter Mittelwert
107 A1 = zeros(MOP); A2 = A1;
109 %------------------------------------------------
111 %------------------------------------------------
114 %H=[y(t-1); H(1:p-1); E ; H(p+1:MOP-1)]
117 if t<=p, H(m+(1:t-1)) = y(t-1:-1:1); %H(p)=mu0; % Autoregressive
118 else H(m+(1:p)) = y(t-1:-1:t-p); %mu0];
121 if t<=q, H(m+p+(1:t-1)) = e(t-1:-1:1); % Moving Average
122 else H(m+p+(1:q)) = e(t-1:-1:t-q);
134 % [zt, t, y(t), E,ESU(t),V(t),H,Z],pause,
142 V(t) = V0*(1-UC)+UC*E2;
145 V(t) = V0; %V(t-1)*(1-UC)+UC*E2;
148 V(t) = V0; %(t-1)*(1-UC)+UC*E2;
150 V0 = V0*(1-UC)+UC*E2;
157 V0=(1-UC)*V0+UC*(E2-ESU(t));
163 V(t) = (1-UC)*V0+UC*(E2-ESU(t));
169 V(t) = V0; % (t-1)*(1-UC)+UC*E2;
174 k = AY / (ESU(t) + V0); % Kalman Gain
179 %W = W; %nop % Schloegl et al. 2003
181 T(t)=(1-UC)*T(t-1)+UC*(E2-Q(t))/(H'*H); % Roberts I 1998
183 if T(t)>0 W=T(t)*eye(MOP); else W=zeros(MOP);end;
185 Q_wo = (H'*C*H + V(t-1)); % Roberts II 1998
186 T(t)=(1-UC)*T(t-1)+UC*(E2-Q_wo)/(H'*H);
187 if T(t)>0 W=T(t)*eye(MOP); else W=zeros(MOP); end;
189 T(t)=(1-UC)*T(t-1)+UC*(E2-Q(t))/(H'*H);
191 if T(t)>0 W=T(t)*eye(MOP); else W=zeros(MOP); end;
198 W = UC*diag(diag(Z));
200 W = (UC*UC)*diag(diag(Z));
204 W = UC*eye(MOP); % Schloegl 1998
213 Z = Z - k*AY'; % Schloegl 1998
219 if any(any(isnan(W))), W=UC*Z; end;
222 Z = Z + W; % Schloegl 1998
228 z(:,1)=M0*z(:,1)./(1-sum(z(:,2:p),2));
231 REV = mean(e.*e)/mean(y.*y);
232 if any(~isfinite(Z(:))), REV=inf; end;