function [ARF,RCF,PE,DC,varargout] = mvar(Y, Pmax, Mode)
% MVAR estimates Multi-Variate AutoRegressive model parameters
% Several estimation algorithms are implemented, all estimators
% can handle data with missing values encoded as NaNs.
%
% [AR,RC,PE] = mvar(Y, p);
% [AR,RC,PE] = mvar(Y, p, Mode);
%
% INPUT:
% Y Multivariate data series
% p Model order
% Mode determines estimation algorithm
%
% OUTPUT:
% AR multivariate autoregressive model parameter
% RC reflection coefficients (= -PARCOR coefficients)
% PE remaining error variance
%
% All input and output parameters are organized in columns, one column
% corresponds to the parameters of one channel.
%
% Mode determines estimation algorithm.
% 1: Correlation Function Estimation method using biased correlation function estimation method
% also called the "multichannel Yule-Walker" [1,2]
% 6: Correlation Function Estimation method using unbiased correlation function estimation method
%
% 2: Partial Correlation Estimation: Vieira-Morf [2] using unbiased covariance estimates.
% In [1] this mode was used and (incorrectly) denominated as Nutall-Strand.
% Its the DEFAULT mode; according to [1] it provides the most accurate estimates.
% 5: Partial Correlation Estimation: Vieira-Morf [2] using biased covariance estimates.
% Yields similar results than Mode=2;
%
% 3: Partial Correlation Estimation: Nutall-Strand [2] (biased correlation function)
% 7: Partial Correlation Estimation: Nutall-Strand [2] (unbiased correlation function)
% 8: Least-Squares
%
% 10: ARFIT [3]
% 11: BURGV [4]
%
% REFERENCES:
% [1] A. Schl\"ogl, Comparison of Multivariate Autoregressive Estimators.
% Signal processing, Elsevier B.V. (in press).
% available at http://dx.doi.org/doi:10.1016/j.sigpro.2005.11.007
% [2] S.L. Marple "Digital Spectral Analysis with Applications" Prentice Hall, 1987.
% [3] Schneider and Neumaier)
% [4] Stijn de Waele, 2003
%
%
% A multivariate inverse filter can be realized with
% [AR,RC,PE] = mvar(Y,P);
% e = mvfilter([eye(size(AR,1)),-AR],eye(size(AR,1)),Y);
%
% see also: MVFILTER, MVFREQZ, COVM, SUMSKIPNAN, ARFIT2
% $Id: mvar.m 5090 2008-06-05 08:12:04Z schloegl $
% Copyright (C) 1996-2006 by Alois Schloegl <a.schloegl@ieee.org>
% This is part of the TSA-toolbox. See also
% http://hci.tugraz.at/schloegl/matlab/tsa/
% http://octave.sourceforge.net/
% http://biosig.sourceforge.net/
%
% Inititialization
[N,M] = size(Y);
if nargin<2,
Pmax=max([N,M])-1;
end;
if iscell(Y)
Pmax = min(max(N ,M ),Pmax);
C = Y;
end;
if nargin<3,
% according to [1], and other internal validations this is in many cases the best algorithm
Mode=2;
end;
[C(:,1:M),n] = covm(Y,'M');
PE(:,1:M) = C(:,1:M)./n;
if 0,
elseif Mode==0; % this method is broken
fprintf('Warning MVAR: Mode=0 is broken.\n')
C(:,1:M) = C(:,1:M)/N;
F = Y;
B = Y;
PEF = C(:,1:M); %?% PEF = PE(:,1:M);
PEB = C(:,1:M);
for K=1:Pmax,
[D,n] = covm(Y(K+1:N,:),Y(1:N-K,:),'M');
D = D/N;
ARF(:,K*M+(1-M:0)) = D/PEB;
ARB(:,K*M+(1-M:0)) = D'/PEF;
tmp = F(K+1:N,:) - B(1:N-K,:)*ARF(:,K*M+(1-M:0))';
B(1:N-K,:) = B(1:N-K,:) - F(K+1:N,:)*ARB(:,K*M+(1-M:0))';
F(K+1:N,:) = tmp;
for L = 1:K-1,
tmp = ARF(:,L*M+(1-M:0)) - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
ARB(:,(K-L)*M+(1-M:0)) = ARB(:,(K-L)*M+(1-M:0)) - ARB(:,K*M+(1-M:0))*ARF(:,L*M+(1-M:0));
ARF(:,L*M+(1-M:0)) = tmp;
end;
RCF(:,K*M+(1-M:0)) = ARF(:,K*M+(1-M:0));
RCB(:,K*M+(1-M:0)) = ARB(:,K*M+(1-M:0));
PEF = [eye(M) - ARF(:,K*M+(1-M:0))*ARB(:,K*M+(1-M:0))]*PEF;
PEB = [eye(M) - ARB(:,K*M+(1-M:0))*ARF(:,K*M+(1-M:0))]*PEB;
PE(:,K*M+(1:M)) = PEF;
end;
elseif Mode==1, %%%%% Levinson-Wiggens-Robinson (LWR) algorithm using biased correlation function
% ===== In [1,2] this algorithm is denoted "Multichannel Yule-Walker" ===== %
C(:,1:M) = C(:,1:M)/N;
PEF = C(:,1:M);
PEB = C(:,1:M);
for K=1:Pmax,
[C(:,K*M+(1:M)),n] = covm(Y(K+1:N,:),Y(1:N-K,:),'M');
C(:,K*M+(1:M)) = C(:,K*M+(1:M))/N;
D = C(:,K*M+(1:M));
for L = 1:K-1,
D = D - ARF(:,L*M+(1-M:0))*C(:,(K-L)*M+(1:M));
end;
ARF(:,K*M+(1-M:0)) = D / PEB;
ARB(:,K*M+(1-M:0)) = D'/ PEF;
for L = 1:K-1,
tmp = ARF(:,L*M+(1-M:0)) - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
ARB(:,(K-L)*M+(1-M:0)) = ARB(:,(K-L)*M+(1-M:0)) - ARB(:,K*M+(1-M:0))*ARF(:,L*M+(1-M:0));
ARF(:,L*M+(1-M:0)) = tmp;
end;
RCF(:,K*M+(1-M:0)) = ARF(:,K*M+(1-M:0));
RCB(:,K*M+(1-M:0)) = ARB(:,K*M+(1-M:0));
PEF = [eye(M) - ARF(:,K*M+(1-M:0))*ARB(:,K*M+(1-M:0))]*PEF;
PEB = [eye(M) - ARB(:,K*M+(1-M:0))*ARF(:,K*M+(1-M:0))]*PEB;
PE(:,K*M+(1:M)) = PEF;
end;
elseif Mode==6, %%%%% Levinson-Wiggens-Robinson (LWR) algorithm using unbiased correlation function
C(:,1:M) = C(:,1:M)/N;
PEF = C(:,1:M); %?% PEF = PE(:,1:M);
PEB = C(:,1:M);
for K = 1:Pmax,
[C(:,K*M+(1:M)),n] = covm(Y(K+1:N,:),Y(1:N-K,:),'M');
C(:,K*M+(1:M)) = C(:,K*M+(1:M))./n;
%C{K+1} = C{K+1}/N;
D = C(:,K*M+(1:M));
for L = 1:K-1,
D = D - ARF(:,L*M+(1-M:0))*C(:,(K-L)*M+(1:M));
end;
ARF(:,K*M+(1-M:0)) = D / PEB;
ARB(:,K*M+(1-M:0)) = D'/ PEF;
for L = 1:K-1,
tmp = ARF(:,L*M+(1-M:0)) - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
ARB(:,(K-L)*M+(1-M:0)) = ARB(:,(K-L)*M+(1-M:0)) - ARB(:,K*M+(1-M:0))*ARF(:,L*M+(1-M:0));
ARF(:,L*M+(1-M:0)) = tmp;
end;
RCF(:,K*M+(1-M:0)) = ARF(:,K*M+(1-M:0));
RCB(:,K*M+(1-M:0)) = ARB(:,K*M+(1-M:0));
PEF = [eye(M) - ARF(:,K*M+(1-M:0))*ARB(:,K*M+(1-M:0))]*PEF;
PEB = [eye(M) - ARB(:,K*M+(1-M:0))*ARF(:,K*M+(1-M:0))]*PEB;
PE(:,K*M+(1:M)) = PEF;
end;
elseif Mode==2 %%%%% Partial Correlation Estimation: Vieira-Morf Method [2] with unbiased covariance estimation
%===== In [1] this algorithm is denoted "Nutall-Strand with unbiased covariance" =====%
%C(:,1:M) = C(:,1:M)/N;
F = Y;
B = Y;
%PEF = C(:,1:M);
%PEB = C(:,1:M);
PEF = PE(:,1:M);
PEB = PE(:,1:M);
for K = 1:Pmax,
[D,n] = covm(F(K+1:N,:),B(1:N-K,:),'M');
D = D./n;
ARF(:,K*M+(1-M:0)) = D / PEB;
ARB(:,K*M+(1-M:0)) = D'/ PEF;
tmp = F(K+1:N,:) - B(1:N-K,:)*ARF(:,K*M+(1-M:0)).';
B(1:N-K,:) = B(1:N-K,:) - F(K+1:N,:)*ARB(:,K*M+(1-M:0)).';
F(K+1:N,:) = tmp;
for L = 1:K-1,
tmp = ARF(:,L*M+(1-M:0)) - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
ARB(:,(K-L)*M+(1-M:0)) = ARB(:,(K-L)*M+(1-M:0)) - ARB(:,K*M+(1-M:0))*ARF(:,L*M+(1-M:0));
ARF(:,L*M+(1-M:0)) = tmp;
end;
RCF(:,K*M+(1-M:0)) = ARF(:,K*M+(1-M:0));
RCB(:,K*M+(1-M:0)) = ARB(:,K*M+(1-M:0));
[PEF,n] = covm(F(K+1:N,:),F(K+1:N,:),'M');
PEF = PEF./n;
[PEB,n] = covm(B(1:N-K,:),B(1:N-K,:),'M');
PEB = PEB./n;
PE(:,K*M+(1:M)) = PEF;
end;
elseif Mode==5 %%%%% Partial Correlation Estimation: Vieira-Morf Method [2] with biased covariance estimation
%===== In [1] this algorithm is denoted "Nutall-Strand with biased covariance" ===== %
F = Y;
B = Y;
PEF = C(:,1:M);
PEB = C(:,1:M);
for K=1:Pmax,
[D,n] = covm(F(K+1:N,:),B(1:N-K,:),'M');
%D=D/N;
ARF(:,K*M+(1-M:0)) = D / PEB;
ARB(:,K*M+(1-M:0)) = D'/ PEF;
tmp = F(K+1:N,:) - B(1:N-K,:)*ARF(:,K*M+(1-M:0)).';
B(1:N-K,:) = B(1:N-K,:) - F(K+1:N,:)*ARB(:,K*M+(1-M:0)).';
F(K+1:N,:) = tmp;
for L = 1:K-1,
tmp = ARF(:,L*M+(1-M:0)) - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
ARB(:,(K-L)*M+(1-M:0)) = ARB(:,(K-L)*M+(1-M:0)) - ARB(:,K*M+(1-M:0))*ARF(:,L*M+(1-M:0));
ARF(:,L*M+(1-M:0)) = tmp;
end;
RCF(:,K*M+(1-M:0)) = ARF(:,K*M+(1-M:0));
RCB(:,K*M+(1-M:0)) = ARB(:,K*M+(1-M:0));
[PEB,n] = covm(B(1:N-K,:),B(1:N-K,:),'M');
%PEB = D/N;
[PEF,n] = covm(F(K+1:N,:),F(K+1:N,:),'M');
%PEF = D/N;
%PE(:,K*M+(1:M)) = PEF;
PE(:,K*M+(1:M)) = PEF./n;
end;
elseif Mode==3 %%%%% Partial Correlation Estimation: Nutall-Strand Method [2] with biased covariance estimation
% C(:,1:M) = C(:,1:M)/N;
F = Y;
B = Y;
PEF = C(:,1:M);
PEB = C(:,1:M);
for K=1:Pmax,
[D, n] = covm(F(K+1:N,:),B(1:N-K,:),'M');
D = D./N;
ARF(:,K*M+(1-M:0)) = 2*D / (PEB+PEF);
ARB(:,K*M+(1-M:0)) = 2*D'/ (PEF+PEB);
tmp = F(K+1:N,:) - B(1:N-K,:)*ARF(:,K*M+(1-M:0)).';
B(1:N-K,:) = B(1:N-K,:) - F(K+1:N,:)*ARB(:,K*M+(1-M:0)).';
F(K+1:N,:) = tmp;
for L = 1:K-1,
tmp = ARF(:,L*M+(1-M:0)) - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
ARB(:,(K-L)*M+(1-M:0)) = ARB(:,(K-L)*M+(1-M:0)) - ARB(:,K*M+(1-M:0))*ARF(:,L*M+(1-M:0));
ARF(:,L*M+(1-M:0)) = tmp;
end;
%RCF{K} = ARF{K};
RCF(:,K*M+(1-M:0)) = ARF(:,K*M+(1-M:0));
[PEF,n] = covm(F(K+1:N,:),F(K+1:N,:),'M');
PEF = PEF./N;
[PEB,n] = covm(B(1:N-K,:),B(1:N-K,:),'M');
PEB = PEB./N;
%PE(:,K*M+(1:M)) = PEF;
PE(:,K*M+(1:M)) = PEF./n;
end;
elseif Mode==7 %%%%% Partial Correlation Estimation: Nutall-Strand Method [2] with unbiased covariance estimation
%C(:,1:M) = C(:,1:M)/N;
F = Y;
B = Y;
%PEF = C(:,1:M);
%PEB = C(:,1:M);
PEF = PE(:,1:M);
PEB = PE(:,1:M);
for K = 1:Pmax,
[D,n] = covm(F(K+1:N,:),B(1:N-K,:),'M');
D = D./n;
ARF(:,K*M+(1-M:0)) = 2*D / (PEB+PEF);
ARB(:,K*M+(1-M:0)) = 2*D'/ (PEF+PEB);
tmp = F(K+1:N,:) - B(1:N-K,:)*ARF(:,K*M+(1-M:0)).';
B(1:N-K,:) = B(1:N-K,:) - F(K+1:N,:)*ARB(:,K*M+(1-M:0)).';
F(K+1:N,:) = tmp;
for L = 1:K-1,
tmp = ARF(:,L*M+(1-M:0)) - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
ARB(:,(K-L)*M+(1-M:0)) = ARB(:,(K-L)*M+(1-M:0)) - ARB(:,K*M+(1-M:0))*ARF(:,L*M+(1-M:0));
ARF(:,L*M+(1-M:0)) = tmp;
end;
%RCF{K} = ARF{K};
RCF(:,K*M+(1-M:0)) = ARF(:,K*M+(1-M:0));
[PEF,n] = covm(F(K+1:N,:),F(K+1:N,:),'M');
PEF = PEF./n;
[PEB,n] = covm(B(1:N-K,:),B(1:N-K,:),'M');
PEB = PEB./n;
%PE{K+1} = PEF;
PE(:,K*M+(1:M)) = PEF;
end;
elseif Mode==4, %%%%% Kay, not fixed yet.
fprintf('Warning MVAR: Mode=4 is broken.\n')
C(:,1:M) = C(:,1:M)/N;
K = 1;
[C(:,M+(1:M)),n] = covm(Y(2:N,:),Y(1:N-1,:));
C(:,M+(1:M)) = C(:,M+(1:M))/N; % biased estimate
D = C(:,M+(1:M));
ARF(:,1:M) = C(:,1:M)\D;
ARB(:,1:M) = C(:,1:M)\D';
RCF(:,1:M) = ARF(:,1:M);
RCB(:,1:M) = ARB(:,1:M);
PEF = C(:,1:M)*[eye(M) - ARB(:,1:M)*ARF(:,1:M)];
PEB = C(:,1:M)*[eye(M) - ARF(:,1:M)*ARB(:,1:M)];
for K=2:Pmax,
[C(:,K*M+(1:M)),n] = covm(Y(K+1:N,:),Y(1:N-K,:),'M');
C(:,K*M+(1:M)) = C(:,K*M+(1:M)) / N; % biased estimate
D = C(:,K*M+(1:M));
for L = 1:K-1,
D = D - C(:,(K-L)*M+(1:M))*ARF(:,L*M+(1-M:0));
end;
ARF(:,K*M+(1-M:0)) = PEB \ D;
ARB(:,K*M+(1-M:0)) = PEF \ D';
for L = 1:K-1,
tmp = ARF(:,L*M+(1-M:0)) - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
ARB(:,(K-L)*M+(1-M:0)) = ARB(:,(K-L)*M+(1-M:0)) - ARB(:,K*M+(1-M:0))*ARF(:,L*M+(1-M:0));
ARF(:,L*M+(1-M:0)) = tmp;
end;
RCF(:,K*M+(1-M:0)) = ARF(:,K*M+(1-M:0)) ;
RCB(:,K*M+(1-M:0)) = ARB(:,K*M+(1-M:0)) ;
PEF = PEF*[eye(M) - ARB(:,K*M+(1-M:0)) *ARF(:,K*M+(1-M:0)) ];
PEB = PEB*[eye(M) - ARF(:,K*M+(1-M:0)) *ARB(:,K*M+(1-M:0)) ];
PE(:,K*M+(1:M)) = PEF;
end;
elseif Mode==8, %%%%% Least Squares
ix = Pmax+1:N;
y = repmat(NaN,N-Pmax,M*Pmax);
for k = 1:Pmax,
y(:,k*M+[1-M:0]) = Y(ix-k,:);
end;
ix2 = ~any(isnan([Y(ix,:),y]),2);
ARF = Y(ix(ix2),:)'/y(ix2,:)';
PE = covm(Y(ix,:) - y*ARF');
RCF = zeros(M,M*Pmax);
elseif Mode==9, %%%%% Least Squares
%% FIXME
for K=1:Pmax,
[C(:,K*M+(1:M)),n] = covm(Y(K+1:N,:),Y(1:N-K,:),'M');
% C(K*M+(1:M),:) = C(K*M+(1:M),:)/N;
end;
ARF = C(:,1:M)'/C(:,M+1:end)';
PE = covm(C(:,1:M)'-ARF*C(:,M+1:end)');
RCF = zeros(M,M*Pmax);
elseif Mode==10, %%%%% ARFIT
try
RCF = [];
[w, ARF, PE] = arfit(Y, Pmax, Pmax, 'sbc', 'zero');
catch
ARF = zeros(M,M*Pmax);
RCF = ARF;
end;
elseif Mode==11, %%%%% de Waele 2003
[pc,R0] = burgv(reshape(Y',[M,1,N]),Pmax);
try,
[ARF,ARB,Pf,Pb,RCF,RCB] = pc2parv(pc,R0);
catch
[RCF,RCB,Pf,Pb] = pc2rcv(pc,R0);
[ARF,ARB,Pf,Pb] = pc2parv(pc,R0);
end;
ARF = -reshape(ARF(:,:,2:end),[M,M*Pmax]);
RCF = -reshape(RCF(:,:,2:end),[M,M*Pmax]);
PE = reshape(Pf,[M,M*(Pmax+1)]);
elseif Mode==12,
% this is equivalent to Mode==11
[pc,R0] = burgv(reshape(Y',[M,1,N]),Pmax);
[rcf,rcb,Pf,Pb] = pc2rcv(pc,R0);
%%%%% Convert reflection coefficients RC to autoregressive parameters
ARF = zeros(M,M*Pmax);
ARB = zeros(M,M*Pmax);
for K = 1:Pmax,
ARF(:,K*M+(1-M:0)) = -rcf(:,:,K+1);
ARB(:,K*M+(1-M:0)) = -rcb(:,:,K+1);
for L = 1:K-1,
tmp = ARF(:,L*M+(1-M:0)) - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
ARB(:,(K-L)*M+(1-M:0)) = ARB(:,(K-L)*M+(1-M:0)) - ARB(:,K*M+(1-M:0))*ARF(:,L*M+(1-M:0));
ARF(:,L*M+(1-M:0)) = tmp;
end;
end;
RCF = -reshape(rcf(:,:,2:end),[M,M*Pmax]);
PE = reshape(Pf,[M,M*(Pmax+1)]);
%%%%% EXPERIMENTAL VERSIONS %%%%%
elseif Mode==83, %%%%% de Waele 2003
[rcf,rcb,pc,R0] = burgv_as1(reshape(Y',[M,1,N]),Pmax);
%[ARF,ARB,Pf,Pb,RCF,RCB] = pc2parv(pc,R0);
RCF = -reshape(rcf(:,:,2:end),[M,M*Pmax]);
% PE = reshape(Pf,[M,M*(Pmax+1)]);
ARF = zeros(M,M*Pmax);
ARB = zeros(M,M*Pmax);
for K = 1:Pmax,
ARF(:,K*M+(1-M:0)) = -rcf(:,:,K+1);
ARB(:,K*M+(1-M:0)) = -rcb(:,:,K+1);
for L = 1:K-1,
tmp = ARF(:,L*M+(1-M:0)) - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
ARB(:,(K-L)*M+(1-M:0)) = ARB(:,(K-L)*M+(1-M:0)) - ARB(:,K*M+(1-M:0))*ARF(:,L*M+(1-M:0));
ARF(:,L*M+(1-M:0)) = tmp;
end;
end;
elseif Mode==84, %%%%% de Waele 2003
[rcf,rcb,pc,R0] = burgv_as2(reshape(Y',[M,1,N]),Pmax);
%[ARF,ARB,Pf,Pb,RCF,RCB] = pc2parv(pc,R0);
RCF = -reshape(rcf(:,:,2:end),[M,M*Pmax]);
% PE = reshape(Pf,[M,M*(Pmax+1)]);
ARF = zeros(M,M*Pmax);
ARB = zeros(M,M*Pmax);
for K = 1:Pmax,
ARF(:,K*M+(1-M:0)) = -rcf(:,:,K+1);
ARB(:,K*M+(1-M:0)) = -rcb(:,:,K+1);
for L = 1:K-1,
tmp = ARF(:,L*M+(1-M:0)) - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
ARB(:,(K-L)*M+(1-M:0)) = ARB(:,(K-L)*M+(1-M:0)) - ARB(:,K*M+(1-M:0))*ARF(:,L*M+(1-M:0));
ARF(:,L*M+(1-M:0)) = tmp;
end;
end;
elseif Mode==85, %%%%% de Waele 2003
[rcf,rcb,pc,R0] = burgv_as1(reshape(Y',[M,1,N]),Pmax);
%[ARF,ARB,Pf,Pb,RCF,RCB] = pc2parv(pc,R0);
RCF = -reshape(rcf(:,:,2:end),[M,M*Pmax]);
% PE = reshape(Pf,[M,M*(Pmax+1)]);
ARF = zeros(M,M*Pmax);
ARB = zeros(M,M*Pmax);
for K = 1:Pmax,
ARF(:,K*M+(1-M:0)) = -rcf(:,:,K+1)';
ARB(:,K*M+(1-M:0)) = -rcb(:,:,K+1)';
for L = 1:K-1,
tmp = ARF(:,L*M+(1-M:0)) - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
ARB(:,(K-L)*M+(1-M:0)) = ARB(:,(K-L)*M+(1-M:0)) - ARB(:,K*M+(1-M:0))*ARF(:,L*M+(1-M:0));
ARF(:,L*M+(1-M:0)) = tmp;
end;
end;
elseif Mode==86, %%%%% de Waele 2003
[rcf,rcb,pc,R0] = burgv_as2(reshape(Y',[M,1,N]),Pmax);
%[ARF,ARB,Pf,Pb,RCF,RCB] = pc2parv(pc,R0);
RCF = -reshape(rcf(:,:,2:end),[M,M*Pmax]);
% PE = reshape(Pf,[M,M*(Pmax+1)]);
ARF = zeros(M,M*Pmax);
ARB = zeros(M,M*Pmax);
for K = 1:Pmax,
ARF(:,K*M+(1-M:0)) = -rcf(:,:,K+1)';
ARB(:,K*M+(1-M:0)) = -rcb(:,:,K+1)';
for L = 1:K-1,
tmp = ARF(:,L*M+(1-M:0)) - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
ARB(:,(K-L)*M+(1-M:0)) = ARB(:,(K-L)*M+(1-M:0)) - ARB(:,K*M+(1-M:0))*ARF(:,L*M+(1-M:0));
ARF(:,L*M+(1-M:0)) = tmp;
end;
end;
elseif Mode==90;
% similar to MODE=0
%% not recommended
C(:,1:M) = C(:,1:M)/N;
F = Y;
B = Y;
PEF = PE(:,1:M); %CHANGED%
PEB = PE(:,1:M); %CHANGED%
for K=1:Pmax,
[D,n] = covm(Y(K+1:N,:),Y(1:N-K,:),'M');
D = D/N;
ARF(:,K*M+(1-M:0)) = D/PEB;
ARB(:,K*M+(1-M:0)) = D'/PEF;
tmp = F(K+1:N,:) - B(1:N-K,:)*ARF(:,K*M+(1-M:0))';
B(1:N-K,:) = B(1:N-K,:) - F(K+1:N,:)*ARB(:,K*M+(1-M:0))';
F(K+1:N,:) = tmp;
for L = 1:K-1,
tmp = ARF(:,L*M+(1-M:0)) - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
ARB(:,(K-L)*M+(1-M:0)) = ARB(:,(K-L)*M+(1-M:0)) - ARB(:,K*M+(1-M:0))*ARF(:,L*M+(1-M:0));
ARF(:,L*M+(1-M:0)) = tmp;
end;
RCF(:,K*M+(1-M:0)) = ARF(:,K*M+(1-M:0));
RCB(:,K*M+(1-M:0)) = ARB(:,K*M+(1-M:0));
PEF = [eye(M) - ARF(:,K*M+(1-M:0))*ARB(:,K*M+(1-M:0))]*PEF;
PEB = [eye(M) - ARB(:,K*M+(1-M:0))*ARF(:,K*M+(1-M:0))]*PEB;
PE(:,K*M+(1:M)) = PEF;
end;
elseif Mode==91, %%%%% Levinson-Wiggens-Robinson (LWR) algorithm using biased correlation function
% ===== In [1,2] this algorithm is denoted "Multichannel Yule-Walker" ===== %
% similar to MODE=1
%% not recommended
C(:,1:M) = C(:,1:M)/N;
PEF = PE(:,1:M); %CHANGED%
PEB = PE(:,1:M); %CHANGED%
for K=1:Pmax,
[C(:,K*M+(1:M)),n] = covm(Y(K+1:N,:),Y(1:N-K,:),'M');
C(:,K*M+(1:M)) = C(:,K*M+(1:M))/N;
D = C(:,K*M+(1:M));
for L = 1:K-1,
D = D - ARF(:,L*M+(1-M:0))*C(:,(K-L)*M+(1:M));
end;
ARF(:,K*M+(1-M:0)) = D / PEB;
ARB(:,K*M+(1-M:0)) = D'/ PEF;
for L = 1:K-1,
tmp = ARF(:,L*M+(1-M:0)) - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
ARB(:,(K-L)*M+(1-M:0)) = ARB(:,(K-L)*M+(1-M:0)) - ARB(:,K*M+(1-M:0))*ARF(:,L*M+(1-M:0));
ARF(:,L*M+(1-M:0)) = tmp;
end;
RCF(:,K*M+(1-M:0)) = ARF(:,K*M+(1-M:0));
RCB(:,K*M+(1-M:0)) = ARB(:,K*M+(1-M:0));
PEF = [eye(M) - ARF(:,K*M+(1-M:0))*ARB(:,K*M+(1-M:0))]*PEF;
PEB = [eye(M) - ARB(:,K*M+(1-M:0))*ARF(:,K*M+(1-M:0))]*PEB;
PE(:,K*M+(1:M)) = PEF;
end;
elseif Mode==96, %%%%% Levinson-Wiggens-Robinson (LWR) algorithm using unbiased correlation function
% similar to MODE=6
%% not recommended
C(:,1:M) = C(:,1:M)/N;
PEF = PE(:,1:M); %CHANGED%
PEB = PE(:,1:M); %CHANGED%
for K = 1:Pmax,
[C(:,K*M+(1:M)),n] = covm(Y(K+1:N,:),Y(1:N-K,:),'M');
C(:,K*M+(1:M)) = C(:,K*M+(1:M))./n;
D = C(:,K*M+(1:M));
for L = 1:K-1,
D = D - ARF(:,L*M+(1-M:0))*C(:,(K-L)*M+(1:M));
end;
ARF(:,K*M+(1-M:0)) = D / PEB;
ARB(:,K*M+(1-M:0)) = D'/ PEF;
for L = 1:K-1,
tmp = ARF(:,L*M+(1-M:0)) - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
ARB(:,(K-L)*M+(1-M:0)) = ARB(:,(K-L)*M+(1-M:0)) - ARB(:,K*M+(1-M:0))*ARF(:,L*M+(1-M:0));
ARF(:,L*M+(1-M:0)) = tmp;
end;
RCF(:,K*M+(1-M:0)) = ARF(:,K*M+(1-M:0));
RCB(:,K*M+(1-M:0)) = ARB(:,K*M+(1-M:0));
PEF = [eye(M) - ARF(:,K*M+(1-M:0))*ARB(:,K*M+(1-M:0))]*PEF;
PEB = [eye(M) - ARB(:,K*M+(1-M:0))*ARF(:,K*M+(1-M:0))]*PEB;
PE(:,K*M+(1:M)) = PEF;
end;
elseif Mode<100,
fprintf('Warning MVAR: Mode=%i not supported\n',Mode);
%%%%% IMPUTATION METHODS %%%%%
else
Mode0 = rem(Mode,100);
if ((Mode0>=10) & (Mode0<20)),
Mode0 = 1;
end;
if 0,
elseif Mode>=2400, % forward and backward
% assuming that past missing values are already IMPUTED with the prediction value + innovation process
% no decaying
[ARF,RCF,PE2] = mvar(Y, Pmax, Mode0);
c = chol( PE2 (:, M*Pmax+(1:M)));
Y1 = Y;
Y1(1,isnan(Y1(1,:))) = 0;
z = [];
for k = 2:size(Y,1)
[z1,z] = mvfilter(ARF,eye(M),Y1(k-1,:)',z);
ix = isnan(Y1(k,:));
z1 = z1 + (randn(1,M)*c)';
Y1(k,ix) = z1(ix);
end;
Y2 = flipud(Y);
[ARB,RCF,PE] = mvar(Y2, Pmax, Mode0);
Y2(1,isnan(Y2(1,:))) = 0;
z = [];
for k = 2:size(Y2,1)
[z2,z] = mvfilter(ARB,eye(M),Y2(k-1,:)',z);
ix = isnan(Y(size(Y,1)-k+1,:));
z2 = z2 + (randn(1,M)*c)';
Y2(k,ix) = (z2(ix)' + Y2(k,ix))/2;
end;
Y2 = flipud(Y2);
Z = (Y2+Y1)/2;
Y(isnan(Y)) = Z(isnan(Y));
[ARF,RCF,PE] = mvar(Y, Pmax, rem(Mode,100));
elseif Mode>=2300, % backward prediction
% assuming that past missing values are already IMPUTED with the prediction value + innovation process
% no decaying
Y = flipud(Y);
[ARB,RCF,PE] = mvar(Y, Pmax, Mode0);
c = chol(PE(:,M*Pmax+(1:M)));
Y1 = Y;
Y1(1,isnan(Y1(1,:))) = 0;
z = [];
for k = 2:size(Y,1)
[z1,z] = mvfilter(ARB,eye(M),Y1(k-1,:)',z);
ix = isnan(Y1(k,:));
z1 = z1 + (randn(1,M)*c)';
Y1(k,ix) = z1(ix);
end;
Y1 = flipud(Y1);
[ARF,RCF,PE] = mvar(Y1, Pmax, rem(Mode,100));
elseif Mode>=2200, % forward predictions,
% assuming that past missing values are already IMPUTED with the prediction value + innovation process
% no decaying
[ARF,RCF,PE] = mvar(Y, Pmax, Mode0);
c = chol(PE(:,M*Pmax+(1:M)));
Y1 = Y;
Y1(1,isnan(Y1(1,:))) = 0;
z = [];
for k = 2:size(Y,1)
[z1,z] = mvfilter(ARF,eye(M),Y1(k-1,:)',z);
ix = isnan(Y1(k,:));
z1 = z1 + (randn(1,M)*c)';
Y1(k,ix) = z1(ix);
end;
[ARF,RCF,PE] = mvar(Y1, Pmax, rem(Mode,100));
elseif Mode>=1400, % forward and backward
%assuming that past missing values are already IMPUTED with the prediction value
[ARF,RCF,PE] = mvar(Y, Pmax, Mode0);
Y1 = Y;
Y1(1,isnan(Y1(1,:))) = 0;
z = [];
for k = 2:size(Y,1)
[z1,z] = mvfilter(ARF,eye(M),Y1(k-1,:)',z);
ix = isnan(Y1(k,:));
Y1(k,ix) = z1(ix);
end;
Y2 = flipud(Y);
[ARB,RCF,PE] = mvar(Y2, Pmax, Mode0);
Y2(1,isnan(Y2(1,:))) = 0;
z = [];
for k = 2:size(Y2,1)
[z2,z] = mvfilter(ARB,eye(M),Y2(k-1,:)',z);
ix = isnan(Y2(k,:));
Y2(k,ix) = z2(ix)';
end;
Y2 = flipud(Y2);
Z = (Y2+Y1)/2;
Y(isnan(Y)) = Z(isnan(Y));
[ARF,RCF,PE] = mvar(Y, Pmax, rem(Mode,100));
elseif Mode>=1300, % backward prediction
Y = flipud(Y);
[ARB,RCF,PE] = mvar(Y, Pmax, Mode0);
Y2 = Y;
Y2(1,isnan(Y2(1,:))) = 0;
z = [];
for k = 2:size(Y,1)
[z2,z] = mvfilter(ARB,eye(M),Y2(k-1,:)',z);
ix = isnan(Y2(k,:));
Y2(k,ix) = z2(ix);
end;
Y2 = flipud(Y2);
[ARF,RCF,PE] = mvar(Y2, Pmax, rem(Mode,100));
elseif Mode>=1200, % forward predictions,
%assuming that past missing values are already IMPUTED with the prediction value
[ARF,RCF,PE] = mvar(Y, Pmax, Mode0);
Y1 = Y;
Y1(1,isnan(Y1(1,:))) = 0;
z = [];
for k = 2:size(Y,1)
[z1,z] = mvfilter(ARF,eye(M),Y1(k-1,:)',z);
ix = isnan(Y1(k,:));
Y1(k,ix) = z1(ix);
end;
[ARF,RCF,PE] = mvar(Y1, Pmax, rem(Mode,100));
elseif Mode>=400, % forward and backward
% assuming that "past" missing values are ZERO
[ARF,RCF,PE] = mvar(Y, Pmax, Mode0);
Y1 = Y;
Y1(isnan(Y)) = 0;
Z1 = mvfilter([zeros(M),ARF],eye(M),Y1')';
Y1(isnan(Y)) = Z1(isnan(Y));
Y = flipud(Y);
[ARB,RCF,PE] = mvar(Y, Pmax, Mode0);
Y2 = Y;
Y2(isnan(Y)) = 0;
Z2 = mvfilter([zeros(M),ARB],eye(M),Y2')';
Y2(isnan(Y)) = Z2(isnan(Y));
Y2 = flipud(Y2);
[ARF,RCF,PE] = mvar((Y1+Y2)/2, Pmax, rem(Mode,100));
elseif Mode>=300, % backward prediction
% assuming that "past" missing values are ZERO
Y = flipud(Y);
[ARB,RCF,PE] = mvar(Y, Pmax, Mode0);
Y2 = Y;
Y2(isnan(Y)) = 0;
Z2 = mvfilter([zeros(M),ARB],eye(M),Y2')';
Y2(isnan(Y)) = Z2(isnan(Y));
Y2 = flipud(Y2);
[ARF,RCF,PE] = mvar(Y2, Pmax, rem(Mode,100));
elseif Mode>=200,
% forward predictions, assuming that past missing values are ZERO
[ARF,RCF,PE] = mvar(Y, Pmax, Mode0);
Y1 = Y;
Y1(isnan(Y)) = 0;
Z1 = mvfilter([zeros(M),ARF],eye(M),Y1')';
Y1(isnan(Y)) = Z1(isnan(Y));
[ARF,RCF,PE] = mvar(Y1, Pmax, rem(Mode,100));
elseif Mode>=100,
[ARF,RCF,PE] = mvar(Y, Pmax, Mode0);
Z1 = mvfilter(ARF,eye(M),Y')';
Z1 = [zeros(1,M); Z1(1:end-1,:)];
Y(isnan(Y)) = Z1(isnan(Y));
[ARF,RCF,PE] = mvar(Y, Pmax, rem(Mode,100));
end;
end;
if any(ARF(:)==inf),
% Test for matrix division bug.
% This bug was observed in LNX86-ML5.3, 6.1 and 6.5, but not in SGI-ML6.5, LNX86-ML6.5, Octave 2.1.35-40; Other platforms unknown.
p = 3;
FLAG_MATRIX_DIVISION_ERROR = ~all(all(isnan(repmat(0,p)/repmat(0,p)))) | ~all(all(isnan(repmat(nan,p)/repmat(nan,p))));
if FLAG_MATRIX_DIVISION_ERROR,
%fprintf(2,'### Warning MVAR: Bug in Matrix-Division 0/0 and NaN/NaN yields INF instead of NAN. Workaround is applied.\n');
warning('MVAR: bug in Matrix-Division 0/0 and NaN/NaN yields INF instead of NAN. Workaround is applied.');
%%%%% Workaround
ARF(ARF==inf)=NaN;
RCF(RCF==inf)=NaN;
end;
end;
%MAR = zeros(M,M*Pmax);
DC = zeros(M);
for K = 1:Pmax,
% VAR{K+1} = -ARF(:,K*M+(1-M:0))';
DC = DC + ARF(:,K*M+(1-M:0))'.^2; %DC meausure [3]
end;
