invert_akaike.m

用于模拟时变非平稳的ARMA过程

function AM= invert_akaike(rPM, BM)
% function AM= invert_akaike(rPM, BM)
%   This file is part of the TFPM toolbox v0.5 (c)
%   [e-mail] and underlies the GPL.
% 
% Implements the Akaike algorithm (SIAMAM 24(2), mar73,
% pp234) for solving the block-Toeplitz system 
%
% btoep(rPM)*AP = BP; M...block size, P...super size,
% size(rPM)= [M, M(2P-1)].

%% Dimensions:
PP= size(rPM);
M= PP(1);
P= (PP(2)/M+1)/2;

% Permute rPM and transpose BM!!
rrPM= [];
for p= 1:2*P-1
   rrPM= [rPM(:, (p-1)*M+1:p*M).' rrPM];
end;
rPM= rrPM;
BM= BM.';

rhs_dim= size(BM, 1);% Dimension of right-hand-side block-vector

%% The output tensor
AM= zeros(P*M, rhs_dim, P);

%% Here is the algorithm for the system [A1..AP] T= [B1..BP]!!
% Init (p= 0):
qinv= rPM(:, (P-1)*M+1:P*M);
sinv= qinv;
am  = BM(:, 1:M)*inv(qinv);
AM(:, :, 1)= [am.'; zeros(M*(P-1), rhs_dim)];

if(P==1)
   return;
end;

Xp= -rPM(:, (P-2)*M+1:(P-1)*M) * inv(qinv);%X1
Yp= -rPM(:, (P  )*M+1:(P+1)*M) * inv(qinv);%Y1

ap= rPM(:, (P-2)*M+1:(P-1)*M);             %a1
rp= rPM(:, (P  )*M+1:(P+1)*M);             %r1

qinv_new= rPM(:, (P-1)*M+1:P*M) - rPM(:, (P-2)*M+1:(P-1)*M)*inv(qinv)*rPM(:, (P  )*M+1:(P+1)*M);
sinv_new= rPM(:, (P-1)*M+1:P*M) - rPM(:, (P  )*M+1:(P+1)*M)*inv(qinv)*rPM(:, (P-2)*M+1:(P-1)*M);

for p= 1:P-2
   Ep= block_exchange(M, p);
   qinv= qinv_new;
   sinv= sinv_new;

   App= zeros(rhs_dim, M);
   for m= 1:p
      App= App + am(:, (m-1)*M+1:m*M)*rPM(:, (P+p-m)*M+1:(P+p+1-m)*M);
   end;
   App= (BM(:, p*M+1:(p+1)*M)-App)*inv(qinv);
   am_new= [];
   for m= 1:p
      am_new= [am_new am(:, (m-1)*M+1:m*M) + App*Xp(:, (m-1)*M+1:m*M)];
   end;
   am= [am_new App];
   AM(:, :, p+1)= [am.'; zeros(M*(P-1-p), rhs_dim)];

   Xp_new= [zeros(M) Xp] - (rPM(:, (P-2-p)*M+1:(P-1-p)*M) + Xp*   ap)*inv(sinv)*[eye(M) Yp];
   Yp_new= [Yp zeros(M)] - (rPM(:, (P  +p)*M+1:(P+1+p)*M) + Yp*Ep*rp)*inv(qinv)*[Xp eye(M)];
   
   ap= [ap; rPM(:, (P-2-p)*M+1:(P-1-p)*M)];
   rp= [rp; rPM(:, (P  +p)*M+1:(P+1+p)*M)];
   
   qinv_new= qinv - Xp_new(:, 1:M)          *Yp_new(:, p*M+1:(p+1)*M) * qinv;
   sinv_new= sinv - Yp_new(:, p*M+1:(p+1)*M)*Xp_new(:, 1:M)           * sinv;

   Xp= Xp_new;
   Yp= Yp_new;
end;

p= P-1;
qinv= qinv_new;
sinv= sinv_new;
App= zeros(rhs_dim, M);
for m= 1:p
   App= App + am(:, (m-1)*M+1:m*M)*rPM(:, (P+p-m)*M+1:(P+p+1-m)*M);
end;
App= (BM(:, p*M+1:(p+1)*M)-App)*inv(qinv);
am_new= [];
for m= 1:p
   am_new= [am_new am(:, (m-1)*M+1:m*M) + App*Xp(:, (m-1)*M+1:m*M)];
end;
am= [am_new App];
AM(:, :, P)= am.';