invert_waxkailath.m

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

function AM= invert_waxkailath(rPM, BM)
% function AM= invert_waxkailath(rPM, BM)
%   This file is part of the TFPM toolbox v0.5 (c)
%   [e-mail] and underlies the GPL.
% 
% Implements the Wax-Kailath Algorithm (IEEEASSP 31(5), oct83,
% pp1218) for solving the Toeplitz/Block-Toeplitz system 
%
% toeptoep(rPM)*AM= BM; P...block size, M...super size,
% size(rPM)= [2M-1, 2P-1].
%
% ATTENTION: (JP*R).' ~= JP*(R.')!!

% Dimensions:
PM= size(rPM);
P= (PM(2)+1)/2;
M= (PM(1)+1)/2;
rhs_dim= size(BM, 2);% Dimension of right-hand-side block-vector

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

% Definitions:
NP = zeros(P);
IP = eye(P);
JP = rot90(IP);
JPn= JP;

% Init (n=0):
Alpha= toep(rPM(M, :));%Alpha_0
Delta= BM(1:P, :);     %Delta_0
if(M==1)
   AM= inv(Alpha)*Delta;
   return;
end;

Beta = JP*toep(rPM(M+1, :));      %Beta_0
Gamma=   (toep(rPM(M-1, :))*JP).';%Gamma_0
am= inv(Alpha)*Delta;             %AM_0
AM(:, :, 1)= [am; zeros(P*(M-1), rhs_dim)];

% Init (n=1):
Vn= -JP*inv(toep(rPM(M, :)))  *Gamma;%V_1
Wn= -JP*inv(toep(rPM(M, :))).'*Beta ;%W_1
Rp= toep(rPM(M+1, :));               %R_{+1}
Tm= toep(rPM(M-1, :));               %R^T_{-1}

for n= 1:M-2
   Alpha= Alpha - Gamma.'*inv(Alpha).'*Beta;
   Beta = Wn.'*JPn*Rp + JP*toep(rPM(M+1+n, :));
   Gamma= Vn.'*(Tm*JPn).'  + (JP*toep(rPM(M-1-n, :)).');
   Delta= Vn.'*BM(1:P*n, :) + BM(P*n+1:P*(n+1), :);

   Vhat_new= [JPn*Vn; NP] - [Wn; IP]*inv(Alpha)  *Gamma;
   What_new= [JPn*Wn; NP] - [Vn; IP]*inv(Alpha).'*Beta ;

   am= [am; zeros(P, rhs_dim)] + [Wn; IP]*inv(Alpha)*Delta;
   AM(:, :, n+1)= [am; zeros(P*(M-1-n), rhs_dim)];

   JPn= rot90(eye((n+1)*P));

   Vn= JPn*Vhat_new;
   Wn= JPn*What_new;

   Rp= [toep(rPM(M+1+n, :)); Rp];
   Tm= [toep(rPM(M-1-n, :))  Tm];
end;
Alpha= Alpha-Gamma.'*inv(Alpha).'*Beta;
Delta= Vn.'*BM(1:P*(M-1), :) + BM(P*(M-1)+1:P*M, :);
am= [am; zeros(P, rhs_dim)] + [Wn; IP]*inv(Alpha)*Delta;
AM(:, :, M)= am;