X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=octave_packages%2Ftsa-4.2.4%2Famarma.m;fp=octave_packages%2Ftsa-4.2.4%2Famarma.m;h=8cf9f65571b3099e2887f065e39595458467be1e;hb=c880e8788dfc484bf23ce13fa2787f2c6bca4863;hp=0000000000000000000000000000000000000000;hpb=1705066eceaaea976f010f669ce8e972f3734b05;p=CreaPhase.git diff --git a/octave_packages/tsa-4.2.4/amarma.m b/octave_packages/tsa-4.2.4/amarma.m new file mode 100644 index 0000000..8cf9f65 --- /dev/null +++ b/octave_packages/tsa-4.2.4/amarma.m @@ -0,0 +1,233 @@ +function [z,e,REV,ESU,V,Z,SPUR] = amarma(y, Mode, MOP, UC, z0, Z0, V0, W); +% Adaptive Mean-AutoRegressive-Moving-Average model estimation +% [z,E,ESU,REV,V,Z,SPUR] = amarma(y, mode, MOP, UC, z0, Z0, V0, W); +% Estimates AAR parameters with Kalman filter algorithm +% y(t) = sum_i(a(i,t)*y(t-i)) + mu(t) + E(t) +% +% State space model: +% z(t)=G*z(t-1) + w(t) w(t)=N(0,W) +% y(t)=H*z(t) + v(t) v(t)=N(0,V) +% +% G = I, +% z = [µ(t)/(1-sum_i(a(i,t))),a_1(t-1),..,a_p(t-p),b_1(t-1),...,b_q(t-q)]; +% H = [1,y(t-1),..,y(t-p),e(t-1),...,e(t-q)]; +% W = E{(z(t)-G*z(t-1))*(z(t)-G*z(t-1))'} +% V = E{(y(t)-H*z(t-1))*(y(t)-H*z(t-1))'} +% +% Input: +% y Signal (AR-Process) +% Mode +% [0,0] uses V0 and W +% +% MOP Model order [m,p,q], default [0,10,0] +% m=1 includes the mean term, m=0 does not. +% p and q must be positive integers +% it is recommended to set q=0. +% UC Update Coefficient, default 0 +% z0 Initial state vector +% Z0 Initial Covariance matrix +% +% Output: +% z AR-Parameter +% E error process (Adaptively filtered process) +% REV relative error variance MSE/MSY +% +% +% see also: AAR +% +% 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. +% [2] Schlögl A, Lee FY, Bischof H, Pfurtscheller G +% Characterization of Four-Class Motor Imagery EEG Data for the BCI-Competition 2005. +% Journal of neural engineering 2 (2005) 4, S. L14-L22 +% +% More references can be found at +% http://www.dpmi.tu-graz.ac.at/~schloegl/publications/ + +% $Id: amarma.m 5376 2008-10-13 15:53:47Z schloegl $ +% Copyright (C) 1998-2002,2005,2006,2007,2008 by Alois Schloegl +% +% This program is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% This program 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 General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with this program. If not, see . + + +[nc,nr]=size(y); + +if nargin<2 Mode=0; +elseif isnan(Mode) return; end; +if nargin<3, MOP=[0,10,0]; end; +if nargin<8, W = nan ; end; +if length(MOP)==0, m=0;p=10; q=0; MOP=p; +elseif length(MOP)==1, m=0;p=MOP(1); q=0; MOP=p; +elseif length(MOP)==2, fprintf(1,'Error AMARMA: MOP is ambiguos\n'); +elseif length(MOP)>2, m=MOP(1); p=MOP(2); q=MOP(3);MOP=m+p+q; +end; + +if prod(size(Mode))>1 + aMode=Mode(1); + eMode=Mode(2); +end; +%fprintf(1,['a' int2str(aMode) 'e' int2str(eMode) ' ']); + + +e = zeros(nc,1); +V = zeros(nc,1);V(1)=V0; +T = zeros(nc,1); +ESU = zeros(nc,1)+nan; +SPUR = zeros(nc,1)+nan; +z = z0(ones(nc,1),:); + +dW = UC/MOP*eye(MOP); % Schloegl + +%------------------------------------------------ +% First Iteration +%------------------------------------------------ + +H = zeros(MOP,1); +if m, + %M0 = z0(1)/(1-sum(z0(2:p+1))); %transformierter Mittelwert + H(1) = 1;%M0; + %z0(1)= 1; +end; + +Z = Z0; +zt= z0; + +A1 = zeros(MOP); A2 = A1; + +%------------------------------------------------ +% Update Equations +%------------------------------------------------ + +for t=1:nc, + %H=[y(t-1); H(1:p-1); E ; H(p+1:MOP-1)] + + + if t<=p, H(m+(1:t-1)) = y(t-1:-1:1); %H(p)=mu0; % Autoregressive + else H(m+(1:p)) = y(t-1:-1:t-p); %mu0]; + end; + + if t<=q, H(m+p+(1:t-1)) = e(t-1:-1:1); % Moving Average + else H(m+p+(1:q)) = e(t-1:-1:t-q); + end; + + % Prediction Error + E = y(t) - zt*H; + + e(t) = E; + + if ~isnan(E), + E2 = E*E; + AY = Z*H; + +% [zt, t, y(t), E,ESU(t),V(t),H,Z],pause, + + ESU(t) = H'*AY; + + if eMode==0 + V(t) = V0; + elseif eMode==1 + V0 = V(t-1); + V(t) = V0*(1-UC)+UC*E2; + elseif eMode==2 + V0 = 1; + V(t) = V0; %V(t-1)*(1-UC)+UC*E2; + elseif eMode==3 + V0 = 1-UC; + V(t) = V0; %(t-1)*(1-UC)+UC*E2; + elseif eMode==4 + V0 = V0*(1-UC)+UC*E2; + V(t) = V0; + elseif eMode==5 + V(t)=V0; + %V0 = V0; + elseif eMode==6 + if E2>ESU(t) + V0=(1-UC)*V0+UC*(E2-ESU(t)); + end; + V(t)=V0; + elseif eMode==7 + V0=V(t); + if E2>ESU(t) + V(t) = (1-UC)*V0+UC*(E2-ESU(t)); + else + V(t) = V0; + end; + elseif eMode==8 + V0=0; + V(t) = V0; % (t-1)*(1-UC)+UC*E2; + end; + +%[t,size(H),size(Z)] + + k = AY / (ESU(t) + V0); % Kalman Gain + zt = zt + k'*E; + %z(t,:) = zt; + + if aMode==0 + %W = W; %nop % Schloegl et al. 2003 + elseif aMode==2 + T(t)=(1-UC)*T(t-1)+UC*(E2-Q(t))/(H'*H); % Roberts I 1998 + Z=Z*V(t-1)/Q(t); + if T(t)>0 W=T(t)*eye(MOP); else W=zeros(MOP);end; + elseif aMode==5 + Q_wo = (H'*C*H + V(t-1)); % Roberts II 1998 + T(t)=(1-UC)*T(t-1)+UC*(E2-Q_wo)/(H'*H); + if T(t)>0 W=T(t)*eye(MOP); else W=zeros(MOP); end; + elseif aMode==6 + T(t)=(1-UC)*T(t-1)+UC*(E2-Q(t))/(H'*H); + Z=Z*V(t)/Q(t); + if T(t)>0 W=T(t)*eye(MOP); else W=zeros(MOP); end; + elseif aMode==11 + %Z = Z - k*AY'; + W = sum(diag(Z))*dW; + elseif aMode==12 + W = UC*UC*eye(MOP); + elseif aMode==13 + W = UC*diag(diag(Z)); + elseif aMode==14 + W = (UC*UC)*diag(diag(Z)); + elseif aMode==15 + W = sum(diag(Z))*dW; + elseif aMode==16 + W = UC*eye(MOP); % Schloegl 1998 + elseif aMode==17 + Z = 0.5*(Z+Z'); + W = UC*Z; + elseif aMode==18 + W = 0.5*UC*(Z+Z'); + %W=W; + end; + + Z = Z - k*AY'; % Schloegl 1998 + else + + V(t) = V0; + + end; + if any(any(isnan(W))), W=UC*Z; end; + + z(t,:) = zt; + Z = Z + W; % Schloegl 1998 + SPUR(t)=trace(Z); +end; + + +if 0,m, + z(:,1)=M0*z(:,1)./(1-sum(z(:,2:p),2)); +end; + +REV = mean(e.*e)/mean(y.*y); +if any(~isfinite(Z(:))), REV=inf; end; +